Feat req: Refactor wrf model_mod and models/wrf directory, unify WRF/WRF-CHEM

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.

DART/models/wrf/model_mod.f90

Lines 323 to 332 in 7ee9ed9

    
           ! 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:

DART/models/wrf/model_mod.f90

Lines 889 to 890 in 7ee9ed9

    
           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.

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

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)

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.

DART/models/wrf/model_mod.f90

Line 3030 in 70e6af8

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

issue for this bug bug: wrf convert_vertical_state does not use which_vert, #1038

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

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)

DART/models/wrf/model_mod.f90

Line 938 in 9729d78

! 10 = polar observation while restrict_polar namelist true

For VERTISLEVEL we have:

DART/models/wrf/model_mod.f90

Lines 895 to 900 in 9729d78

    
           ! 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

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

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.

Mask variables: XLAND, QTY_SURFACE_TYPE

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

DART/observations/forward_operators/obs_def_rttov13_mod.f90

Line 2065 in 9729d78

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

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

DART/models/wrf/model_mod.f90

Line 7092 in 9729d78

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

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]

Do we like the results?

DART/models/wrf/model_mod.f90

Line 5144 in 9729d78

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

DART/models/wrf/model_mod.f90

Line 5215 in 9729d78

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

Note on extrapolation 2D vs 3D: see discussion on #773, #768
2D fields (e.g. 2m temp) seem to be temperature rather than potential temperature, so be aware that testing that the code is (or is not) doing the pot_temp -> temp conversion appropriately.

NSF NCAR MMM statement "Given the extensive use of WRF in the U.S. and international atmospheric modeling communities, we are not planning to retire WRF or step away from its support and maintenance in the foreseeable future."
https://www.mmm.ucar.edu/about/wrf-mpas-support

In the namelist section, not actually in the namelist:

DART/models/wrf/model_mod.f90

Lines 229 to 231 in d908222

    
           logical :: log_vert_interp  = .true. 
        
           logical :: log_horz_interpM = .false. 
        
           logical :: log_horz_interpQ = .false.

support mete and chem output in different files?

415GB mete and chem:
/glade/derecho/scratch/hkershaw/DART/WRF-CHEM/29Aug2025

Cloud water units:
https://github.com/NCAR/DART/pull/896/files#diff-e12d14ed6b47c2ee908ccf4ad30e45a718ebd9e26297e77254a978e8804d7403R1430-R1435

Separate mete and chem files: 9091f62

This is just a note, on QTY_POTENTIAL_TEMPERATURE interpolate

This will give different answers at O2:

fld(1, e) = ts0 + dym*( dxmx_ill(e) + dxx_ilr(e)) + dy*( dxmx_iul(e) + dxx_iur(e) )

fld(1, e) = dym*( dxmx_ill(e) + dxx_ilr(e)) + dy*( dxmx_iul(e) + dxx_iur(e) ) + ts0

DART/models/wrf/model_mod.f90

Line 1824 in a32ab03

    
           fld(1, e) = ts0 + dym*( dxm*x_ill(e) + dx*x_ilr(e)) + dy*( dxm*x_iul(e) + dx*x_iur(e) )

This is just a note, on QTY_POTENTIAL_TEMPERATURE interpolate

This will give different answers at O2:

fld(1, e) = ts0 + dym*( dxm_x_ill(e) + dx_x_ilr(e)) + dy*( dxm_x_iul(e) + dx_x_iur(e) )

fld(1, e) = dym*( dxm_x_ill(e) + dx_x_ilr(e)) + dy*( dxm_x_iul(e) + dx_x_iur(e) ) + ts0

DART/models/wrf/model_mod.f90

Line 1824 in a32ab03

fld(1, e) = ts0 + dym*( dxmx_ill(e) + dxx_ilr(e)) + dy*( dxmx_iul(e) + dxx_iur(e) )

refactored uses a function for the interpolation, so you do not get bitwise differences between 00, 01, 02
t0 + func()
func() + t0
I think these interpolation lines are (one reason) why dart-wrf has need -fp model precise for bitwise testing. With the function the compiler has to materialize the function result before using it in a calculation.

Note wrf_unified currently not doing the case for observation in scaleheight, convert to a different vertical.
I don't think observations come in scaleheight, but it is in the wrf model mod.

WRF
model interpolate does:
horizontal interpolation on the quad, then vertical interpolation between two levels

vert_convert does:
vertical interpolation at each quad point, then horizontal interpolation on the quad

model_interpolate:

DART/models/wrf/model_mod.f90

Lines 1343 to 1346 in 4c89e64

    
           ! The big horizontal interp loop below computes the data values in the level 
        
           ! below and above the actual location, and then does a separate vertical 
        
           ! interpolation (if the obs is not a 2d surface obs).  The two values are 
        
           ! stored in fld(1:2).  Set them to missing here, and if the code below cannot

vert_convert:

DART/models/wrf/model_mod.f90

Lines 3445 to 3458 in 4c89e64

    
           ! compute height at all neighboring mass points and interpolate 
        
           hgta = model_height_w_distrib(ll(1), ll(2), k  ,id, state_handle) 
        
           hgtb = model_height_w_distrib(ll(1), ll(2), k+1,id, state_handle) 
        
           hgt1 = dzm*hgta + dz*hgtb 
        
           hgta = model_height_w_distrib(lr(1), lr(2), k  ,id, state_handle) 
        
           hgtb = model_height_w_distrib(lr(1), lr(2), k+1,id, state_handle) 
        
           hgt2 = dzm*hgta + dz*hgtb 
        
           hgta = model_height_w_distrib(ul(1), ul(2), k  ,id, state_handle) 
        
           hgtb = model_height_w_distrib(ul(1), ul(2), k+1,id, state_handle) 
        
           hgt3 = dzm*hgta + dz*hgtb 
        
           hgta = model_height_w_distrib(ur(1), ur(2), k  ,id, state_handle) 
        
           hgtb = model_height_w_distrib(ur(1), ur(2), k+1,id, state_handle) 
        
           hgt4 = dzm*hgta + dz*hgtb 
        
           zout = dym*( dxm*hgt1 + dx*hgt2 ) + dy*( dxm*hgt3 + dx*hgt4 )

fix-bitwise: vertical conversion of obs vertical then horiz interp

WRF has biliear interpolation all over the place: #860 (comment)
and then clamps to zero in model interpolate.
Refactor - not put bilinear interpolation as a function. Maybe worth doing this to switch out bilinear to idw in one step if needed.

Navigation Menu

Repository navigation

Feat req: Refactor wrf model_mod and models/wrf directory, unify WRF/WRF-CHEM #404

Use case

Is your feature request related to a problem?

Describe your preferred solution

Describe any alternatives you have considered

Vertical Convertion

39 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	! 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

	var_type = wrf%dom(id)%var_type(var_id)
	dart_type = wrf%dom(id)%dart_kind(var_id)

	! 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

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

	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

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

	! 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

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

	!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

	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( dxmreal(wrf%dom(id)%land(ll(1), ll(2))) + &
	dx*real(wrf%dom(id)%land(lr(1), lr(2))) ) + &
	dy( dxmreal(wrf%dom(id)%land(ul(1), ul(2))) + &
	dx*real(wrf%dom(id)%land(ur(1), ur(2))) ) - 1

	endif

	! 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( dxmreal(wrf%dom(id)%land(ll(1), ll(2))) + &
	dx*real(wrf%dom(id)%land(lr(1), lr(2))) ) + &
	dy( dxmreal(wrf%dom(id)%land(ul(1), ul(2))) + &
	dx*real(wrf%dom(id)%land(ur(1), ur(2))) )

	endif

	! 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

	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

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

	logical :: log_vert_interp = .true.
	logical :: log_horz_interpM = .false.
	logical :: log_horz_interpQ = .false.

	! The big horizontal interp loop below computes the data values in the level
	! below and above the actual location, and then does a separate vertical
	! interpolation (if the obs is not a 2d surface obs). The two values are
	! stored in fld(1:2). Set them to missing here, and if the code below cannot

	! compute height at all neighboring mass points and interpolate
	hgta = model_height_w_distrib(ll(1), ll(2), k ,id, state_handle)
	hgtb = model_height_w_distrib(ll(1), ll(2), k+1,id, state_handle)
	hgt1 = dzmhgta + dzhgtb
	hgta = model_height_w_distrib(lr(1), lr(2), k ,id, state_handle)
	hgtb = model_height_w_distrib(lr(1), lr(2), k+1,id, state_handle)
	hgt2 = dzmhgta + dzhgtb
	hgta = model_height_w_distrib(ul(1), ul(2), k ,id, state_handle)
	hgtb = model_height_w_distrib(ul(1), ul(2), k+1,id, state_handle)
	hgt3 = dzmhgta + dzhgtb
	hgta = model_height_w_distrib(ur(1), ur(2), k ,id, state_handle)
	hgtb = model_height_w_distrib(ur(1), ur(2), k+1,id, state_handle)
	hgt4 = dzmhgta + dzhgtb
	zout = dym( dxmhgt1 + dxhgt2 ) + dy( dxmhgt3 + dxhgt4 )

Feat req: Refactor wrf model_mod and models/wrf directory, unify WRF/WRF-CHEM #404

Description

Use case

Is your feature request related to a problem?

Describe your preferred solution

Describe any alternatives you have considered

Activity

hkershaw-brown commented on Oct 21, 2022

hkershaw-brown commented on May 26, 2023

hkershaw-brown commented on May 30, 2023

Vertical Convertion

hkershaw-brown commented on Aug 2, 2023

hkershaw-brown commented on Aug 9, 2023

hkershaw-brown commented on Aug 10, 2023

nancycollins commented on Aug 10, 2023

hkershaw-brown commented on Aug 14, 2023

hkershaw-brown commented on Aug 15, 2023

hkershaw-brown commented on Aug 15, 2023

hkershaw-brown commented on Aug 18, 2023

39 remaining items

hkershaw-brown commented on Nov 11, 2024

hkershaw-brown commented on Nov 17, 2024

hkershaw-brown commented on Jun 10, 2025

hkershaw-brown commented on Aug 29, 2025

hkershaw-brown commented on Sep 2, 2025

hkershaw-brown commented on Nov 26, 2025

hkershaw-brown commented on Dec 1, 2025

hkershaw-brown commented on Jan 9, 2026

hkershaw-brown commented on Jan 9, 2026

hkershaw-brown commented on Jan 12, 2026

hkershaw-brown commented on Jan 30, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions