Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature request: Refactor wrf model_mod and models/wrf directory #404

Open
hkershaw-brown opened this issue Oct 4, 2022 · 29 comments · May be fixed by #683
Open

Feature request: Refactor wrf model_mod and models/wrf directory #404

hkershaw-brown opened this issue Oct 4, 2022 · 29 comments · May be fixed by #683
Assignees
Labels
missing_state missing_r8 in the state wrf Weather Research & Forecasting Model

Comments

@hkershaw-brown
Copy link
Member

Use case

Refactor the wrf model_mod to make it easier/possible to:

  • add state variables
  • swap between wrf v3 and v4
  • not assume the structure of the state variable in model_mod.

Is your feature request related to a problem?

The WRF model_mod has wrf_static_data_for_dart which is replicating functionality available in state_structure_mod.

  • Adding variable X to the state means altering the model_mod code to have a wrf%dom(id)%type_X, rather than simply adding a line in the model_nml::wrf_state_variables.

  • WRF has its own integer types (qtys), replicating (but not quite) the DART qty functionality.

    ! JPH local variables to hold type indices
    integer :: type_u, type_v, type_w, type_t, type_qv, type_qr, type_hdiab, &
    type_qndrp, type_qnsnow, type_qnrain, type_qngraupel, type_qnice, &
    type_qc, type_qg, type_qi, type_qs, type_gz, type_refl, type_fall_spd, &
    type_dref, type_spdp, type_qh, type_qnhail, type_qhvol, type_qgvol, &
    type_cldfra
    integer :: type_u10, type_v10, type_t2, type_th2, type_q2, &
    type_ps, type_mu, type_tsk, type_tslb, type_sh2o, &
    type_smois, type_2dflash

    so you see code like this:
    var_type = wrf%dom(id)%var_type(var_id)
    dart_type = wrf%dom(id)%dart_kind(var_id)

  • WRF v3 vs. v4 - it is not clear for users what to change when moving from WRF v3 to v4. If possible, it would be good for the model_mod to be able to check the WRF version from the netcdf file.

  • Wider issue Model_mods which use their own structure to keep track of state variables #389. Moving control of the state vector to core dart modules allows us to have general solutions for things like compression, removing missing values (so dart core code can assume every element of the state is valid), and other cool stuff.

Describe your preferred solution

Remove the state vector information from wrf_static_data_for_dart

Describe any alternatives you have considered

I think the wrf_static_data_for_dart is a road block for core dart development, however wrf-dart is widely used and has a lot of additional user developed code in circulation. There is 8000 lines of code in for wrf model_mod (+ unknown user modifications) so the refactoring needs to be approached with this in mind.
I think it is possible to simplify the bounds checks for each grid (staggered and various flavors of periodic), but I do not think the bounds check affects development of core dart.

@hkershaw-brown hkershaw-brown changed the title Feature request: Refactor wrf model_mod Feature request: Refactor wrf model_mod and models/wrf directory Oct 11, 2022
@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Oct 21, 2022

edit:
https://github.com/hkershaw-brown/DART/tree/wrf-refactor is the initial attempt, having a go at updating the wrf model_mod

https://github.com/NCAR/DART/tree/wrf-fresh is starting from scratch redoing the wrf model_mod

@hkershaw-brown hkershaw-brown added the missing_state missing_r8 in the state label Mar 17, 2023
@hkershaw-brown
Copy link
Member Author

wrf model_mod get_close has a test for missing_r8

https://github.com/NCAR/DART/blob/72e5740b71ef6c7cde4c6f946e2c1df2019365fb/models/wrf/model_mod.f90#L6419C3-L6428

should check for verti_is_undefined before doing this (or not have the missing_r8 check)

@hkershaw-brown hkershaw-brown self-assigned this May 26, 2023
@hkershaw-brown hkershaw-brown added the wrf Weather Research & Forecasting Model label May 26, 2023
@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented May 30, 2023

Vertical Convertion

There are several references to the vertical location of state being converted in get_state_meta_data. This is incorrect, get_state_meta_data does not do vertical conversion. For example:

DART/models/wrf/model_mod.f90

Lines 6406 to 6417 in 72e5740

! Convert local vertical coordinate to requested vertical coordinate if necessary.
! This should only be necessary for obs priors, as state location information already
! contains the correct vertical coordinate (filter_assim's call to get_state_meta_data).
if (vertical_localization_on()) then
if (local_which /= wrf%dom(1)%localization_coord) then
call vert_convert(state_handle, local_loc, loc_qtys(t_ind), istatus2)
! Store the "new" location into the original full local array
locs(t_ind) = local_loc !HK Overwritting the location
else
istatus2 = 0
endif
endif

get_close_obs and get_close_state are identical - they both call get_close.
get_close calls vert_convert

convert_vertical_obs calls vert_convert
convert_vertical_state - builds columns from the state

  • unneeded call to get horizontal location, the location is passed into convert_vertical_state

    DART/models/wrf/model_mod.f90

    Lines 3026 to 3027 in 70e6af8

    ! first obtain lat/lon from (ip,jp)
    call get_wrf_horizontal_location( ip, jp, var_type, id, lon, lat )

  • 🐛 Bug: convert_vertical_state is using the module localization coord, not the which_vert passed into convert_vertical_state.

    if (wrf%dom(id)%localization_coord == VERTISLEVEL) then

  • check all surface variables, model_pressure_distrib

    DART/models/wrf/model_mod.f90

    Lines 4807 to 4812 in 70e6af8

    elseif( var_type == wrf%dom(id)%type_mu .or. var_type == wrf%dom(id)%type_tslb .or. &
    var_type == wrf%dom(id)%type_ps .or. var_type == wrf%dom(id)%type_u10 .or. &
    var_type == wrf%dom(id)%type_v10 .or. var_type == wrf%dom(id)%type_t2 .or. &
    var_type == wrf%dom(id)%type_th2 .or. &
    var_type == wrf%dom(id)%type_q2 .or. var_type == wrf%dom(id)%type_tsk .or. &
    var_type == wrf%dom(id)%type_smois .or. var_type == wrf%dom(id)%type_sh2o) then

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Aug 2, 2023

To do:

  • wrf has a namelist option for calendar type. Are there observations not in GREGORIAN (2d?)

  • variable_is_on_domain - assuming <10 domains (i1):

    DART/models/wrf/model_mod.f90

    Lines 7745 to 7746 in f996f26

    read(domain_id_string(i:i),'(i1)') domain_int
    if ( domain_int == id ) variable_is_on_domain = .true.

  • check whether map utils varies between WRF versions module_map_utils.F misc_definitions_module.F

  • restrict_polar is hardcoded in two places (and is not a namelist option)

    ! 10 = polar observation while restrict_polar namelist true

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Aug 9, 2023

For VERTISLEVEL we have:

! return staggered levels for vertical types? check this later.
if( (var_type == wrf%dom(id)%type_w ) .or. (var_type == wrf%dom(id)%type_gz) ) then
lev = real(kp) - 0.5_r8
else
lev = real(kp)
endif

bottom_top = 29 ;
bottom_top_stag = 30 ;
float W(Time, bottom_top_stag, south_north, west_east) ;

Edit: then this:

DART/models/wrf/model_mod.f90

Lines 1857 to 1858 in 9729d78

! Adjust zloc for staggered ZNW grid (or W-grid, as compared to ZNU or M-grid)
zloc = zloc + 0.5_r8

@hkershaw-brown
Copy link
Member Author

note in the code about latlon project settings only valid for global domains

DART/models/wrf/model_mod.f90

Lines 7307 to 7316 in 9729d78

!JPH --- this latinc/loninc calculations are only valid for global domains
if ( wrf%dom(id)%scm ) then
! JPH -- set to zero which should cause the map utils to return NaNs if called
latinc = 0.0_r8
loninc = 0.0_r8
else
latinc = 180.0_r8/wrf%dom(id)%sn
loninc = 360.0_r8/wrf%dom(id)%we
endif

@nancycollins
Copy link
Collaborator

nancycollins commented Aug 10, 2023

for help with the acronyms:

  • code marked 'nc' is actually code done by greg lawson with minimal input from me
  • 'scm' is single column model and was done by josh hacker (JPH)
  • the vortex forward operator code was done by ryan torn (and eventually should be pulled out into an obs_def module

comments on these particular few lines:
'global' doesn't seem to be the opposite of single column, but i'm also sure i don't understand exactly what projection settings are needed to use the wrf code that converts a location into an index into the wrf grid (which is how this model_mod works). most wrf cases are NOT global - they are regional grids. perhaps latinc and loninc aren't used in the regional case.

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Aug 14, 2023

Mask variables: XLAND, QTY_SURFACE_TYPE

Should you be interpolating a mask?
edit: rounded to the nearest integer in rttov obs_def:

surftype = nint(atmos%surftype(imem))

wrf netcdf file XLAND:description = "LAND MASK (1 FOR LAND, 2 FOR WATER)" ;
comment and QTY_SURFACE_TYPE ! 0 = land, 1 = water, 2 = sea ice (not yet supported)

edit:
wrf 1 for land, rttov forward operator expecting 0 for land.
QTY_SURFACE_TYPE has the -1, QTY_LAND_MASK does not.

DART/models/wrf/model_mod.f90

Lines 2825 to 2843 in 9729d78

else if( obs_kind == QTY_SURFACE_TYPE ) then ! 0 = land, 1 = water, 2 = sea ice (not yet supported)
if ( debug ) print*,'Getting land mask'
! 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 ) ) 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 XLAND rc = ', rc
! Interpolation for the XLAND field -- XLAND is NOT part of state vector x, but rather
! in the associated domain meta data
fld(1, :) = dym*( dxm*real(wrf%dom(id)%land(ll(1), ll(2))) + &
dx*real(wrf%dom(id)%land(lr(1), lr(2))) ) + &
dy*( dxm*real(wrf%dom(id)%land(ul(1), ul(2))) + &
dx*real(wrf%dom(id)%land(ur(1), ur(2))) ) - 1
endif

DART/models/wrf/model_mod.f90

Lines 2801 to 2823 in 9729d78

! Land Mask has been added to accommodate satellite observations.
! XLAND is not in the dart_ind vector, so get it from wrf%dom(id)%land
else if( obs_kind == QTY_LANDMASK ) then
if( my_task_id() == 0 ) print*, '*** Land mask forward operator not tested'
if ( debug ) print*,'Getting land mask'
! 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 ) ) 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 XLAND rc = ', rc
! Interpolation for the XLAND field -- XLAND is NOT part of state vector x, but rather
! in the associated domain meta data
fld(1, :) = dym*( dxm*real(wrf%dom(id)%land(ll(1), ll(2))) + &
dx*real(wrf%dom(id)%land(lr(1), lr(2))) ) + &
dy*( dxm*real(wrf%dom(id)%land(ul(1), ul(2))) + &
dx*real(wrf%dom(id)%land(ur(1), ur(2))) )
endif

@hkershaw-brown
Copy link
Member Author

I'm not sure what is not allowed here

DART/models/wrf/model_mod.f90

Lines 6738 to 6740 in 9729d78

! As of 22 Oct 2007, this option is not allowed!
! Note that j = 0 can only happen if we are on the M (or U) wrt to latitude
if ( wrf%dom(id)%polar .and. j == 0 .and. .not. restrict_polar ) then

There are some notes in the code about periodic_x,y and polar

subroutine assign_boundary_conditions(id)

It is looking like these are still wrf namelist options and not in the wrf output netcdf file:
namelist_variables
wrfoutput

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Aug 15, 2023

Boundary condition options:

The boundary conditions periodic_x,y, polar in model_mod.f90 only apply to domain 1. I guess this is why it is a namelist option rather than part of the wrf netcdf metadata(?). There is a lot of checking for in model_mod for (id)%periodic but periodic is only allowed to be true for domain1.

DART/models/wrf/model_mod.f90

Lines 7114 to 7124 in 9729d78

if ( id == 1 ) then
wrf%dom(id)%periodic_x = periodic_x
wrf%dom(id)%periodic_y = periodic_y
wrf%dom(id)%polar = polar
wrf%dom(id)%scm = scm
else
wrf%dom(id)%periodic_x = .false.
wrf%dom(id)%periodic_y = .false.
wrf%dom(id)%polar = .false.
wrf%dom(id)%scm = .false.
endif

Edit: the wrf namelist has periodic_x,y polar for each domain - namelist (max_dom)

The wrf docs has polar as use polar boundary condition (see table). The model_mod.f90 has logical :: polar = .false. ! wrap over the poles Are these the same thing?

polar .false. polar boundary condition (v=0 at polarward-most v-point)

Edit: What is the difference between polar and periodic_y?
Can you have periodic_y and not polar?

if ( wrf%dom(id)%periodic_x .and. .not. wrf%dom(id)%periodic_y  ) then
elseif ( wrf%dom(id)%periodic_x .and. wrf%dom(id)%periodic_y ) then
else
   if ( wrf%dom(id)%polar ) then
   else
     !   NOT Periodic X & M_grid ==> [1 we]
     !   NOT Periodic Y & M_grid ==> [1 sn]
   endif
endif

no .not. wrf%dom(id)%periodic_x .and. wrf%dom(id)%periodic_y option
Just polar:

DART/models/wrf/model_mod.f90

Lines 6275 to 6277 in 70e6af8

if ( wrf%dom(id)%polar ) then
! NOT Periodic X & M_grid ==> [1 we]
! Periodic Y & M_grid ==> [0.5 sn+0.5]

@hkershaw-brown
Copy link
Member Author

Do we like the results?

! once we like the results, remove the log_horz_interpQ test.

! once we like the results, remove the log_horz_interpQ test.

@hkershaw-brown
Copy link
Member Author

  • do you need the check for 0 surface pressure?

    DART/models/wrf/model_mod.f90

    Lines 4600 to 4606 in 9729d78

    ! I'm not quite sure where this comes from, but I will trust them on it....
    ! Do you have to do this per ensemble?
    !>@todo This is messy
    do e = 1,ens_size
    if ( x_ill(e) /= 0.0_r8 .and. x_ilr(e) /= 0.0_r8 .and. x_iul(e) /= 0.0_r8 .and. &
    x_iur(e) /= 0.0_r8 ) then

@hkershaw-brown
Copy link
Member Author

  • qc needed for obs rejected above model surface but below the lowest sigma level?

    DART/models/wrf/model_mod.f90

    Lines 1158 to 1178 in 9729d78

    ! If location is above model surface but below the lowest sigma level,
    ! the default is to reject it. But if the namelist value is true, then
    ! accept the observation and later on extrapolate the values from levels
    ! 1 and 2 downward.
    !HK ensemble loop, however do you reject the obs for all ensembles later anyway?
    do e = 1, ens_size
    if (is_lev0(e)) then
    ! the pres_to_zk() routine has returned a valid zloc in case we
    ! want to use it. the default is to reject the observation and so
    ! we overwrite it with missing -- but, if the namelist value is set
    ! to true, leave zloc alone.
    if (.not. allow_obs_below_vol) zloc(e) = missing_r8
    if (debug .and. .not. allow_obs_below_vol) print*, 'setting zloc missing'
    ! else need to set a qc here?
    endif
    enddo

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Aug 18, 2023

Edit also a set of constants in wrf model_mod:

DART/models/wrf/model_mod.f90

Lines 7997 to 8006 in 9729d78

! Parameters below from WGS-84 model software inside GPS receivers.
parameter(semi_major_axis = 6378.1370d3) ! (m)
parameter(semi_minor_axis = 6356.7523142d3) ! (m)
parameter(grav_polar = 9.8321849378) ! (m/s2)
parameter(grav_equator = 9.7803253359) ! (m/s2)
parameter(earth_omega = 7.292115d-5) ! (rad/s)
parameter(grav_constant = 3.986004418d14) ! (m3/s2)
parameter(grav = 9.80665d0) ! (m/s2) WMO std g at 45 deg lat
parameter(eccentricity = 0.081819d0) ! unitless
parameter(pi2 = 3.14159265358979d0/180.d0)

Edit2: function compute_geometric_height uses gravity = 9.80665d0 vs. the calling code model_height_distrib using gravity = 9.81_r8

hkershaw-brown added a commit that referenced this issue Aug 18, 2023
Needs logic to deal with obs below lowest level
Needs static data to calculate rho
Using two different values for gravity see #404 (comment)
I'm not sure about the check for 0.0 for surface pressure. Does this happen?
For model_rho_t checking speficically for MU variable name rather than qty_pressure
@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Aug 21, 2023

  • geopotential height is using the first ensemble member only. Is this the same across the ensemble?

    ill = get_dart_vector_index(ll(1), ll(2), k(1), domain_id(id), wrf%dom(id)%type_gz)

  • The force non-negative values, are these complete? e.g. there are clamped variables not in this list

    fld = max(0.0_r8, fld) ! Don't accept negative

  • localization coordinate is a namelist option, same across all domains, but is put into the wrf%dom structure.

    if (vert_localization_coord == VERTISLEVEL) then
    wrf%dom(:)%localization_coord = VERTISLEVEL
    Would it change depending on the domain (no?)

  • casting int to real should be real(j, KIND=r8)

    dx = x - real (j)

  • shortest time between assimilations. Current code you can get 0 window size

    model_dt = nint(wrf%dom(1)%dt)
    ! The integer arithmetic does its magic.
    assim_dt = (assimilation_period_seconds / model_dt) * model_dt
    shortest_time_between_assimilations = set_time(assim_dt)

