diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 0000000..ed8f062 --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,32 @@ +# .readthedocs.yaml +# Read the Docs configuration file +# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details + +# Required +version: 2 + +# Set the OS, Python version and other tools you might need +build: + os: ubuntu-22.04 + tools: + python: "3.12" + # You can also specify other tool versions: + # nodejs: "19" + # rust: "1.64" + # golang: "1.19" + +# Build documentation in the "doc/" directory with Sphinx +sphinx: + configuration: doc/conf.py + +# Optionally build your docs in additional formats such as PDF and ePub +# formats: +# - pdf +# - epub + +# Optional but recommended, declare the Python requirements required +# to build your documentation +# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html +python: + install: + - requirements: doc/requirements.txt diff --git a/CMakeLists.txt b/CMakeLists.txt index bfb2a60..e6d19ec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.18) project( ufs_fire_behavior VERSION 0.0.1 - LANGUAGES Fortran CXX) + LANGUAGES Fortran CXX C ) # user defined cmake modules list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake/modules/") @@ -37,7 +37,11 @@ set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/mod") install(DIRECTORY ${CMAKE_Fortran_MODULE_DIRECTORY} DESTINATION ${CMAKE_INSTALL_PREFIX}) # third party libraries -find_package(NetCDF REQUIRED Fortran) +if(${CMAKE_PROJECT_NAME} STREQUAL "WRF") + find_package(netCDF-Fortran REQUIRED) +else() + find_package(NetCDF REQUIRED Fortran) +endif() # turn on DM_PARALLEL preprocessor directive if(DM_PARALLEL) @@ -45,6 +49,11 @@ if(DM_PARALLEL) add_compile_definitions(firelib PUBLIC DM_PARALLEL=1) endif() +if(OPENMP) + find_package(OpenMP REQUIRED) + message(STATUS "OpenMP enabled") +endif() + #------------------------------------------------------------------------------ # Add source files @@ -55,9 +64,11 @@ list(APPEND _state_files state/state_mod.F90 list(APPEND _share_files share/proj_lc_mod.F90 share/datetime_mod.F90 + share/emis_mod.F90 share/fuel_mod.F90 share/ros_mod.F90 share/fmc_mod.F90 + share/interp_mod.F90 share/constants_mod.F90) # IO @@ -88,14 +99,31 @@ list(APPEND _driver_files driver/advance_mod.F90 # configure fire behavior library add_library(firelib STATIC ${_state_files} ${_share_files} ${_io_files} ${_wrffire_physics} ${_driver_files}) -target_include_directories(firelib PUBLIC ${NetCDF_INCLUDE_DIRS}) -target_link_libraries(firelib PUBLIC $) +if(${CMAKE_PROJECT_NAME} STREQUAL "WRF") + target_link_libraries(firelib PUBLIC netCDF::netcdff ) +else() + target_include_directories(firelib PUBLIC ${NetCDF_INCLUDE_DIRS}) + target_link_libraries(firelib PUBLIC $) +endif() +if (DM_PARALLEL) + target_link_libraries(firelib PUBLIC $) +endif() +if(OPENMP) + target_link_libraries(firelib PUBLIC $) +endif() +set_target_properties( + firelib + PROPERTIES + Fortran_FORMAT FREE + ) # configure fire behavior executable add_executable(fire_behavior.exe driver/fire_behavior.F90) add_dependencies(fire_behavior.exe firelib) target_link_libraries(fire_behavior.exe PUBLIC firelib) -target_link_libraries(fire_behavior.exe PUBLIC NetCDF::NetCDF_Fortran) +if(NOT ${CMAKE_PROJECT_NAME} STREQUAL "WRF") + target_link_libraries(fire_behavior.exe PUBLIC NetCDF::NetCDF_Fortran) +endif() install(TARGETS fire_behavior.exe DESTINATION bin) if(NUOPC OR ESMX) diff --git a/cmake/modules/FireBuildSettings.cmake b/cmake/modules/FireBuildSettings.cmake index ed897b8..1c05de1 100644 --- a/cmake/modules/FireBuildSettings.cmake +++ b/cmake/modules/FireBuildSettings.cmake @@ -1,11 +1,11 @@ macro(set_fire_defaults) - if(DM_PARALLEL STREQUAL "") + if("${DM_PARALLEL}" STREQUAL "") set(DM_PARALLEL ON) endif() - if(NUOPC STREQUAL "") + if("${NUOPC}" STREQUAL "") set(NUOPC OFF) endif() - if(ESMX STREQUAL "") + if("${ESMX}" STREQUAL "") set(ESMX OFF) endif() endmacro() @@ -14,10 +14,12 @@ function(set_fire_cache) set (DM_PARALLEL ${DM_PARALLEL} CACHE BOOL "Turn on MPI build" FORCE) set (NUOPC ${NUOPC} CACHE BOOL "Build NUOPC cap" FORCE) set (ESMX ${ESMX} CACHE BOOL "Build ESMX application" FORCE) + set (OPENMP ${OPENMP} CACHE BOOL "Enable OpenMP parallelization" FORCE) endfunction() function(print_fire_settings) message (STATUS "FIRE BEHAVIOR BUILD SETTINGS\n" + "\tOPENMP: ${OPENMP}\n" "\tDM_PARALLEL: ${DM_PARALLEL}\n" "\tNUOPC: ${NUOPC}\n" "\tESMX: ${ESMX}") diff --git a/compile.sh b/compile.sh index 0ad3459..4e7fac3 100755 --- a/compile.sh +++ b/compile.sh @@ -23,6 +23,8 @@ usage () { printf " build NUOPC library and module\n" printf " --esmx, -x\n" printf " build ESMX application (includes NUOPC)\n" + printf " --openmp-on\n" + printf " enable OpenMP parallelization\n" printf " --prefix=INSTALL_PREFIX\n" printf " installation prefix\n" printf " --verbose, -v\n" @@ -81,6 +83,7 @@ VERBOSE=false TEST=false TEST_NAME="" CLEAN=false +OPENMP=false #------------------------------------------------------------------------------ @@ -129,6 +132,7 @@ while :; do --test=?*|-t=?*) TEST=true; TEST_NAME=${1#*=} ;; --test=) TEST=true ;; --clean) CLEAN=true ;; + --openmp-on) OPENMP=true ;; --clean=?*) printf "ERROR: $1 argument ignored.\n"; usage; exit 1 ;; --clean=) printf "ERROR: $1 argument ignored.\n"; usage; exit 1 ;; -?*) printf "ERROR: Unknown option $1\n"; usage; exit 1 ;; @@ -213,6 +217,11 @@ if [ "${ESMX}" = true ]; then else CMAKE_SETTINGS+=("-DESMX=OFF") fi +if [ "${OPENMP}" = true ]; then + CMAKE_SETTINGS+=("-DOPENMP=ON") +else + CMAKE_SETTINGS+=("-DOPENMP=OFF") +fi cmake -S${FIRE_DIR} -B${BUILD_DIR} ${CMAKE_SETTINGS[@]} if [ "$?" != "0" ]; then echo "$0 Failed: (cmake)" diff --git a/doc/CFBM-NUOPC.jpeg b/doc/CFBM-NUOPC.jpeg new file mode 100644 index 0000000..df0a450 Binary files /dev/null and b/doc/CFBM-NUOPC.jpeg differ diff --git a/doc/Configuration.rst b/doc/Configuration.rst new file mode 100644 index 0000000..473e2fe --- /dev/null +++ b/doc/Configuration.rst @@ -0,0 +1,252 @@ +.. _Configuration: + +============================================== +Configuring the Community Fire Behavior Module +============================================== + + +.. _domain_config: + +Configuring a domain with the WRF Pre-processing System (WPS) +============================================================= + +Because the CFBM was originally developed as part of the WRF model, creating a domain must be done using the WRF Pre-processing System (WPS). These instructions can be found in `the WRF Users Guide `_. To run the CFBM with the UFS or with WRF data, users need to provide a geo_em.d01.nc file containting the interpolated static data fields. + +Future releases will include a method for creating domains without needing to compile WPS. + +.. _namelist: + +Namelist Configuration +====================== + +The options specific to the CFBM are controlled by a :term:`namelist` file ``namelist.fire``. This namelist file consists of three sections: ``&time``, ``&atm``, and ``&fire``. The available options in each section are described below. + +Example namelists can be found in the various test subdirectories under the ``tests/`` directory. + + +&time +--------------------------------- + +``start_year``: *integer* (**Required**) + Start year of the simulation. + +``start_month``: *integer* (**Required**) + Start month of the simulation. + +``start_day``: *integer* (**Required**) + Start day of the simulation. + +``start_hour``: *integer* (**Required**) + Start hour of the simulation. + +``start_minute``: *integer* (**Required**) + Start minute of the simulation. + +``start_second``: *integer* (**Required**) + Start second of the simulation. + +``end_year``: *integer* (**Required**) + End year of the simulation. + +``end_month``: *integer* (**Required**) + End month of the simulation. + +``end_day``: *integer* (**Required**) + End day of the simulation. + +``end_hour``: *integer* (**Required**) + End hour of the simulation. + +``end_minute``: *integer* (**Required**) + End minute of the simulation. + +``end_second``: *integer* (**Required**) + End second of the simulation. + +``dt``: *real* (Default: ``2.0``) + Atmospheric time step in seconds. + +``interval_output``: *integer* (**Required**) + [Units: s] + Specifies the time interval (in seconds) for writing to the history output files + +``num_tiles``: *integer* (Default: ``1``) + Number of tiles for MPI domain decomposition. Not yet implemented. + + +&atm +---- +``kde``: *integer* (Default: ``2``) + Number of vertical levels for the atmospheric simulation + +``interval_atm``: *integer* (Default: ``0``) + [Units: s] + Time interval (in seconds) for incoming atmospheric data. When running a coupled model, this value represents the atmospheric timestep. In offline mode, it determines the frequency of reading atmospheric data from the input file. + + +&fire +----- + +``fire_print_msg``: *integer* (Default: ``0``) + Debug print level for the fire module. + Levels greater than 0 will print extra messages at run time. + + 0: no extra prints + + 1: Extra prints + + 2: More extra prints + + 3: Even more extra prints + +``fire_atm_feedback``: *real* (Default: ``1.0``) + Multiplier for heat fluxes from the fire to the atmosphere. + 0.0: one-way (atmosphere --> fire) coupling. + + 1.0: normal two-way coupling. + + Intermediate values will vary the amount of forcing provided from the fire to the dynamical core. + +``fire_upwinding``: *integer* (Default: ``9``) + This option controls the type of upwinding scheme used for calculating the normal spread of the fire front. The choice of upwinding scheme significantly impacts the accuracy of fire spread simulations. Higher-order schemes, like WENO3 and WENO5, generally offer better accuracy but can be more computationally expensive. + 0: Central Difference: Uses central differences for calculating gradients, combining left- and right-sided differences for both x- and y-directions to compute a central gradient approximation. + + 1: Standard: Employs an upwind scheme, selecting between left- and right-sided differences based on flow direction. + + 2: Godunov: The Godunov scheme is a first-order upwind scheme based on Osher & Fedkiw + + 3: ENO1: The First-Order Essentially Non-Oscillatory (ENO1) scheme uses the smoothest stencil to avoid sharp gradients, which can lead to underestimations of fire area and errors in the rate of spread. + + 4: Sethian scheme :cite:`SethianMethod` + + 5: 2nd-order: Calculates gradients using a second-order central difference. + + 6: WENO3: Third-Order Weighted Essentially Non-Oscillatory (WENO3) scheme. + + 7: WENO5: Fifth-Order Weighted Essentially Non-Oscillatory (WENO5) scheme. + + 8: Hybrid WENO3/ENO1: A hybrid scheme that combines WENO3 in a band surrounding the fire front interface with ENO1 for regions further away. This approach reduces computational cost while maintaining accuracy near the front. + + 9: Hybrid WENO5/ENO1 (default): Similar to option 8, but uses WENO5 instead of WENO3 in the band surrounding the fire front. This approach reduces computational cost while maintaining accuracy near the front. + +``fire_viscosity``: *real* (Default: ``0.4``) + Artificial viscocity in :term:`level-set method` away from the near-front region. + +``fire_lsm_reinit``: *logical* (Default: ``.true.``) + Flag to activate reinitialization of the :term:`level-set method` + +``fire_lsm_reinit_iter``: *integer* (Default: ``1``) + Number of iterations for reinitialization :term:`PDE` + +``fire_upwinding_reinit``: *integer* (Default: ``4``) + Numerical scheme (space) for reinitialization :term:`PDE`. + 1: WENO3 + + 2: WENO5 + + 3: hybrid WENO3-ENO1 + + 4: hybrid WENO5-ENO1 + +``fire_lsm_band_ngp``: *integer* (Default: ``4``) + When using ``fire_upwinding_reinit=3,4`` and ``fire_upwinding=8/9``, the number of grid points around lfn=0 that WENO5/3 is used + +``fire_lsm_zcoupling``: *logical* (Default: ``1``) + When true, uses ``fire_lsm_zcoupling_ref`` instead of ``fire_wind_height`` as a reference height to calculate the logarithmic surface layer wind profile + +``fire_lsm_zcoupling_ref``: *real* (Default: ``50.0``) + [Units: m] + Reference height from which the velocity at ``fire_wind_height`` is calculated using a logarithmic profile + +``fire_viscosity_bg``: *real* (Default: ``0.4``) + Artificial viscosity in the near-front region + +``fire_viscosity_band``: *real* (Default: ``0.5``) + Number of times the hybrid advection band to transition from ``fire_viscosity_bg`` to ``fire_viscosity`` + +``fire_viscosity_ngp``: *integer* (Default: ``2``) + Number of grid points around lfn=0 where ``fire_viscosity_bg`` is used + +``fmoist_run``: *logical* (Default: ``.false.``) + Runs moisture model on the atmospheric grid, outputting the result as a variable named ``fmc_gc`` + +``fmoist_freq``: *integer* (Default: ``0``) + Frequency to run moisture model. + 0: use ``fmoist_dt`` + + k>0: every "k" timesteps + +``fmoist_dt``: *real* (Default: ``600.0``) + [Units: s] + Time step of moisture model (only used if ``fmoist_freq=0``) + +``fire_wind_height``: *integer* (Default: ``6.096``) + [Units: m] + Height of uah,vah wind in fire spread formula + +``fire_is_real_perim``: *logical* (Default: ``.false.``) + Determines if perimeter represents a real fire boundary. + .true. = observed perimeter + + .false. = point/line ignition + +``frac_fburnt_to_smoke``: *real* (Default: ``0.02``) + [Units: g/kg] + Parts per unit of burned fuel converted to smoke, represented as grams of smoke per kilogram of air. + +``fuelmc_g``: *real* (Default: ``0.08``) + Fuel moisture content ground (Dead :term:`FMC`) + +``fuelmc_g_live``: *real* (Default: ``0.30``) + Fuel moisture content ground (Live :term:`FMC`). 30% Completely cured, treat as dead fuel + +``fuelmc_c``: *real* (Default: ``1.00``) + Fuel moisture content of the canopy + +``fuel_opt``: *integer* (Default: ``1``) + Fuel type model. + 1: Anderson fuel model (only option currently implemented) + +``ros_opt``: *integer* (Default: ``0``) + Rate of Spread (ROS) parameterization option. + 0: Rothermel model (only option currently implemented) + +``fmc_opt``: *integer* (Default: ``1``) + :term:`FMC` model + -1 = Constant fuel moisture (only option currently implemented) + +``fire_num_ignitions``: *integer* (Default: ``1``) + Number of ignitions for fire initiation. Maximum of 5. + +.. note:: + For each additional fire ignition, you must specify an additional set of ignition parameters below, with increasing numerical suffixes ( *e.g.* ``fire_ignition_start_lon2``, ``fire_ignition_start_lon3``, etc. ) + +``fire_ignition_start_lon1``: *real* (Default: ``0.0``) + Longitude of first ignition start point. + +``fire_ignition_start_lat1``: *real* (Default: ``0.0``) + Latitude of first ignition start point. + +``fire_ignition_end_lon1``: *real* (Default: ``0.0``) + Longitude of first ignition end point. + +``fire_ignition_end_lat1``: *real* (Default: ``0.0``) + Latitude of first ignition end point. + +``fire_ignition_ros1``: *real* (Default: ``0.01``) + [Units: m/s] + Rate of spread for first ignition (Rothermel parameterization). + +``fire_ignition_start_time1``: *real* (Default: ``0.0``) + [Units: s] + Start time of first ignition in seconds (counting from the beginning of the simulation) + +``fire_ignition_end_time1``: *real* (Default: ``1``) + [Units: s] + End time of first ignition in seconds (counting from the beginning of the simulation) + +``fire_ignition_radius1``: *real* (Default: ``0.0``) + [Units: m] + Radius of the ignition area for first ignition. + + diff --git a/doc/Glossary.rst b/doc/Glossary.rst new file mode 100644 index 0000000..ce867d4 --- /dev/null +++ b/doc/Glossary.rst @@ -0,0 +1,73 @@ +.. _Glossary: + +************************* +Glossary +************************* + +.. glossary:: + + advect + advection + According to the American Meteorological Society (AMS) definition, `advection `_ is "The process of transport of an atmospheric property solely by the mass motion (velocity field) of the atmosphere." In common parlance, advection is movement of atmospheric substances that are carried around by the wind. + + CCPP + The `Common Community Physics Package `_ is a forecast-model agnostic, vetted collection of code containing atmospheric physical parameterizations and suites of parameterizations for use in Numerical Weather Prediction (NWP) along with a framework that connects the physics to the host forecast model. + + ESMF + `Earth System Modeling Framework `__. The ESMF defines itself as “a suite of software tools for developing high-performance, multi-component Earth science modeling applications.” + + FMC + Fuel moisture content + + FV3 + The Finite-Volume Cubed-Sphere :term:`dynamical core` (dycore). Developed at NOAA's `Geophysical + Fluid Dynamics Laboratory `__ (GFDL), it is a scalable and flexible dycore capable of both hydrostatic and non-hydrostatic atmospheric simulations. It is the dycore used in the UFS Weather Model. + + HRRR + `High Resolution Rapid Refresh `__. The HRRR is a NOAA real-time 3-km resolution, hourly updated, cloud-resolving, convection-allowing atmospheric model initialized by 3km grids with 3km radar assimilation. Radar data is assimilated in the HRRR every 15 min over a 1-h period adding further detail to that provided by the hourly data assimilation from the 13km radar-enhanced Rapid Refresh. + + ICs + Initial conditions + + level-set method + A commonly used method in wildland fire modeling used to track and propogate the fire perimiter. More details available in :cite:alp:`LevelSetMethod`. + + LBCs + Lateral boundary conditions + + namelist + A namelist defines a group of variables or arrays. Namelists are an I/O feature for format-free input and output of variables by key-value assignments in Fortran compilers. Fortran variables can be read from and written to plain-text files in a standardised format, usually with a ``.nml`` file ending. + + NCAR + The National Science Foundation's `National Center for Atmospheric Research `__. + + NCEP + National Centers for Environmental Prediction (NCEP) is an arm of the National Weather Service + consisting of nine centers. More information can be found at https://www.ncep.noaa.gov. + + netCDF + NetCDF (`Network Common Data Form `__) is a file format and community standard for storing multidimensional scientific data. It includes a set of software libraries and machine-independent data formats that support the creation, access, and sharing of array-oriented scientific data. + + NUOPC + The `National Unified Operational Prediction Capability `__ Layer "defines conventions and a set of generic components for building coupled models using the Earth System Modeling Framework (:term:`ESMF`)." + + NWS + The `National Weather Service `__ (NWS) is an agency of the United States government that is tasked with providing weather forecasts, warnings of hazardous weather, and other weather-related products to organizations and the public for the purposes of protection, safety, and general information. It is a part of the National Oceanic and Atmospheric Administration (NOAA) branch of the Department of Commerce. + + Parameterizations + Simplified functions that approximate the effects of small-scale processes (e.g., microphysics, gravity wave drag) that cannot be explicitly resolved by a model grid’s representation of the earth. + + PDE + Partial differential equation + + SDF + Suite Definition File. An external file containing information about the construction of a physics suite. It describes the schemes that are called, in which order they are called, whether they are subcycled, and whether they are assembled into groups to be called together. + + tracer + According to the American Meteorological Society (AMS) definition, a `tracer `_ is "Any substance in the atmosphere that can be used to track the history [i.e., movement] of an air mass." Tracers are carried around by the motion of the atmosphere (i.e., by :term:`advection`). These substances are usually gases (e.g., water vapor, CO2), but they can also be non-gaseous (e.g., rain drops in microphysics parameterizations). In weather models, temperature (or potential temperature), absolute humidity, and radioactivity are also usually treated as tracers. According to AMS, "The main requirement for a tracer is that its lifetime be substantially longer than the transport process under study." + + UFS + The Unified Forecast System is a community-based, coupled, comprehensive Earth modeling + system consisting of several applications (apps). These apps span regional to global + domains and sub-hourly to seasonal time scales. The UFS is designed to support the :term:`Weather Enterprise` and to be the source system for NOAA's operational numerical weather prediction applications. For more information, visit https://ufs.epic.noaa.gov/. + diff --git a/doc/Idealized.rst b/doc/Idealized.rst new file mode 100644 index 0000000..4a2f5da --- /dev/null +++ b/doc/Idealized.rst @@ -0,0 +1,13 @@ +.. _Idealized: + +===================== +CFBM: Idealized cases +===================== + +Under construction + +.. image:: https://media.tenor.com/4fu8LKc2vZ4AAAAi/under-construction-wip.gif + :height: 400 + :alt: An animated construction sign + :align: center + diff --git a/doc/Introduction.rst b/doc/Introduction.rst new file mode 100644 index 0000000..6d7fc31 --- /dev/null +++ b/doc/Introduction.rst @@ -0,0 +1,16 @@ +.. _Introduction: + +============== +Introduction +============== + +Welcome to the Community Fire Behavior Model (CFBM) User’s Guide! + +CFBM is a cutting-edge tool designed to simulate the evolution of wildland fires. It leverages the Earth System Modeling Framework (ESMF) and supports both idealized or offline fire spread simulations, as well as coupling with atmospheric models such as WRF or the Univied Forecast System (UFS) atmospheric (FV3) components.The fire model relies on wind and other near-surface variables from the atmospheric component to propagate fires. In turn, it provides feedback to the atmospheric component by delivering heat and moisture fluxes generated by the fire. Additionally, the model releases smoke, which can be advected and diffused as a tracer by the atmospheric component. + +CFBM draws from the proven capabilities of Weather Research and Forecasting Model (WRF-Fire) version 4.3.3, yet continues to expand and evolve to meet the demands of both research and operational communities. This documentation serves as a practical guide for setting up and running simulations using CFBM. + +References +================= + +.. bibliography:: references.bib diff --git a/doc/Makefile b/doc/Makefile new file mode 100644 index 0000000..c91f2f1 --- /dev/null +++ b/doc/Makefile @@ -0,0 +1,27 @@ +# Makefile for Sphinx documentation + +SPHINXOPTS = -a -n #-W +SPHINXBUILD = sphinx-build +SOURCEDIR = . +BUILDDIR = build +LINKCHECKDIR = $(BUILDDIR)/linkcheck + +.PHONY: help Makefile linkcheck + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) + +doc: + make clean + $(MAKE) linkcheck + $(MAKE) html + +linkcheck: + $(SPHINXBUILD) -b linkcheck $(SPHINXOPTS) $(SOURCEDIR) $(LINKCHECKDIR) + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). + +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) -w $(BUILDDIR)/warnings.log \ No newline at end of file diff --git a/doc/SRW.rst b/doc/SRW.rst new file mode 100644 index 0000000..c283cdf --- /dev/null +++ b/doc/SRW.rst @@ -0,0 +1,25 @@ +.. _SRW: + +================================================ +Coupling to the UFS: Running the CFBM in the SRW +================================================ + +This schematic figure illustrates the integration of the Community Fire Behavior Model (CFBM) as a component within the Unified Forecast System (UFS) using the NUOPC (National Unified Operational Prediction Capability) framework. + +The blue box at the top represents the UFS atmosphere component, which provides atmospheric variables (e.g., wind, temperature) to drive the fire model. The orange box represents CFBM that simulates fire spread and calculates heat, moisture, and fire smoke, which feed back into the atmospheric component. The NUOPC cap acts as an interface, facilitating the exchange of data between the UFS atmosphere and the CFBM. The CFBM requires static input in the form of geo_em.d01, which contains refined grids, fire fuels, and detailed topography to accurately simulate fire spread. + +.. !.. image:: https://github.com/NCAR/fire_behavior/blob/develop/doc/CFBM-NUOPC.jpeg +.. ! :width: 400 +.. ! :alt: CFBM-NUOPC coupling schematic +.. ! :align: center + +UFS Short-Range Weather Application (SRW) +========================================= + +The CFBM has been coupled to the UFS Weather Model for both one-way (atmosphere -> fire) and two-way coupled simulations. +Simulations can be run using the UFS Short-Range Weather Application (SRW), a community-supported application for running +numerical weather prediction simulations on limited-area domains. For information on using this capability, see +`the SRW Users Guide `_. + +A preprint on scientific results using this new UFS Fire capability is available (:cite:`CFBM`). + diff --git a/doc/WRFData.rst b/doc/WRFData.rst new file mode 100644 index 0000000..d4faeaf --- /dev/null +++ b/doc/WRFData.rst @@ -0,0 +1,24 @@ +.. _WRF_data: + +================================ +Running real cases with WRF data +================================ + +The CFBM can run in offline mode using WRF atmospheric fields as input. To run CFBM in the offline mode, users need to provide WRF atmospheric data, including all available timestamps, in a file named ```wrf.nc```, along with static inpupts in ```geo_em.d01.nc```. The atmospheric data must include wind components (U, V), geopotential heights (PH, PHB), surface variables used with the fuel moisture model (RAINC, RAINNC, T2, Q2, PSFC), and roughness length (ZNT). + +The simulation configuration is set up in the ```namelist.fire``` file, as described in Section 2.2. Users must specify the number of vertical levels (```kde```) and the time interval for incoming atmospheric data (```interval_atm```) to match the data provided in the ```wrf.nc``` file. + +Users can obtain the model from the github repository: + +```git clone https://github.com/NCAR/fire_behavior.git``` + +To compile the code on derecho, run: + +```./compile.sh --env-auto``` + +If the compilation is successful, the model can be run using ```fire_behavior.exe``` located in the ```build``` directory. + +An example is provided in ```test/test7/``` directory with the ```namelist.fire```. + +If the simulation is successful, the model outputs are written to files named 'fire_output_*', and diagnostic messages are recorded in the 'log' file. + diff --git a/doc/conf.py b/doc/conf.py new file mode 100644 index 0000000..0d2e167 --- /dev/null +++ b/doc/conf.py @@ -0,0 +1,214 @@ +# -*- coding: utf-8 -*- +# +# Configuration file for the Sphinx documentation builder. +# +# This file does only contain a selection of the most common options. For a +# full list see the documentation: +# http://www.sphinx-doc.org/en/master/config + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +import sphinx +from sphinx.util import logging +sys.path.insert(0, os.path.abspath('../physics')) +sys.path.insert(0, os.path.abspath('../tests')) + + +# -- Project information ----------------------------------------------------- + +project = 'Community Fire Behavior Model User\'s Guide' +copyright = '2024, ' +author = ' ' + +# The short X.Y version +version = 'develop' +# The full version, including alpha/beta/rc tags +release = f'version {version}' +html_logo = "https://www.archives.ucar.edu/sites/default/files/images/NSF-NCAR_Lockup-UCAR-Dark_102523%20%282%29.png" + +numfig = True + +nitpick_ignore = [('py:class', 'obj'),('py:class', + 'yaml.dumper.Dumper'),('py:class', + 'xml.etree.ElementTree'),] + +# -- General configuration --------------------------------------------------- + +# Sphinx extension module names: +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.doctest', + 'sphinx.ext.intersphinx', + 'sphinx.ext.extlinks', + 'sphinx.ext.mathjax', + 'sphinxcontrib.bibtex', +] + +bibtex_bibfiles = ['references.bib'] +bibtex_reference_style = 'author_year' + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# +# source_suffix = ['.rst', '.md'] +source_suffix = '.rst' + +# The master toctree document. +master_doc = 'index' + +# The language for content autogenerated by Sphinx. +# Not set because default is 'en'. +# language = 'en' + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ['_build', + '.DS_Store',] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# Documentation-wide substitutions + +rst_prolog = """ +.. |wflow_env| replace:: ``srw_app`` +.. |graphics_env| replace:: ``regional_workflow`` +.. |cmaq_env| replace:: ``regional_workflow_cmaq`` +.. |activate| replace:: ``conda activate srw_app`` +.. |prompt| replace:: ``(srw_app)`` +.. |latestr| replace:: v2.2.0 +.. |branch| replace:: ``develop`` +.. |data| replace:: develop +""" + +# Linkcheck options + +# Avoid a 403 Forbidden error when accessing certain links (e.g., noaa.gov) +user_agent = "Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/121.0.0.0 Safari/537.36" + +# Ignore working links that cause a linkcheck 403 error. +#linkcheck_ignore = [r'https://www\.intel\.com/content/www/us/en/docs/cpp\-compiler/developer\-guide\-reference/2021\-10/thread\-affinity\-interface\.html', +# r'https://www\.intel\.com/content/www/us/en/developer/tools/oneapi/hpc\-toolkit\-download\.html', +# r'https://glossary.ametsoc.org/.*', +# ] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'classic' +html_logo= "https://www.archives.ucar.edu/sites/default/files/images/NSF-NCAR_Lockup-UCAR-Dark_102523%20%282%29.png" + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +# html_theme_options = {} +html_theme_options = { + "body_max_width": "none", + "navigation_depth": 8, + } + + +# -- Options for HTMLHelp output --------------------------------------------- + +# Output file base name for HTML help builder. +htmlhelp_basename = 'UFS-SRWeather-App' + + +# -- Options for LaTeX output ------------------------------------------------ + +latex_engine = 'pdflatex' +latex_elements = { + # The paper size ('letterpaper' or 'a4paper'). + 'papersize': 'letterpaper', + + # The font size ('10pt', '11pt' or '12pt'). + 'pointsize': '11pt', + + # Additional stuff for the LaTeX preamble. + 'preamble': r''' + \usepackage{charter} + \usepackage[defaultsans]{lato} + \usepackage{inconsolata} + ''', + # Release name prefix + 'releasename': ' ', + + # Latex figure (float) alignment + # + # 'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'CFBM-UsersGuide.tex', 'Community Fire Behavior Module Documentation', + ' ', 'manual'), +] + + + +# -- Options for Epub output ------------------------------------------------- + +# Bibliographic Dublin Core info. +epub_title = project + +# The unique identifier of the text. This can be a ISBN number +# or the project homepage. +# +# epub_identifier = '' + +# A unique identification for the text. +# +# epub_uid = '' + +# A list of files that should not be packed into the epub file. +epub_exclude_files = ['search.html'] + + +# -- Extension configuration ------------------------------------------------- + + +# -- Options for intersphinx extension --------------------------------------- + +# Example configuration for intersphinx: refer to the Python standard library. +intersphinx_mapping = { + 'python': ('https://docs.python.org/3', None), + 'hpc-stack': ('https://hpc-stack-epic.readthedocs.io/en/develop/', None), + 'spack-stack': ('https://spack-stack.readthedocs.io/en/develop/', None), + 'met': ('https://met.readthedocs.io/en/develop/', None), + 'metplus': ('https://metplus.readthedocs.io/en/develop/', None), + 'ufs-wm': ('https://ufs-weather-model.readthedocs.io/en/develop/', None), + 'upp': ('https://upp.readthedocs.io/en/develop/', None), + 'ufs-utils': ('https://noaa-emcufs-utils.readthedocs.io/en/latest/', None), + 'ccpp-techdoc': ('https://ccpp-techdoc.readthedocs.io/en/ufs_srw_app_v2.2.0/', None), + 'stochphys': ('https://stochastic-physics.readthedocs.io/en/latest/', None), + 'srw_v2.2.0': ('https://ufs-srweather-app.readthedocs.io/en/release-public-v2.2.0/', None), +} + +# -- Options for extlinks extension --------------------------------------- + +extlinks_detect_hardcoded_links = True +extlinks = {'github-docs': ('https://docs.github.com/en/%s', '%s'), + 'nco': ('https://www.nco.ncep.noaa.gov/idsb/implementation_standards/%s', '%s'), + "rst": ("https://www.sphinx-doc.org/en/master/usage/restructuredtext/%s", "%s"), + "rtd": ("https://readthedocs.org/projects/ufs-srweather-app/%s", "%s"), + 'srw-repo': ('https://github.com/ufs-community/ufs-srweather-app/%s', '%s'), + 'srw-wiki': ('https://github.com/ufs-community/ufs-srweather-app/wiki/%s','%s'), + 'uw': ('https://uwtools.readthedocs.io/en/main/%s', '%s'), + } + diff --git a/doc/index.rst b/doc/index.rst new file mode 100644 index 0000000..3aad61a --- /dev/null +++ b/doc/index.rst @@ -0,0 +1,13 @@ +Community Fire Behavior Module (CFBM) Documentation (|version|) +=============================================================== + +.. toctree:: + :numbered: + :maxdepth: 3 + + Introduction + Configuration + SRW + WRFData + Idealized + Glossary diff --git a/doc/make.bat b/doc/make.bat new file mode 100644 index 0000000..295c7e3 --- /dev/null +++ b/doc/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=doc +set BUILDDIR=build + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% + +:end +popd diff --git a/doc/references.bib b/doc/references.bib new file mode 100644 index 0000000..ce4fd0c --- /dev/null +++ b/doc/references.bib @@ -0,0 +1,24 @@ +@article{CFBM, + title={The Community Fire Behavior Model for coupled fire-atmosphere modeling: Implementation in the Unified Forecast System}, + author={Pedro Angel {Jimenez y Munoz} and Maria Frediani and Masih Eghdami and Daniel Rosen and Michael Kavulich and and Timothy W. Juliano}, + journal={Geosci. Model Dev.}, + year={2024}, + doi={10.5194/gmd-2024-124}, + } + +@article{LevelSetMethod, + title={An accurate fire-spread algorithm in the Weather Research and Forecasting model using the level-set method}, + author={Domingo Muñoz-Esparza and Branko Kosović and Pedro A. Jiménez and Janice L. Coen}, + journal={J. Adv. Model. Earth Syst.}, + year={2018}, + doi={10.1002/2017MS001108}, + } + +@article{SethianMethod, + title={A fast marching level set method for monotonically advancing fronts}, + author={J A Sethian}, + journal={P. Natl. Acad. Sci. USA}, + year={1996}, + doi={10.1073/pnas.93.4.159}, + } + diff --git a/doc/requirements.in b/doc/requirements.in new file mode 100644 index 0000000..26c778f --- /dev/null +++ b/doc/requirements.in @@ -0,0 +1,3 @@ +sphinx>=6.0.0 +sphinx_rtd_theme +sphinxcontrib-bibtex diff --git a/doc/requirements.txt b/doc/requirements.txt new file mode 100644 index 0000000..15617eb --- /dev/null +++ b/doc/requirements.txt @@ -0,0 +1,3 @@ +sphinx-rtd-theme==2.0.0 +sphinxcontrib-bibtex==2.6.2 + diff --git a/driver/fire_behavior.F90 b/driver/fire_behavior.F90 index 22b8a9e..78480e0 100644 --- a/driver/fire_behavior.F90 +++ b/driver/fire_behavior.F90 @@ -17,6 +17,7 @@ program fire_behavior type (wrf_t) :: atm_state type (namelist_t) :: config_flags + #ifdef DM_PARALLEL call mpi_init (ierr) if (ierr /= MPI_SUCCESS) then @@ -28,14 +29,25 @@ program fire_behavior ! Read namelist call config_flags%Initialization (file_name = 'namelist.fire') - call Init_atm_state (atm_state, config_flags) - call Init_fire_state (grid, config_flags, atm_state) + select case (config_flags%ideal_opt) + case (0) + call Init_atm_state (atm_state, config_flags) + call Init_fire_state (grid, config_flags, atm_state) + + case (1) + call Init_fire_state (grid, config_flags) + + case default + write (ERROR_UNIT, *) 'ERROR: ideal_opt option not supported: ', config_flags%ideal_opt + stop + + end select call grid%Save_state () do while (grid%datetime_now < grid%datetime_end) call Advance_state (grid, config_flags) call grid%Handle_output (config_flags) - call grid%Handle_wrfdata_update (atm_state, config_flags) + if (config_flags%ideal_opt == 0) call grid%Handle_wrfdata_update (atm_state, config_flags) end do #ifdef DM_PARALLEL diff --git a/driver/initialize_mod.F90 b/driver/initialize_mod.F90 index 0c500a1..d14d270 100644 --- a/driver/initialize_mod.F90 +++ b/driver/initialize_mod.F90 @@ -9,7 +9,7 @@ module initialize_mod private - public :: Init_fire_state, Init_atm_state + public :: Init_fire_state, Init_atm_state, Init_fire_state_within_wrf contains @@ -46,20 +46,25 @@ subroutine Init_fire_state (grid, config_flags, wrf) if (DEBUG_LOCAL) call Print_message (' Entering subroutine Init_state') - ! Fire state initialization - if (DEBUG_LOCAL) call Print_message (' Reading geogrid file') - geogrid = geogrid_t (file_name = 'geo_em.d01.nc') + if (DEBUG_LOCAL) call Print_message (' Initialization...') + if (config_flags%ideal_opt == 0) then + ! Real world + if (DEBUG_LOCAL) call Print_message (' Reading geogrid file') + geogrid = geogrid_t (file_name = 'geo_em.d01.nc') - if (DEBUG_LOCAL) call Print_message (' Initializing state') - call grid%Initialization (config_flags, geogrid) + if (DEBUG_LOCAL) call Print_message (' Initializing fire state') + call grid%Initialization (config_flags, geogrid) - ! Atmosphere to Fire - if (present (wrf)) then - if (DEBUG_LOCAL) call Print_message (' Initializing atmospheric state') - call grid%Handle_wrfdata_update (wrf, config_flags) + if (present (wrf)) then + if (DEBUG_LOCAL) call Print_message (' Initializing atmospheric state') + call grid%Handle_wrfdata_update (wrf, config_flags) + end if + else + ! Ideal + if (DEBUG_LOCAL) call Print_message (' Initializing fire state') + call grid%Initialization (config_flags) end if - ! Fire init call Init_fire_components (grid, config_flags) if (DEBUG_LOCAL) then @@ -90,4 +95,37 @@ subroutine Init_fire_state (grid, config_flags, wrf) end subroutine Init_fire_state + subroutine Init_fire_state_within_wrf (state, config_flags, & + ifds, ifde, ifms, ifme, ifps, ifpe, & + jfds, jfde, jfms, jfme, jfps, jfpe, & + kfds, kfde, kfms, kfme, kfps, kfpe, & + kfts, kfte, ide, jde, dx, dy, sr_x, sr_y, & + map_proj, cen_lat, cen_lon, truelat1, truelat2, stand_lon, & + nfuel_cat, zsf, dzdxf, dzdyf) + + implicit none + + type (state_fire_t), intent (in out) :: state + type (namelist_t), intent (in) :: config_flags + integer, intent (in) :: ifds, ifde, ifms, ifme, ifps, ifpe, & + jfds, jfde, jfms, jfme, jfps, jfpe, & + kfds, kfde, kfms, kfme, kfps, kfpe, & + kfts, kfte, map_proj, sr_x, sr_y, ide, jde + real :: dx, dy, cen_lat, cen_lon, truelat1, truelat2, stand_lon + real, dimension(ifms:ifme, jfms:jfme), intent (in) :: nfuel_cat, zsf, dzdxf, dzdyf + + + call state%Initialization (config_flags, & + ifds = ifds, ifde = ifde, ifms = ifms, ifme = ifme, ifps = ifps, ifpe = ifpe, & + jfds = jfds, jfde = jfde, jfms = jfms, jfme = jfme, jfps = jfps, jfpe = jfpe, & + kfds = kfds, kfde = kfde, kfms = kfms, kfme = kfme, kfps = kfps, kfpe = kfpe, & + kfts = kfts, kfte = kfte, ide = ide, jde = jde, & + cen_lat = cen_lat, cen_lon = cen_lon, truelat1 = truelat1, & + truelat2 = truelat2, stand_lon = stand_lon, dx = dx, dy = dy, sr_x = sr_x, sr_y = sr_y, & + nfuel_cat = nfuel_cat, zsf = zsf, dzdxf = dzdxf, dzdyf = dzdyf) + + call Init_fire_components (state, config_flags) + + end subroutine Init_fire_state_within_wrf + end module initialize_mod diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index cc1cdda..04a32f0 100644 --- a/io/namelist_mod.F90 +++ b/io/namelist_mod.F90 @@ -9,6 +9,12 @@ module namelist_mod public :: namelist_t, FIRE_MAX_IGNITIONS_IN_NAMELIST integer, parameter :: FIRE_MAX_IGNITIONS_IN_NAMELIST = 5 + ! Default ideal block + real, parameter :: DX_DEFAULT = 100.0, U_DEFAULT = 5.0, V_DEFAULT = 0.0, LAT_DEFAULT = 40.3636, LON_DEFAULT = -4.4035, & + DZ_DX_DEFAULT = 0.0, ELEVATION_DEFAULT = 0.0 + integer, parameter :: IDEAL_OPT_DEFAULT = 0, NX_DEFAULT = 100, FUEL_CAT_DEFAULT = 1 + ! Default options + integer, parameter :: NO_FMC_MODEL = -1 type :: namelist_t integer :: start_year = -1, start_month = -1, start_day = -1, start_hour = -1, start_minute = -1, start_second = -1, & @@ -17,6 +23,7 @@ module namelist_mod real :: dt = 2.0 integer :: num_tiles = 1 + integer :: tile_strategy = 0 integer :: fire_print_msg = 0 ! "write fire statistics, 0 no writes, 1+ for more" "" real :: fire_atm_feedback = 1.0 ! "the heat fluxes to the atmosphere are multiplied by this" "1" @@ -36,6 +43,7 @@ module namelist_mod ! for fire_upwinding_reinit=4,5 and fire_upwinding=8,9 options" real :: fire_wind_height = 6.096 ! "height of uah,vah wind in fire spread formula" "m" + integer :: wind_vinterp_opt = 0 ! "wind (adjustment factor) interpolation option" logical :: fire_lsm_zcoupling = .false. ! "flag to activate reference velocity at a different height from fire_wind_height" real :: fire_lsm_zcoupling_ref = 50.0 ! "reference height from wich u at fire_wind_hegiht is calculated using a logarithmic profile" "m" @@ -48,10 +56,13 @@ module namelist_mod integer :: fmoist_freq = 0 ! frequency to run moisture model 0: use fmoist_dt, k>0: every k timesteps real :: fmoist_dt = 600.0 ! moisture model time step [s] + integer :: ideal_opt = IDEAL_OPT_DEFAULT ! 0) real world, 1) ideal + ! Objects integer :: fuel_opt = 1 ! Fuel model integer :: ros_opt = 0 ! ROS parameterization - integer :: fmc_opt = -1 ! FMC model + integer :: fmc_opt = NO_FMC_MODEL ! FMC model + integer :: emis_opt = 0 ! Object to be added. 0) WRF-Fire emiss, 1) PM2.5 as a function of FMC ! Ignitions integer :: fire_num_ignitions = 0 ! "number of ignition lines" @@ -101,17 +112,49 @@ module namelist_mod real :: fire_ignition_end_time5 = 0.0 real :: fire_ignition_radius5 = 0.0 + ! Ideal block + real :: dx = DX_DEFAULT + real :: dy = DX_DEFAULT + integer :: nx = NX_DEFAULT + integer :: ny = NX_DEFAULT + + real :: zonal_wind = U_DEFAULT + real :: meridional_wind = V_DEFAULT + integer :: fuel_cat = FUEL_CAT_DEFAULT + real :: dz_dx = DZ_DX_DEFAULT + real :: dz_dy = DZ_DX_DEFAULT + real :: elevation = ELEVATION_DEFAULT + + real :: cen_lat = LAT_DEFAULT + real :: cen_lon = LON_DEFAULT + real :: stand_lon = LON_DEFAULT + real :: true_lat_1 = LAT_DEFAULT + real :: true_lat_2 = LAT_DEFAULT + ! Atmosphere integer :: kds = 1, kde = 1 contains + procedure, public :: Check_nml => Check_nml procedure, public :: Initialization => Init_namelist procedure, public :: Init_fire_block => Init_fire_block + procedure, public :: Init_ideal_block => Init_ideal_block procedure, public :: Init_time_block => Init_time_block procedure, public :: Init_atm_block => Init_atm_block_legacy end type namelist_t contains + subroutine Check_nml (this) + + implicit none + + class (namelist_t), intent (in out) :: this + + if (this%ideal_opt /= 0 .and. this%fmoist_run) & + call Stop_simulation ('ideal runs do not support a FMC model') + + end subroutine Check_nml + subroutine Init_atm_block_legacy (this, file_name) implicit none @@ -172,16 +215,20 @@ subroutine Init_fire_block (this, file_name) integer :: fmoist_freq = 0 ! "frequency to run moisture model 0: use fmoist_dt, k>0: every k timesteps" "1" real :: fmoist_dt = 600 ! "moisture model time step" "s" real :: fire_wind_height = 6.096 ! "height of uah,vah wind in fire spread formula" "m" + integer :: wind_vinterp_opt = 0 ! "wind (adjustment factor) interpolation option" logical :: fire_is_real_perim = .false. ! .false. = point/line ignition, .true. = observed perimeter" real :: frac_fburnt_to_smoke = 0.02 ! "parts per unit of burned fuel becoming smoke " "g_smoke/kg_air" real :: fuelmc_g = 0.08 ! Fuel moisture content ground (Dead FMC) real :: fuelmc_g_live = 0.30 ! Fuel moisture content ground (Live FMC). 30% Completely cured, treat as dead fuel real :: fuelmc_c = 1.00 ! Fuel moisture content canopy + integer :: ideal_opt = IDEAL_OPT_DEFAULT + ! Objects integer :: fuel_opt = 1 ! Fuel model integer :: ros_opt = 0 ! ROS parameterization integer :: fmc_opt = -1 ! FMC model + integer :: emis_opt = 0 ! Smoke emissions ! ignitions integer :: fire_num_ignitions = 0 @@ -218,8 +265,10 @@ subroutine Init_fire_block (this, file_name) fmoist_freq, fmoist_dt, & fire_wind_height, fire_is_real_perim, frac_fburnt_to_smoke, fuelmc_g, & fuelmc_g_live, fuelmc_c, & + ideal_opt, & ! objects - fuel_opt, ros_opt, fmc_opt, & + fuel_opt, ros_opt, fmc_opt, emis_opt, & + wind_vinterp_opt, & ! Ignitions fire_num_ignitions, & ! Ignition 1 @@ -276,9 +325,13 @@ subroutine Init_fire_block (this, file_name) this%fuelmc_g_live = fuelmc_g_live this%fuelmc_c = fuelmc_c + this%ideal_opt = ideal_opt + this%fuel_opt = fuel_opt this%ros_opt = ros_opt this%fmc_opt = fmc_opt + this%emis_opt = emis_opt + this%wind_vinterp_opt = wind_vinterp_opt this%fire_num_ignitions = fire_num_ignitions @@ -329,6 +382,74 @@ subroutine Init_fire_block (this, file_name) end subroutine Init_fire_block + subroutine Init_ideal_block (this, file_name) + + implicit none + + class (namelist_t), intent (in out) :: this + character (len = *), intent (in) :: file_name + + real :: dx, dy, zonal_wind, meridional_wind, cen_lat, cen_lon, stand_lon, true_lat_1, true_lat_2, & + dz_dx, dz_dy, elevation + integer :: nx, ny, fuel_cat + + character (len = :), allocatable :: msg + integer :: unit_nml, io_stat + + namelist /ideal/ dx, dy, nx, ny, zonal_wind, meridional_wind, fuel_cat, & + dz_dx, dz_dy, elevation, cen_lat, cen_lon, stand_lon, true_lat_1, true_lat_2 + + + ! Set default values + dx = DX_DEFAULT + dy = DX_DEFAULT + nx = NX_DEFAULT + ny = NX_DEFAULT + + zonal_wind = U_DEFAULT + meridional_wind = V_DEFAULT + fuel_cat = FUEL_CAT_DEFAULT + dz_dx = DZ_DX_DEFAULT + dz_dy = DZ_DX_DEFAULT + elevation = ELEVATION_DEFAULT + + cen_lat = LAT_DEFAULT + cen_lon = LON_DEFAULT + stand_lon = LON_DEFAULT + true_lat_1 = LAT_DEFAULT + true_lat_2 = LAT_DEFAULT + + open (newunit = unit_nml, file = trim (file_name), action = 'read', iostat = io_stat) + if (io_stat /= 0) then + msg = 'Problems opening namelist file ' // trim (file_name) + call Stop_simulation (msg) + end if + + read (unit_nml, nml = ideal, iostat = io_stat) + if (io_stat /= 0) call Stop_simulation ('Problems reading namelist ideal block') + close (unit_nml) + + this%dx = dx + this%dy = dy + + this%nx = nx + this%ny = ny + + this%zonal_wind = zonal_wind + this%meridional_wind = meridional_wind + this%fuel_cat = fuel_cat + this%dz_dx = dz_dx + this%dz_dy = dz_dy + this%elevation = elevation + + this%cen_lat = cen_lat + this%cen_lon = cen_lon + this%stand_lon = stand_lon + this%true_lat_1 = true_lat_1 + this%true_lat_2 = true_lat_2 + + end subroutine Init_ideal_block + subroutine Init_time_block (this, file_name) implicit none @@ -338,7 +459,7 @@ subroutine Init_time_block (this, file_name) integer :: start_year, start_month, start_day, start_hour, start_minute, start_second, & end_year, end_month, end_day, end_hour, end_minute, end_second, interval_output, & - num_tiles + num_tiles, tile_strategy real :: dt character (len = :), allocatable :: msg @@ -364,6 +485,7 @@ subroutine Init_time_block (this, file_name) dt = 2.0 interval_output = 0 num_tiles = 1 + tile_strategy = 0 open (newunit = unit_nml, file = trim (file_name), action = 'read', iostat = io_stat) if (io_stat /= 0) then @@ -391,6 +513,7 @@ subroutine Init_time_block (this, file_name) this%interval_output = interval_output this%num_tiles = num_tiles + this%tile_strategy = tile_strategy end subroutine Init_time_block @@ -409,6 +532,9 @@ subroutine Init_namelist (this, file_name) call this%Init_time_block (file_name = trim (file_name)) call this%Init_fire_block (file_name = trim (file_name)) call this%Init_atm_block (file_name = trim (file_name)) + if (this%ideal_opt > 0) call this%Init_ideal_block (file_name = trim (file_name)) + + call this%Check_nml () if (DEBUG_LOCAL) call Print_message (' Leaving subroutine Read_namelist') diff --git a/io/wrf_mod.F90 b/io/wrf_mod.F90 index a246e26..757db4e 100644 --- a/io/wrf_mod.F90 +++ b/io/wrf_mod.F90 @@ -5,13 +5,14 @@ module wrf_mod use namelist_mod, only : namelist_t use netcdf_mod, only : Get_netcdf_var, Get_netcdf_att, Get_netcdf_dim, Is_netcdf_file_present use proj_lc_mod, only : proj_lc_t - use stderrout_mod, only : Print_message + use stderrout_mod, only : Print_message, Stop_simulation + use interp_mod, only : Interp_profile implicit none private - public :: wrf_t, G, RERADIUS + public :: wrf_t, G, RERADIUS, Interp_wrf2dvar_to_cfbm, Interp_wrfwinds_to_cfbm, Provide_atm_feedback real, parameter :: G = 9.81 ! acceleration due to gravity [m s-2] real, parameter :: RERADIUS = 1.0 / 6370.0e03 ! reciprocal of earth radius (m^-1) @@ -57,6 +58,82 @@ module wrf_mod contains + pure function Calc_rh (p, t, qv) result (rh) + + implicit none + + real, intent(in) :: p, t, qv + real :: rh + + real, parameter :: PQ0 = 379.90516, A2 = 17.2693882, A3= 273.16, A4 = 35.86, RHMIN = 1.0, RHMAX = 100.0 + real :: q, qs + integer :: i, j, k + + + q = qv / (1.0 + qv) + qs = PQ0 / p * exp (A2 * (t - A3) / (t - A4)) + rh = 100.0 * q / qs + rh = max (min (rh, RHMAX), RHMIN) + + end function Calc_rh + + subroutine Calc_smoke_aod (dz8w, p_phy, t_phy, qv, rho, smoke_tracer, aod5502d_smoke, & + ids, ide, kds, kde, jds, jde, & + ims, ime, kms, kme, jms, jme, & + its, ite, kts, kte, jts, jte) + + implicit none + + integer, intent (in) :: ids, ide, kds, kde, jds, jde, & + ims, ime, kms, kme, jms, jme, & + its, ite, kts, kte, jts, jte + real, dimension(ims:ime, kms:kme, jms:jme), intent (in) :: dz8w, p_phy, t_phy, qv, rho, smoke_tracer + real, dimension(ims:ime, jms:jme), intent (out) :: aod5502d_smoke + + ! 4.0 abs + 0.5 scattering [m2 g-1] + real, parameter :: MASS_EXT_COEF = 4.5, RH_CRIT = 0.3, RH_MAX = 0.95, CONVERT_PERCENT_TO_UNITLESS = 0.01 + real :: rh, augm_ext_coef + integer :: i, k, j + logical, parameter :: DEBUG_LOCAL = .false. + + + Loop_j_aod : do j = jts, min (jte, jde - 1) + Loop_i_aod : do i = its, min (ite, ide - 1) + aod5502d_smoke(i, j) = 0.0 + Loop_k_aod : do k = kts, min (kte, kde - 1) + rh = CONVERT_PERCENT_TO_UNITLESS * Calc_rh (p_phy(i, k, j), t_phy(i, k, j), qv(i, k, j)) + rh = min (rh, RH_MAX) + if (rh > RH_CRIT) then + augm_ext_coef = MASS_EXT_COEF * ((1.0 - RH_CRIT) / (1.0 - rh)) ** 0.18 + else + augm_ext_coef = MASS_EXT_COEF + end if ! [m2 g-1] [g smoke kg-1 air] [kg air m-3] [m] + aod5502d_smoke(i, j) = aod5502d_smoke(i, j) + augm_ext_coef * smoke_tracer(i, k, j) * rho(i, k, j) * dz8w(i, k, j) + if (DEBUG_LOCAL) call Print_profile (i, j) + end do Loop_k_aod + end do Loop_i_aod + end do Loop_j_aod + + contains + + subroutine Print_profile (i, j) + + use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT + + implicit none + + integer, intent (in) :: i, j + integer, parameter :: I_TO_PRINT = 31, J_TO_PRINT = 34 + + + if (i == I_TO_PRINT .and. j == J_TO_PRINT) then + write (OUTPUT_UNIT, *) k, dz8w(i, k, j), rh, augm_ext_coef, rho(i, k, j), smoke_tracer(i, k, j), aod5502d_smoke(i, j) + end if + + end subroutine Print_profile + + end subroutine Calc_smoke_aod + subroutine Destroy_geopotential_levels (this) implicit none @@ -643,6 +720,69 @@ subroutine Interp_var2grid_nearest (this, lats_in, lons_in, var_name, vals_out) end subroutine Interp_var2grid_nearest + subroutine Interp_wrf2dvar_to_cfbm (wrfatm2dvar, ims, ime, jms, jme, ifms, ifme, jfms, jfme, ifps, ifpe, jfps, jfpe, & + lats_in, lons_in, proj, vals_out) + + implicit none + + integer, intent (in) :: ims, ime, jms, jme, & + ifms, ifme, jfms, jfme, & + ifps, ifpe, jfps, jfpe + real, dimension(ims:ime, jms:jme), intent (in) :: wrfatm2dvar + real, dimension(ifms:ifme, jfms:jfme), intent (in) :: lats_in, lons_in + type (proj_lc_t), intent (in) :: proj + real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: vals_out + + integer :: i, j, i_wrf, j_wrf + real :: i_real, j_real + + + do j = jfps, jfpe + do i = ifps, ifpe + call proj%Calc_ij (lats_in(i, j), lons_in(i, j), i_real, j_real) + i_wrf = min (max (ims, nint (i_real)), ime) + j_wrf = min (max (jms, nint (j_real)), jme) + vals_out(i, j) = wrfatm2dvar(i_wrf, j_wrf) + end do + end do + + end subroutine Interp_wrf2dvar_to_cfbm + + subroutine Interp_wrfwinds_to_cfbm (u_phy, v_phy, z_at_w, ims, ime, kms, kme, jms, jme, ifms, ifme, jfms, jfme, ifps, ifpe, jfps, jfpe, & + kfds, kfde, lats_in, lons_in, proj, z0f, fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, u_out, v_out) + + implicit none + + integer, intent (in) :: ims, ime, kms, kme, jms, jme, & + ifms, ifme, jfms, jfme, & + ifps, ifpe, jfps, jfpe, & + kfds, kfde + logical, intent (in) :: fire_lsm_zcoupling + real, dimension(ims:ime, kms:kme, jms:jme), intent (in) :: u_phy, v_phy, z_at_w + real, dimension(ifms:ifme, jfms:jfme), intent (in) :: lats_in, lons_in, z0f + type (proj_lc_t), intent (in) :: proj + real, intent (in) :: fire_wind_height, fire_lsm_zcoupling_ref + real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: u_out, v_out + + integer :: i, j, i_wrf, j_wrf + real :: i_real, j_real, uout, vout + + + do j = jfps, jfpe + do i = ifps, ifpe + call proj%Calc_ij (lats_in(i, j), lons_in(i, j), i_real, j_real) + i_wrf = min (max (ims, nint (i_real)), ime) + j_wrf = min (max (jms, nint (j_real)), jme) + call Interp_profile (fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, kfds, kfde, & + u_phy(i_wrf, :, j_wrf), v_phy(i_wrf, :, j_wrf), z_at_w(i_wrf, :, j_wrf), z0f(i, j), & + uout, vout) + u_out(i, j) = uout + v_out(i, j) = vout + end do + end do + + end subroutine Interp_wrfwinds_to_cfbm + subroutine Print_domain (this) use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT @@ -667,6 +807,281 @@ subroutine Print_domain (this) end subroutine Print_domain + subroutine Provide_atm_feedback (config_flags, & + ifms, ifme, jfms, jfme, & + ifts, ifte, jfts, jfte, & + ifps, ifpe, jfps, jfpe, & + ids, ide, kds, kde, jds, jde, & + ims, ime, kms, kme, jms, jme, & + its, ite, kts, kte, jts, jte, & + sr_x, sr_y, & + emis_smoke, smoke_tracer, tracer_opt, & + p_phy, t_phy, qv, & + aod5502d_smoke, & + fgrnhfx, fgrnqfx, & + grnhfx, grnqfx, canhfx, canqfx, & + grnsmk, & + alfg, alfc, z1can, & + rho, dz8w, z_at_w, & + mu, c1h, c2h, & + rthfrten, rqvfrten) + + implicit none + + type (namelist_t), intent (in) :: config_flags + integer, intent (in) :: ifms, ifme, jfms, jfme, & + ifps, ifpe, jfps, jfpe, & + ids, ide, kds, kde, jds, jde, & + ims, ime, kms, kme, jms, jme, & + its, ite, kts, kte, jts, jte, & + ifts, ifte, jfts, jfte, sr_x, sr_y + + real, dimension(ifms:ifme, jfms:jfme), intent (in) :: emis_smoke + real, dimension(ims:ime, kms:kme, jms:jme), intent (in out), optional :: smoke_tracer + integer, intent (in) :: tracer_opt + real, dimension(ifms:ifme, jfms:jfme), intent (in) :: fgrnhfx, fgrnqfx + real, dimension(ims:ime, kms:kme, jms:jme), intent (in) :: rho, dz8w, z_at_w, p_phy, t_phy, qv + real, intent(in), dimension(ims:ime, jms:jme) :: mu ! dry air mass (pa) + real, intent(in), dimension(kms:kme) :: c1h, c2h ! hybrid coordinate weights + real, intent(in) :: alfg ! extinction depth surface fire heat (m) + real, intent(in) :: alfc ! extinction depth crown fire heat (m) + real, intent(in) :: z1can ! height of crown fire heat release (m) + real, dimension(ims:ime, jms:jme), intent (out) :: grnhfx, grnqfx, canhfx, canqfx, grnsmk, aod5502d_smoke + real, intent(out), dimension(ims:ime, kms:kme, jms:jme) :: & + rthfrten, & ! theta tendency from fire (in mass units) + rqvfrten ! Qv tendency from fire (in mass units) + + logical, parameter :: DEBUG_LOCAL = .true. + integer :: i, j, ibase, jbase, i_f, j_f, ioff, joff + real :: avgw, convert_kg_m2_to_g_kg + + + if (DEBUG_LOCAL) call Check_dims (its, ite, jts, jte, ifts, ifte, jfts, jfte, sr_x, sr_y) + + avgw = 1.0 / (sr_x * sr_y) + do j = max (jds + 1, jts), min (jte, jde - 2) + jbase = jfts + sr_y * (j - jts) + do i = max (ids + 1, its), min (ite, ide - 2) + ibase = ifts + sr_x * (i - its) + canqfx(i, j) = 0.0 + canhfx(i, j) = 0.0 + grnsmk(i, j) = 0.0 + grnhfx(i, j) = 0.0 + grnqfx(i, j) = 0.0 + convert_kg_m2_to_g_kg = 1000.0 / (rho(i, kts, j) * dz8w(i, kts, j)) + do joff = 0, sr_y - 1 + j_f = joff + jbase + do ioff = 0, sr_x - 1 + i_f = ioff + ibase + grnsmk(i, j) = grnsmk(i, j) + emis_smoke(i_f, j_f) + grnhfx(i, j) = grnhfx(i, j) + fgrnhfx(i_f, j_f) ! * config_flags%fire_atm_feedback + grnqfx(i, j) = grnqfx(i, j) + fgrnqfx(i_f, j_f) ! * config_flags%fire_atm_feedback + end do + end do + grnhfx(i, j) = grnhfx(i, j) * avgw + grnqfx(i, j) = grnqfx(i, j) * avgw + grnsmk(i, j) = grnsmk(i, j) * convert_kg_m2_to_g_kg * avgw + if (tracer_opt == 3) smoke_tracer(i, kts, j) = smoke_tracer(i, kts, j) + grnsmk(i, j) + end do + end do + + call Fire_tendency ( & + ids,ide - 1,kds,kde,jds,jde - 1, & ! dimensions + ims,ime,kms,kme,jms,jme, & + its,min (ite, ide-1),kts,kte,jts,min (jte, jde - 1), & + grnhfx,grnqfx,canhfx,canqfx, & ! heat fluxes summed up to atm grid + alfg,alfc,z1can, & ! coeffients, properties, geometry + z_at_w,dz8w,mu,c1h,c2h,rho, & + config_flags%fire_atm_feedback, & + rthfrten,rqvfrten) ! theta and Qv tendencies + + if (tracer_opt == 3) call Calc_smoke_aod (dz8w, p_phy, t_phy, qv, rho, smoke_tracer, aod5502d_smoke, & + ids, ide, kds, kde, jds, jde, & + ims, ime, kms, kme, jms, jme, & + its, ite, kts, kte, jts, jte) + + contains + + subroutine Check_dims (its, ite, jts, jte, ifts, ifte, jfts, jfte, sr_x, sr_y) + + use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT + + implicit none + + integer, intent (in) :: its, ite, jts, jte, ifts, ifte, jfts, jfte, sr_x, sr_y + + integer :: isz1, jsz1, isz2, jsz2, ir, jr + logical, parameter :: DEBUG_LOCAL = .false. + + + isz1 = ite - its + 1 + jsz1 = jte - jts + 1 + isz2 = ifte - ifts + 1 + jsz2 = jfte - jfts + 1 + ir = isz2 / isz1 + jr = jsz2 / jsz1 + + if (DEBUG_LOCAL) write (OUTPUT_UNIT, *) 'its, ite, jts, jte =', its, ite, jts, jte + if (DEBUG_LOCAL) write (OUTPUT_UNIT, *) 'ifts, ifte, jfts, jfte =', ifts, ifte, jfts, jfte + if (DEBUG_LOCAL) write (OUTPUT_UNIT, *) 'isz1, jsz1, isz2, jsz2 =', isz1, jsz1, isz2, jsz2 + if (DEBUG_LOCAL) write (OUTPUT_UNIT, *) 'ir, jz =', ir, jr + + if (ir /= sr_x .or. jr /= sr_y) call Stop_simulation ('Tile dims do not preserve fire/atm ratio') + + end subroutine Check_dims + + end subroutine Provide_atm_feedback + + subroutine Fire_tendency( & + ids,ide, kds,kde, jds,jde, & ! dimensions + ims,ime, kms,kme, jms,jme, & + its,ite, kts,kte, jts,jte, & + grnhfx,grnqfx,canhfx,canqfx, & ! heat fluxes summed up to atm grid + alfg,alfc,z1can, & ! coeffients, properties, geometry + z_at_w,dz8w,mu,c1h,c2h,rho, & + fire_atm_feedback, & + rthfrten,rqvfrten) ! theta and Qv tendencies + + ! This routine is atmospheric physics + + ! --- this routine takes fire generated heat and moisture fluxes and + ! calculates their influence on the theta and water vapor + ! --- note that these tendencies are valid at the Arakawa-A location + + use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT + implicit none + + ! --- incoming variables + + integer, intent (in) :: ids, ide, kds, kde, jds, jde, & + ims, ime, kms, kme, jms, jme, & + its, ite, kts, kte, jts, jte + + real, intent(in), dimension(ims:ime, jms:jme) :: grnhfx,grnqfx ! w/m^2 + real, intent(in), dimension(ims:ime, jms:jme) :: canhfx,canqfx ! w/m^2 + real, intent(in), dimension(ims:ime, jms:jme) :: mu ! dry air mass (pa) + real, intent(in), dimension(kms:kme) :: c1h, c2h ! hybrid coordinate weights + + real, intent(in), dimension(ims:ime, kms:kme, jms:jme) :: z_at_w ! m abv sealvl + real, intent(in), dimension(ims:ime, kms:kme, jms:jme) :: dz8w ! dz across w-lvl + real, intent(in), dimension(ims:ime, kms:kme, jms:jme) :: rho ! density + + real, intent(in) :: alfg ! extinction depth surface fire heat (m) + real, intent(in) :: alfc ! extinction depth crown fire heat (m) + real, intent(in) :: z1can ! height of crown fire heat release (m) + + real, intent(in) :: fire_atm_feedback + + ! --- outgoing variables + + real, intent(out), dimension(ims:ime, kms:kme, jms:jme) :: & + rthfrten, & ! theta tendency from fire (in mass units) + rqvfrten ! Qv tendency from fire (in mass units) + ! --- local variables + + integer :: i,j,k + integer :: i_st,i_en, j_st,j_en, k_st,k_en + + real :: cp_i + real :: rho_i + real :: xlv_i + real :: z_w + real :: fact_g, fact_c + real :: alfg_i, alfc_i + + real, dimension( its:ite,kts:kte,jts:jte ) :: hfx,qfx + + +! write (OUTPUT_UNIT, *) 'pajm: its, ite, jts, jte, kts, kte = ', its, ite, jts, jte, kts, kte +! write (OUTPUT_UNIT, *) 'pajm: ids, ide, jds, jde, kds, kde = ', ids, ide, jds, jde, kds, kde +! write (OUTPUT_UNIT, *) 'pajm: alfg,alfc,z1can = ', alfg,alfc,z1can + + do j=jts,jte + do k=kts,min(kte+1,kde) + do i=its,ite + rthfrten(i,k,j)=0. + rqvfrten(i,k,j)=0. + enddo + enddo + enddo + + if (fire_atm_feedback <= 0.0) return + + ! --- set some local constants + + cp_i = 1./cp ! inverse of specific heat + xlv_i = 1./xlv ! inverse of latent heat + alfg_i = 1./alfg + alfc_i = 1./alfc + + ! --- set loop indicies : note that + + i_st = MAX(its,ids+1) + i_en = MIN(ite,ide-1) + k_st = kts + k_en = MIN(kte,kde-1) + j_st = MAX(jts,jds+1) + j_en = MIN(jte,jde-1) + ! --- distribute fluxes + + do j = j_st,j_en + do k = k_st,k_en + do i = i_st,i_en + + ! --- set z (in meters above ground) + + z_w = z_at_w(i,k,j) - z_at_w(i, 1, j) + + ! --- heat flux + + fact_g = cp_i * EXP( - alfg_i * z_w ) + if ( z_w < z1can ) then + fact_c = cp_i + else + fact_c = cp_i * EXP( - alfc_i * (z_w - z1can) ) + end if + hfx(i,k,j) = fact_g * grnhfx(i,j) * fire_atm_feedback+ fact_c * canhfx(i,j) + + ! --- vapor flux + + fact_g = xlv_i * EXP( - alfg_i * z_w ) + if (z_w < z1can) then + fact_c = xlv_i + else + fact_c = xlv_i * EXP( - alfc_i * (z_w - z1can) ) + end if + qfx(i,k,j) = fact_g * grnqfx(i,j) * fire_atm_feedback + fact_c * canqfx(i,j) + +! if ((grnhfx(i,j) * fire_atm_feedback >0.) .and. (k == 1)) then +! write (OUTPUT_UNIT, *) 'masih: grnhfx, grnqfx', grnhfx(i,j), grnqfx(i,j) +! write (OUTPUT_UNIT, *) 'masih: hfx, qfx', hfx(i,1,j), qfx(i,1,j), hfx(i,1,j), qfx(i,1,j) +! end if + + end do + end do + end do + ! --- add flux divergence to tendencies + ! + ! multiply by dry air mass (mu) to eliminate the need to + ! call sr. calculate_phy_tend (in dyn_em/module_em.F) + + do j = j_st,j_en + do k = k_st,k_en-1 + do i = i_st,i_en + + rho_i = 1./rho(i,k,j) + + rthfrten(i,k,j) = - (c1h(k)*mu(i,j)+c2h(k)) * rho_i * (hfx(i,k+1,j)-hfx(i,k,j)) / dz8w(i,k,j) + rqvfrten(i,k,j) = - (c1h(k)*mu(i,j)+c2h(k)) * rho_i * (qfx(i,k+1,j)-qfx(i,k,j)) / dz8w(i,k,j) + + end do + end do + end do + + return + + end subroutine Fire_tendency + subroutine Update_atm_state (this, datetime_now) implicit none diff --git a/nuopc/fire_behavior_nuopc.F90 b/nuopc/fire_behavior_nuopc.F90 index cc60bf7..2c5365f 100644 --- a/nuopc/fire_behavior_nuopc.F90 +++ b/nuopc/fire_behavior_nuopc.F90 @@ -13,6 +13,7 @@ module fire_behavior_nuopc use initialize_mod, only : Init_fire_state use advance_mod, only : Advance_state use constants_mod, only : G, XLV, CP, FVIRT, R_D + use stderrout_mod, only : Stop_simulation implicit none @@ -37,6 +38,9 @@ module fire_behavior_nuopc real(ESMF_KIND_R8), pointer :: ptr_hflx_fire(:,:) real(ESMF_KIND_R8), pointer :: ptr_evap_fire(:,:) real(ESMF_KIND_R8), pointer :: ptr_smoke_fire(:,:) + real(ESMF_KIND_R8), pointer :: ptr_u10(:,:) + real(ESMF_KIND_R8), pointer :: ptr_v10(:,:) + integer :: clb(2), cub(2), clb3(3), cub3(3) logical :: imp_rainrte = .FALSE. logical :: imp_rainacc = .FALSE. @@ -215,6 +219,23 @@ subroutine Advertise(model, rc) file=__FILE__)) & return ! bail out + ! importable field: inst_zonal_wind_height10m + call NUOPC_Advertise(importState, & + StandardName="inst_zonal_wind_height10m", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + + ! importable field: inst_merid_wind_height10m + call NUOPC_Advertise(importState, & + StandardName="inst_merid_wind_height10m", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + + #endif !#define WITHEXPORTFIELDS_disable @@ -574,6 +595,44 @@ subroutine Realize(model, rc) file=__FILE__)) & return ! bail out + ! importable field on Grid: inst_zonal_wind_height10m + field = ESMF_FieldCreate(name="inst_zonal_wind_height10m", grid=fire_grid, & + typekind=ESMF_TYPEKIND_R8, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + call NUOPC_Realize(importState, field=field, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + ! Get Field memory + call ESMF_FieldGet(field, localDe=0, farrayPtr=ptr_u10, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + + ! importable field on Grid: inst_merid_wind_height10m + field = ESMF_FieldCreate(name="inst_merid_wind_height10m", grid=fire_grid, & + typekind=ESMF_TYPEKIND_R8, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + call NUOPC_Realize(importState, field=field, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + ! Get Field memory + call ESMF_FieldGet(field, localDe=0, farrayPtr=ptr_v10, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + #endif ! #ifdef WITHEXPORTFIELDS @@ -763,24 +822,36 @@ subroutine Advance(model, rc) #endif - do j = 1, grid%jfde - do i = 1, grid%ifde - call grid%Interpolate_profile (config_flags, & ! for debug output, <= 0 no output - config_flags%fire_wind_height, & ! interpolation height - grid%kfds, grid%kfde, & ! fire grid dimensions - atm_u3d(i,j,:),atm_v3d(i,j,:), & ! atm grid arrays in - atm_ph(i,j,:), & - grid%uf(i,j),grid%vf(i,j),grid%fz0(i,j)) - - ! avoid arithmatic error - wspd = (grid%uf(i,j) ** 2. + grid%vf(i,j) ** 2.) ** .5 - if (wspd < 0.001) then - grid%uf(i,j) = sign(0.001, grid%uf(i,j)) - grid%vf(i,j) = sign(0.001, grid%vf(i,j)) - endif - - enddo - enddo + select case (config_flags%wind_vinterp_opt) + case (0) + do j = 1, grid%jfde + do i = 1, grid%ifde + call grid%Interpolate_profile (config_flags, & ! for debug output, <= 0 no output + config_flags%fire_wind_height, & ! interpolation height + grid%kfds, grid%kfde, & ! fire grid dimensions + atm_u3d(i,j,:),atm_v3d(i,j,:), & ! atm grid arrays in + atm_ph(i,j,:), & + grid%uf(i,j),grid%vf(i,j),grid%fz0(i,j)) + + ! avoid arithmatic error + wspd = (grid%uf(i,j) ** 2. + grid%vf(i,j) ** 2.) ** .5 + if (wspd < 0.001) then + grid%uf(i,j) = sign(0.001, grid%uf(i,j)) + grid%vf(i,j) = sign(0.001, grid%vf(i,j)) + endif + + enddo + enddo + case (1) + do j = 1, grid%jfde + do i = 1, grid%ifde + grid%uf(i,j) = grid%fuels%waf(int(grid%nfuel_cat(i,j))) * ptr_u10(i,j) + grid%vf(i,j) = grid%fuels%waf(int(grid%nfuel_cat(i,j))) * ptr_v10(i,j) + end do + end do + case default + call Stop_simulation ('Error: wrong wind_vinterp_opt') + end select if (grid%datetime_now == grid%datetime_start) call grid%Save_state () diff --git a/nuopc/wrf_nuopc.F90 b/nuopc/wrf_nuopc.F90 index 3eab1f1..651315d 100644 --- a/nuopc/wrf_nuopc.F90 +++ b/nuopc/wrf_nuopc.F90 @@ -161,6 +161,22 @@ subroutine Advertise(model, rc) ! return ! bail out ! 2D + ! exportable field: inst_zonal_wind_levels + call NUOPC_Advertise(exportState, & + StandardName="inst_zonal_wind_height10m", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + + ! exportable field: inst_merid_wind_levels + call NUOPC_Advertise(exportState, & + StandardName="inst_merid_wind_height10m", rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + ! exportable field: inst_surface_roughness call NUOPC_Advertise(exportState, & StandardName="inst_surface_roughness", rc=rc) @@ -293,7 +309,7 @@ subroutine Realize(model, rc) rc=rc) if(ESMF_STDERRORCHECK(rc)) return - grid = ESMF_GridCreate(name='ATM', & + grid = ESMF_GridCreate(name='ATM', & distgrid=distgrid, coordSys = ESMF_COORDSYS_SPH_DEG, & rc = rc) if(ESMF_STDERRORCHECK(rc)) return @@ -617,6 +633,31 @@ subroutine Realize(model, rc) file=__FILE__)) & return ! bail out + ! exportable field on Grid: inst_zonal_wind_height10m + field = ESMF_FieldCreate(name="inst_zonal_wind_height10m", grid=grid, & + typekind=ESMF_TYPEKIND_R8, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + call NUOPC_Realize(exportState, field=field, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + + ! exportable field on Grid: inst_merid_wind_height10m + field = ESMF_FieldCreate(name="inst_merid_wind_height10m", grid=grid, & + typekind=ESMF_TYPEKIND_R8, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out + call NUOPC_Realize(exportState, field=field, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail out datetime_now = datetime_t (config_flags%start_year, config_flags%start_month, & config_flags%start_day, config_flags%start_hour, config_flags%start_minute, & @@ -759,7 +800,7 @@ subroutine Advance(model, rc) end subroutine subroutine Update_atm_state (datetime) - + type (datetime_t), intent (in) :: datetime @@ -773,7 +814,7 @@ subroutine Update_atm_state (datetime) call state%Get_phl(datetime) ! call state%Get_pres(datetime) - ! convert m to cm + ! convert m to cm ptr_z0(clb(1):cub(1),clb(2):cub(2))= & state%z0(1:size(state%lats, dim=1),1:size(state%lats, dim=2)) * 100.0 ptr_q2(clb(1):cub(1),clb(2):cub(2))= & @@ -793,7 +834,7 @@ subroutine Update_atm_state (datetime) state%phl(1:size(state%lats, dim=1),1:size(state%lats, dim=2), 1:state%kde - 1) ! ptr_pres(clb3(1):cub3(1),clb3(2):cub3(2),clb3(3):cub3(3))= & ! state%pres(1:size(state%lats, dim=1),1:size(state%lats, dim=2), 1:state%kde - 1) - + end subroutine end module wrf_nuopc diff --git a/physics/fire_driver_mod.F90 b/physics/fire_driver_mod.F90 index c30e1c7..ca57a53 100644 --- a/physics/fire_driver_mod.F90 +++ b/physics/fire_driver_mod.F90 @@ -70,6 +70,8 @@ subroutine Init_fire_components (grid, config_flags) end select call grid%ros_param%Init (grid%ifms, grid%ifme, grid%jfms, grid%jfme) + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij) do ij = 1, grid%num_tiles call Extrapol_var_at_bdys (grid%ifms, grid%ifme, grid%jfms, grid%jfme, grid%ifds, grid%ifde, & grid%jfds, grid%jfde, grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), & @@ -82,6 +84,7 @@ subroutine Init_fire_components (grid, config_flags) call grid%ros_param%Set_params (grid%ifms, grid%ifme, grid%jfms, grid%jfme, grid%i_start(ij), grid%i_end(ij), & grid%j_start(ij), grid%j_end(ij), grid%fuels, grid%nfuel_cat, grid%fmc_g) end do + !$OMP END PARALLEL DO end subroutine Init_fire_components @@ -94,8 +97,11 @@ subroutine Advance_fire_components (grid, config_flags) integer, parameter :: PRINT_LEVEL = 1 integer :: ij + logical, parameter :: DEBUG_LOCAL = .false. + if (DEBUG_LOCAL) call Print_message ('Entering Advance_fire_components...') + if (config_flags%fmoist_run) call grid%fmc_param%Advance_fmc_model (config_flags%fmoist_freq, config_flags%fmoist_dt, & grid%itimestep, grid%dt, grid%ifms, grid%ifme, grid%jfms, grid%jfme, & grid%i_start, grid%i_end, grid%j_start, & @@ -103,13 +109,12 @@ subroutine Advance_fire_components (grid, config_flags) grid%fire_rain_old, grid%fire_t2_old, grid%fire_q2_old, grid%fire_psfc_old, grid%fire_rh_fire, config_flags%fuelmc_g, & grid%fmc_g, grid%nfuel_cat, grid%fuels, grid%ros_param) - do ij = 1, grid%num_tiles - call Advance_fire_model (config_flags, grid, & - grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij)) - end do + call Advance_fire_model (config_flags, grid) if (config_flags%fire_print_msg >= PRINT_LEVEL) call Print_summary (config_flags, grid) + if (DEBUG_LOCAL) call Print_message ('Leaving Advance_fire_components...') + end subroutine Advance_fire_components function Calc_domain_stats (fun, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, a, b) result (return_value) diff --git a/physics/fire_model_mod.F90 b/physics/fire_model_mod.F90 index 0332f67..08a4523 100644 --- a/physics/fire_model_mod.F90 +++ b/physics/fire_model_mod.F90 @@ -5,6 +5,7 @@ module fire_model_mod use namelist_mod, only : namelist_t use ros_mod, only : ros_t use state_mod, only: state_fire_t + use stderrout_mod, only : Print_message private @@ -31,7 +32,7 @@ pure subroutine Copy_lfnout_to_lfn (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jf end subroutine Copy_lfnout_to_lfn - subroutine Advance_fire_model (config_flags, grid, i_start, i_end, j_start, j_end) + subroutine Advance_fire_model (config_flags, grid) ! Purpose advance the fire from time_start to time_start + dt @@ -39,22 +40,19 @@ subroutine Advance_fire_model (config_flags, grid, i_start, i_end, j_start, j_en type (namelist_t), intent (in) :: config_flags type (state_fire_t), intent (in out) :: grid - integer, intent (in) :: i_start, i_end, j_start, j_end - integer :: ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme + integer :: ij, ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme real :: tbound, time_start + logical, parameter :: DEBUG_LOCAL = .false. + if (DEBUG_LOCAL) call Print_message ('Entering Advance_fire_model...') + ifds = grid%ifds ifde = grid%ifde jfds = grid%jfds jfde = grid%jfde - ifts = i_start - ifte = i_end - jfts = j_start - jfte = j_end - ifms = grid%ifms ifme = grid%ifme jfms = grid%jfms @@ -62,36 +60,129 @@ subroutine Advance_fire_model (config_flags, grid, i_start, i_end, j_start, j_en time_start = grid%itimestep * grid%dt - call Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, time_start, grid%dt, grid%dx, grid%dy, & + if (DEBUG_LOCAL) call Print_message ('calling Prop_level_set...') + call Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & + grid%num_tiles, grid%i_start, grid%i_end, grid%j_start, grid%j_end, time_start, grid%dt, grid%dx, grid%dy, & config_flags%fire_upwinding, config_flags%fire_viscosity, config_flags%fire_viscosity_bg, config_flags%fire_viscosity_band, & config_flags%fire_viscosity_ngp, config_flags%fire_lsm_band_ngp, tbound, grid%lfn, grid%lfn_0, grid%lfn_1, grid%lfn_2, & grid%lfn_out, grid%tign_g, grid%ros, grid%uf, grid%vf, grid%dzdxf, grid%dzdyf, grid%ros_param) - call Stop_if_close_to_bdy (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, ifds, jfds, ifde, jfde, grid%lfn_out) + if (DEBUG_LOCAL) call Print_message ('calling Stop_if_close_to_bdy...') + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, grid%num_tiles + ifts = grid%i_start(ij) + ifte = grid%i_end(ij) + jfts = grid%j_start(ij) + jfte = grid%j_end(ij) - call Update_ignition_times (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, ifds, jfds, ifde, jfde, & - time_start, grid%dt, grid%lfn, grid%lfn_out, grid%tign_g) + call Stop_if_close_to_bdy (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, ifds, jfds, ifde, jfde, grid%lfn_out) + end do + !$OMP END PARALLEL DO + + if (DEBUG_LOCAL) call Print_message ('calling Update_ignition_times...') + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, grid%num_tiles + ifts = grid%i_start(ij) + ifte = grid%i_end(ij) + jfts = grid%j_start(ij) + jfte = grid%j_end(ij) + + call Update_ignition_times (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, ifds, jfds, ifde, jfde, & + time_start, grid%dt, grid%lfn, grid%lfn_out, grid%tign_g) + end do + !$OMP END PARALLEL DO - call Calc_flame_length (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, & - grid%ros, grid%ros_param%iboros, grid%flame_length, grid%ros_front, grid%fire_area) - if (config_flags%fire_lsm_reinit) call Reinit_level_set (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, & + if (DEBUG_LOCAL) call Print_message ('calling Calc_flame_length...') + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, grid%num_tiles + ifts = grid%i_start(ij) + ifte = grid%i_end(ij) + jfts = grid%j_start(ij) + jfte = grid%j_end(ij) + + call Calc_flame_length (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, & + grid%ros, grid%ros_param%iboros, grid%flame_length, grid%ros_front, grid%fire_area) + end do + !$OMP END PARALLEL DO + + if (DEBUG_LOCAL) call Print_message ('calling Reinit_level_set...') + if (config_flags%fire_lsm_reinit) call Reinit_level_set (grid%num_tiles, grid%i_start, grid%i_end, grid%j_start, grid%j_end, & + ifms, ifme, jfms, jfme, & ifds, ifde, jfds, jfde, time_start, grid%dt, grid%dx, grid%dy, config_flags%fire_upwinding_reinit, & config_flags%fire_lsm_reinit_iter, config_flags%fire_lsm_band_ngp, grid%lfn, grid%lfn_2, grid%lfn_s0, & grid%lfn_s1, grid%lfn_s2, grid%lfn_s3, grid%lfn_out, grid%tign_g) - call Copy_lfnout_to_lfn (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, grid%lfn_out, grid%lfn) + if (DEBUG_LOCAL) call Print_message ('calling Copy_lfnout_to_lfn...') + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, grid%num_tiles + ifts = grid%i_start(ij) + ifte = grid%i_end(ij) + jfts = grid%j_start(ij) + jfte = grid%j_end(ij) - call Ignite_prescribed_fires (grid, config_flags, time_start, ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde) - - call Calc_fuel_left (ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, ifts, ifte, jfts, jfte, & - grid%lfn,grid%tign_g,grid%fuel_time, time_start + grid%dt, grid%fuel_frac, grid%fire_area, & - grid%fuel_frac_burnt_dt) - - call Calc_fire_fluxes (grid%dt, grid, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, & - ifts, ifte, jfts, jfte, grid%fuel_load_g, grid%fuel_frac_burnt_dt, grid%fgrnhfx, grid%fgrnqfx) + call Copy_lfnout_to_lfn (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, grid%lfn_out, grid%lfn) + end do + !$OMP END PARALLEL DO + + if (DEBUG_LOCAL) call Print_message ('calling Ignite_prescribed_fires...') + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, grid%num_tiles + ifts = grid%i_start(ij) + ifte = grid%i_end(ij) + jfts = grid%j_start(ij) + jfte = grid%j_end(ij) + + call Ignite_prescribed_fires (grid, config_flags, time_start, ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde) + end do + !$OMP END PARALLEL DO + + if (DEBUG_LOCAL) call Print_message ('calling Calc_fuel_left...') + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, grid%num_tiles + ifts = grid%i_start(ij) + ifte = grid%i_end(ij) + jfts = grid%j_start(ij) + jfte = grid%j_end(ij) + call Calc_fuel_left (ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, ifts, ifte, jfts, jfte, & + grid%lfn,grid%tign_g,grid%fuel_time, time_start + grid%dt, grid%fuel_frac, grid%fire_area, & + grid%fuel_frac_burnt_dt) + end do + !$OMP END PARALLEL DO + + if (DEBUG_LOCAL) call Print_message ('calling Calc_fire_fluxes...') + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, grid%num_tiles + ifts = grid%i_start(ij) + ifte = grid%i_end(ij) + jfts = grid%j_start(ij) + jfte = grid%j_end(ij) + call Calc_fire_fluxes (grid%dt, grid, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, & + ifts, ifte, jfts, jfte, grid%fuel_load_g, grid%fuel_frac_burnt_dt, grid%fgrnhfx, grid%fgrnqfx) + end do + !$OMP END PARALLEL DO + + if (DEBUG_LOCAL) call Print_message ('calling Calc_smoke_emissions...') + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, grid%num_tiles + ifts = grid%i_start(ij) + ifte = grid%i_end(ij) + jfts = grid%j_start(ij) + jfte = grid%j_end(ij) + + call Calc_smoke_emissions (grid, config_flags, ifts, ifte, jfts, jfte) + end do + !$OMP END PARALLEL DO - call Calc_smoke_emissions (grid, config_flags, ifts, ifte, jfts, jfte) + if (DEBUG_LOCAL) call Print_message ('Leaving Advance_fire_model...') end subroutine Advance_fire_model diff --git a/physics/fire_physics_mod.F90 b/physics/fire_physics_mod.F90 index f441f6f..b0b6ab9 100644 --- a/physics/fire_physics_mod.F90 +++ b/physics/fire_physics_mod.F90 @@ -4,6 +4,8 @@ module fire_physics_mod use constants_mod, only : XLV, CMBCNST use state_mod, only: state_fire_t use namelist_mod, only: namelist_t + use emis_mod, only : EMIS_WRFFIRE, EMIS_FMC_PM2P5 + use stderrout_mod, only : Stop_simulation private @@ -81,13 +83,44 @@ subroutine Calc_smoke_emissions (grid, config_flags, ifts, ifte, jfts, jfte) integer, intent(in) :: ifts, ifte, jfts, jfte integer :: i, j - - - do j = jfts, jfte - do i = ifts, ifte - grid%emis_smoke(i, j) = config_flags%frac_fburnt_to_smoke * grid%fuel_frac_burnt_dt(i, j) * grid%fuel_load_g(i, j) ! kg/m^2 - end do - end do + real :: emis_pm2p5 + real, parameter :: PERCEN_CARBON_LITTER = 0.5 + + + select case (config_flags%emis_opt) + case (EMIS_WRFFIRE) +! write (*, *) 'Smoke WRF-Fire', config_flags%frac_fburnt_to_smoke + do j = jfts, jfte + do i = ifts, ifte + grid%emis_smoke(i, j) = config_flags%frac_fburnt_to_smoke * grid%fuel_frac_burnt_dt(i, j) * grid%fuel_load_g(i, j) ! kg/m^2 + end do + end do + + case (EMIS_FMC_PM2P5) + ! Assuminmg the fuel is litter for the moment +! write (*, *) 'Smoke CFBM' + do j = jfts, jfte + do i = ifts, ifte + if (grid%fmc_g(i, j) <= 0.035) then + emis_pm2p5 = 8.1 + else if (grid%fmc_g(i, j) <= 0.1) then + emis_pm2p5 = 8.1 + (160.5 - 8.1) / (0.1 - 0.035) * (grid%fmc_g(i, j) - 0.035) + else if (grid%fmc_g(i, j) <= 0.2) then + emis_pm2p5 = 160.5 + (179.9 - 160.5) / (0.2 - 0.1) * (grid%fmc_g(i, j) - 0.1) + else + emis_pm2p5 = 179.9 + end if + ! 0.001 to convert to Kg which would be later coverted to g, until fixed + emis_pm2p5 = emis_pm2p5 * PERCEN_CARBON_LITTER * 0.001 + grid%emis_smoke(i, j) = emis_pm2p5 * grid%fuel_frac_burnt_dt(i, j) * grid%fuel_load_g(i, j) ! kg/m^2 +! if (grid%fuel_load_g(i, j)>0.0) write (*, *) i, j, grid%fmc_g(i, j), emis_pm2p5, grid%emis_smoke(i, j) + end do + end do + + case default + call Stop_simulation ('The emission option selected does not exist.') + + end select end subroutine Calc_smoke_emissions diff --git a/physics/fmc_wrffire_mod.F90 b/physics/fmc_wrffire_mod.F90 index 01d7420..7110ac9 100644 --- a/physics/fmc_wrffire_mod.F90 +++ b/physics/fmc_wrffire_mod.F90 @@ -176,19 +176,28 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms, end if if (this%run_advance_moisture) then + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij) do ij = 1, num_tiles call this%Advance_moisture_classes (itimestep == 1, ifms, ifme, jfms, jfme, i_start(ij), i_end(ij), j_start(ij), j_end(ij), & fire_rain, fire_t2, fire_q2, fire_psfc, fire_rain_old, fire_t2_old, fire_q2_old, fire_psfc_old, fire_rh_fire, fuelmc_g) end do + !$OMP END PARALLEL DO + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij) do ij = 1, num_tiles call this%Average_moisture_classes (ifms, ifme, jfms, jfme, i_start(ij), i_end(ij), j_start(ij), j_end(ij), nfuel_cat, fmc_g) end do + !$OMP END PARALLEL DO + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij) do ij = 1, num_tiles call ros_param%Set_params (ifms, ifme, jfms, jfme, i_start(ij), i_end(ij), j_start(ij), j_end(ij), & fuels, nfuel_cat, fmc_g) end do + !$OMP END PARALLEL DO end if end subroutine Advance_fmc_model diff --git a/physics/fuel_anderson_mod.F90 b/physics/fuel_anderson_mod.F90 index 605cc4a..e9ec224 100644 --- a/physics/fuel_anderson_mod.F90 +++ b/physics/fuel_anderson_mod.F90 @@ -117,6 +117,9 @@ subroutine Init_anderson_fuel_model (this, fuelmc_c) end if end do + ! Initialize wind adjustment factor + call this%Calc_wind_adjustment_factor() + end subroutine Init_anderson_fuel_model end module fuel_anderson_mod diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index cf558e0..a7bdb29 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -14,7 +14,7 @@ module level_set_mod ! https://doi.org/10.1002/2017MS001108 use ros_wrffire_mod, only: ros_wrffire_t - use stderrout_mod, only: Stop_simulation + use stderrout_mod, only: Stop_simulation, Print_message use state_mod, only: state_fire_t use ignition_line_mod, only : ignition_line_t use ros_mod, only : ros_t @@ -267,7 +267,7 @@ subroutine Calc_fuel_left_at_grid_point (lfn00, lfn01, lfn10, & end subroutine Calc_fuel_left_at_grid_point subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & - ifts, ifte, jfts, jfte, ts, dt, dx, dy, fire_upwinding, fire_viscosity, & + num_tiles, i_start, i_end, j_start, j_end, ts, dt, dx, dy, fire_upwinding, fire_viscosity, & fire_viscosity_bg, fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, & tbound, lfn_in, lfn_0, lfn_1, lfn_2, lfn_out, tign, ros, uf, vf, dzdxf, dzdyf, ros_model) @@ -275,8 +275,9 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & implicit none - integer, intent(in) :: ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & + integer, intent(in) :: ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, num_tiles, & fire_upwinding, fire_viscosity_ngp, fire_lsm_band_ngp + integer, dimension (num_tiles), intent (in) :: i_start, i_end, j_start, j_end real, intent(in) :: fire_viscosity, fire_viscosity_bg, fire_viscosity_band real, dimension(ifms:ifme, jfms:jfme), intent (in) :: uf, vf, dzdxf, dzdyf real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: lfn_in, tign, lfn_1, lfn_2, lfn_0 @@ -287,66 +288,157 @@ subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & ! to store tendency (rhs of the level set pde) real, dimension(ifms:ifme, jfms:jfme) :: tend - real :: tbound2, tbound3 - integer :: i, j - character (len = :), allocatable :: msg + real :: tbound2, tbound3, tbound_thread, tbound_min + integer :: i, j, ij, ifts, ifte, jfts, jfte + character (len = 128) :: msg + logical, parameter :: DEBUG_LOCAL = .false. - ! Runge-Kutta step 1 - do j = jfts, jfte - do i = ifts, ifte - lfn_0(i, j) = lfn_in(i, j) + if (DEBUG_LOCAL) call Print_message ('Entering sub Prop_level_set...') + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + do j = jfts, jfte + do i = ifts, ifte + lfn_0(i, j) = lfn_in(i, j) + end do end do end do + !$OMP END PARALLEL DO - call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & - ifms, ifme, jfms, jfme, ts, dt, dx, dy, fire_upwinding, & - fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_0, tbound, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) - - do j = jfts, jfte - do i = ifts, ifte - lfn_1(i, j) = lfn_0(i, j) + (dt / 3.0) * tend(i, j) + ! Runge-Kutta step 1 + if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 1...') + + tbound_min = huge(tbound_min) + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & + ifms, ifme, jfms, jfme, ts, dt, dx, dy, fire_upwinding, & + fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_0, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + + tbound_min = min(tbound_min, tbound_thread) + end do + !$OMP END PARALLEL DO + tbound = tbound_min + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + do j = jfts, jfte + do i = ifts, ifte + lfn_1(i, j) = lfn_0(i, j) + (dt / 3.0) * tend(i, j) + end do end do end do + !$OMP END PARALLEL DO ! Runge-Kutta step 2 - call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & - ifms,ifme,jfms,jfme, ts + dt, dt, dx, dy, fire_upwinding, & - fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_1, tbound2, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) - - do j = jfts, jfte - do i = ifts, ifte - lfn_2(i, j) = lfn_0(i, j) + (dt / 2.0) * tend(i, j) + if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 2...') + + tbound_min = huge(tbound_min) + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Calc_tend_ls (ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, & + ifms,ifme,jfms,jfme, ts + dt, dt, dx, dy, fire_upwinding, & + fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_1, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + + tbound_min = min(tbound_min, tbound_thread) + end do + !$OMP END PARALLEL DO + tbound2 = tbound_min + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + do j = jfts, jfte + do i = ifts, ifte + lfn_2(i, j) = lfn_0(i, j) + (dt / 2.0) * tend(i, j) + end do end do end do + !$OMP END PARALLEL DO ! Runge-Kutta step 3 - call Calc_tend_ls (ifds,ifde,jfds,jfde, ifts, ifte, jfts, jfte, & - ifms, ifme, jfms, jfme, ts + dt, dt, dx, dy, fire_upwinding, & - fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & - fire_viscosity_ngp, fire_lsm_band_ngp, lfn_2, tbound3, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) - - do j = jfts, jfte - do i = ifts, ifte - lfn_out(i, j) = lfn_0(i, j) + dt * tend(i, j) + if (DEBUG_LOCAL) call Print_message ('call Calc_tend_ls 3...') + + tbound_min = huge(tbound_min) + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte, tbound_thread) SHARED(tbound_min) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Calc_tend_ls (ifds,ifde,jfds,jfde, ifts, ifte, jfts, jfte, & + ifms, ifme, jfms, jfme, ts + dt, dt, dx, dy, fire_upwinding, & + fire_viscosity, fire_viscosity_bg, fire_viscosity_band, & + fire_viscosity_ngp, fire_lsm_band_ngp, lfn_2, tbound_thread, tend, ros, uf, vf, dzdxf, dzdyf, ros_model) + + tbound_min = min(tbound_min, tbound_thread) + end do + !$OMP END PARALLEL DO + tbound3 = tbound_min + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + do j = jfts, jfte + do i = ifts, ifte + lfn_out(i, j) = lfn_0(i, j) + dt * tend(i, j) + end do end do - end do + end do + !$OMP END PARALLEL DO ! CFL check, tbound is the max allowed time step - tbound = min (tbound, tbound2, tbound3) + tbound = min(min(tbound, tbound2), tbound3) if (dt > tbound) then !$omp critical - write (msg, '(2(a, f10.2))') 'CFL violation: time step ', dt, ' > bound =', tbound + write (msg, '(a, f10.2, a, f10.2)') 'CFL violation: time step ', dt, ' > bound =', tbound call Stop_simulation (msg) !$omp end critical end if + if (DEBUG_LOCAL) call Print_message ('leaving sub Prop_level_set...') + end subroutine Prop_level_set - subroutine Reinit_level_set (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, & + subroutine Reinit_level_set (num_tiles, i_start, i_end, j_start, j_end, ifms, ifme, jfms, jfme, & ifds, ifde, jfds, jfde, ts, dt, dx, dy, fire_upwinding_reinit, & fire_lsm_reinit_iter, fire_lsm_band_ngp, lfn_in, lfn_2, lfn_s0, & lfn_s1, lfn_s2, lfn_s3, lfn_out, tign) @@ -363,7 +455,9 @@ subroutine Reinit_level_set (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, & implicit none - integer, intent (in) :: ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme + integer, intent (in) :: num_tiles + integer, dimension (num_tiles), intent (in) :: i_start, i_end, j_start, j_end + integer, intent (in) :: ifms, ifme, jfms, jfme integer, intent (in) :: ifds, ifde, jfds, jfde integer, intent (in) :: fire_upwinding_reinit, fire_lsm_reinit_iter, fire_lsm_band_ngp real, dimension (ifms:ifme, jfms:jfme), intent (in out) :: lfn_in, tign @@ -372,62 +466,154 @@ subroutine Reinit_level_set (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme, & real, intent (in) :: dx, dy, ts, dt real :: dt_s, threshold_hlu - integer :: nts, i, j + integer :: nts, i, j, ij, ifts, ifte, jfts, jfte threshold_hlu = fire_lsm_band_ngp * dx ! Define S0 based on current lfn values - do j = jfts, jfte - do i = ifts, ifte - lfn_s0(i, j) = lfn_out(i,j) / sqrt (lfn_out(i, j) ** 2.0 + dx ** 2.0) - lfn_s3(i, j) = lfn_out(i,j) + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + do j = jfts, jfte + do i = ifts, ifte + lfn_s0(i, j) = lfn_out(i,j) / sqrt (lfn_out(i, j) ** 2.0 + dx ** 2.0) + lfn_s3(i, j) = lfn_out(i,j) + end do end do end do + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ifds, ifde, & + jfds, jfde, ifts, ifte, jfts, jfte, lfn_s3) + end do + !$OMP END PARALLEL DO - call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ifds, ifde, & - jfds, jfde, ifts, ifte, jfts, jfte, lfn_s3) + dt_s = 0.01 * dx + dt_s = 0.0001 * dx - dt_s = 0.0001 * dx ! iterate to solve to steady state reinit PDE ! 1 iter each time step is enoguh - do nts = 1, fire_lsm_reinit_iter + Loop_iter: do nts = 1, fire_lsm_reinit_iter ! Runge-Kutta step 1 - call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & - ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & - lfn_s0, lfn_s3, lfn_s3, lfn_s1, 1.0 / 3.0, & ! sign funcition, initial ls, current stage ls, next stage advanced ls, RK coefficient - fire_upwinding_reinit) - - call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ifds, ifde, & - jfds, jfde, ifts, ifte, jfts, jfte, lfn_s1) + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & + ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & + lfn_s0, lfn_s3, lfn_s3, lfn_s1, 1.0 / 3.0, & ! sign funcition, initial ls, current stage ls, next stage advanced ls, RK coefficient + fire_upwinding_reinit) + end do + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ifds, ifde, & + jfds, jfde, ifts, ifte, jfts, jfte, lfn_s1) + end do + !$OMP END PARALLEL DO ! Runge-Kutta step 2 - call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & - ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & - lfn_s0, lfn_s3, lfn_s1, lfn_s2, 1.0 / 2.0, & - fire_upwinding_reinit) - - call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ifds, ifde, & - jfds, jfde, ifts, ifte, jfts, jfte, lfn_s2) + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & + ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & + lfn_s0, lfn_s3, lfn_s1, lfn_s2, 1.0 / 2.0, & + fire_upwinding_reinit) + end do + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ifds, ifde, & + jfds, jfde, ifts, ifte, jfts, jfte, lfn_s2) + end do + !$OMP END PARALLEL DO ! Runge-Kutta step 3 - call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & - ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & - lfn_s0, lfn_s3, lfn_s2, lfn_s3, 1.0, & - fire_upwinding_reinit) - - call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ifds, ifde, & - jfds,jfde, ifts, ifte, jfts, jfte, lfn_s3) - end do - - do j = jfts, jfte - do i = ifts, ifte - ! assing to lfn_out the reinitialized level-set function - lfn_out(i, j) = lfn_s3(i, j) - ! fire area can only increase - lfn_out(i, j) = min (lfn_out(i, j), lfn_in(i, j)) + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Advance_ls_reinit (ifms, ifme, jfms, jfme, ifds, ifde, jfds, jfde, & + ifts, ifte, jfts, jfte, dx, dy, dt_s, threshold_hlu, & + lfn_s0, lfn_s3, lfn_s2, lfn_s3, 1.0, & + fire_upwinding_reinit) + end do + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ifds, ifde, & + jfds,jfde, ifts, ifte, jfts, jfte, lfn_s3) + end do + !$OMP END PARALLEL DO + end do Loop_iter + + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) + do ij = 1, num_tiles + ifts = i_start(ij) + ifte = i_end(ij) + jfts = j_start(ij) + jfte = j_end(ij) + + do j = jfts, jfte + do i = ifts, ifte + ! assing to lfn_out the reinitialized level-set function + lfn_out(i, j) = lfn_s3(i, j) + ! fire area can only increase + lfn_out(i, j) = min (lfn_out(i, j), lfn_in(i, j)) + end do end do end do + !$OMP END PARALLEL DO end subroutine Reinit_level_set @@ -613,15 +799,20 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm threshold_hlu, threshold_av, fire_viscosity_var integer :: i, j character (len = :), allocatable :: msg + logical, parameter :: DEBUG_LOCAL = .false. + + if (DEBUG_LOCAL) call Print_message ('Entering subroutine Calc_tend_ls') threshold_hll = -fire_lsm_band_ngp * dx threshold_hlu = fire_lsm_band_ngp * dx threshold_av = fire_viscosity_ngp * dx + if (DEBUG_LOCAL) call Print_message ('calling Extrapol_var_at_bdys') call Extrapol_var_at_bdys (ifms, ifme, jfms, jfme, ids, ide, jds, jde, & its, ite, jts, jte, lfn) + if (DEBUG_LOCAL) call Print_message ('starting ij loops') tbound = 0.0 do j = jts, jte do i = its, ite @@ -783,6 +974,8 @@ subroutine Calc_tend_ls (ids, ide, jds, jde, its, ite, jts, jte, ifms, ifme, jfm ! final CFL bound tbound = 1.0 / (tbound + TOL) + if (DEBUG_LOCAL) call Print_message ('Leaving subroutine Calc_tend_ls') + end subroutine Calc_tend_ls pure function Select_upwind (diff_lx, diff_rx) result (return_value) diff --git a/share/emis_mod.F90 b/share/emis_mod.F90 new file mode 100644 index 0000000..48292de --- /dev/null +++ b/share/emis_mod.F90 @@ -0,0 +1,11 @@ + module emis_mod + + implicit none + + private + + public :: EMIS_WRFFIRE, EMIS_FMC_PM2P5 + + integer, parameter :: EMIS_WRFFIRE = 0, EMIS_FMC_PM2P5 = 1 + + end module emis_mod diff --git a/share/fuel_mod.F90 b/share/fuel_mod.F90 index 86519dd..529da46 100644 --- a/share/fuel_mod.F90 +++ b/share/fuel_mod.F90 @@ -3,7 +3,7 @@ module fuel_mod implicit none private - + public :: fuel_t, FUEL_ANDERSON, Crosswalk_from_scottburgan_to_anderson integer, parameter :: FUEL_ANDERSON = 1, UNKNOWN_FUEL_CAT = 0 @@ -36,12 +36,15 @@ module fuel_mod real, dimension(:), allocatable :: fuelmce ! fuel loading 1-h, 10-h, 100-h, 1000-h, and live [ton/acre] real, dimension(:), allocatable :: fgi_1h, fgi_10h, fgi_100h, fgi_1000h, fgi_live + ! wind adjustment factor + real, dimension(:), allocatable :: waf contains procedure (Initialization), deferred :: Initialization + procedure, public :: Calc_wind_adjustment_factor => Calc_wind_adjustment_factor end type fuel_t abstract interface - subroutine Initialization (this, fuelmc_c) + subroutine Initialization (this, fuelmc_c) import :: fuel_t class (fuel_t), intent (in out) :: this real, intent(in) :: fuelmc_c @@ -124,4 +127,32 @@ pure function Crosswalk_from_scottburgan_to_anderson (fuel_cat) result (return_v end function Crosswalk_from_scottburgan_to_anderson + subroutine Calc_wind_adjustment_factor (this) + + implicit none + + class (fuel_t), intent (in out) :: this + + ! local + integer :: i + real :: wh, flamelength, h + real :: HF ! extent of the flame above the vegetation + + + allocate (this%waf(this%n_fuel_cat)) + wh = 10.0 ! input wind height in meter + + do i = 1,this%n_fuel_cat + h = this%fueldepthm(i) + flamelength = 2.0 * h ! assume flamelength is double the fuel bed height + if (flamelength > h) then + hf = max(flamelength - h, 0.13 * h) + this%waf(i) = (1.0 + 0.36 * h / hf) / log(((wh-h) + 0.36 * h) / (0.13 * h)) * (log( (hf/h + 0.36)/0.13) - 1.0) + else + this%waf(i) = 0.15 + endif + end do + + end subroutine Calc_wind_adjustment_factor + end module fuel_mod diff --git a/share/interp_mod.F90 b/share/interp_mod.F90 new file mode 100644 index 0000000..57932b3 --- /dev/null +++ b/share/interp_mod.F90 @@ -0,0 +1,99 @@ + module interp_mod + + implicit none + + private + + public :: Interp_profile + + contains + + subroutine Interp_profile (fire_lsm_zcoupling, fire_lsm_zcoupling_ref, fire_wind_height, kfds, kfde, & + uin, vin, z_at_w, z0f, uout, vout) + + implicit none + + real, intent (in) :: fire_wind_height, fire_lsm_zcoupling_ref + integer, intent (in) :: kfds, kfde + real, dimension(:), intent (in) :: uin, vin, z_at_w + real, intent (in) :: z0f + logical, intent (in) :: fire_lsm_zcoupling + real, intent (out) :: uout, vout + + + real, parameter :: VK_KAPPA = 0.4 + real, dimension (kfds:kfde - 1) :: altw, hgt + integer :: k, kdmax + real :: loght, loglast, logz0, logfwh, ht, r_nan, fire_wind_height_local, z0fc, & + ust_d, wsf, wsf1, uf_temp, vf_temp + + + ! max layer to interpolate from, can be less + kdmax = kfde - 2 + do k = kfds, kdmax + 1 + ! altitude of the bottom w-point +! altw(k) = phl(k) / G + altw(k) = z_at_w(k) + end do + + do k = kfds, kdmax + ! height of the mass point above the ground + hgt(k) = 0.5 * (altw(k) + altw(k + 1)) - altw(kfds) + end do + + ! extrapolate mid-flame height from fire_lsm_zcoupling_ref? + if (fire_lsm_zcoupling) then + logfwh = log (fire_lsm_zcoupling_ref) + fire_wind_height_local = fire_lsm_zcoupling_ref + else + logfwh = log (fire_wind_height) + fire_wind_height_local = fire_wind_height + end if + + ! interpolate u + if (fire_wind_height_local > z0f)then + do k = kfds, kdmax + ht = hgt(k) + if (ht >= fire_wind_height_local) then + ! found layer k this point is in + loght = log(ht) + if (k == kfds) then + ! first layer, log linear interpolation from 0 at zr + logz0 = log(z0f) + uout = uin(k) * (logfwh - logz0) / (loght - logz0) + vout = vin(k) * (logfwh - logz0) / (loght - logz0) + else + ! log linear interpolation + loglast = log (hgt(k - 1)) + uout = uin(k - 1) + (uin(k) - uin(k - 1)) * (logfwh - loglast) / (loght - loglast) + vout = vin(k - 1) + (vin(k) - vin(k - 1)) * (logfwh - loglast) / (loght - loglast) + end if + exit + end if + if (k == kdmax) then + ! last layer, still not high enough + uout = uin(k) + vout = vin(k) + end if + end do + else + ! roughness higher than the fire wind height + uout = 0.0 + vout = 0.0 + end if + + ! Extrapol wind to target height + if (fire_lsm_zcoupling) then + uf_temp = uout + vf_temp = vout + wsf = max (sqrt (uf_temp ** 2.0 + vf_temp ** 2.0), 0.1) + z0fc = z0f + ust_d = wsf * VK_KAPPA / log(fire_lsm_zcoupling_ref / z0fc) + wsf1 = (ust_d / VK_KAPPA) * log((fire_wind_height + z0fc) / z0fc) + uout = wsf1 * uf_temp / wsf + vout = wsf1 * vf_temp / wsf + end if + + end subroutine Interp_profile + + end module interp_mod diff --git a/state/state_mod.F90 b/state/state_mod.F90 index 64d333c..ce9fb9a 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -72,6 +72,7 @@ module state_mod class (ros_t), allocatable :: ros_param type (ignition_line_t) :: ignition_lines class (fmc_t), allocatable :: fmc_param + type (proj_lc_t) :: proj ! New vars defined on fire grid for NUOPC coupling real, dimension(:, :), allocatable :: fire_psfc ! "Surface Pressure" "Pa" @@ -91,6 +92,7 @@ module state_mod integer :: ny ! "number of latitudinal grid points" "1" real :: cen_lat, cen_lon contains + procedure, public :: Allocate_vars => Allocate_vars procedure, public :: Convert_sb_to_ander => Convert_scottburgan_to_anderson procedure, public :: Handle_output => Handle_output procedure, public :: Handle_wrfdata_update => Handle_wrfdata_update @@ -99,14 +101,70 @@ module state_mod procedure, public :: Init_ignition_lines => Init_ignition_lines procedure :: Init_latlons => Init_latlons procedure :: Init_tiles => Init_tiles + procedure :: Init_tiles_in_wrf => Init_tiles_in_wrf procedure :: Interpolate_vars_atm_to_fire => Interpolate_vars_atm_to_fire procedure, public :: Interpolate_profile => Interpolate_profile procedure, public :: Print => Print_domain ! private + procedure, public :: Print_tiles => Print_tiles procedure, public :: Save_state => Save_state + procedure, public :: Set_vars_to_default => Set_vars_to_default + procedure, public :: Set_time_stamps => Set_time_stamps end type state_fire_t contains + subroutine Allocate_vars (this, ifms, ifme, jfms, jfme) + + implicit none + + class (state_fire_t), intent(in out) :: this + integer, intent (in) :: ifms, ifme, jfms, jfme + + + allocate (this%uf(ifms:ifme, jfms:jfme)) + allocate (this%vf(ifms:ifme, jfms:jfme)) + allocate (this%fmc_g(ifms:ifme, jfms:jfme)) + allocate (this%lfn(ifms:ifme, jfms:jfme)) + allocate (this%lfn_hist(ifms:ifme, jfms:jfme)) + allocate (this%lfn_0(ifms:ifme, jfms:jfme)) + allocate (this%lfn_1(ifms:ifme, jfms:jfme)) + allocate (this%lfn_2(ifms:ifme, jfms:jfme)) + allocate (this%lfn_s0(ifms:ifme, jfms:jfme)) + allocate (this%lfn_s1(ifms:ifme, jfms:jfme)) + allocate (this%lfn_s2(ifms:ifme, jfms:jfme)) + allocate (this%lfn_s3(ifms:ifme, jfms:jfme)) + allocate (this%lfn_out(ifms:ifme, jfms:jfme)) + allocate (this%fuel_load_g(ifms:ifme, jfms:jfme)) + allocate (this%flame_length(ifms:ifme, jfms:jfme)) + allocate (this%ros_front(ifms:ifme, jfms:jfme)) + allocate (this%tign_g(ifms:ifme, jfms:jfme)) + allocate (this%fuel_frac(ifms:ifme, jfms:jfme)) + allocate (this%fire_area(ifms:ifme, jfms:jfme)) + allocate (this%fuel_frac_burnt_dt(ifms:ifme, jfms:jfme)) + allocate (this%fgrnhfx(ifms:ifme, jfms:jfme)) + allocate (this%fgrnqfx(ifms:ifme, jfms:jfme)) + allocate (this%fcanhfx(ifms:ifme, jfms:jfme)) + allocate (this%fcanqfx(ifms:ifme, jfms:jfme)) + allocate (this%ros(ifms:ifme, jfms:jfme)) + allocate (this%fz0(ifms:ifme, jfms:jfme)) + allocate (this%fuel_time(ifms:ifme, jfms:jfme)) + allocate (this%fire_psfc(ifms:ifme, jfms:jfme)) + allocate (this%fire_rain(ifms:ifme, jfms:jfme)) + allocate (this%fire_t2(ifms:ifme, jfms:jfme)) + allocate (this%fire_q2(ifms:ifme, jfms:jfme)) + allocate (this%fire_rh_fire(ifms:ifme, jfms:jfme)) + allocate (this%fire_psfc_old(ifms:ifme, jfms:jfme)) + allocate (this%fire_rain_old(ifms:ifme, jfms:jfme)) + allocate (this%fire_t2_old(ifms:ifme, jfms:jfme)) + allocate (this%fire_q2_old(ifms:ifme, jfms:jfme)) + allocate (this%zsf(ifms:ifme, jfms:jfme)) + allocate (this%dzdxf(ifms:ifme, jfms:jfme)) + allocate (this%dzdyf(ifms:ifme, jfms:jfme)) + allocate (this%nfuel_cat(ifms:ifme, jfms:jfme)) + allocate (this%emis_smoke(ifms:ifme, jfms:jfme)) + + end subroutine Allocate_vars + subroutine Convert_scottburgan_to_anderson (this) implicit none @@ -116,6 +174,8 @@ subroutine Convert_scottburgan_to_anderson (this) integer :: i, j, ij, ifts, ifte, jfts, jfte + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, ifts, ifte, jfts, jfte) do ij = 1, this%num_tiles ifts = this%i_start(ij) ifte = this%i_end(ij) @@ -128,6 +188,7 @@ subroutine Convert_scottburgan_to_anderson (this) end do end do end do + !$OMP END PARALLEL DO end subroutine Convert_scottburgan_to_anderson @@ -174,154 +235,225 @@ subroutine Handle_wrfdata_update (this, wrf, config_flags) end subroutine Handle_wrfdata_update - subroutine Init_domain (this, config_flags, geogrid) + subroutine Init_domain (this, config_flags, geogrid, & + ifds, ifde, ifms, ifme, ifps, ifpe, & + jfds, jfde, jfms, jfme, jfps, jfpe, & + kfds, kfde, kfms, kfme, kfps, kfpe, & + kfts, kfte, ide, jde, & + cen_lat, cen_lon, truelat1, truelat2, stand_lon, & + dx, dy, sr_x, sr_y, nfuel_cat, zsf, dzdxf, dzdyf) + implicit none class (state_fire_t), intent(in out) :: this type (namelist_t), intent (in) :: config_flags - type (geogrid_t), intent (in) :: geogrid - + type (geogrid_t), intent (in), optional :: geogrid + integer, intent (in), optional :: ifds, ifde, ifms, ifme, ifps, ifpe, & + jfds, jfde, jfms, jfme, jfps, jfpe, & + kfds, kfde, kfms, kfme, kfps, kfpe, & + kfts, kfte, sr_x, sr_y, ide, jde + real, intent (in), optional :: cen_lat, cen_lon, truelat1, truelat2, stand_lon, dx, dy + real, dimension(:, :), intent (in), optional :: nfuel_cat, zsf, dzdxf, dzdyf + + integer, parameter :: INIT_MODE_NONE = 0, INIT_MODE_GEOGRID = 1, INIT_MODE_WRF = 2, INIT_MODE_IDEAL = 3 + type (proj_lc_t) :: proj logical, parameter :: DEBUG_LOCAL = .false. - integer :: ids0, ide0, jds0, jde0 - - - ! Domain dimensions - ids0 = geogrid%ifds - ide0 = geogrid%ifde - jds0 = geogrid%jfds - jde0 = geogrid%jfde - - this%ifds = ids0 - this%ifde = ide0 - this%ifms = ids0 - N_POINTS_IN_HALO - this%ifme = ide0 + N_POINTS_IN_HALO - this%ifps = ids0 - this%ifpe = ide0 - - this%jfds = jds0 - this%jfde = jde0 - this%jfms = jds0 - N_POINTS_IN_HALO - this%jfme = jde0 + N_POINTS_IN_HALO - this%jfps = jds0 - this%jfpe = jde0 - - this%kfds = config_flags%kds - this%kfde = config_flags%kde - this%kfms = config_flags%kds - this%kfme = config_flags%kde - this%kfps = config_flags%kds - this%kfpe = config_flags%kde - this%kfts = config_flags%kds - this%kfte = config_flags%kde - - ! Datetimes - this%datetime_start = datetime_t (config_flags%start_year, config_flags%start_month, config_flags%start_day, & - config_flags%start_hour, config_flags%start_minute, config_flags%start_second) - this%datetime_end = datetime_t (config_flags%end_year, config_flags%end_month, config_flags%end_day, & - config_flags%end_hour, config_flags%end_minute, config_flags%end_second) - this%datetime_now = this%datetime_start + integer :: ids0, ide0, jds0, jde0, i, j, init_mode + + + init_mode = INIT_MODE_NONE + if (config_flags%ideal_opt == 1) init_mode = INIT_MODE_IDEAL + if (present (geogrid)) init_mode = INIT_MODE_GEOGRID + if (present (ifds) .and. present (ifde) .and. present (ifms) .and. present (ifme) .and. present (ifps) .and. present (ifpe) .and. & + present (jfds) .and. present (jfde) .and. present (jfms) .and. present (jfme) .and. present (jfps) .and. present (jfpe) .and. & + present (kfds) .and. present (kfde) .and. present (kfms) .and. present (kfme) .and. present (kfps) .and. present (kfpe) .and. & + present (kfts) .and. present (kfte) .and. present (ide) .and. present (jde) .and. & + present (cen_lat) .and. present (cen_lon) .and. present (truelat1) .and. present (truelat2) .and. present (stand_lon) .and. & + present (dx) .and. present (dy) .and. present (sr_x) .and. present (sr_y) .and. present (nfuel_cat) .and. present (zsf) .and. & + present (dzdxf) .and. present (dzdyf)) & + init_mode = INIT_MODE_WRF + + if (init_mode == INIT_MODE_NONE) & + call Stop_simulation ('Not enough information to initialize domain') + + ! Set dimensions + Set_dims: select case (init_mode) + case (INIT_MODE_GEOGRID, INIT_MODE_IDEAL) + if (init_mode == INIT_MODE_GEOGRID) then + ids0 = geogrid%ifds + ide0 = geogrid%ifde + jds0 = geogrid%jfds + jde0 = geogrid%jfde + else if (init_mode == INIT_MODE_IDEAL) then + ids0 = 1 + ide0 = config_flags%nx + jds0 = 1 + jde0 = config_flags%ny + end if + + this%ifds = ids0 + this%ifde = ide0 + this%ifms = ids0 - N_POINTS_IN_HALO + this%ifme = ide0 + N_POINTS_IN_HALO + this%ifps = ids0 + this%ifpe = ide0 + + this%jfds = jds0 + this%jfde = jde0 + this%jfms = jds0 - N_POINTS_IN_HALO + this%jfme = jde0 + N_POINTS_IN_HALO + this%jfps = jds0 + this%jfpe = jde0 + + this%kfds = config_flags%kds + this%kfde = config_flags%kde + this%kfms = config_flags%kds + this%kfme = config_flags%kde + this%kfps = config_flags%kds + this%kfpe = config_flags%kde + this%kfts = config_flags%kds + this%kfte = config_flags%kde + + call this%Init_tiles (config_flags) + + case (INIT_MODE_WRF) + this%ifds = ifds + this%ifde = ifde + this%ifms = ifms + this%ifme = ifme + this%ifps = ifps + this%ifpe = ifpe + + this%jfds = jfds + this%jfde = jfde + this%jfms = jfms + this%jfme = jfme + this%jfps = jfps + this%jfpe = jfpe + + this%kfds = kfds + this%kfde = kfde + this%kfms = kfms + this%kfme = kfme + this%kfps = kfps + this%kfpe = kfpe + this%kfts = kfts + this%kfte = kfte + + call this%Init_tiles_in_wrf (config_flags, sr_x, sr_y) + + case default + + call Stop_simulation ('Not ready to complete fire state initialization 1') + + end select Set_dims + + call this%Print_tiles () - this%datetime_next_output = this%datetime_start - call this%datetime_next_output%Add_seconds (config_flags%interval_output) + this%nx = this%ifde + this%ny = this%jfde + this%dt = config_flags%dt - this%datetime_next_atm_update = this%datetime_start + ! Init memory + call this%Allocate_vars (this%ifms, this%ifme, this%jfms, this%jfme) - this%cen_lat = geogrid%cen_lat - this%cen_lon = geogrid%cen_lon + ! Set projection + Set_proj: select case (init_mode) + case (INIT_MODE_GEOGRID) + proj = geogrid%Get_atm_proj () + call this%Init_latlons (proj, srx = geogrid%sr_x, sry = geogrid%sr_y) - this%dx = geogrid%dx / geogrid%sr_x - this%dy = geogrid%dy / geogrid%sr_y + this%cen_lat = geogrid%cen_lat + this%cen_lon = geogrid%cen_lon - if (DEBUG_LOCAL) call this%Print() + this%dx = geogrid%dx / geogrid%sr_x + this%dy = geogrid%dy / geogrid%sr_y - this%nx = this%ifde - this%ny = this%jfde + case (INIT_MODE_WRF) + proj = proj_lc_t (cen_lat = cen_lat , cen_lon = cen_lon, dx = dx, dy = dy, & + standard_lon = stand_lon, true_lat_1 = truelat1, true_lat_2 = truelat2, nx = ide - 1, ny = jde - 1) + call this%Init_latlons (proj, srx = sr_x, sry = sr_y) - allocate (this%uf(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%vf(this%ifms:this%ifme, this%jfms:this%jfme)) - this%uf = 0. - this%vf = 0. - allocate (this%fmc_g(this%ifms:this%ifme, this%jfms:this%jfme)) - this%fmc_g = config_flags%fuelmc_g + this%cen_lat = cen_lat + this%cen_lon = cen_lon - ! Init lfn more than the largest domain side - allocate (this%lfn(this%ifms:this%ifme, this%jfms:this%jfme)) - this%lfn(this%ifds:this%ifde, this%jfds:this%jfde) = 2.0 * & - max ((this%ifde - this%ifds + 1) * this%dx, (this%jfde - this%jfds + 1) * this%dy) + this%dx = dx / sr_x + this%dy = dy / sr_y - allocate (this%lfn_hist(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%lfn_0(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%lfn_1(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%lfn_2(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%lfn_s0(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%lfn_s1(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%lfn_s2(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%lfn_s3(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%lfn_out(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fuel_load_g(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%flame_length(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%ros_front(this%ifms:this%ifme, this%jfms:this%jfme)) + case (INIT_MODE_IDEAL) + this%dx = config_flags%dx + this%dy = config_flags%dy - ! Init tign_g a bit into the future - allocate (this%tign_g(this%ifms:this%ifme, this%jfms:this%jfme)) - this%tign_g(this%ifps:this%ifpe, this%jfps:this%jfpe) = epsilon (this%tign_g) + this%cen_lat = config_flags%cen_lat + this%cen_lon = config_flags%cen_lon - allocate (this%fuel_frac(this%ifms:this%ifme, this%jfms:this%jfme)) - this%fuel_frac(this%ifds:this%ifde, this%jfds:this%jfde) = 1.0 + proj = proj_lc_t (cen_lat = this%cen_lat , cen_lon = this%cen_lon, dx = this%dx, dy = this%dy, & + standard_lon = config_flags%stand_lon, true_lat_1 = config_flags%true_lat_1, & + true_lat_2 = config_flags%true_lat_2, nx = config_flags%nx, ny = config_flags%ny) - allocate (this%fire_area(this%ifms:this%ifme, this%jfms:this%jfme)) - this%fire_area(this%ifds:this%ifde, this%jfds:this%jfde) = 0.0 + call this%Init_latlons (proj) - allocate (this%fuel_frac_burnt_dt(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fgrnhfx(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fgrnqfx(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fcanhfx(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fcanqfx(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%ros(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fz0(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fuel_time(this%ifms:this%ifme, this%jfms:this%jfme)) - - allocate (this%fire_psfc(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fire_rain(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fire_t2(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fire_q2(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fire_rh_fire(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fire_psfc_old(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fire_rain_old(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fire_t2_old(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%fire_q2_old(this%ifms:this%ifme, this%jfms:this%jfme)) + case default + call Stop_simulation ('Not ready to complete fire state initialization 2') - this%dt = config_flags%dt + end select Set_proj + this%proj = proj - allocate (this%zsf(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%dzdxf(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%dzdyf(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%nfuel_cat(this%ifms:this%ifme, this%jfms:this%jfme)) - allocate (this%emis_smoke(this%ifms:this%ifme, this%jfms:this%jfme)) - this%emis_smoke = 0.0 + ! Init vars + call this%Set_vars_to_default (config_flags) - this%zsf(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%elevations - this%dzdxf(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%dz_dxs - this%dzdyf(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%dz_dys - this%nfuel_cat(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%fuel_cats - - if (config_flags%fire_is_real_perim) then - if (allocated (geogrid%lfn_init)) then - this%lfn_hist(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%lfn_init - else - Call Stop_simulation ('Attenting to initialize fire from given perimeter but no initialization data present') - end if - end if + Set_topo_fuels: select case (init_mode) + case (INIT_MODE_GEOGRID) + this%zsf(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%elevations + this%dzdxf(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%dz_dxs + this%dzdyf(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%dz_dys + this%nfuel_cat(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%fuel_cats - this%unit_fxlat = 2.0 * PI / (360.0 * RERADIUS) ! earth circumference in m / 360 degrees - this%unit_fxlong = cos (this%cen_lat * 2.0 * PI / 360.0) * this%unit_fxlat ! latitude - call this%Init_latlons (geogrid) + if (config_flags%fire_is_real_perim) then + if (allocated (geogrid%lfn_init)) then + this%lfn_hist(this%ifds:this%ifde, this%jfds:this%jfde) = geogrid%lfn_init + else + Call Stop_simulation ('Attenting to initialize fire from given perimeter but no initialization data present') + end if + end if - call this%Init_tiles (config_flags) + case (INIT_MODE_WRF) + this%zsf(this%ifms:this%ifme, this%jfms:this%jfme) = zsf + this%dzdxf(this%ifms:this%ifme, this%jfms:this%jfme) = dzdxf + this%dzdyf(this%ifms:this%ifme, this%jfms:this%jfme) = dzdyf + this%nfuel_cat(this%ifms:this%ifme, this%jfms:this%jfme) = nfuel_cat + if (config_flags%fire_is_real_perim) & + !this%lfn_hist(this%ifms:this%ifme, this%jfms:this%jfme) = lfn_hist + call Stop_simulation ('Not ready to initialize from fire perimeter inside WRF') + + case (INIT_MODE_IDEAL) + do j = this%jfds, this%jfde + do i = this%ifds, this%ifde + this%zsf(i, j) = config_flags%elevation + & + (i - this%ifds) * config_flags%dz_dx * config_flags%dx + & + (j - this%jfds) * config_flags%dz_dy * config_flags%dy + end do + end do + this%dzdxf(this%ifds:this%ifde, this%jfds:this%jfde) = config_flags%dz_dx + this%dzdyf(this%ifds:this%ifde, this%jfds:this%jfde) = config_flags%dz_dy + this%nfuel_cat(this%ifds:this%ifde, this%jfds:this%jfde) = config_flags%fuel_cat + + if (config_flags%fire_is_real_perim) & + call Stop_simulation ('Not ready to initialize from fire perimeter in idealized mode') + + case default + call Stop_simulation ('Not ready to complete fire state initialization 3') + + end select Set_topo_fuels if (config_flags%fuel_opt == FUEL_ANDERSON) call this%Convert_sb_to_ander () + ! Set clock + call this%Set_time_stamps (config_flags) + + if (DEBUG_LOCAL) call this%Print() + end subroutine Init_domain subroutine Init_fuel_vars (this) @@ -333,6 +465,8 @@ subroutine Init_fuel_vars (this) integer :: ij, i, j, ifts, ifte, jfts, jfte, k + !$OMP PARALLEL DO & + !$OMP PRIVATE (ij, i, j, k, ifts, ifte, jfts, jfte) do ij = 1, this%num_tiles ifts = this%i_start(ij) ifte = this%i_end(ij) @@ -355,6 +489,7 @@ subroutine Init_fuel_vars (this) end do end do end do + !$OMP END PARALLEL DO end subroutine Init_fuel_vars @@ -371,33 +506,39 @@ subroutine Init_ignition_lines (this, config_flags) end subroutine Init_ignition_lines - subroutine Init_latlons (this, geogrid) + subroutine Init_latlons (this, proj, srx, sry) implicit none class (state_fire_t), intent (in out) :: this - type (geogrid_t), intent(in) :: geogrid + type (proj_lc_t), intent(in) :: proj + integer, optional :: srx, sry real, parameter :: OFFSET = 0.5 - type (proj_lc_t) :: proj - integer :: i, j + integer :: i, j, sr_x, sr_y real :: i_atm, j_atm, offset_corners_x, offset_corners_y + if (present (srx) .and. present (sry)) then + sr_x = srx + sr_y = sry + else + sr_x = 1 + sr_y = 1 + end if + allocate (this%lons(this%ifms:this%ifme, this%jfms:this%jfme)) allocate (this%lats(this%ifms:this%ifme, this%jfms:this%jfme)) allocate (this%lons_c(this%nx + 1, this%ny + 1)) allocate (this%lats_c(this%nx + 1, this%ny + 1)) - proj = geogrid%Get_atm_proj () - - offset_corners_x = (1.0 / real (geogrid%sr_x)) / 2.0 - offset_corners_y = (1.0 / real (geogrid%sr_y)) / 2.0 + offset_corners_x = (1.0 / real (sr_x)) / 2.0 + offset_corners_y = (1.0 / real (sr_y)) / 2.0 do j = 1, this%ny do i = 1, this%nx - i_atm = (i - OFFSET) / geogrid%sr_x + OFFSET - j_atm = (j - OFFSET) / geogrid%sr_y + OFFSET + i_atm = (i - OFFSET) / sr_x + OFFSET + j_atm = (j - OFFSET) / sr_y + OFFSET call proj%Calc_latlon (i = i_atm, j = j_atm, lat = this%lats(i, j), lon = this%lons(i, j)) call proj%Calc_latlon (i = i_atm - offset_corners_x, j = j_atm - offset_corners_y, & lat = this%lats_c(i, j), lon = this%lons_c(i, j)) @@ -405,21 +546,21 @@ subroutine Init_latlons (this, geogrid) end do do j = 1, this%ny - i_atm = (this%nx - OFFSET) / geogrid%sr_x + OFFSET - j_atm = (j - OFFSET) / geogrid%sr_y + OFFSET + i_atm = (this%nx - OFFSET) / sr_x + OFFSET + j_atm = (j - OFFSET) / sr_y + OFFSET call proj%Calc_latlon (i = i_atm + offset_corners_x, j = j_atm - offset_corners_y, & lat = this%lats_c(this%nx + 1, j), lon = this%lons_c(this%nx + 1, j)) end do do i = 1, this%nx - i_atm = (i - OFFSET) / geogrid%sr_x + OFFSET - j_atm = (this%ny - OFFSET) / geogrid%sr_y + OFFSET + i_atm = (i - OFFSET) / sr_x + OFFSET + j_atm = (this%ny - OFFSET) / sr_y + OFFSET call proj%Calc_latlon (i = i_atm - offset_corners_x, j = j_atm + offset_corners_y, & lat = this%lats_c(i, this%ny + 1), lon = this%lons_c(i, this%ny + 1)) end do - i_atm = (this%nx - OFFSET) / geogrid%sr_x + OFFSET - j_atm = (this%ny - OFFSET) / geogrid%sr_y + OFFSET + i_atm = (this%nx - OFFSET) / sr_x + OFFSET + j_atm = (this%ny - OFFSET) / sr_y + OFFSET call proj%Calc_latlon (i = i_atm + offset_corners_x, j = j_atm + offset_corners_y, & lat = this%lats_c(this%nx + 1, this%ny + 1), lon = this%lons_c(this%nx + 1, this%ny + 1)) @@ -432,18 +573,49 @@ subroutine Init_tiles (this, config_flags) class (state_fire_t), intent(in out) :: this type (namelist_t), intent (in) :: config_flags - integer :: num_tiles - num_tiles = config_flags%num_tiles - call Calc_tiles_dims (this%ifps, this%ifpe, this%jfps, this%jfpe, num_tiles, & + this%num_tiles = config_flags%num_tiles + call Calc_tiles_dims (this%ifps, this%ifpe, this%jfps, this%jfpe, this%num_tiles, config_flags%tile_strategy, & this%i_start, this%i_end, this%j_start, this%j_end) - if (num_tiles /= config_flags%num_tiles) then + if (this%num_tiles /= config_flags%num_tiles) then call Stop_simulation ('Not able to use the number of tiles specified') end if end subroutine Init_tiles + subroutine Init_tiles_in_wrf (this, config_flags, sr_x, sr_y) + + implicit none + + class (state_fire_t), intent(in out) :: this + type (namelist_t), intent (in) :: config_flags + integer, intent (in) :: sr_x, sr_y + + integer :: ips, ipe, jps, jpe, ij + + + this%num_tiles = config_flags%num_tiles + ips = (this%ifps - 1) / sr_x + 1 + ipe = this%ifpe / sr_x + jps = (this%jfps - 1) / sr_y + 1 + jpe = this%jfpe / sr_y + call Calc_tiles_dims (ips, ipe, jps, jpe, this%num_tiles, config_flags%tile_strategy, & + this%i_start, this%i_end, this%j_start, this%j_end) + + if (this%num_tiles /= config_flags%num_tiles) then + call Stop_simulation ('Not able to use the number of tiles specified') + end if + + do ij = 1, this%num_tiles + this%i_start(ij) = this%i_start(ij) * sr_x - sr_x + 1 + this%i_end(ij) = this%i_end(ij) * sr_x + this%j_start(ij) = this%j_start(ij) * sr_y - sr_y + 1 + this%j_end(ij) = this%j_end(ij) * sr_y + end do + + end subroutine Init_tiles_in_wrf + subroutine Interpolate_vars_atm_to_fire (this, wrf, config_flags) implicit none @@ -617,6 +789,24 @@ subroutine Print_domain (this) end subroutine Print_domain + subroutine Print_tiles (this) + + implicit none + + class (state_fire_t), intent(in) :: this + + integer :: ij + character (len = 300) :: msg + + + do ij = 1, this%num_tiles + write (msg, '(a10, 1x, i3, a4, i7, a4, i7, a4, i7, a4, i7)') & + 'CFBM TILE', ij, ' IS', this%i_start(ij), ' IE', this%i_end(ij), ' JS', this%j_start(ij), ' JE', this%j_end(ij) + call Print_message (trim (msg)) + end do + + end subroutine Print_tiles + subroutine Save_state (this) implicit none @@ -655,5 +845,58 @@ subroutine Save_state (this) end subroutine Save_state + subroutine Set_time_stamps (this, config_flags) + + implicit none + + class (state_fire_t), intent (in out) :: this + type (namelist_t), intent (in) :: config_flags + + + this%datetime_start = datetime_t (config_flags%start_year, config_flags%start_month, config_flags%start_day, & + config_flags%start_hour, config_flags%start_minute, config_flags%start_second) + this%datetime_end = datetime_t (config_flags%end_year, config_flags%end_month, config_flags%end_day, & + config_flags%end_hour, config_flags%end_minute, config_flags%end_second) + this%datetime_now = this%datetime_start + + this%datetime_next_output = this%datetime_start + call this%datetime_next_output%Add_seconds (config_flags%interval_output) + + this%datetime_next_atm_update = this%datetime_start + + end subroutine Set_time_stamps + + subroutine Set_vars_to_default (this, config_flags) + + implicit none + + class (state_fire_t), intent (in out) :: this + type (namelist_t), intent (in) :: config_flags + + + if (config_flags%ideal_opt == 1) then + this%uf(this%ifds:this%ifde, this%jfds:this%jfde) = config_flags%zonal_wind + this%vf(this%ifds:this%ifde, this%jfds:this%jfde) = config_flags%meridional_wind + else + this%uf = 0.0 + this%vf = 0.0 + end if + this%fmc_g = config_flags%fuelmc_g + ! Init lfn more than the largest domain side + this%lfn(this%ifds:this%ifde, this%jfds:this%jfde) = 2.0 * & + max ((this%ifde - this%ifds + 1) * this%dx, (this%jfde - this%jfds + 1) * this%dy) + ! Init tign_g a bit into the future + this%tign_g(this%ifps:this%ifpe, this%jfps:this%jfpe) = epsilon (this%tign_g) + + this%fuel_frac(this%ifds:this%ifde, this%jfds:this%jfde) = 1.0 + this%fire_area(this%ifds:this%ifde, this%jfds:this%jfde) = 0.0 + + this%emis_smoke = 0.0 + + this%unit_fxlat = 2.0 * PI / (360.0 * RERADIUS) ! earth circumference in m / 360 degrees + this%unit_fxlong = cos (this%cen_lat * 2.0 * PI / 360.0) * this%unit_fxlat ! latitude + + end subroutine Set_vars_to_default + end module state_mod diff --git a/state/tiles_mod.F90 b/state/tiles_mod.F90 index ada8f29..24cf4a9 100644 --- a/state/tiles_mod.F90 +++ b/state/tiles_mod.F90 @@ -1,20 +1,24 @@ module tiles_mod + use stderrout_mod, only : Stop_simulation + implicit none private + integer, parameter :: TILE_NONE = 0, TILE_X = 1, TILE_Y = 2, TILE_XY = 3 + public :: Calc_tiles_dims contains - subroutine Calc_tiles_dims (spx, epx, spy, epy, num_tiles, i_start, i_end, j_start, j_end) + subroutine Calc_tiles_dims (spx, epx, spy, epy, num_tiles, tile_strategy, i_start, i_end, j_start, j_end) ! This code is borrowed from module_tiles.F in WRF, specifically SUBROUTINE set_tiles2 implicit none - integer, intent (in) :: spx, epx, spy, epy + integer, intent (in) :: spx, epx, spy, epy, tile_strategy integer, intent (in out) :: num_tiles integer, parameter :: MIN_TILES_IN_X = 1, MIN_TILES_IN_Y = 1 @@ -37,7 +41,23 @@ subroutine Calc_tiles_dims (spx, epx, spy, epy, num_tiles, i_start, i_end, j_sta allocate (j_end(num_tiles)) ! Calc number of tiles in x and y based on total number of tiles - call least_aspect (num_tiles, MIN_TILES_IN_Y, MIN_TILES_IN_X, num_tiles_y, num_tiles_x) + select case (tile_strategy) + case (TILE_NONE, TILE_Y) + if (num_tiles > (epy - spy + 1)) call Stop_simulation ('Number of tiles is too high for TILE_Y strategy') + num_tiles_x = 1 + num_tiles_y = num_tiles + + case (TILE_X) + num_tiles_x = num_tiles + num_tiles_y = 1 + + case (TILE_XY) + call least_aspect (num_tiles, MIN_TILES_IN_Y, MIN_TILES_IN_X, num_tiles_y, num_tiles_x) + + case default + call Stop_simulation ('The tile strategy selected is not valid.') + + end select ! Calc start and end tile indices nt = 1 diff --git a/tests/fd_fire.yaml b/tests/fd_fire.yaml index 3480859..148c451 100644 --- a/tests/fd_fire.yaml +++ b/tests/fd_fire.yaml @@ -73,3 +73,13 @@ alias: inst_temp_height_lowest_from_phys canonical_units: K description: atmosphere export - bottom layer temperature + # + - standard_name: Sa_u10m + alias: inst_zonal_wind_height10m + canonical_units: m s-1 + description: atmosphere export - zonal wind height 10m + # + - standard_name: Sa_v10m + alias: inst_merid_wind_height10m + canonical_units: m s-1 + description: atmosphere export - meridional wind height 10m diff --git a/tests/testx/namelist.fire b/tests/testx/namelist.fire index 58ce971..fffaea8 100644 --- a/tests/testx/namelist.fire +++ b/tests/testx/namelist.fire @@ -37,3 +37,5 @@ fire_upwinding = 9, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian / + &ideal + / diff --git a/tests/testx/testx.yaml b/tests/testx/testx.yaml index 6f0e829..7752071 100644 --- a/tests/testx/testx.yaml +++ b/tests/testx/testx.yaml @@ -60,3 +60,5 @@ ATMD: inst_pres_height_lowest_from_phys: {dim: 2, val: 85000} inst_spec_humid_height_lowest_from_phys: {dim: 2, val: 0.005} inst_temp_height_lowest_from_phys: {dim: 2, val: 310} + inst_zonal_wind_height10m: {dim: 2, val: 1} + inst_merid_wind_height10m: {dim: 2, val: 3}