diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 5e2c1706f..db39b1910 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,15 @@ individual files. The changes are now listed with the most recent at the top. +**August 15 2024 :: WRF fwd operator bug fixes. Tag v11.6.1** + +WRF-DART bug-fixes: + + - Bug fix for surface temperature observations to use QTY_2M_TEMPERATURE + - Bug fix for conversion of vapor mixing ratio to specific humidity + - Bug fix for diagnostics_obs.csh + - Improved documentation for WRF model_mod and WRF-DART Tutorial + **July 26 2024 :: Library build tools for DART. Tag v11.6.0** - Buildtools for compiling DART as a shared or a static library. diff --git a/conf.py b/conf.py index deb053a0b..fdfe8997c 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '11.6.0' +release = '11.6.1' root_doc = 'index' # -- General configuration --------------------------------------------------- diff --git a/models/wrf/model_mod.f90 b/models/wrf/model_mod.f90 index b29c0dc79..9ca5ddbe7 100644 --- a/models/wrf/model_mod.f90 +++ b/models/wrf/model_mod.f90 @@ -1648,86 +1648,96 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte elseif ( obs_kind == QTY_TEMPERATURE ) then ! This is for 3D temperature field -- surface temps later - !print*, 'k ', k + if(.not. surf_var) then - if ( wrf%dom(id)%type_t >= 0 ) then + if ( wrf%dom(id)%type_t >= 0 ) then - do uk = 1, count ! for the different ks + do uk = 1, count ! for the different ks - ! Check to make sure retrieved integer gridpoints are in valid range - if ( boundsCheck( i, wrf%dom(id)%periodic_x, id, dim=1, type=wrf%dom(id)%type_t ) .and. & - boundsCheck( j, wrf%dom(id)%polar, id, dim=2, type=wrf%dom(id)%type_t ) .and. & - boundsCheck( uniquek(uk), .false., id, dim=3, type=wrf%dom(id)%type_t ) ) then + ! Check to make sure retrieved integer gridpoints are in valid range + if ( boundsCheck( i, wrf%dom(id)%periodic_x, id, dim=1, type=wrf%dom(id)%type_t ) .and. & + boundsCheck( j, wrf%dom(id)%polar, id, dim=2, type=wrf%dom(id)%type_t ) .and. & + boundsCheck( uniquek(uk), .false., id, dim=3, type=wrf%dom(id)%type_t ) ) then - call getCorners(i, j, id, wrf%dom(id)%type_t, ll, ul, lr, ur, rc ) - if ( rc .ne. 0 ) & - print*, 'model_mod.f90 :: model_interpolate :: getCorners T rc = ', rc + call getCorners(i, j, id, wrf%dom(id)%type_t, ll, ul, lr, ur, rc ) + if ( rc .ne. 0 ) & + print*, 'model_mod.f90 :: model_interpolate :: getCorners T rc = ', rc - ! Interpolation for T field at level k - ill = get_dart_vector_index(ll(1), ll(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_t) - iul = get_dart_vector_index(ul(1), ul(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_t) - ilr = get_dart_vector_index(lr(1), lr(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_t) - iur = get_dart_vector_index(ur(1), ur(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_t) - - x_iul = get_state(iul, state_handle) - x_ill = get_state(ill, state_handle) - x_ilr = get_state(ilr, state_handle) - x_iur = get_state(iur, state_handle) - - ! In terms of perturbation potential temperature - a1 = dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) - - pres1 = model_pressure_t_distrib(ll(1), ll(2), uniquek(uk), id, state_handle, ens_size) - pres2 = model_pressure_t_distrib(lr(1), lr(2), uniquek(uk), id, state_handle, ens_size) - pres3 = model_pressure_t_distrib(ul(1), ul(2), uniquek(uk), id, state_handle, ens_size) - pres4 = model_pressure_t_distrib(ur(1), ur(2), uniquek(uk), id, state_handle, ens_size) - - ! Pressure at location - pres = dym*( dxm*pres1 + dx*pres2 ) + dy*( dxm*pres3 + dx*pres4 ) - - do e = 1, ens_size - if ( k(e) == uniquek(uk) ) then - ! Full sensible temperature field - fld(1, e) = (ts0 + a1(e))*(pres(e)/ps0)**kappa - endif - enddo - - ! Interpolation for T field at level k+1 - ill = get_dart_vector_index(ll(1), ll(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_t) - iul = get_dart_vector_index(ul(1), ul(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_t) - ilr = get_dart_vector_index(lr(1), lr(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_t) - iur = get_dart_vector_index(ur(1), ur(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_t) - - x_ill = get_state(ill, state_handle) - x_iul = get_state(iul, state_handle) - x_iur = get_state(iur, state_handle) - x_ilr = get_state(ilr, state_handle) - - ! In terms of perturbation potential temperature - a1 = dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) - - pres1 = model_pressure_t_distrib(ll(1), ll(2), uniquek(uk)+1, id, state_handle, ens_size) - pres2 = model_pressure_t_distrib(lr(1), lr(2), uniquek(uk)+1, id, state_handle, ens_size) - pres3 = model_pressure_t_distrib(ul(1), ul(2), uniquek(uk)+1, id, state_handle, ens_size) - pres4 = model_pressure_t_distrib(ur(1), ur(2), uniquek(uk)+1, id, state_handle, ens_size) - - ! Pressure at location - pres = dym*( dxm*pres1 + dx*pres2 ) + dy*( dxm*pres3 + dx*pres4 ) - - do e = 1, ens_size - if ( k(e) == uniquek(uk) ) then - ! Full sensible temperature field - fld(2, e) = (ts0 + a1(e))*(pres(e)/ps0)**kappa - endif - enddo - endif - enddo + ! Interpolation for T field at level k + ill = get_dart_vector_index(ll(1), ll(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_t) + iul = get_dart_vector_index(ul(1), ul(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_t) + ilr = get_dart_vector_index(lr(1), lr(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_t) + iur = get_dart_vector_index(ur(1), ur(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_t) + + x_iul = get_state(iul, state_handle) + x_ill = get_state(ill, state_handle) + x_ilr = get_state(ilr, state_handle) + x_iur = get_state(iur, state_handle) + + ! In terms of perturbation potential temperature + a1 = dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) + + pres1 = model_pressure_t_distrib(ll(1), ll(2), uniquek(uk), id, state_handle, ens_size) + pres2 = model_pressure_t_distrib(lr(1), lr(2), uniquek(uk), id, state_handle, ens_size) + pres3 = model_pressure_t_distrib(ul(1), ul(2), uniquek(uk), id, state_handle, ens_size) + pres4 = model_pressure_t_distrib(ur(1), ur(2), uniquek(uk), id, state_handle, ens_size) + + ! Pressure at location + pres = dym*( dxm*pres1 + dx*pres2 ) + dy*( dxm*pres3 + dx*pres4 ) + + do e = 1, ens_size + if ( k(e) == uniquek(uk) ) then + ! Full sensible temperature field + fld(1, e) = (ts0 + a1(e))*(pres(e)/ps0)**kappa + endif + enddo + + ! Interpolation for T field at level k+1 + ill = get_dart_vector_index(ll(1), ll(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_t) + iul = get_dart_vector_index(ul(1), ul(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_t) + ilr = get_dart_vector_index(lr(1), lr(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_t) + iur = get_dart_vector_index(ur(1), ur(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_t) + + x_ill = get_state(ill, state_handle) + x_iul = get_state(iul, state_handle) + x_iur = get_state(iur, state_handle) + x_ilr = get_state(ilr, state_handle) + + ! In terms of perturbation potential temperature + a1 = dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) + + pres1 = model_pressure_t_distrib(ll(1), ll(2), uniquek(uk)+1, id, state_handle, ens_size) + pres2 = model_pressure_t_distrib(lr(1), lr(2), uniquek(uk)+1, id, state_handle, ens_size) + pres3 = model_pressure_t_distrib(ul(1), ul(2), uniquek(uk)+1, id, state_handle, ens_size) + pres4 = model_pressure_t_distrib(ur(1), ur(2), uniquek(uk)+1, id, state_handle, ens_size) + + ! Pressure at location + pres = dym*( dxm*pres1 + dx*pres2 ) + dy*( dxm*pres3 + dx*pres4 ) + + do e = 1, ens_size + if ( k(e) == uniquek(uk) ) then + ! Full sensible temperature field + fld(2, e) = (ts0 + a1(e))*(pres(e)/ps0)**kappa + endif + enddo + endif + enddo + endif else - fld = missing_r8 - end if + + ! This is for surface temperature (T2) + if ( wrf%dom(id)%type_t2 >= 0 ) then + call surface_interp_distrib(fld, wrf, id, i, j, obs_kind, wrf%dom(id)%type_t2, dxm, dx, dy, dym, ens_size, state_handle) + if (all(fld == missing_r8)) goto 200 + else + call error_handler(E_MSG, 'model_mod section 1.b Sensible Temperature:', & + 'WARNING: Surface temperature variable not found in &model_nml') + fld = missing_r8 + end if + end if elseif (obs_kind == QTY_2M_TEMPERATURE) then ! This is for 2-meter temperature - if ( wrf%dom(id)%type_t2 >= 0 ) then ! HK is there a better way to do this? + if ( wrf%dom(id)%type_t2 >= 0 ) then call surface_interp_distrib(fld, wrf, id, i, j, obs_kind, wrf%dom(id)%type_t2, dxm, dx, dy, dym, ens_size, state_handle) if (all(fld == missing_r8)) goto 200 else @@ -1801,7 +1811,7 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte call surface_interp_distrib(fld, wrf, id, i, j, obs_kind, wrf%dom(id)%type_th2, dxm, dx, dy, dym, ens_size, state_handle) if (all(fld == missing_r8)) goto 200 - endif + endif endif !----------------------------------------------------- @@ -1890,9 +1900,8 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte if (all(fld == missing_r8)) goto 200 !----------------------------------------------------- - ! 1.f Specific Humidity (SH, SH2) - ! Look at me - ! Convert water vapor mixing ratio to specific humidity: + ! 1.f Specific Humidity (QV,Q2) + ! Water vapor mixing ratio is used to calculate specific humidity: else if( obs_kind == QTY_SPECIFIC_HUMIDITY ) then ! This is for 3D specific humidity -- surface spec humidity later @@ -1910,9 +1919,9 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte call getCorners(i, j, id, wrf%dom(id)%type_t, ll, ul, lr, ur, rc ) ! HK why is this type_t if ( rc .ne. 0 ) & - print*, 'model_mod.f90 :: model_interpolate :: getCorners SH rc = ', rc + print*, 'model_mod.f90 :: model_interpolate :: getCorners QV rc = ', rc - ! Interpolation for SH field at level k + ! Interpolation for QV field at level k ill = get_dart_vector_index(ll(1), ll(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_qv) iul = get_dart_vector_index(ul(1), ul(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_qv) ilr = get_dart_vector_index(lr(1), lr(2), uniquek(uk), domain_id(id), wrf%dom(id)%type_qv) @@ -1926,11 +1935,13 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte do e = 1, ens_size if ( k(e) == uniquek(uk) ) then ! interpolate only if it is the correct k a1 = dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) ! think this can go outside the k loop + + ! Vapor mixing ratio (QV) to specific humidity conversion fld(1,e) = a1(e) /(1.0_r8 + a1(e)) endif enddo - ! Interpolation for SH field at level k+1 + ! Interpolation for QV field at level k+1 ill = get_dart_vector_index(ll(1), ll(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_qv) iul = get_dart_vector_index(ul(1), ul(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_qv) ilr = get_dart_vector_index(lr(1), lr(2), uniquek(uk)+1, domain_id(id), wrf%dom(id)%type_qv) @@ -1942,8 +1953,10 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte x_iur = get_state(iur, state_handle) do e = 1, ens_size - if ( k(e) == uniquek(uk) ) then ! interpolate only if it is the correct + if ( k(e) == uniquek(uk) ) then ! interpolate only if it is the correct k a1 = dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) + + ! Vapor mixing ratio (QV) to specific humidity conversion fld(2,e) = a1(e) /(1.0_r8 + a1(e)) endif enddo @@ -1952,7 +1965,7 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte endif - ! This is for surface specific humidity (calculated from Q2) + ! This is for surface specific humidity calculated from vapor mixing ratio (Q2) else ! confirm that field is in the DART state vector @@ -1966,7 +1979,7 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte if ( rc .ne. 0 ) & print*, 'model_mod.f90 :: model_interpolate :: getCorners Q2 rc = ', rc - ! Interpolation for the SH2 field + ! Interpolation for the Q2 field ill = get_dart_vector_index(ll(1), ll(2), 1, domain_id(id), wrf%dom(id)%type_q2) iul = get_dart_vector_index(ul(1), ul(2), 1, domain_id(id), wrf%dom(id)%type_q2) ilr = get_dart_vector_index(lr(1), lr(2), 1, domain_id(id), wrf%dom(id)%type_q2) @@ -1978,6 +1991,8 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte x_ilr = get_state(ilr, state_handle) a1 = dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) + + ! Vapor mixing ratio (Q2) to specific humidity conversion fld(1,:) = a1 /(1.0_r8 + a1) endif @@ -2002,11 +2017,10 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte endif endif - ! Don't accept negative water vapor amounts (?) + ! Don't accept negative water vapor amounts fld = max(0.0_r8, fld) else if (obs_kind == QTY_2M_SPECIFIC_HUMIDITY ) then - ! FIXME: Q2 is actually a mixing ratio, not a specific humidity if ( wrf%dom(id)%type_q2 >= 0 ) then ! Check to make sure retrieved integer gridpoints are in valid range if ( ( boundsCheck( i, wrf%dom(id)%periodic_x, id, dim=1, type=wrf%dom(id)%type_t ) .and. & @@ -2017,7 +2031,7 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte if ( rc .ne. 0 ) & print*, 'model_mod.f90 :: model_interpolate :: getCorners Q2 rc = ', rc - ! Interpolation for the SH2 field + ! Interpolation for the Q2 field ill = get_dart_vector_index(ll(1), ll(2), 1, domain_id(id), wrf%dom(id)%type_q2) iul = get_dart_vector_index(ul(1), ul(2), 1, domain_id(id), wrf%dom(id)%type_q2) ilr = get_dart_vector_index(lr(1), lr(2), 1, domain_id(id), wrf%dom(id)%type_q2) @@ -2029,7 +2043,10 @@ subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expecte x_ilr = get_state(ilr, state_handle) a1 = dym*( dxm*x_ill + dx*x_ilr ) + dy*( dxm*x_iul + dx*x_iur ) - fld(1,:) = a1 + + ! Vapor mixing ratio (Q2) to specific humidity conversion + fld(1,:) = a1/(1.0_r8 + a1) + endif endif diff --git a/models/wrf/readme.rst b/models/wrf/readme.rst index 0f37ad09f..ebb322751 100644 --- a/models/wrf/readme.rst +++ b/models/wrf/readme.rst @@ -4,34 +4,42 @@ WRF Overview -------- - -DART interface module for the Weather Research and Forecasting +The following is a description of the DART interface module for the +Weather Research and Forecasting model `(WRF) `__ -model. This page documents the details of the -module compiled into DART that interfaces with the WRF data in the state vector. +model. This page provides an overview of the module compiled into DART +that interfaces with the WRF data in the state vector. **The WRF-DART interface is compatible with WRF versions 4 and later, and is no longer backwards compatible with WRFv3.9 and earlier.** For more information on the interface changes required between -different WRF versions see the WRF tutorial link in the next section. +different WRF versions, read through this documentation *and* the +WRF-DART tutorial link in the next section. -WRF+DART Tutorial ------------------ +There have been several important updates to the WRF-DART interface starting +with `DARTv11.5.0. `__ +Some important WRF-DART updates include: + +- Version 11.4.1: Detects use of the Hybrid Vertical Coordinate system + (terrain following at surface) and accounts for this in the forward + operator calculations. + +- Version 11.5.0: Improves compatibility with WRFv4+ versions where + the prognostic 3D temperature variable is THM. + +It is always recommended that you update your DART version to the +`latest release `__ before beginning new research. -**There is additional overview and tutorial documentation for running a WRF/DART -assimilation in** :doc:`./tutorial/README` +WRF-DART Tutorial +----------------- -Please work through the tutorial in order to learn how to run WRF and DART. +This tutorial provides a real-world example of assimilating a wide variety of atmospheric +observations during an extreme storm event for the United States during April 2017. +**It is strongly recommended that you also review and perform the tutorial for +running a WRF-DART assimilation** `here. `__ -Items of Note -~~~~~~~~~~~~~ -- The ``model_mod`` reads WRF netCDF files directly to acquire the model state - data. The ``wrf_to_dart`` and ``dart_to_wrf`` programs are no longer - necessary. -- A netCDF file named ``wrfinput_d01`` is required and must be at the same - resolution and have the same surface elevation data as the files converted to - create the DART initial conditions. No data will be read from this file, but - the grid information must match exactly. +WRF Interface Overview +---------------------- The model interface code supports WRF configurations with multiple domains. Data for all domains is read into the DART state vector. During the computation of @@ -44,21 +52,31 @@ domain case the data in the state vector that came from ``wrfinput_d04`` is searched first, then ``wrfinput_d03``, ``wrfinput_d02``, and finally ``wrfinput_d01``. -The forward operator is computed from the first domain grid that contains the -lat/lon of the observation. During the assimilation phase, when the state values -are adjusted based on the correlations and assimilation increments, all points -in all domains that are within the localization radius are adjusted, regardless -of domain. The impact of an observation on the state depends only on the -distance between the observation and the state vector point, and the regression -coefficient based on the correlation between the distributions of the ensemble -of state vector points and the ensemble of observation forward operator values. +.. Important:: + + Although the model interface code is compatible with multiple domains, the + supporting `shell scripts `__ + and WRF-DART tutorial are currently designed for a single domain and will + require some modifications for multiple (nested) domain functionality. If you + need help with these modifications please contact DART support. + + +In summary, the forward operator is computed from the first domain grid (searching from +finest grid to coarsest grid) that contains the lat/lon of the observation. During the +assimilation phase, however, when the state values are adjusted based on the correlations +and assimilation increments, all points in all domains that are within the +localization radius are adjusted, regardless of domain. The impact of an observation +on the state depends only on the distance between the observation and the state +vector point, and the regression coefficient based on the correlation between the +distributions of the ensemble of state vector points and the ensemble of observation +forward operator values. The fields from WRF that are copied into the DART state vector are controlled by -namelist. See below for the documentation on the &model_nml entries. The state -vector should include all fields needed to restart a WRF run. There may be -additional fields needed depending on the microphysics scheme selected. See the -ascii file ``wrf_state_variables_table`` in the ``models/wrf`` directory for a -list of fields that are often included in the DART state. +the namelist within ``input.nml``. See below for the documentation on the ``&model_nml`` entries within +``input.nml``. The state vector should include all fields needed to restart a WRF run. +There may be additional fields needed depending on the microphysics scheme selected. See the +ascii file `wrf_state_variables_table `__ +that includes a list of fields that may be added to the DART state. Namelist -------- @@ -68,7 +86,7 @@ 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 namelist. -.. code-block:: +:: &model_nml default_state_variables = .true. @@ -89,279 +107,290 @@ prematurely terminating the namelist. periodic_y = .false. scm = .false. allow_perturbed_ics = .false. # testing purposes only - / - - # Notes for model_nml: - # (1) vert_localization_coord must be one of: - # 1 = model level - # 2 = pressure - # 3 = height - # 4 = scale height - # (2) see bottom of this file for explanations of polar, periodic_x, - # periodic_y, and scm - # (3) calendar = 3 is GREGORIAN, which is what WRF uses. - # (4) if 'default_state_variables' is .true. the model_mod.f90 code will - # fill the state variable table with the following wrf vars: - # U, V, W, PH, T, MU - # you must set it to false before you change the value - # of 'wrf_state_variables' and have it take effect. - # (5) the format for 'wrf_state_variables' is an array of 5 strings: - # wrf netcdf variable name, dart QTY_xxx string, type string (must be - # unique, will soon be obsolete, we hope), 'UPDATE', and '999' if the - # array is part of all domains. otherwise, it is a string with the domain - # numbers (e.g. '12' for domains 1 and 2, '13' for domains 1 and 3). - # example: - # wrf_state_variables='U','QTY_U_WIND_COMPONENT','TYPE_U','UPDATE','999', - # 'V','QTY_V_WIND_COMPONENT','TYPE_V','UPDATE','999', - # 'W','QTY_VERTICAL_VELOCITY','TYPE_W','UPDATE','999', - # 'T','QTY_POTENTIAL_TEMPERATURE','TYPE_T','UPDATE','999', - # 'PH','QTY_GEOPOTENTIAL_HEIGHT','TYPE_GZ','UPDATE','999', - # 'MU','QTY_PRESSURE','TYPE_MU','UPDATE','999', - # 'QVAPOR','QTY_VAPOR_MIXING_RATIO','TYPE_QV','UPDATE','999', - # 'QCLOUD','QTY_CLOUD_LIQUID_WATER','TYPE_QC','UPDATE','999', - # 'QRAIN','QTY_RAINWATER_MIXING_RATIO','TYPE_QR','UPDATE','999', - # 'U10','QTY_U_WIND_COMPONENT','TYPE_U10','UPDATE','999', - # 'V10','QTY_V_WIND_COMPONENT','TYPE_V10','UPDATE','999', - # 'T2','QTY_TEMPERATURE','TYPE_T2','UPDATE','999', - # 'TH2','QTY_POTENTIAL_TEMPERATURE','TYPE_TH2','UPDATE','999', - # 'Q2','QTY_SPECIFIC_HUMIDITY','TYPE_Q2','UPDATE','999', - # 'PSFC','QTY_PRESSURE','TYPE_PS','UPDATE','999', - # (6) the format for 'wrf_state_bounds' is an array of 4 strings: - # wrf netcdf variable name, minimum value, maximum value, and either - # FAIL or CLAMP. FAIL will halt the program if an out of range value - # is detected. CLAMP will set out of range values to the min or max. - # The special string 'NULL' will map to plus or minus infinity and will - # not change the values. arrays not listed in this table will not - # be changed as they are read or written. - # - # - # polar and periodic_x are used in global wrf. if polar is true, the - # grid interpolation routines will wrap over the north and south poles. - # if periodic_x is true, when the east and west edges of the grid are - # reached the interpolation will wrap. note this is a separate issue - # from regional models which cross the GMT line; those grids are marked - # as having a negative offset and do not need to wrap; this flag controls - # what happens when the edges of the grid are reached. - - # the scm flag is used for the 'single column model' version of WRF. - # it needs the periodic_x and periodic_y flags set to true, in which - # case the X and Y directions are periodic; no collapsing of the grid - # into a single location like the 3d-spherical polar flag implies. - -Description of each namelist entry -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -+---------------------------------------+-------------------+---------------------------------------+ -| Item | Type | Description | -+=======================================+===================+=======================================+ -| default_state_variables | logical | If *.true.*, the dart state vector | -| | | contains the fields U, V, W, PH, T, | -| | | MU, in that order, and only those. | -| | | Any values listed in the | -| | | *wrf_state_variables* namelist item | -| | | will be ignored. | -+---------------------------------------+-------------------+---------------------------------------+ -| wrf_state_variables | character(:, 5) | A 2D array of strings, 5 per wrf | -| | | array to be added to the dart state | -| | | vector. If *default_state_variables* | -| | | is *.true.*, this is ignored. When | -| | | *.false.*, this list of array names | -| | | controls which arrays and the order | -| | | that they are added to the state | -| | | vector. The 5 strings are: | -| | | | -| | | #. WRF field name - must match netcdf | -| | | name exactly | -| | | #. DART KIND name - must match a | -| | | valid DART QTY_xxx exactly | -| | | #. TYPE_NN - will hopefully be | -| | | obsolete, but for now NN should | -| | | match the field name. | -| | | #. the string UPDATE. at some future | -| | | point, non-updatable fields may | -| | | become part of the state vector. | -| | | #. A numeric string listing the | -| | | domain numbers this array is part | -| | | of. The specical string 999 means | -| | | all domains. For example, '12' | -| | | means domains 1 and 2, '13' means | -| | | 1 and 3. | -+---------------------------------------+-------------------+---------------------------------------+ -| wrf_state_bounds | character(:, 4) | A 2D array of strings, 4 per wrf | -| | | array. During the copy of data to and | -| | | from the wrf netcdf file, variables | -| | | listed here will have minimum and | -| | | maximum values enforced. The 4 | -| | | strings are: | -| | | | -| | | #. WRF field name - must match netcdf | -| | | name exactly | -| | | #. Minimum -- specified as a string | -| | | but must be a numeric value (e.g. | -| | | '0.1') Can be 'NULL' to allow any | -| | | minimum value. | -| | | #. Maximum -- specified as a string | -| | | but must be a numeric value (e.g. | -| | | '0.1') Can be 'NULL' to allow any | -| | | maximum value. | -| | | #. Action -- valid strings are | -| | | 'CLAMP', 'FAIL'. 'FAIL' means if a | -| | | value is found outside the range, | -| | | the code fails with an error. | -| | | 'CLAMP' simply sets the out of | -| | | range values to the given minimum | -| | | or maximum without error. | -+---------------------------------------+-------------------+---------------------------------------+ -| num_domains | integer | Total number of WRF domains, | -| | | including nested domains. | -+---------------------------------------+-------------------+---------------------------------------+ -| calendar_type | integer | Calendar type. Should be 3 | -| | | (GREGORIAN) for WRF. | -+---------------------------------------+-------------------+---------------------------------------+ -| assimilation_period_seconds | integer | The time (in seconds) between | -| | | assimilations. This is modified if | -| | | necessary to be an integer multiple | -| | | of the underlying model timestep. | -+---------------------------------------+-------------------+---------------------------------------+ -| periodic_x | logical | If *.true.*, the grid is periodic in | -| | | longitude, and points above the last | -| | | grid cell and points below the first | -| | | grid cell are wrapped. Note this is | -| | | not the same as a grid which crosses | -| | | the prime meridian. WRF handles that | -| | | with an offset in longitude and | -| | | points beyond the last grid index are | -| | | outside the domain. | -+---------------------------------------+-------------------+---------------------------------------+ -| periodic_y | logical | Used for the Single Column Model to | -| | | make the grid wrap in Y (see scm | -| | | below). This is NOT the same as | -| | | wrapping in latitude (see polar | -| | | below). | -+---------------------------------------+-------------------+---------------------------------------+ -| polar | logical | If *.true.*, points at the poles are | -| | | wrapped across the grid. It is not | -| | | clear this is a good idea since the | -| | | grid is degnerate here. | -+---------------------------------------+-------------------+---------------------------------------+ -| scm | logical | If *.true.* the Single Column Model | -| | | is assumed. The grid is a single | -| | | vertical column, and there are 9 | -| | | cells arranged in a 3x3 grid. See the | -| | | WRF documentation for more | -| | | information on this configuration. | -| | | *periodic_x* and *periodic_y* should | -| | | also be *.true.* in this case. | -+---------------------------------------+-------------------+---------------------------------------+ -| sfc_elev_max_diff | real(r8) | If > 0, the maximum difference, in | -| | | meters, between an observation marked | -| | | as a 'surface obs' as the vertical | -| | | type (with the surface elevation, in | -| | | meters, as the numerical vertical | -| | | location), and the surface elevation | -| | | as defined by the model. Observations | -| | | further away from the surface than | -| | | this threshold are rejected and not | -| | | assimilated. If the value is | -| | | negative, this test is skipped. | -+---------------------------------------+-------------------+---------------------------------------+ -| allow_obs_below_vol | logical | If *.false.* then if an observation | -| | | with a vertical coordinate of | -| | | pressure or height (i.e. not a | -| | | surface observation) is below the | -| | | lowest 3d sigma level, it is outside | -| | | the field volume and the | -| | | interpolation routine rejects it. If | -| | | this is set to *.true.* and the | -| | | observation is above the surface | -| | | elevation but below the lowest field | -| | | volume level, the code will | -| | | extrapolate downward from data values | -| | | at levels 1 and 2. | -+---------------------------------------+-------------------+---------------------------------------+ -| center_search_half_length | real(r8) | The model_mod now contains two | -| | | schemes for searching for a vortex | -| | | center location. If the **old** | -| | | scheme is compiled in, then this and | -| | | the center_spline_grid_scale namelist | -| | | items are used. (Search code for | -| | | 'use_old_vortex'.) Half length (in | -| | | meters) of a square box for searching | -| | | the vortex center. | -+---------------------------------------+-------------------+---------------------------------------+ -| center_spline_grid_scale | integer | The model_mod now contains two | -| | | schemes for searching for a vortex | -| | | center location. If the **old** | -| | | scheme is compiled in, then this and | -| | | the center_search_half_length | -| | | namelist items are used. (Search code | -| | | for 'use_old_vortex'.) Ratio of | -| | | refining grid for | -| | | spline-interpolation in determining | -| | | the vortex center. | -+---------------------------------------+-------------------+---------------------------------------+ -| circulation_pres_level | real(r8) | The model_mod now contains two | -| | | schemes for searching for a vortex | -| | | center location. If the **new** | -| | | scheme is compiled in, then this and | -| | | the circulation_radius namelist items | -| | | are used. (Search code for | -| | | 'use_old_vortex'.) Pressure, in | -| | | pascals, of the level at which the | -| | | circulation is computed when | -| | | searching for the vortex center. | -+---------------------------------------+-------------------+---------------------------------------+ -| circulation_radius | real(r8) | The model_mod now contains two | -| | | schemes for searching for a vortex | -| | | center location. If the **new** | -| | | scheme is compiled in, then this and | -| | | the circulation_pres_level namelist | -| | | items are used. (Search code for | -| | | 'use_old_vortex'.) Radius, in meters, | -| | | of the circle over which the | -| | | circulation calculation is done when | -| | | searching for the vortex center. | -+---------------------------------------+-------------------+---------------------------------------+ -| vert_localization_coord | integer | Vertical coordinate for vertical | -| | | localization. | -| | | | -| | | - 1 = model level | -| | | - 2 = pressure (in pascals) | -| | | - 3 = height (in meters) | -| | | - 4 = scale height (unitless) | -+---------------------------------------+-------------------+---------------------------------------+ -| allow_perturbed_ics | logical | *allow_perturbed_ics* should not be | -| | | used in most cases. It is provided | -| | | only as a means to create a tiny | -| | | ensemble for non-advancing tests. | -| | | Creating an initial ensemble is | -| | | covered in :doc:`./tutorial/README` | -+---------------------------------------+-------------------+---------------------------------------+ - - -The following items used to be in the WRF namelist but have been removed. The -first 4 are no longer needed, and the last one was moved to the -``&dart_to_wrf_nml`` namelist in 2010. In the Lanai release having these values -in the namelist does not cause a fatal error, but more recent versions of the -code will fail if any of these values are specified. Remove them from your -namelist to avoid errors. - -=================== ================= ========================================= -Item Type Description -=================== ================= ========================================= -``surf_obs`` logical OBSOLETE -- now an error to specify this. -``soil_data`` logical OBSOLETE -- now an error to specify this. -``h_diab`` logical OBSOLETE -- now an error to specify this. -``num_moist_vars`` integer OBSOLETE -- now an error to specify this. -``adv_mod_command`` character(len=32) OBSOLETE -- now an error to specify this. -=================== ================= ========================================= - -Files ------ - -- model_nml in input.nml -- wrfinput_d01, wrfinput_d02, ... (one file for each domain) -- netCDF output state diagnostics files + / + + +Namelist Description: +~~~~~~~~~~~~~~~~~~~~~ + ++-------------------------------+-------------------+---------------------------------------+ +| Item | Type | Description | ++===============================+===================+=======================================+ +| default_state_variables | logical | If *.true.*, the dart state vector | +| | | contains the fields U, V, W, PH, THM, | +| | | MU, in that order, and only those. | +| | | Any values listed in the | +| | | *wrf_state_variables* namelist item | +| | | will be ignored. | ++-------------------------------+-------------------+---------------------------------------+ +| wrf_state_variables | character(:,5) | A 2D array of strings, 5 per wrf | +| | | array to be added to the dart state | +| | | vector. If *default_state_variables* | +| | | is *.true.*, this is ignored. When | +| | | *.false.*, this list of array names | +| | | controls which arrays and the order | +| | | that they are added to the state | +| | | vector. The 5 strings are: | +| | | | +| | | #. WRF field name - must match netcdf | +| | | name exactly | +| | | #. DART Quantity name - must match a | +| | | valid DART QTY_xxx exactly | +| | | #. WRF Type - supplements the quantity| +| | | name to control the operation of | +| | | forward operator. | +| | | #. The string UPDATE or NO_COPY_BACK | +| | | Determines whether WRF state | +| | | is updated by DART | +| | | #. A numeric string listing the | +| | | domain(s) that include the WRF | +| | | state variable. | +| | | The special string '999' means | +| | | all domains. For example, '12' | +| | | means domains 1 and 2, '13' means | +| | | 1 and 3. | ++-------------------------------+-------------------+---------------------------------------+ +| wrf_state_bounds | character(:,4) | A 2D array of strings, 4 per wrf | +| | | array. During the copy of data to and | +| | | from the WRF (wrfinput*) file, | +| | | variables listed here will have | +| | | minimum and maximum values enforced. | +| | | The 4 strings are: | +| | | | +| | | #. WRF field name - must match | +| | | WRF variable name exactly | +| | | #. Minimum -- specified as a string | +| | | but must be a numeric value (e.g. | +| | | '0.1') Can be 'NULL' to allow any | +| | | minimum value. | +| | | #. Maximum -- specified as a string | +| | | but must be a numeric value (e.g. | +| | | '0.1') Can be 'NULL' to allow any | +| | | maximum value. | +| | | #. Action -- valid strings are | +| | | 'CLAMP' or 'FAIL'. Ignored by | +| | | filter. Filter will always clamp | +| | | if min and/or max is set. | ++-------------------------------+-------------------+---------------------------------------+ +| num_domains | integer | Total number of WRF domains, | +| | | including nested domains. | ++-------------------------------+-------------------+---------------------------------------+ +| calendar_type | integer | Calendar type. Should be 3 | +| | | (GREGORIAN) for WRF. | ++-------------------------------+-------------------+---------------------------------------+ +| assimilation_period_seconds | integer | The time (in seconds) between | +| | | assimilations. This is modified if | +| | | necessary to be an integer multiple | +| | | of the underlying model timestep. | ++-------------------------------+-------------------+---------------------------------------+ +| periodic_x | logical | If *.true.*, the grid is periodic in | +| | | longitude, and points above the last | +| | | grid cell and points below the first | +| | | grid cell are wrapped. Note this is | +| | | not the same as a grid which crosses | +| | | the prime meridian. WRF handles that | +| | | with an offset in longitude and | +| | | points beyond the last grid index are | +| | | outside the domain. | ++-------------------------------+-------------------+---------------------------------------+ +| periodic_y | logical | Used for the WRF single column model | +| | | to make the grid wrap in Y (see scm | +| | | below). This is NOT the same as | +| | | wrapping in latitude (see polar | +| | | below). | ++-------------------------------+-------------------+---------------------------------------+ +| polar | logical | If *.true.*, points at the poles are | +| | | wrapped across the grid. It is not | +| | | clear this is a good idea because the | +| | | grid is degnerate here. | ++-------------------------------+-------------------+---------------------------------------+ +| scm | logical | If *.true.* the single column model | +| | | is assumed. The grid is a single | +| | | vertical column, and there are 9 | +| | | cells arranged in a 3x3 grid. See the | +| | | WRF documentation for more | +| | | information on this configuration. | +| | | *periodic_x* and *periodic_y* should | +| | | also be *.true.* in this case. | ++-------------------------------+-------------------+---------------------------------------+ +| sfc_elev_max_diff | real(r8) | The maximum elevation difference | +| | | (in meters) between a 'surface' | +| | | observation and the land surface | +| | | elevation defined in WRF. | +| | | If the value is > 0, that value is | +| | | the threshold at which the surface | +| | | observations are rejected. If the | +| | | value is negative the test is skipped.| ++-------------------------------+-------------------+---------------------------------------+ +| allow_obs_below_vol | logical | If *.false.* then if an observation | +| | | with a vertical coordinate of | +| | | pressure or height (i.e. not a | +| | | surface observation) is below the | +| | | lowest 3d sigma level, it is outside | +| | | the field volume and the | +| | | interpolation routine rejects it. If | +| | | this is set to *.true.* and the | +| | | observation is above the surface | +| | | elevation but below the lowest field | +| | | volume level, the code will | +| | | extrapolate downward from data values | +| | | at levels 1 and 2. | ++-------------------------------+-------------------+---------------------------------------+ +| center_search_half_length | real(r8) | A parameter in the 'use_old_vortex' | +| | | scheme used to search for a vortex | +| | | center location. It is the half-length| +| | | (meters) of a square box used during | +| | | the vortex search. This value and the | +| | | 'center_spline_grid_scale' namelist | +| | | items are required. To implement, set | +| | | ``use_old_vortex = .true.`` in | +| | | ``model_mod.f90`` prior to compiling | +| | | DART. | ++-------------------------------+-------------------+---------------------------------------+ +| center_spline_grid_scale | integer | A parameter in the 'use_old_vortex' | +| | | scheme used to search for a vortex | +| | | center location. It is the fine grid | +| | | ratio for the spline interpolation | +| | | used during the vortex search. This | +| | | value and the | +| | | 'center_search_half_length' namelist | +| | | items are required. To implement, set | +| | | ``use_old_vortex = .true.`` in | +| | | ``model_mod.f90`` prior to compiling | +| | | DART. | ++-------------------------------+-------------------+---------------------------------------+ +| circulation_pres_level | real(r8) | A parameter in the 'circulation' | +| | | scheme used to search for a vortex | +| | | center location. It is the pressure | +| | | (Pascals) at which the circulation is | +| | | computed during the vortex search. | +| | | This value and the | +| | | 'circulation_radius' namelist items | +| | | are required. To implement, set | +| | | ``use_old_vortex = .false.`` in | +| | | ``model_mod.f90`` prior to compiling | +| | | DART. | ++-------------------------------+-------------------+---------------------------------------+ +| circulation_radius | real(r8) | A parameter in the 'circulation' | +| | | scheme used to search for a vortex | +| | | center location. It is the radius | +| | | (meters) of the circle over which the | +| | | search for the vortex center is | +| | | performed. This value and the | +| | | 'circulation_pres_level' namelist | +| | | items are required. To implement, | +| | | set ``use_old_vortex = .false.`` in | +| | | ``model_mod.f90`` prior to compiling | +| | | DART. | ++-------------------------------+-------------------+---------------------------------------+ +| vert_localization_coord | integer | Vertical coordinate for vertical | +| | | localization. | +| | | | +| | | - 1 = model level | +| | | - 2 = pressure (in pascals) | +| | | - 3 = height (in meters) | +| | | - 4 = scale height (unitless) | ++-------------------------------+-------------------+---------------------------------------+ +| allow_perturbed_ics | logical | *allow_perturbed_ics* should not be | +| | | used in most cases. It is provided | +| | | only as a means to create a tiny | +| | | ensemble for non-advancing tests. | +| | | Creating an initial ensemble is | +| | | covered in :doc:`./tutorial/README` | ++-------------------------------+-------------------+---------------------------------------+ + + +Additional Namelist Information +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +- default_state_variables + +You must set ``default_state_variables = .false.`` before changing the value +of ``wrf_state_variables`` to have it take effect. + + +- wrf_state_variables + +The format for ``wrf_state_variables`` is an array of 5 strings: +WRF output field, DART Quantity, WRF TYPE, 'UPDATE' or 'NO_COPY_BACK', and a numerical +string 'XXX'. If XXX=999 the variable is part of all domains, otherwise it is limited +to specific domains (e.g. '12' for domains 1 and 2, '13' for domains 1 and 3). +For example: + +:: + + wrf_state_variables='U','QTY_U_WIND_COMPONENT','TYPE_U','UPDATE','999', + 'V','QTY_V_WIND_COMPONENT','TYPE_V','UPDATE','999', + 'W','QTY_VERTICAL_VELOCITY','TYPE_W','UPDATE','999', + 'THM','QTY_POTENTIAL_TEMPERATURE','TYPE_T','UPDATE','999', + 'PH','QTY_GEOPOTENTIAL_HEIGHT','TYPE_GZ','UPDATE','999', + 'MU','QTY_PRESSURE','TYPE_MU','UPDATE','999', + 'QVAPOR','QTY_VAPOR_MIXING_RATIO','TYPE_QV','UPDATE','999', + 'QCLOUD','QTY_CLOUD_LIQUID_WATER','TYPE_QC','UPDATE','999', + 'QRAIN','QTY_RAINWATER_MIXING_RATIO','TYPE_QR','UPDATE','999', + 'U10','QTY_U_WIND_COMPONENT','TYPE_U10','UPDATE','999', + 'V10','QTY_V_WIND_COMPONENT','TYPE_V10','UPDATE','999', + 'T2','QTY_TEMPERATURE','TYPE_T2','UPDATE','999', + 'TH2','QTY_POTENTIAL_TEMPERATURE','TYPE_TH2','UPDATE','999', + 'Q2','QTY_SPECIFIC_HUMIDITY','TYPE_Q2','UPDATE','999', + 'PSFC','QTY_PRESSURE','TYPE_PS','UPDATE','999', + + +- polar, periodic_x + +The ``Polar`` and ``periodic_x`` namelist values are used in global WRF simulations. +If ``polar`` is true, the grid interpolation routines will wrap over the north and south poles. +If ``periodic_x`` is true, when the east and west edges of the grid are +reached the interpolation will wrap. Note this is a separate issue +from regional models which cross the GMT line. Those grids are marked +as having a negative offset and do not need to wrap. This flag controls +what happens when the edges of the grid are reached. + + +- Single Column Model (scm) + +The ``scm`` flag is used for the single column model version of WRF. +It needs the periodic_x and periodic_y flags set to true, in which +case the X and Y directions are periodic. There is no collapsing of the grid +into a single location like the 3d-spherical polar flag implies. + + +- sfc_elev_max_diff + +The intent of the ``sfc_elev_max_diff`` quality control check is to eliminate +surface observations that are mismatched from the WRF model's surface elevation. +Mismatch can occur if the WRF land surface elevation is not finely resolved (coarse grid) +thus there is a significant representation mismatch between a point observation +and the WRF model. Assimilating surface observations with large mismatch can +deprecate assimilation forecast skill. +This check can only be applied to **surface observations** which are automatically +assigned to observations that use the ``VERTISSURFACE`` vertical coordinate +defined in the ``obs_seq.out`` file. + + +- allow_obs_below_vol + +The ``allow_obs_below_vol`` enables vertical extrapolation in cases where the +observation vertical location is below the lowest WRF model vertical layer, thus +used as an alternative for the standard vertical interpolation routine. +The bottom WRF layer can vary based on total vertical levels, however, in general, +descends to (roughly) 10-50 meters above the surface and does not encompass common +surface observations at 2 and 10 meters. This is not recommended given +(linear) extrapolation is a poor approximation of surface observations at the +land-atmosphere boundary where energy and vapor exchange are controlled by +similarity theory. When using surface observations it is preferred +(and the default of the WRF ``model_mod.f90``) to operate on the WRF 2D +surface output (e.g. T2, U10) instead of WRF 3D output (e.g. T, THM) to +avoid the need for extrapolation. + + +- Vortex option + +The vortex searching namelist options are only required during WRF simulations +where the spatial domain of interest is dynamic such as with a hurricane. + + + References ---------- diff --git a/models/wrf/shell_scripts/driver.csh b/models/wrf/shell_scripts/driver.csh index dbe9b6c11..f0573607d 100755 --- a/models/wrf/shell_scripts/driver.csh +++ b/models/wrf/shell_scripts/driver.csh @@ -473,19 +473,18 @@ while ( 1 == 1 ) else if ( $SUPER_PLATFORM == 'derecho' ) then - # Prevent double submission for member 1 only - if ( $n == 1) then - sleep 5 - endif if ( `qstat -wa | grep assim_advance_${n} | wc -l` == 0 ) then - echo "assim_advance_${n} is missing from the queue" - qsub assim_advance_mem${n}.csh + echo "Warning, detected that assim_advance_${n} is missing from the queue" + echo "If this warning leads to missing output from ensemble ${n}" + echo "consider enabling the qsub command within keep_trying while statement in driver.csh" + + #qsub assim_advance_mem${n}.csh endif endif - sleep 15 + sleep 5 end set start_time = `head -1 start_member_${n}` diff --git a/models/wrf/shell_scripts/mean_increment.ncl b/models/wrf/shell_scripts/mean_increment.ncl index d6e230e41..c12b8e7af 100644 --- a/models/wrf/shell_scripts/mean_increment.ncl +++ b/models/wrf/shell_scripts/mean_increment.ncl @@ -1,7 +1,7 @@ ; find the mean state space increment, output the fields to a single mean file ; that can be used to make plots ; G. Romine 2011-12 - +; Updating for 1 domain only B. Raczka 2024-08 begin ; get the list of files to read in @@ -23,9 +23,7 @@ begin ListSetType(fil, "join") pull_2D_field_names = (/"T2", "Q2", "U10", "V10", "PSFC"/) - pull_2D_field_names(:) = pull_2D_field_names(:)+"_d01" - pull_3D_field_names = (/"U", "V", "T", "QVAPOR"/) - pull_3D_field_names(:) = pull_3D_field_names(:)+"_d01" + pull_3D_field_names = (/"U", "V", "THM", "QVAPOR"/) npulls = dimsizes(pull_2D_field_names) ; Below will dump out the data to a file for replotting later @@ -34,16 +32,17 @@ begin do i=0,npulls-1 print(" Extracting 2d variable "+pull_2D_field_names(i)) do fil_num=0,nfils-1 -; print(" reading file "+flist(fil_num)) -; dimensions are ncljoin, copy, Time, south_north, west_east +; print(" reading file "+flist(fil_num)) +; dimensions are ncljoin, Time, south_north, west_east ; copy zero is the ensemble mean - pull_var = fil[fil_num]->$pull_2D_field_names(i)$(:,0,:,:,:) + pull_var = fil[fil_num]->$pull_2D_field_names(i)$(:,:,:,:) dims = dimsizes(pull_var) if (fil_num .eq. 0) then ; first iteration, make var alltimes_var = new ( (/nfils,dims(2),dims(3)/), typeof(pull_var) ) end if ; printVarSummary(pull_var) alltimes_var(fil_num,:,:) = pull_var(0,0,:,:) +; printVarSummary(alltimes_var) delete(pull_var) end do ; average over time (first dimension) @@ -69,9 +68,9 @@ begin print(" Extracting 3d variable "+pull_3D_field_names(i)) do fil_num=0,nfils-1 ; print(" reading file "+flist(fil_num)) -; dimensions are ncljoin, copy, Time, level, south_north, west_east +; dimensions are ncljoin, Time, level, south_north, west_east ; copy zero is the ensemble mean - pull_var = fil[fil_num]->$pull_3D_field_names(i)$(:,0,:,:,:,:) + pull_var = fil[fil_num]->$pull_3D_field_names(i)$(:,:,:,:,:) dims = dimsizes(pull_var) if (fil_num .eq. 0) then ; first iteration, make var alltimes_var = new ( (/nfils,dims(2),dims(3),dims(4)/), typeof(pull_var) ) diff --git a/models/wrf/tutorial/README.rst b/models/wrf/tutorial/README.rst index 80cd82c34..2dd00180d 100644 --- a/models/wrf/tutorial/README.rst +++ b/models/wrf/tutorial/README.rst @@ -600,7 +600,7 @@ also want to modify this script to test running a single member first — just in case you have some debugging to do. However, be warned that to successfully complete the tutorial, including -running the *driver.csh* script in Step 5, using a smaller ensemble +running the *driver.csh* script in Step 6, using a smaller ensemble (e.g. < 20 members) can lead to spurious updates during the analysis step, causing the WRF simulation to fail. @@ -617,7 +617,7 @@ Step 3: Prepare observations [OPTIONAL] The observation sequence files to run this tutorial are already provided for you. If you want to run with the provided tutorial observations, you - can skip to Step 4 right now. If you are interested in using custom + can skip to Step 5 right now. If you are interested in using custom observations for a WRF experiment other than the tutorial you should read on. The remaining instructions provided below in Step 3 are meant as a guideline to converting raw PREPBUFR data files into the required ``obs_seq`` format @@ -807,8 +807,228 @@ you would do the following: namelist options to consider, and you must provide a *wrfinput* file for the program to access the analysis domain information. +Step 4: Overview of forward (observation) operators [OPTIONAL] +-------------------------------------------------------------- -Step 4: Creating the first set of adaptive inflation files +This section is for informational purposes only and does not include any +instructions to complete the tutorial. It provides a description of +the DART settings that control the forward operator which +calculates the prior and posterior model estimates for the observations. +An introduction to important namelist variables that control the operation of the forward +operator are located in the `WRF namelist documentation. +<../../../models/wrf/readme.html#namelist>`__ + + +The ``obs_seq.out`` file generated as described in Step 3 includes a total +of 30 observation types. Here we examine an exerpt of that file, focusing +on a single temperature observation to describe the process: + +:: + + obs_sequence + obs_kind_definitions + 30 + 41 METAR_TEMPERATURE_2_METER + .. + .. + num_copies: 1 num_qc: 1 + num_obs: 70585 max_num_obs: 70585 + NCEP BUFR observation + NCEP QC index + first: 1 last: 70585 + OBS 1 + 288.750000000000 + 1.00000000000000 + -1 2 -1 + obdef + loc3d + 4.819552185804497 0.6141813398083548 518.0000000000000 -1 + kind + 41 + 43200 152057 + 3.06250000000000 + .. + .. + .. + + +A critical piece of observation metadata includes the observation type +(``METAR_TEMPERATURE_2_METER``) which is linked to the quantity +(``QTY_2M_TEMPERATURE``) through the observation definition file +(``obs_def_metar_mod.f90``). This file is included within the +``&preprocess_nml`` section of the namelist file as: + +:: + + &preprocess_nml + overwrite_output = .true. + input_obs_qty_mod_file = '../../../assimilation_code/modules/observations/DEFAULT_obs_kind_mod.F90' + output_obs_qty_mod_file = '../../../assimilation_code/modules/observations/obs_kind_mod.f90' + input_obs_def_mod_file = '../../../observations/forward_operators/DEFAULT_obs_def_mod.F90' + output_obs_def_mod_file = '../../../observations/forward_operators/obs_def_mod.f90' + quantity_files = '../../../assimilation_code/modules/observations/atmosphere_quantities_mod.f90' + obs_type_files = '../../../observations/forward_operators/obs_def_reanalysis_bufr_mod.f90', + '../../../observations/forward_operators/obs_def_altimeter_mod.f90', + '../../../observations/forward_operators/obs_def_radar_mod.f90', + '../../../observations/forward_operators/obs_def_metar_mod.f90', + .. + .. + .. + +During the DART compilation described within Step 1 this information is +included within the ``obs_def_mod.f90``. + +The vertical coordinate type is the 4th column beneath the loc3d header within ``obs_seq.out``. +In this example the value -1 indicates the vertical coordinate is ``VERTISSURFACE``. It defines the +vertical units of the observation (e.g. pressure, meters above sea level, model levels etc). +This serves two purposes -- foremost it is required during the vertical spatial interpolation +to calculate the precise location of the expected observation. +A second crtical function is that it defines whether it is a surface observation. +Observations with a vertical coordinate of ``VERTISSURFACE`` are defined as surface +observations. All other coordinates are considered non-surface observations +(e.g. profile observations). Of note is that the vertical coordinate ``VERTISSURFACE`` and +``VERTISHEIGHT`` are functionally identical (i.e. meters above sea level), however +only the ``VERTISSURFACE`` is a surface observation. + +For more information on the vertical coordinate metadata see the detailed structure of +an `obs_seq file. <../../../guide/creating-obs-seq-real.html#observation-location>`__ + +In order to connect this observation to the appropriate WRF output variables +the ``wrf_state_variables`` within ``&model_nml`` defines the *WRF field name* and +the *WRF TYPE* in the 1st and 3rd columns as shown in the tutorial example below: + +:: + + &model_nml + wrf_state_variables = 'T2','QTY_TEMPERATURE','TYPE_T2','UPDATE','999' + + .. + .. + +For more information on the ``&model_nml`` variables see the `WRF documentation page +<../../../models/wrf/readme.html#namelist>`__ + + +As described above, the linkage between the observation type and the WRF output field +is defined through the physical quantity, surface variable designation (observation +vertical coordinate), and WRF TYPE. The current design of the WRF ``model_mod.f90`` +is such that the quantity is a general classification (e.g. temperature, wind +specific humidity), whereas the WRF TYPE classification is more precisely +mapped to the WRF output field. The table below summarizes the dependency between +the observation type and the WRF output field for a select number of observation types +within the tutorial. + +.. Note:: + + The number of WRF output fields required to support an observation type can vary. For + observation types where there is a direct analog to a WRF output field, the forward + operator consists of only spatial interpolation, thus requires only a single output + variable (e.g. METAR_TEMPERATURE_2_METER). For observation types that require multiple + WRF output fields, the forward operator is more complex than a simple spatial interpolation. + For more information see the notes below the table. A rule of thumb is a surface + observation should use a surface output field (e.g. T2, U10). WRF surface output fields + are appended by a numeric value indicating surface height in meters. It is possible to use + a non-surface WRF output field (3D field) to estimate a surface observation, however, this + requires a vertical interpolation of the 3D WRF field where the observed surface height does + not coincide with the model levels. This either requires a vertical interpolation or an + extrapolation which can be **inaccurate and is not recommended**. + + + + ++----------------------------------+---------+-------------------------------+--------------+------------+ +| DART Observation Type | Surface | DART Quantity | WRF Type | WRF output | +| | Obs ? | | | field | ++==================================+=========+===============================+==============+============+ +| ``METAR_TEMPERATURE_2_METER`` | Yes | ``QTY_2M_TEMPERATURE`` | ``TYPE_T2`` | ``T2`` | +| | | | | | ++----------------------------------+---------+-------------------------------+--------------+------------+ +| ``RADIOSONDE_TEMPERATURE`` | No | ``QTY_POTENTIAL_TEMPERATURE`` | ``TYPE_T`` | ``THM`` | +| | | ``QTY_VAPOR_MIXING_RATIO`` | ``TYPE_QV`` | ``QVAPOR`` | +| | | ``QTY_PRESSURE`` | ``TYPE_MU`` | ``MU PH`` | +| | | ``QTY_GEOPOTENTIAL_HEIGHT`` | ``TYPE_GZ`` | | ++----------------------------------+---------+-------------------------------+--------------+------------+ +| ``METAR_U_10_METER_WIND`` | Yes | ``QTY_U_WIND_COMPONENT`` | ``TYPE_U10`` | ``U10`` | +| | | ``QTY_V_WIND_COMPONENT`` | ``TYPE_V10`` | ``V10`` | ++----------------------------------+---------+-------------------------------+--------------+------------+ +| ``ACARS_U_WIND_COMPONENT`` | No | ``QTY_U_WIND_COMPONENT`` | ``TYPE_U`` | ``U`` | +| | | ``QTY_V_WIND_COMPONENT`` | ``TYPE_V`` | ``V`` | ++----------------------------------+---------+-------------------------------+--------------+------------+ +| ``METAR_DEWPOINT_2_METER`` | Yes | ``QTY_DEWPOINT`` | | | +| | | ``QTY_SPECIFIC_HUMIDITY`` | ``TYPE_Q2`` | ``Q2`` | +| | | ``QTY_PRESSURE`` | ``TYPE_PS`` | ``PSFC`` | ++----------------------------------+---------+-------------------------------+--------------+------------+ +| ``RADIOSONDE_SPECIFIC_HUMIDITY`` | No | ``QTY_SPECIFIC_HUMIDITY`` | ``TYPE_QV`` | ``QVAPOR`` | +| | | | | | ++----------------------------------+---------+-------------------------------+--------------+------------+ + + + +Surface Temperature (e.g. METAR_TEMPERATURE_2_METER) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +WRF output includes a direct analog for sensible temperature surface observations (e.g. T2), thus +the forward operator requires only 1 variable to calculate the expected observation. +The calculation includes a horizontal interpolation of the 2D temperature variable (e.g. T2). + + +Non-Surface Temperature (e.g. RADIOSONDE_TEMPERATURE) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +In contrast to surface temperature observations, non-surface temperature observations require 4 WRF +output fields. This is because observations are sensible temperature, whereas the 3D WRF +temperature field is provided in perturbation potential temperature. Thus, the forward +operator first requires a physical conversion between perturbation potential temperature to +sensible temperature, followed by a spatial interpolation (this includes horizontal interpolation +on WRF levels k and k+1, followed by vertical interpolation). + +.. Important:: + + There are two different 3D temperature WRF output fields that can work to calculate non- + surface temperature observations (e.g. T or THM, T=THM when use_theta_m=0). However, and **of + utmost importance** is the variable THM is required to be within the ``&model_nml`` if the + 3D temperature field is to be updated in the ``filter`` step. **This is because the WRF field *T* + is a diagnostic variable with no impact on the forecast step, whereas the WRF field *THM* is + a prognostic field which will impact the forecast.** + + +Surface Wind (e.g. METAR_U_10_METER_WIND) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Surface winds have a direct WRF output analog (e.g. U10) +and requires horizontal interpolation of the 2D zonal wind field. However, the +meridional wind (e.g. V10) is also required in order to convert from modeled *gridded* winds to +*true* wind observations. This requirement is an artifact of winds measured on a sphere being +mapped on a 2D grid. + + +Non-Surface Wind (e.g. ACARS_U_WIND_COMPONENT) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This is identical to surface winds as described above, except the spatial interpolation requires +horizontal interpolation on the k and k+1 WRF levels, followed by vertical interpolation. + + +Surface Dewpoint (e.g. METAR_DEWPOINT_2_METER) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The calculation of surface dewpoint requires a physical conversion using both surface +pressure (PSFC) and surface vapor mixing ratio (Q2), follwed by horizontal interpolation. + + +Non-Surface Specific Humidity (e.g. RADIOSONDE_SPECIFIC_HUMIDITY) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Specific humidity observations require the (water) vapor mixing ratio (QVAPOR) for the forward operator. +Although specific humidity and vapor mixing ratio are nearly identical, especially in dry +air, they are actually two distinct physical properties -- the ratio of water mass to total air mass +versus ratio of water vapor mass to dry air mass respectively. Therefore the forward operator +includes this physical conversion followed by a spatial interpolation (i.e. horizontal interpolation of k and +k+1 WRF vertical levels followed by vertical interpolation). + + + +Step 5: Creating the first set of adaptive inflation files ---------------------------------------------------------- In this section we describe how to create initial adaptive inflation @@ -874,7 +1094,7 @@ cycle. -Step 5: Cycled analysis system +Step 6: Cycled analysis system ------------------------------ While the DART system provides executables to perform individual tasks @@ -924,10 +1144,10 @@ continue to cycle until the final analysis time has been reached. -Step 6: Diagnosing the assimilation results +Step 7: Diagnosing the assimilation results ------------------------------------------- -Once you have successfully completed steps 1-5, it is important to +Once you have successfully completed steps 1-6, it is important to check the quality of the assimilation. In order to do this, DART provides analysis system diagnostics in both state and observation space. @@ -989,7 +1209,7 @@ The tools below provide methods to visualize the spatial patterns, statistics an failure mode for all observations. The observation diagnostics use the **obs_epoch*.nc** file as input. This file is -automatically generated by the **obs_diagnostic.csh** script within Step 5 of this +automatically generated by the **obs_diagnostic.csh** script within Step 6 of this tutorial. The **obs_epoch*.nc** file is located in the output directory of each time step. @@ -1211,7 +1431,7 @@ quite high (>90%). This high acceptance percentage is typical of a high-quality assimilation and consistent with the strong reduction in RMSE. -The same plot as above is given below except for the observation type: +The same plot as above except for the observation type: ``RADIOSONE_SPECIFIC_HUMIDITY``. +-------------------------------------------------------------+ diff --git a/models/wrf/wrf_state_variables_table b/models/wrf/wrf_state_variables_table index 73ee3b537..a224cfb8a 100644 --- a/models/wrf/wrf_state_variables_table +++ b/models/wrf/wrf_state_variables_table @@ -4,14 +4,13 @@ # # DART $Id$ -This table contains the description of a value, followed by the line(s) +This table contains the description of an atmospheric property, followed by the line(s) that should be added to the DART input.nml namelist &model_nml, 'wrf_state_variables' list. All items must be strings. The 5 columns are: -Exact variable name in WRF netcdf file, DART Kind String, a Type -column that has a suggested value but is unused, a flag to say -whether the variable should be updated or not (currently only 'UPDATE' is -supported), and a flag to say which of multiple domains should have this -variable, where '999' means all domains. +exact variable name in WRF output file, DART Quantity , a WRF Type, +a flag to determine whether the variable is updated during the analysis +('UPDATE' 'NO_COPY_BACK'), and a flag to identify which domains include + this variable, where '999' means all domains. Horizontal Winds: @@ -21,7 +20,7 @@ Horizontal Winds: 'V10', 'QTY_V_WIND_COMPONENT', 'TYPE_V10', 'UPDATE', '999', Sensible Temperature: - 'T', 'QTY_TEMPERATURE', 'TYPE_T', 'UPDATE', '999', + 'THM', 'QTY_TEMPERATURE', 'TYPE_T', 'UPDATE', '999', 'T2', 'QTY_TEMPERATURE', 'TYPE_T2', 'UPDATE', '999', Potential Temperature: