MPAS refactor notes #753

hkershaw-brown · 2024-10-11T19:52:08Z

Unresolved issues:

Surface elevation check
! Reject obs if the station height is far way from the model terrain.

but fwd operators (rttov) need these surface values
see discussion #332

Code:

branch (be aware I may force push to this)
https://github.com/NCAR/DART/tree/mpas-refactor

Spec (note 7 v 8)
https://docs.google.com/document/d/1wD4MKM5PquOUbxK3YxIjcbIaFThJfa318heDn_bFUbs/edit

The text was updated successfully, but these errors were encountered:

hkershaw-brown · 2024-10-25T19:15:00Z

In compute_u_with_rbf there is a loop around ens_size, for get_reconstruct but uval(1:ens_size) is set to set to the single value ureconstructzonal/ureconstructmeridional so you are getting the whole ensemble set each step of the loop

whole ensemble == last value from the loop.

DART/models/mpas_atm/model_mod.f90

Lines 6214 to 6226 in e218864

    
           do e = 1, ens_size 
        
              call get_reconstruct(nedges, lat*deg2rad, lon*deg2rad, & 
        
                          coeffs_reconstruct, on_a_sphere, veldata(:, e), & 
        
                          ureconstructx, ureconstructy, ureconstructz, & 
        
                          ureconstructzonal, ureconstructmeridional) 
        
              if (zonal) then 
        
                 uval = ureconstructzonal 
        
              else 
        
                 uval = ureconstructmeridional 
        
              endif 
        
           enddo

hkershaw-brown · 2024-10-25T19:49:04Z

todo can be different across the ensemble

DART/models/mpas_atm/model_mod.f90

Lines 6183 to 6192 in e218864

    
           !>@todo lower (upper) could be different levels in pressure 
        
           !lowval = x(index1 + (edgelist(i)-1) * nvert + lower(i)-1) 
        
           lowindx = int(index1,i8) + int((edgelist(i)-1) * nvert,i8) + int(lower(i,1)-1,i8) 
        
           lowval =  get_state(lowindx, state_handle) 
        
           !uppval = x(index1 + (edgelist(i)-1) * nvert + upper(i)-1) 
        
           upindx = int(index1,i8) + int((edgelist(i)-1) * nvert,i8) + int(upper(i,1)-1,i8) 
        
           uppval =  get_state(upindx, state_handle) 
        
           veldata(i, :) = lowval*(1.0_r8 - fract(i, :)) + uppval*fract(i, :)

hkershaw-brown · 2024-10-28T20:28:49Z

on_bound? not used

DART/models/mpas_atm/model_mod.f90

Line 4871 in e218864

logical :: at_surf, do_norm, on_bound

DART/models/mpas_atm/model_mod.f90

Lines 4930 to 4938 in e218864

    
           ! if the entire ensemble has missing vertical values we can return now. 
        
           ! otherwise we need to continue to convert the members with good vertical values. 
        
           ! boundary cells will be updated by the assimilation. 
        
           ! if all the vertical localization coord values are missing,  
        
           ! we don't call this routine again, and return. 
        
           if (all(zin == missing_r8)) then ! .or. on_bound) then 
        
              istatus(:) = 0 
        
              return 
        
           endif

hkershaw-brown · 2024-10-28T21:55:41Z

find_vert_indices is a wrap around find_vert_level. Remove

DART/models/mpas_atm/model_mod.f90

Line 6030 in e218864

    
           subroutine find_vert_indices (state_handle, ens_size, loc, nc, c, lower, upper, fract, ier)

hkershaw-brown · 2024-10-28T22:13:10Z

convert_vert_distrib_state (state)
convert_vert_distrib (obs)

super repetitive & according to this comment state is calling the interpolation anyway:

DART/models/mpas_atm/model_mod.f90

Lines 4968 to 4972 in e218864

    
           !>@todo FIXME - this is the original code from the 
        
           !> observation version which does horizontal interpolation. 
        
           !> in this code we know we are on a state vector location 
        
           !> so no interp is needed.  we should be able to make this 
        
           !> more computationally efficient.

hkershaw-brown · 2024-10-28T22:19:25Z

If rbf is depricated 11 years ago, possibly it is time to remove the code?

DART/models/mpas_atm/model_mod.f90

Line 119 in e218864

! RBF (radial basis function) modules, donated by LANL. currently deprecated

note see comment on incorrect rbf calulcaiton.

Is the rbf code depricated #753 (comment)

hkershaw-brown · 2024-10-29T12:46:52Z

Shells scripts are out of date, user complaints about cycling. Lots of csh (poor).

Documentation is poor, also points people to this MPAS page which sends people to none existent documentation (https://svn-dares-dart.cgd.ucar.edu/)

hkershaw-brown · 2024-10-29T13:04:15Z

Nope:

DART/models/mpas_atm/model_mod.f90

Line 188 in e218864

! try adjusting what static_init_model does, before it is called.

DART/models/mpas_atm/model_mod.f90

Line 285 in e218864

! if set by: call force_u_into_state() *BEFORE* calling static_init_model(),

DART/models/mpas_atm/model_mod.f90

Lines 551 to 552 in e218864

    
           subroutine static_init_model() 
        
           !>@todo FIXME - can we add an optional arg here for update_bc use?

Note, beware of whatever the regional stuff is doing. Seems like some bizarre choices here.

DART/models/mpas_atm/model_mod.f90

Lines 2097 to 2098 in e218864

    
           ! modify what static_init_model does.  this *must* be called before 
        
           ! calling static_init_model().

more nope:

DART/models/mpas_atm/update_bc.f90

Lines 100 to 102 in e218864

    
           ! Note that force_u_into_state should be called before static_init_model, which is unusual. 
        
           call force_u_into_state() 
        
           call set_lbc_variables(bdy_template_filename)

hkershaw-brown · 2024-10-29T13:09:37Z

update_bc.f90 is calling use direct_netcdf_mod,only : read_variables directly. 🙃

hkershaw-brown · 2024-10-29T13:12:15Z

why is this compile time, rather than runtime:

DART/models/mpas_atm/model_mod.f90

Lines 238 to 243 in e218864

    
           ! compile-time control over whether grid information is written to the 
        
           ! diagnostic files or not.  if it is, the files are self-contained (e.g. for 
        
           ! ease of plotting), but are also much larger than they would be otherwise. 
        
           ! change this private variable to control whether the grid information is 
        
           ! written or not. 
        
           logical :: add_static_data_to_diags = .false.

hkershaw-brown · 2024-10-29T13:17:23Z

why is this compile time, rather than runtime:

DART/models/mpas_atm/model_mod.f90

Lines 238 to 243 in e218864

! compile-time control over whether grid information is written to the

! diagnostic files or not. if it is, the files are self-contained (e.g. for

! ease of plotting), but are also much larger than they would be otherwise.

! change this private variable to control whether the grid information is

! written or not.

logical :: add_static_data_to_diags = .false.

oh add_static_data_to_diags is not used and there is also a compile time (namelist) option write_grid_to_diag_file

hkershaw-brown · 2024-10-29T13:20:12Z

! FIXME: it may be desirable to read in xCell(:), yCell(:), zCell(:)
! to keep from having to compute them on demand, especially since we
! have converted the radian lat/lon of the cell centers into degrees.
! we have to convert back, then take a few sin and cos to get xyz.
! time/space/accuracy tradeoff here.

https://github.com/NCAR/DART/blame/e2188646b97573f198c66f7ea63482ebc34afa11/models/mpas_atm/model_mod.f90#L416

hkershaw-brown · 2024-10-29T13:21:20Z

model size needs to be i8:

DART/models/mpas_atm/model_mod.f90

Line 455 in e218864

integer :: model_size ! the state vector length

hkershaw-brown · 2024-10-29T13:30:10Z

currently unused (11 years). Not needed for regional? Or should be used?

DART/models/mpas_atm/model_mod.f90

Lines 477 to 480 in e218864

    
           ! currently unused; for a regional model it is going to be necessary to know 
        
           ! if the grid is continuous around longitudes (wraps in east-west) or not, 
        
           ! and if it covers either of the poles. 
        
           character(len= 64) :: ew_boundary_type, ns_boundary_type

hkershaw-brown · 2024-10-29T13:34:22Z

This data on edge logic seems overcomplicated

DART/models/mpas_atm/model_mod.f90

Lines 464 to 475 in e218864

    
           ! Do we have any state vector items located on the cell edges? 
        
           ! If not, avoid reading in or using the edge arrays to save space. 
        
           ! FIXME: this should be set after looking at the fields listed in the 
        
           ! namelist which are to be read into the state vector - if any of them 
        
           ! are located on the edges then this flag should be changed to .true. 
        
           ! however, the way the code is structured these arrays are allocated 
        
           ! before the details of field list is examined.  since right now the 
        
           ! only possible field array that is on the edges is the 'u' edge normal 
        
           ! winds, search specifically for that in the state field list and set 
        
           ! this based on that.  if any other data might be on edges, this routine 
        
           ! will need to be updated: is_edgedata_in_state_vector() 
        
           logical :: data_on_edges = .false.

variable dimension tells you edge vs cell, .e.g

float u(Time, nEdges, nVertLevels) ;
                u:units = "m s^{-1}" ;
                u:long_name = "Horizontal normal velocity at edges" ;

float qs(Time, nCells, nVertLevels) ;
                qs:units = "kg kg^{-1}" ;
                qs:long_name = "Snow mixing ratio" ;

hkershaw-brown · 2024-10-29T13:37:59Z

"causes mpas namelists to be read" - does not appear to be true

DART/models/mpas_atm/model_mod.f90

Lines 603 to 606 in e218864

    
           ! Set the time step ... causes mpas namelists to be read. 
        
           ! Ensures model_timestep is multiple of 'dynamics_timestep' 
        
           call set_calendar_type( calendar )   ! comes from model_mod_nml

hkershaw-brown · 2024-10-29T13:43:05Z

This comment, not sure what it is getting at:

DART/models/mpas_atm/model_mod.f90

Lines 843 to 845 in e218864

    
           !>@todo FIXME: if we want to do increments, we could also add a 
        
           ! third domain which is the original forecast fields before 
        
           ! the assimilation (so we can compute increments)

hkershaw-brown · 2024-10-29T13:48:41Z

boundary domain variables do not get bounds for clamping:

DART/models/mpas_atm/model_mod.f90

Lines 856 to 858 in e218864

    
           lbc_domid = add_domain(bdy_template_filename, lbc_nfields,    & 
        
                                  var_names  = lbc_variables) 
        
                                ! FIXME clamp_vals = variable_bounds(1:nfields,:) )

