Merge pull request #650 from NCAR/WRF_vert_hybrid

Add WRF hybrid vertical coordinate

CHANGELOG.rst

@@ -22,22 +22,32 @@ individual files.
 
 The changes are now listed with the most recent at the top.
 
+**April 23 2024 :: Bug-fix: WRF hybrid vertical coordinate. Tag v11.4.1**
+
+- DART now detects whether WRF is using Hybrid Vertical Coordinate (HVC) introduced in WRFv3.9 or terrain following (TF) system.
+  This fix is also compatible with pre WRFv3.9 versions which did not include explicit attribute information for vertical coordinate system.
+- Improved obs_impact_tool documentation.
+
 **March 27 2024 :: WRF-Hydro Developments; AIRS converter documentation update; Add citation.cff file. Tag v11.4.0**
 
 - WRF-Hydro:
+
   - Added a new perfect model obs experimental capability to HydroDART
-  - Modified the Streamflow obs converter to allow for better diagnostics: allows DART to 
+  - Modified the Streamflow obs converter to allow for better diagnostics: allows DART to
     compute obs space diagnostics on all gauges from the Routelink
   - Enhanced performance in the model_mod and noah_hydro_mod when running a full CONUS domain
   - Improved HydroDART Diagnostics with new capabilities (saves the hydrographs in a high-resolution
     pdf, handles hybrid DA components, separate plots for the hybrid statistics, allows the openloop
     to have different ens size and gauges than the DA runs)
+
 - AIRS and AMSU-A observation converters:
+
   - Updated the documentation to use up-to-date build suggestions for the HDFEOS library
   - Updated the AIRS converter code to be able to use version 7 of the AIRS data formats
   - Removed unused and non-functional code: AIRS/BUILD_HDF-EOS.sh, AIRS/L1_AMSUA_to_netcdf.f90,
     AIRS/shell_scripts/Build_HDF_to_netCDF.sh, AIRS/shell_scripts/Convert_HDF_to_netCDF.csh
   - Removed the unnecessary entries from obs_def_rttov_nml in the input.nml
+
 - Added a citation.cff file to help users correctly cite DART software - creates a link to cite
   the repository on the landing page sidebar on GitHub.
 

CITATION.cff

@@ -5,5 +5,5 @@ version: "X.Y.Z"
 date-released: "2024-03-13"
 doi: "10.5065/D6WQ0202"
 authors:
-  - name: "UCAR/NCAR/CISL/DAReS"
+  - name: "UCAR/NSF NCAR/CISL/DAReS"
     city: "Boulder, Colorado"

README.md

@@ -28,7 +28,7 @@ git clone https://github.com/NCAR/DART.git
 
 To cite DART, please use the following text:
 
-> The Data Assimilation Research Testbed (Version X.Y.Z) [Software]. (2021). Boulder, Colorado: UCAR/NCAR/CISL/DAReS.
+> The Data Assimilation Research Testbed (Version X.Y.Z) [Software]. (2021). Boulder, Colorado: UCAR/NSF NCAR/CISL/DAReS.
 > http://doi.org/10.5065/D6WQ0202
 
 and update the DART version and year as appropriate.

assimilation_code/modules/assimilation/assim_tools_mod.rst

@@ -1,3 +1,5 @@
+.. _assim_tools:
+
 MODULE assim_tools_mod
 ======================
 
@@ -289,7 +291,7 @@ Description of each namelist entry
    *type:* character(len=256)
 
    If adjust_obs_impact is true, the name of the file with the observation types and quantities and state quantities
-   that should have have an additional factor applied to the correlations during assimilation.
+   that should have an additional factor applied to the correlations during assimilation.
 
 ``allow_any_impact_values``
    *type:* logical

assimilation_code/programs/obs_impact_tool/obs_impact_tool.rst

@@ -4,47 +4,33 @@ PROGRAM ``obs_impact_tool``
 Overview
 --------
 