hkershaw-brown added a commit that referenced this issue Aug 21, 2023
Original code is using one ensemble member see:
#404 (comment)
hkershaw-brown added a commit that referenced this issue Aug 28, 2023
see not on periodic_y vs. polar
#404 (comment)

remove mass grid comment from getCorners since this checks for U,V
grid
@NCAR NCAR deleted a comment from nancycollins Sep 1, 2023
hkershaw-brown added a commit that referenced this issue Sep 8, 2023
The vertical conversion is difference for state vs obs.
see #404 (comment)
for note on the (no longer true) assumption in the old code that state elements
have already been vertically converted in get_state_meta data.
@hkershaw-brown
Copy link
Member Author

  • shortest time between assimilations. Current code you can get 0 window size
    model_dt = nint(wrf%dom(1)%dt)
    ! The integer arithmetic does its magic.
    assim_dt = (assimilation_period_seconds / model_dt) * model_dt
    shortest_time_between_assimilations = set_time(assim_dt)

hkershaw-brown added a commit that referenced this issue Sep 12, 2023
…el dt

Previous code would let you set assimilation_period_seconds = 0
which gave you a 0 sized assimilation window
#404 (comment)
@hkershaw-brown
Copy link
Member Author

For convert_vertical_obs, it looks like there is a biggish difference when you change order of calculating geometric from geopotential height. Maybe this is just a bug on my end, but making a note of it here:

  1. 4 pointsx2 levels:
    geopotential height at point, compute_geometric_height(geop, j) at that point
    ! compute height at all neighboring mass points and interpolate

    model_height_w_distrib = compute_geometric_height(geop, wrf%dom(id)%latitude(i, j))

Then interpolate to geometric_height observation location.

  1. geopotential height at location of obs (from 4 points x2 levels), then compute_geometric height(geop, lat of obs)

Example of differences:

geop both methods = 6268.30551946850 

method 1: geometric =  6284.95504384888
method 2: geometric  = 6284.93919703895 

hkershaw-brown added a commit that referenced this issue Oct 4, 2023
…ight

This is quite different from main. I'm not sure if this is a bug on my
end or the order of caluculation makes a big difference.
See #404 (comment)
for details
hkershaw-brown added a commit that referenced this issue Oct 4, 2023
…main)

This is method 1 in #404 (comment) for discussion
zout same value as main
@nancycollins
Copy link
Collaborator

i don't know if this is a related issue, but this reminds me that we had many heated discussions about interpolation order in a cube. i believe if everything is linear the order doesn't matter (within numerical error limits). but if the values are not - like with pressure - then it makes a difference in the results. i'll try to describe this without images.

if you label the 8 corners of a cube (A,B,C,D) for the top, clockwise, and (E,F,G,H) for the bottom, also clockwise, then the two ways you can do interpolation to a point in the interior are:

  1. linearly interpolate in the vertical along each of (A,E), (B,F), (C,G) and (D,H) to get 4 new corner points of a 2d quad. then do bilinear interpolation in that quad for the horizontal location.

  2. bilinearly interpolate in quad (A,B,C,D), and bilinearly interpolate in (E,F,G,H) to get 2 values at the correct horizontal location. then do linear interpolation in the vertical.

we have (i think) tended to do method 2, but honestly right now i can't reconstruct the arguments for and against. it has to do with non-linear fields in the vertical near orography.

@hkershaw-brown
Copy link
Member Author

Comments from Jeff and Moha from standup:
This difference is expected and not big.

For easier comparison (bitwise crutch) with existing wrf code I'm leaving in function interpolate_geometric_height
even though you could have less code and do convert_to_geometric(function geopotential_height_interpolate)

@hkershaw-brown
Copy link
Member Author

@mjs2369