see comment #753 (comment) Remove not used variables: add_static_data_to_diags - not used ew_boundary_type, ns_boundary_type -not used Shortened some very verbose comments.

hkershaw-brown · 2024-10-29T14:22:40Z

query_vert_localization_coord not used & not requried.

DART/models/mpas_atm/model_mod.f90

Lines 7481 to 7488 in e218864

    
           !> pass the vertical localization coordinate to assim_tools_mod 
        
           function query_vert_localization_coord() 
        
           integer :: query_vert_localization_coord 
        
           query_vert_localization_coord = vert_localization_coord 
        
           end function query_vert_localization_coord

#753 (comment) removed done todo comment

hkershaw-brown · 2024-10-29T16:32:41Z

get_init_template_filename

DART/models/mpas_atm/model_mod.f90

Line 169 in e218864

public :: get_init_template_filename, &

called but value unused in update_bc ?

DART/models/mpas_atm/update_bc.f90

Line 105 in e218864

call get_init_template_filename(static_filename)

hkershaw-brown · 2024-10-29T16:36:47Z

Hardcoded domain id:

DART/models/mpas_atm/update_bc.f90

Lines 110 to 113 in e218864

    
           ! use get_num_variables() after setting the domains 
        
           ! separate number for analysis file, boundary file 
        
           nanlvars = get_num_variables(1) 
        
           nbdyvars = get_num_variables(2)

hkershaw-brown · 2024-10-29T16:37:41Z

serial loop around files for update bc:

DART/models/mpas_atm/update_bc.f90

Line 122 in e218864

fileloop: do ! until out of files

hkershaw-brown · 2024-10-29T18:09:58Z

write model time overloaded

DART/models/mpas_atm/model_mod.f90

Lines 506 to 509 in e218864

    
           interface write_model_time 
        
              module procedure write_model_time_file 
        
              module procedure write_model_time_restart 
        
           end interface

hkershaw-brown · 2024-10-29T18:11:31Z

is_global_grid

Only used by mpas_dart_obs_preprocess, then global does not appear to be used

DART/models/mpas_atm/mpas_dart_obs_preprocess.f90

Line 1013 in e218864

global = is_global_grid()

hkershaw-brown · 2024-10-29T18:17:07Z

Hardcoded dom_id in update_mpas_states:

DART/models/mpas_atm/update_mpas_states.f90

Line 54 in e218864

integer :: dom_id = 1

hkershaw-brown · 2024-10-29T18:18:56Z

serial loop around files for update_mpas_states

DART/models/mpas_atm/update_mpas_states.f90

Lines 92 to 96 in e218864

    
           !---------------------------------------------------------------------- 
        
           ! Reads lists of input mpas (prior) and filter (analysis) files  
        
           !---------------------------------------------------------------------- 
        
           filenum = 1 
        
           fileloop: do        ! until out of files

no dependencies in the loop? - run multiple update_mpas_states in parallel.

hkershaw-brown · 2024-10-29T18:49:02Z

progvar is removed and replaced with state structure. Notes on done/not done routines

I'm concerned about update_bc and update_mpas_states hard coding domains and using read_variables directory.
Also,

memory wise - reading the whole state in
- clamping - this is already done on writing out the state from filter. Is this additional? (or unneeded) Would you not clamp from filter, but clamp here? Do you need the fail if so?
- updating winds, you only need to read in the winds, not the whole state
performance wise, looping around each restart file in update_mpas_states

hkershaw-brown · 2024-10-29T20:07:44Z

todo docs:

remove progvar from docs.
whatever is going on with clamping.
It would be good to explain "The wind options during a DART assimilation are controlled by combinations of 4 different namelist values." what the 4 are and what the combinations are.

hkershaw-brown · 2024-10-29T20:44:45Z

oncenters is not binary (3 choices)

find_vert_indices -> find_vert_level( ..., oncenters=.true., ...)

the logic for oncenters is a bit weird to be passing all the way down the call tree.

Here it is hardcoded .true. for pressure otherwise .false.

DART/models/mpas_atm/model_mod.f90

Lines 6136 to 6157 in e218864

    
           if (verttype == VERTISPRESSURE) then 
        
              ! get all cells which share any edges with the 3 cells which share 
        
              ! the closest vertex. 
        
              call make_cell_list(vertexid, 3, ncells, celllist) 
        
              call find_vert_level(state_handle, ens_size, loc, ncells, celllist, .true., & 
        
                                   lower, upper, fract, ier) 
        
              if (all(ier /= 0)) return 
        
              ! now have pressure at all cell centers - need to interp to get pressure 
        
              ! at edge centers. 
        
              call move_pressure_to_edges(ncells, celllist, lower, upper, fract, & 
        
                                          nedges, edgelist, ier) 
        
              if (all(ier /= 0)) return 
        
           else 
        
              ! need vert index for the vertical level 
        
              call find_vert_level(state_handle, ens_size, loc, nedges, edgelist, .false., & 
        
                                   lower, upper, fract, ier) 
        
              if (all(ier /= 0)) return 
        
           endif

but oncenters ignored for pressure:

DART/models/mpas_atm/model_mod.f90

Lines 5374 to 5406 in e218864

    
           if(verttype == VERTISPRESSURE ) then 
        
              track_ier = 0 
        
              vert_array = vert 
        
              ! Need to get base offsets for the potential temperature, density, and water 
        
              ! vapor mixing fields in the state vector 
        
              call get_index_range(QTY_POTENTIAL_TEMPERATURE, pt_base_offset) 
        
              call get_index_range(QTY_DENSITY, density_base_offset) 
        
              call get_index_range(QTY_VAPOR_MIXING_RATIO, qv_base_offset) 
        
              do i=1, nc 
        
                 call find_pressure_bounds(state_handle, ens_size, vert_array, ids(i), nVertLevels, & 
        
                       pt_base_offset, density_base_offset, qv_base_offset,  & 
        
                       lower(i, :), upper(i, :), fract(i, :), ier) 
        
                 !if(debug > 9) print '(A,5I5,F10.4)', & 
        
                 !   '    after find_pressure_bounds: ier, i, cellid, lower, upper, fract = ', & 
        
                 !      ier, i, ids(i), lower(i, :), upper(i, :), fract(i, :) 
        
                 !if(debug > 5 .and. ier(e) /= 0) print '(A,4I5,F10.4,I3,F10.4)', & 
        
                 !'fail in find_pressure_bounds: ier, nc, i, id, vert, lower, fract: ', & 
        
                 !   ier, nc, i, ids(i), vert_array, lower(i, :), fract(i, :) 
        
                 ! we are inside a loop over each corner. consolidate error codes 
        
                 ! so that we return an error for that ensemble member if any 
        
                 ! of the corners fails the pressure bounds test. 
        
                 where (ier /= 0 .and. track_ier == 0) track_ier = ier 
        
              enddo 
        
              ier = track_ier 
        
              return 
        
           endif

oncenters(/faces/edges) only relevant for HEIGHT?
oncenters/faces/edges Can be queried from the obs_kind.

see

hkershaw-brown · 2024-12-04T20:12:16Z

related to #757 "VERTISSURFACE describes variables on A surface (not necessarily THE surface)"

DART/models/mpas_atm/model_mod.f90

Line 5336 in 2d5882c

! VERTISSURFACE describes variables on A surface (not necessarily THE surface)

hkershaw-brown self-assigned this Oct 11, 2024

hkershaw-brown added a commit that referenced this issue Oct 28, 2024

uncommented find_surrounding_edges

8d718c4

Is the rbf code depricated #753 (comment)

hkershaw-brown added the mpas Model for Prediction Across Scales (MPAS) label Oct 29, 2024

hkershaw-brown added a commit that referenced this issue Oct 29, 2024

chore: removed unused routine query_vert_localization_coord

268f5af

#753 (comment) removed done todo comment

hkershaw-brown mentioned this issue Dec 9, 2024

mpas_atm refactor to use state structure module #780

Draft

15 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

MPAS refactor notes #753

MPAS refactor notes #753

hkershaw-brown commented Oct 11, 2024

hkershaw-brown commented Oct 25, 2024

hkershaw-brown commented Oct 25, 2024

hkershaw-brown commented Oct 28, 2024

hkershaw-brown commented Oct 28, 2024

hkershaw-brown commented Oct 28, 2024

hkershaw-brown commented Oct 28, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024 •

edited

Loading

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024 •

edited

Loading

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024 •

edited

Loading

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Dec 4, 2024

MPAS refactor notes #753

MPAS refactor notes #753

Comments

hkershaw-brown commented Oct 11, 2024

Unresolved issues:

Code:

hkershaw-brown commented Oct 25, 2024

hkershaw-brown commented Oct 25, 2024

hkershaw-brown commented Oct 28, 2024

hkershaw-brown commented Oct 28, 2024

hkershaw-brown commented Oct 28, 2024

hkershaw-brown commented Oct 28, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024 • edited Loading

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024 • edited Loading

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Oct 29, 2024 • edited Loading

hkershaw-brown commented Oct 29, 2024

hkershaw-brown commented Dec 4, 2024

hkershaw-brown commented Oct 29, 2024 •

edited

Loading

hkershaw-brown commented Oct 29, 2024 •

edited

Loading

hkershaw-brown commented Oct 29, 2024 •

edited

Loading