-The standard DART algorithms compute increments for an observation and then compute corresponding increments for each
-model state variable due to that observation. To do this, DART computes a sample regression coefficient using the prior
-ensemble distributions of a state variable and the observation. The increments for each member of the observation are
-multiplied by this regression coefficient and then added to the corresponding prior ensemble member for the state
-variable. However, in many cases, it is appropriate to reduce the impact of an observation on a state variable; this is
-called localization. The standard DART algorithms allow users to specify a localization that is a function of the
-horizontal (and optionally vertical) distance between the observation and the state variable. The localization is a
-value between 0 and 1 and multiplies the regression coefficient when updating state ensemble members.
-
-Sometimes, it may be desirable to do an additional localization that is a function of the 
-type of observation and the
-state vector quantity. This program allows users to construct a table that is read by 
-filter at run-time to localize the
-impact of sets of observation types on sets of state vector quantities. Users can create 
-named sets of observation types
-and sets of state vector quantities and specify a localization for the impact of the 
-specified observation types on the state vector quantities.
-
-An example would be to create a subset of observations of tracer concentration for a variety of tracers, and a subset of
-dynamic state variable quantities like temperatures and wind components. It has been common to set this localization
-value to 0 so that tracer observations have no impact on dynamic state quantities, however, the tool allows values
-between 0 and 1 to be specified.
-
-This tool allows related collections of observation types and state vector quantities to be named and then express the
-relationship of the named groups to each other in a concise way. It can also define relationships by exceptions.
-
-All the listed observation types and state vector quantities must be known by the system.
-If they are not, look at the
-&preprocess_nml :: input_items namelist which specifies which *obs_def_xxx_mod.f90* files 
-are included, which is where observation types are defined.
-Quantities for different regimes (atmosphere, ocean, land, etc.) are defined in
-``assimilation_code/modules/observations/xxx_quantities_mod.f90`` and explained in
-:doc:`../../modules/observations/obs_kind_mod`
-
-Format of the input file can be any combination of these types of sections:
+The standard DART algorithms work by calculating increments for an observation and then determining corresponding
+increments for each variable in the state due to that observation. This is done by computing a sample regression 
+coefficient using the prior ensemble distributions of a state variable and the observation. The increments for each member 
+of the ensemble are multiplied by this coefficient and then added to the corresponding prior ensemble member for the variable.
 
-.. container::
+However, in many cases it is necessary to limit the influence of an observation on a variable; this is known as localization.
+DART provides a way to specify a localization, known as cutoff, based on the horizontal and vertical distance between the observation 
+and the state variable.
+
+In some situations, you may want additional localization based on the type of observation and the state quanity. 
+``obs_impact_tool`` allows you to create a table that filter reads during runtime to localize the impact of certain types of 
+observations on specific state vector quantities. You can define sets of observation types and state vector quantities, and 
+specify localization for the impact of those observation types on the state vector quantities.
+
+For example, you can create a subset of observations related to tracer concentration for various tracers, and a subset of 
+dynamic state variables like temperatures and wind components. Typically, it is common practice to set this localization value 
+to 0 to prevent tracer observations from affecting dynamic state quantities. However, ``obs_impact_tool`` allows you to specify values
+between 0 and 1.
 
-   ::
 
+#. Build ``obs_sequence_tool`` by adding ``obs_impact_tool`` to the list of serial_programs in the quickbuild.sh script for the model you are using.
+   Run ./quickbuild.sh to build all the DART programs.
+#. Create an input file for ``obs_sequence_tool`` to define the impacts of observations. In the examples on this page, the input file
+   is called `cross_correlations.txt`.  
+   The format of the input file can be any combination of the following types of sections:
 
+   .. code:: bash
 
       # hash mark starts a comment.
 
@@ -99,10 +85,48 @@ Format of the input file can be any combination of these types of sections:
        groupname1  groupname1   0.0
       END IMPACT
 
-Namelist interface ``&obs_impact_tool_nml`` must be read from file ``input.nml``.
+   The following is an example of an input file to prevent chemistry species from impacting the meterological variables in the model state, and vice versa:
 
-Namelist
---------
+   .. code:: bash
+
+       GROUP chem
+        QTY_CO QTY_NO QTY_C2H4
+       END GROUP
+
+       GROUP met
+        ALLQTYS EXCEPT chem
+       END GROUP
+
+       IMPACT
+       chem   met    0.0
+       met    chem   0.0
+       END IMPACT
+
+
+#. Run ``obs_impact_tool`` using your `cross_correlations.txt` as input. ``obs_impact_tool`` will create an output file,
+   named `control_impact_runtime.txt` in this example.
+
+   .. code:: text
+   
+      &obs_impact_tool_nml
+        input_filename          = 'cross_correlations.txt'
+        output_filename         = 'control_impact_runtime.txt'
+        /
+   
+
+#. Set the following namelist options in :ref:`&assim_tools_nml<assim_tools>` to use `control_impact_runtime.txt` in filter. 
+   Filter will apply your selected observation impacts during assimilation.
+
+   .. code:: text
+   
+      &assim_tools_nml
+        adjust_obs_impact               = .true.
+        obs_impact_filename             = 'control_impact_runtime.txt'
+        /
+
+
+obs_impact_tool Namelist
+------------------------
 
 This namelist is read from the file ``input.nml``. Namelists start with an ampersand '&' and terminate with a slash '/'.
 Character strings that contain a '/' must be enclosed in quotes to prevent them from prematurely terminating the
@@ -135,47 +159,3 @@ namelist.
    +-----------------+--------------------+-----------------------------------------------------------------------------+
    | debug           | logical            | If true print out debugging info.                                           |
    +-----------------+--------------------+-----------------------------------------------------------------------------+
-
-| 
-
-Examples
---------
-
-To prevent chemistry species from impacting the meterological variables in the model state, and vice versa:
-
-.. container::
-
-   ::
-
-      GROUP chem
-       QTY_CO QTY_NO QTY_C2H4
-      END GROUP
-
-      GROUP met
-       ALLQTYS EXCEPT chem
-      END GROUP
-
-      IMPACT
-       chem   met    0.0
-       met    chem   0.0
-      END IMPACT
-
-Modules used
-------------
-
-::
-
-   types_mod
-   utilities_mod
-   parse_args_mod
-
-Files
------
-
--  two text files, one input and one output.
--  obs_impact_tool.nml
-
-References
-----------
-
--  none

conf.py

@@ -21,7 +21,7 @@
 author = 'Data Assimilation Research Section'
 
 # The full version, including alpha/beta/rc tags
-release = '11.4.0'
+release = '11.4.1'
 root_doc = 'index'
 
 # -- General configuration ---------------------------------------------------

index.rst

@@ -246,8 +246,8 @@ Citing DART
 
 Cite DART using the following text:
 
-   The Data Assimilation Research Testbed (Version X.Y.Z) [Software]. (2019).
-   Boulder, Colorado: UCAR/NCAR/CISL/DAReS. http://doi.org/10.5065/D6WQ0202
+   The Data Assimilation Research Testbed (Version X.Y.Z) [Software]. (2024).
+   Boulder, Colorado: UCAR/NSF NCAR/CISL/DAReS. http://doi.org/10.5065/D6WQ0202
 
 Update the DART version and year as appropriate.
 

models/wrf/model_mod.f90

@@ -284,6 +284,10 @@ module model_mod
 real (kind=r8), PARAMETER    :: kappa = 2.0_r8/7.0_r8 ! gas_constant / cp
 real (kind=r8), PARAMETER    :: ts0 = 300.0_r8        ! Base potential temperature for all levels.
 
+! hybrid_opt 2 is a hybrid coordinate - terrain following at the
+! surface, and straight pressure levels at the top.
+integer, parameter :: VERT_HYBRID = 2
+
 !---- private data ----
 
 ! Got rid of surf_var as global private variable for model_mod and just defined it locally
@@ -309,7 +313,8 @@ module model_mod
 
    integer  :: domain_size
    integer  :: localization_coord
-   real(r8), dimension(:),     pointer :: znu, dn, dnw, zs, znw
+   integer  :: hybrid_opt
+   real(r8), dimension(:),     pointer :: znu, dn, dnw, zs, znw, c1h, c2h
    real(r8), dimension(:,:),   pointer :: mub, hgt
    real(r8), dimension(:,:),   pointer :: latitude, latitude_u, latitude_v
    real(r8), dimension(:,:),   pointer :: longitude, longitude_u, longitude_v
@@ -3766,7 +3771,7 @@ subroutine nc_write_model_atts( ncid, id )
 
 integer :: DNVarID, ZNUVarID, DNWVarID, phbVarID, &
      MubVarID, LonVarID, LatVarID, ilevVarID, XlandVarID, hgtVarID , LatuVarID, &
-     LatvVarID, LonuVarID, LonvVarID, ZNWVarID
+     LatvVarID, LonuVarID, LonvVarID, ZNWVarID, C1HVarID, C2HVarID
 
 integer :: TimeDimID
 
@@ -4042,6 +4047,34 @@ subroutine nc_write_model_atts( ncid, id )
                  ''), &
                  'nc_write_model_atts','def_var DNW'//' units')
 
+if (wrf%dom(id)%hybrid_opt == VERT_HYBRID) then
+   call nc_check(nf90_def_var(ncid, name='C1H', xtype=nf90_real, &
+                    dimids= btDimID, varid=C1HVarID), &
+                    'nc_write_model_atts','def_var C1H')
+   call nc_check(nf90_put_att(ncid, C1HVarID, 'long_name', &
+                    'dn values on full (w) levels'), &
+                    'nc_write_model_atts','def_var C1H'//' long_name')
+   call nc_check(nf90_put_att(ncid, C1HVarID, 'description', &
+                    'dn values on full (w) levels'), &
+                    'nc_write_model_atts','def_var C1H'//' description')
+   call nc_check(nf90_put_att(ncid, C1HVarID, 'units', &
+                    ''), &
+                    'nc_write_model_atts','def_var C1H'//' units')
+
+   call nc_check(nf90_def_var(ncid, name='C2H', xtype=nf90_real, &
+                    dimids= btDimID, varid=C2HVarID), &
+                    'nc_write_model_atts','def_var C2H')
+   call nc_check(nf90_put_att(ncid, C2HVarID, 'long_name', &
+                    'dn values on full (w) levels'), &
+                    'nc_write_model_atts','def_var C2H'//' long_name')
+   call nc_check(nf90_put_att(ncid, C2HVarID, 'description', &
+                    'dn values on full (w) levels'), &
+                    'nc_write_model_atts','def_var C2H'//' description')
+   call nc_check(nf90_put_att(ncid, C2HVarID, 'units', &
+                    ''), &
+                    'nc_write_model_atts','def_var C2H'//' units')
+endif
+
 !
 !    float MUB(Time, south_north, west_east) ;
 !            MUB:FieldType = 104 ;
@@ -4307,6 +4340,12 @@ subroutine nc_write_model_atts( ncid, id )
               'nc_write_model_atts','put_var znw')
 call nc_check(nf90_put_var(ncid,      DNWVarID, wrf%dom(id)%dnw), &
               'nc_write_model_atts','put_var dnw')
+if (wrf%dom(id)%hybrid_opt == VERT_HYBRID) then
+   call nc_check(nf90_put_var(ncid,      C1HVarID, wrf%dom(id)%c1h), &
+                 'nc_write_model_atts','put_var c1h')
+   call nc_check(nf90_put_var(ncid,      C2HVarID, wrf%dom(id)%c2h), &
+                 'nc_write_model_atts','put_var c2h')
+endif
 
 ! defining horizontal
 call nc_check(nf90_put_var(ncid,      mubVarID, wrf%dom(id)%mub), &
@@ -5276,7 +5315,11 @@ function model_rho_t_distrib(i,j,k,id,state_handle, ens_size)
 
 ! now calculate rho = - mu / dphi/deta
 
-model_rho_t_distrib(:) = - (wrf%dom(id)%mub(i,j)+x_imu) / ph_e
+if (wrf%dom(id)%hybrid_opt == VERT_HYBRID) then
+   model_rho_t_distrib(:) = - (wrf%dom(id)%c1h(k)*(wrf%dom(id)%mub(i,j)+x_imu) + wrf%dom(id)%c2h(k)) / ph_e
+else
+   model_rho_t_distrib(:) = - (wrf%dom(id)%mub(i,j)+x_imu) / ph_e
+endif
 
 end function model_rho_t_distrib
 
@@ -7044,6 +7087,7 @@ subroutine read_wrf_file_attributes(ncid,id)
 
 integer, intent(in)   :: ncid, id
 logical, parameter    :: debug = .false.
+integer :: ret
 
 ! get meta data and static data we need
 
@@ -7082,6 +7126,13 @@ subroutine read_wrf_file_attributes(ncid,id)
                      'static_init_model', 'get_att STAND_LON')
    if(debug) write(*,*) ' stdlon is ',stdlon
 
+   ! this attribute is not present in older wrf files, which means it is not
+   ! using a hybrid vertical coordinate.  not having this attribute should
+   ! not cause an error.
+   ret =  nf90_get_att(ncid, nf90_global, 'HYBRID_OPT', wrf%dom(id)%hybrid_opt)
+   if (ret /= NF90_NOERR) wrf%dom(id)%hybrid_opt = 0
+   if(debug) write(*,*) 'hybrid_opt is ', wrf%dom(id)%hybrid_opt
+
    RETURN
 
 end subroutine read_wrf_file_attributes
@@ -7176,6 +7227,22 @@ subroutine read_wrf_static_data(ncid,id)
                      'read_wrf_static_data','get_var DNW')
    if(debug) write(*,*) ' dnw is ',wrf%dom(id)%dnw
 
+   if (wrf%dom(id)%hybrid_opt == VERT_HYBRID) then
+      allocate(wrf%dom(id)%c1h(1:wrf%dom(id)%bt))
+      call nc_check( nf90_inq_varid(ncid, "C1H", var_id), &
+                        'read_wrf_static_data','inq_varid C1H')
+      call nc_check( nf90_get_var(ncid, var_id, wrf%dom(id)%c1h), &
+                        'read_wrf_static_data','get_var C1H')
+      if(debug) write(*,*) ' c1h is ',wrf%dom(id)%c1h
+
+      allocate(wrf%dom(id)%c2h(1:wrf%dom(id)%bt))
+      call nc_check( nf90_inq_varid(ncid, "C2H", var_id), &
+                        'read_wrf_static_data','inq_varid C2H')
+      call nc_check( nf90_get_var(ncid, var_id, wrf%dom(id)%c2h), &
+                        'read_wrf_static_data','get_var C2H')
+      if(debug) write(*,*) ' c2h is ',wrf%dom(id)%c2h
+   endif
+
    allocate(wrf%dom(id)%zs(1:wrf%dom(id)%sls))
    call nc_check( nf90_inq_varid(ncid, "ZS", var_id), &
                      'read_wrf_static_data','inq_varid ZS')