https://github.com/hkershaw-brown/DART/tree/wrf-refactor is the initial attempt, having a go at updating the wrf model_mod. It is essentially my notes about the code. I plopped it on my fork, just to have the notes. It is dead.

https://github.com/NCAR/DART/tree/wrf-fresh is starting from scratch redoing the wrf model_mod.f90.
wrf-fresh is the branch to use.

I'm getting bitwise (last digit) differences in the distance calculations when doing vertical conversion.
I have a tiny (run on a laptop) test case. I'll put this on derecho so you can grab it -stay tuned.

@hkershaw-brown
Copy link
Member Author

tiny test case:
/glade/work/hkershaw/tiny_wrf

hkershaw@derecho8:/glade/work/hkershaw/tiny_wrf$ ls
check_diffs.sh main  restarts  wrf-refactor

restarts - contains the restarts
main - drop in the main branch version of filter
wrf-refactor - drop in the wrf-refactor version of filter

This case is really tiny, you can run on the login node.

checkdiffs.sh is just for comparing netcdf files and obs_seq.final

$ ./checkdiff.sh main wrf-refactor

@hkershaw-brown
Copy link
Member Author

same or different qty?

   wrf_state_variables         = 'U',     'QTY_U_WIND_COMPONENT',     'TYPE_U',   'UPDATE','999',
                                 'V',     'QTY_V_WIND_COMPONENT',      'TYPE_V',  'UPDATE','999',
                                 'U10',   'QTY_U_WIND_COMPONENT',      'TYPE_U10','UPDATE','999',
                                 'V10',   'QTY_V_WIND_COMPONENT',      'TYPE_V10','UPDATE','999',

@hkershaw-brown
Copy link
Member Author

same or different qty?

                'T',     'QTY_POTENTIAL_TEMPERATURE', 'TYPE_T',  'UPDATE','999',
                'T2',    'QTY_TEMPERATURE',           'TYPE_T2', 'UPDATE','999',
                'TH2',   'QTY_POTENTIAL_TEMPERATURE', 'TYPE_TH2','UPDATE','999',

@hkershaw-brown
Copy link
Member Author

DART/models/wrf/model_mod.f90

Lines 7485 to 7491 in 52e6e45

! NOTE: PSFC can be labeled either QTY_PRESSURE or QTY_SURFACE_PRESSURE
! in the namelist, but make sure however it is labeled that for now we
! allow surface pressure interpolation. this may go away once we work out
! QTY_FOO vs QTY_SURFACE_FOO - are they fundamentally different things
! or should the decision be made based on a QTY_FOO and the vertical
! location type -- if it is VERTISSURFACE, then you do the 2d calc in the
! surface field, otherwise you do the full-up 3d interpolation.

hkershaw-brown added a commit that referenced this issue Jan 16, 2024
based on the latest 3d model_mod template

* state variables added by namelist
* removes wrf_static_data_for_data, uses state_stucture_mod

Note for bitwise comparison with the existing wrf model_mod,
using function interpolate_geometric_height
#404 (comment)
@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Jan 16, 2024

@mjs2369

https://github.com/hkershaw-brown/DART/tree/wrf-refactor is the initial attempt, having a go at updating the wrf model_mod. It is essentially my notes about the code. I plopped it on my fork, just to have the notes. It is dead.

https://github.com/NCAR/DART/tree/wrf-fresh is starting from scratch redoing the wrf model_mod.f90. wrf-fresh is the branch to use.

I'm getting bitwise (last digit) differences in the distance calculations when doing vertical conversion. I have a tiny (run on a laptop) test case. I'll put this on derecho so you can grab it -stay tuned.

The bitwise differences were due to this issue #621.
Moved development to https://github.com/NCAR/DART/tree/wrf-wrf-chem (branched from v11 and squished commits)
Todo:

  • Force 2D quantities before add_domain
  • Test with v3, v4, wrf-chem
  • clamping

If you need to look at wrf-fresh:
https://github.com/hkershaw-brown/DART/tree/wrf-fresh

@hkershaw-brown
Copy link
Member Author

I don't think clamp_or_fail is used:

character(len=10), dimension(:),pointer :: clamp_or_fail

https://github.com/search?q=repo%3ANCAR%2FDART%20clamp_or_fail&type=code

clamp_or_fail not passed to state structure
No fail in clamp_variable

!-------------------------------------------------------------------------------
!> Check a variable for out of bounds and clamp or fail if needed.
!> If the variable has clamping limits, this routine returns .TRUE.
!> If the variable is unbounded, this routine returns .FALSE.
!> The return value is not an indication of whether or not the values have
!> actually been modified.
!-------------------------------------------------------------------------------
subroutine clamp_variable(dom_id, var_index, variable)
integer, intent(in) :: dom_id ! domain id
integer, intent(in) :: var_index ! variable index
real(r8), intent(inout) :: variable(:) ! variable
real(r8) :: minclamp, maxclamp, my_minmax(2)
character(len=NF90_MAX_NAME) :: varname ! for informational log messages
logical :: allow_missing ! used in CLM for state variables
! if neither bound is set, return early
minclamp = get_io_clamping_minval(dom_id, var_index)
maxclamp = get_io_clamping_maxval(dom_id, var_index)
if (minclamp == missing_r8 .and. maxclamp == missing_r8) return
! if we get here, either the min, max or both have a clamping value.
!>@todo this is what the code needs to be for CLM and any other
! model that allows missing values in the state. right now that
! is defined in assim_tools_mod but i don't think we can use it
! because of circular module dependencies. it should be defined
! maybe in filter? and set into some low level module (like types
! or constants or options_mod so anyone can query it).
!
! if we allow missing values in the state (which jeff has never
! liked because it makes the statistics funny), then these next
! two lines need to be:
allow_missing = get_missing_ok_status()
if (allow_missing) then
my_minmax(1) = minval(variable, mask=(variable /= missing_r8))
my_minmax(2) = maxval(variable, mask=(variable /= missing_r8))
else
! get the min/max for this variable before we start
my_minmax(1) = minval(variable)
my_minmax(2) = maxval(variable)
endif
varname = get_variable_name(dom_id, var_index)
! is lower bound set?
if ( minclamp /= missing_r8 ) then ! missing_r8 is flag for no clamping
if ( my_minmax(1) < minclamp ) then
!>@todo again, if we're allowing missing in state, this has to be masked:
if (allow_missing) then
where(variable /= missing_r8) variable = max(minclamp, variable)
else
variable = max(minclamp, variable)
endif
! TJH TOO VERBOSE write(msgstring, *) trim(varname)// ' lower bound ', minclamp, ' min value ', my_minmax(1)
! TJH TOO VERBOSE call error_handler(E_ALLMSG, 'clamp_variable', msgstring, &
! TJH TOO VERBOSE source)
endif
endif ! min range set
! is upper bound set?
if ( maxclamp /= missing_r8 ) then ! missing_r8 is flag for no clamping
if ( my_minmax(2) > maxclamp ) then
!>@todo again, if we're allowing missing in state, this has to be masked:
if (allow_missing) then
where(variable /= missing_r8) variable = min(maxclamp, variable)
else
variable = min(maxclamp, variable)
endif
! TJH TOO VERBOSE write(msgstring, *) trim(varname)// ' upper bound ', maxclamp, ' max value ', my_minmax(2)
! TJH TOO VERBOSE call error_handler(E_ALLMSG, 'clamp_variable', msgstring, &
! TJH TOO VERBOSE source)
endif
endif ! max range set
end subroutine clamp_variable
!-------------------------------------------------------------------------------

@hkershaw-brown
Copy link
Member Author

Todo Follow up on wrf 4 hybrid vertical coordinate

https://github.com/nancycollins/dart_dev_archive/commit/c28a6e82f0a7775f2075555ab84ea32668676adf

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)+wrf%dom(id)%c2h(k) + &
        wrf%dom(id)%c1h(k)*x_imu) / ph_e
else
   model_rho_t_distrib(:) = - (wrf%dom(id)%mub(i,j)+x_imu) / ph_e
endif

@nancycollins
Copy link
Collaborator

i've made a wrf_hybrid branch on my fork of the DART repo with the code changes for the hybrid coords in the wrf model_mod. it compiles. i have no input files to test with.

@hkershaw-brown
Copy link
Member Author

Hi @braczka I'm moving your comment to a new issue, running DART with wrf 4+

@NCAR NCAR deleted a comment from braczka Apr 5, 2024
@hkershaw-brown hkershaw-brown linked a pull request May 29, 2024 that will close this issue
15 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
missing_state missing_r8 in the state wrf Weather Research & Forecasting Model
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants