search

NCAR / DART

Data Assimilation Research Testbed
https://dart.ucar.edu/
Apache License 2.0
197 stars 145 forks source link

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

Open hkershaw-brown opened 2 years ago

hkershaw-brown commented 2 years ago

Use case

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

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.

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 commented 2 years ago

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 commented 1 year ago

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 commented 1 year ago

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:

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

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

hkershaw-brown commented 1 year ago

To do:

hkershaw-brown commented 1 year ago

For VERTISLEVEL we have:

https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L895-L900

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

Edit: then this:

https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L1857-L1858

hkershaw-brown commented 1 year ago

note in the code about latlon project settings only valid for global domains https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L7307-L7316

nancycollins commented 1 year ago

for help with the acronyms:

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 commented 1 year ago

Mask variables: XLAND, QTY_SURFACE_TYPE

Should you be interpolating a mask?
edit: rounded to the nearest integer in rttov obs_def: https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/observations/forward_operators/obs_def_rttov13_mod.f90#L2065

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.

https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L2825-L2843

https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L2801-L2823

hkershaw-brown commented 1 year ago

I'm not sure what is not allowed here

https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L6738-L6740

There are some notes in the code about periodic_x,y and polar https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L7092

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

hkershaw-brown commented 1 year ago

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.

https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L7114-L7124

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: https://github.com/NCAR/DART/blob/70e6af803a52d14b9f77f872c94b1fe11d5dc2d9/models/wrf/model_mod.f90#L6275-L6277

hkershaw-brown commented 1 year ago

Do we like the results?

https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L5144

https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L5215

hkershaw-brown commented 1 year ago
hkershaw-brown commented 1 year ago
hkershaw-brown commented 1 year ago

Edit also a set of constants in wrf model_mod: https://github.com/NCAR/DART/blob/9729d784226295a197ca3bf00c917e4aaab5003b/models/wrf/model_mod.f90#L7997-L8006

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

hkershaw-brown commented 1 year ago
hkershaw-brown commented 1 year ago
hkershaw-brown commented 1 year ago

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 https://github.com/NCAR/DART/blob/6eeb4539f1d140c3184c2f9521a1b3bdf9d3396f/models/wrf/model_mod.f90#L3381 https://github.com/NCAR/DART/blob/6eeb4539f1d140c3184c2f9521a1b3bdf9d3396f/models/wrf/model_mod.f90#L5768

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
nancycollins commented 1 year ago

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 commented 1 year ago

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 commented 11 months ago

@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 commented 11 months ago

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 commented 10 months ago

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 commented 10 months ago

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 commented 10 months ago

https://github.com/NCAR/DART/blob/52e6e45f8f6ff07a90227766c8de72f1474162f2/models/wrf/model_mod.f90#L7485-L7491

hkershaw-brown commented 10 months ago

@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 https://github.com/NCAR/DART/issues/621.
Moved development to https://github.com/NCAR/DART/tree/wrf-wrf-chem (branched from v11 and squished commits) Todo:

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

hkershaw-brown commented 10 months ago

I don't think clamp_or_fail is used: https://github.com/NCAR/DART/blob/52e6e45f8f6ff07a90227766c8de72f1474162f2/models/wrf/model_mod.f90#L343

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

https://github.com/NCAR/DART/blob/52e6e45f8f6ff07a90227766c8de72f1474162f2/assimilation_code/modules/io/direct_netcdf_mod.f90#L1448-L1534

hkershaw-brown commented 10 months ago

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 commented 10 months ago

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 commented 7 months ago

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

hkershaw-brown commented 4 months ago
nancycollins commented 4 months ago
  • [ ] interpolate from 3D fields if 2M field is not in the state? 3D(:,:,1) == 2M ? , 3D(:,:,1) ~= 2M ?, 3D(:,:,1) != 2M ?

it can't be an interpolation, because the 3D fields are located at the vertical cell centers, not the bottom cell face. from what i remember from discussions at old wrf/dart meetings, they never recommend extrapolating but i can't justify why.

braczka commented 4 months ago

I don't know for sure but can imagine two reasons extrapolation would be poor: 1) extrapolating is generating a simple model outside the spatial range of calibration where it likely is not valid and also the temperature profile near the surface can change rapidly/non-linearly with height depending upon atmospheric stability, land/surface energy balance etc.

I have to imagine the way the WRF outputs the 2M and 10M met variables directly is much more sophisticated way to calculate the expected surface observation as compared to a simple extrapolation from the 3D field - but would have to look at WRF documentation to verify this.

hkershaw-brown commented 4 months ago

extrapolate is an option to all the vertical interpolation in the current code, e.g. https://github.com/NCAR/DART/blob/924182f29ce23f2a4e4c53aa57cc7d00fdccf374/models/wrf/model_mod.f90#L5082

braczka commented 4 months ago

extrapolate is an option to all the vertical interpolation in the current code, e.g.

https://github.com/NCAR/DART/blob/924182f29ce23f2a4e4c53aa57cc7d00fdccf374/models/wrf/model_mod.f90#L5082

I think extrapolation could/should be used for any observation that falls in between the bottom level of the model domain, but above 10 meters (surface observation). For a true surface observation (e.g. 10 m, 2 m), I feel like most atmospheric models would have this as standard output. I don't have good intuition yet of how large this gap is (i.e. height between bottom level of model and 10 meters.

nancycollins commented 4 months ago

this is a science question, not a software question and you're the atmospheric scientist, so i'll defer to you. after thinking about it for a while, i remember people talking about the atmospheric model "boundary layer" where below that height surface effects and turbulence are dominant. the models are not as accurate predicting what happens to the air there and below. there are lots of surface obs so the surface conditions can be more easily modeled than air slightly higher up.

extrapolation is there as a software option, but before recommending it be used by default, someone with a lot of low obs should look at the results with and without them. it may be that they are fine to use, or it may be they are different enough that the outlier threshold discards them, or they may make the assimilation results worse. i'm not sure how many "upper air" obs below the lowest level there actually are - maybe this is a moot point if there aren't any/many.

hkershaw-brown commented 1 week ago

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.

hkershaw-brown commented 1 week ago

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