Feature request: CAM SE Manhattan

[hkershaw-brown]

Use case

Manhattan release for CAM-SE

Is your feature request related to a problem?

• cam-se not released for Manhattan
  • 🐛 watch out for indices swap here: bug: cam-fv quad_vals  #246
• Duplicated code across cam-fv and cam-se

Describe your preferred solution

I believe this is the code that people have been using for cam-se using the cam-common-code:
https://github.com/jlaucar/DART/tree/CAM_SE

Brought up to date with main here:
https://github.com/NCAR/DART/tree/CAM_SE
Be aware: index swap #246

edit: 9/22/22 Just in case anyone wants to resurrect the CAM_SE branch here is copy:
https://github.com/hkershaw-brown/DART/tree/CAM_SE

Tidy (not changed the code): https://github.com/NCAR/DART/tree/cam-fv-se-shared f93ba7d

• model directories
  • cam-fv
  • cam-se
  • cam-common-code
• svn junk removal
• path_names
• obsolete file removal

CAM_SE...cam-fv-se-shared

[hkershaw-brown]

This shared module seems like it should be a submodule. I think then you do not have to faff about having copies of the namelist options in the shared module. Will finish the split, then try out a submodule version.

[hkershaw-brown]

precise_dry_mass_get_close is a namelist option (module global) but is passed as a subroutine argument

namelist /model_nml/  
   precise_dry_mass_get_close

There are two different cases when the use of precise vertical coordinates is relevant. One if for doing get_close computations, and the namelist flag controls this. Generally, it is okay to be less that precise for get_close because it's just doing a broad modulation of the impact. However, if there were a need for very tight vertical localization, it could have more impact. JLA

What is the goal of this hard coded precise?

subroutine model_interpolate:

! Should dry mass vertical coordinate be used to build vertical pressure columns? 
! This is expensive but probably necessary to get unbiased forward operators
logical :: precise = .true.

This second 'precise' is for doing the interpolation for forward operators. Here, precision can often be much more important because a lack of precision results in a biased vertical coordinate computation which then leads to a biased forward operator. I hard coded this here because I don't think it should be an option (unless there is enormous cost associated with it in some particular application). JLA

It is passed all the way down to here, find_se_vertical_levels. Hard coded switch off for dry_mass_vertical_coordinate?

The calling tree for get_close and for model_interpolated intersect down at low-levels, so the precise argument is needed. There are clearly some intermediate levels in the get_close computation where the namelist could still be used rather than the argument. While I was developing, I found it easier to just convert the namelist variable into an argument so I didn't have to worry about the details of where the hard-coded precise for interpolation 'merged' into the calling tree. JLA

find_se_vertical_levels:

if(dry_mass_vertical_coordinate .and. precise) then
   call build_dry_mass_pressure_columns(ens_handle, ens_size, ref_nlevels, column, surf_pressure, &
            pressure_array, status1)

[hkershaw-brown]

SE test cases:

Test case from svn trunk:

/glade/scratch/hkershaw/DART/CAM/test-cam-se-old-files

The netcdf files are from /glade/p/cisl/dares/DART_test_cases/cam/cam_se/test_trunk/nc_restarts/ (~2015)

--

Run of CAM SE with cesm 2.2.0:

/glade/scratch/hkershaw/DART/CAM/test-cam-se

The netcdf files are from /glade/scratch/hkershaw/cam-se.FHIST.ne30_g17/run/ (cesm 2.2.0)
Looking at a similar run but with cam-fv /glade/scratch/hkershaw/cam-fv.fhist.f19_g17/run (cesm 2.2.0). I think I'm missing outputting initial files (vs. restart files). The 2015 case is initial files I think.

Not sure about the cesm2.2.0 nenpnp vs ncol yet.

cam-se ne30_g17  cesm 2.2.0


netcdf caminput {
dimensions:
	lev = 32 ;
	ilev = 33 ;
	time = UNLIMITED ; // (1 currently)
	nenpnp = 86400 ;
	ncol = 48602 ;
	pbuf_00032 = 32 ;
	pbuf_00033 = 33 ;
	pbuf_00128 = 128 ;

        double T(time, lev, nenpnp) ;

old file:


netcdf caminput {
dimensions:
	ncol = 48602 ;
	time = UNLIMITED ; // (1 currently)
	nbnd = 2 ;
	chars = 8 ;
	lev = 30 ;
	ilev = 31 ;

        double T(time, lev, ncol) ;
                T:mdims = 1 ;
                T:units = "K" ;
                T:long_name = "Temperature" ;

Edit: running cam-se cesm2.2.0 but outputting initial files gives the familiar:

netcdf cam-se.FHIST.ne30_g17.cam.i.1979-01-06-00000 {
dimensions:
        ncol = 48602 ;
        time = UNLIMITED ; // (1 currently)
        nbnd = 2 ;
        chars = 8 ;
        lev = 32 ;
        ilev = 33 ;

 double T(time, lev, ncol) ;

To output initial files add inithist = 'ENDOFRUN' to user_nl_cam:

cat user_nl_cam 
! Users should add all user specific namelist changes below in the form of 
! namelist_var = new_namelist_value 
 inithist               = 'ENDOFRUN' 

edit 2:

To output the SEMapping.nc file

se_write_gll_corners = .true.

[hkershaw-brown]

model_interpolate( obs_qty )
   -> coord_ind_cs( obs_qty )
       -> get_close( obs_qty )  

In get_close there is a check on the type:

      if (base_type < 1 .or. base_type > size(gc%type_to_cutoff_map)) then

but in model_interpolate, you don't have type you have only quantity

If number of qts > number of types you get:

 Before computing prior observation values TIME: 2022/04/15 09:47:45
 ERROR FROM:
  source : threed_sphere/location_mod.f90
  routine: get_close
  message: base_type out of range, is      328
  message: ... must be between        1      95
  message: ... Lon/Lat(deg): 243.03000000  33.74700000  Vert:  280.9000000 hPa

If number of qts < number of types, then you will get through this if statement without exiting, but then I think your distance calculation can be affected by per_type_vert_norm(loc1%which_vert, type1) in get_dist.

[nancycollins]

there is a get_close_obs() that includes an obs type as one of the args. there is also a get_close_state() that includes a quantity and an index into the model state as part of the arg list. not sure which one is wanted in this context.

[hkershaw-brown]

@nancycollins so it doesn't matter, get_close is using the type to cuttoff map which is initialized in get_close_init

DART/assimilation_code/location/threed_sphere/location_mod.f90

Line 1458 in 2d3f153

if (base_type < 1 .or. base_type > size(gc%type_to_cutoff_map)) then

There is a note (FIXME Nancy) in the model code, about using XYZ location_mod:

DART/models/cam-se/model_mod.f90

Line 2846 in 24814a5

! FIXME: Nancy has a location_xyz:find_closest_???? which will return the N closest points,

Is using the XYZ module the way to go? The goal is to find the corners of the cell which enclose a location.

I intentionally removed use of the xyz get_close from threed_cartesian in the update. This simplifies code and handles the concerns about precise computation. The comment is left for a future exploration of computational cost. JLA

[hkershaw-brown]

Note on approximate_distance = .true. vs. approximate_distance = .false.

approximate_distance = .true.
one close points, distance = 0. Can it really be zero? (un)lucky observation location?

approximate_distance = .false.
five close points

Edit: separating this out into a new issue #346 Can test this independently of CAM

[kdraeder]

Identity obs might have 0 distance. Kevin

…

On Thu, Apr 21, 2022 at 9:59 AM hkershaw-brown ***@***.***> wrote: Note on approximate_distance = .true. vs. approximate_distance = .false. approximate_distance = .true. one close points, distance = 0. Can it really be zero? (un)lucky observation location? [image: Screen Shot 2022-04-21 at 9 13 10 AM] <https://user-images.githubusercontent.com/20047007/164501696-07074e67-b0da-43b6-9b61-5ebfb9ce1532.png> approximate_distance = .false. five close points [image: Screen Shot 2022-04-21 at 9 50 22 AM] <https://user-images.githubusercontent.com/20047007/164501783-be3f4c52-e548-483a-b042-1037459241d4.png> — Reply to this email directly, view it on GitHub <#343 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADFFGEGCOTITZQFMDQVTQODVGF3MHANCNFSM5SWL7HWQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[hkershaw-brown]

good suggestion Kevin, the killer is the two runs above are with the same obs_seq.out (one obs).

[kdraeder]

Ah, sorry that I didn't read all the context first!

…

On Thu, Apr 21, 2022 at 10:47 AM hkershaw-brown ***@***.***> wrote: good suggestion Kevin, the killer is the two runs above are with the same obs_seq.out (one obs). — Reply to this email directly, view it on GitHub <#343 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADFFGEEY4NTI3V5JMRPQECTVGGBDFANCNFSM5SWL7HWQ> . You are receiving this because you commented.Message ID: ***@***.***>

[hkershaw-brown]

Never be sorry about chiming in. Sometimes without context it the best way to view things.

[kdraeder]

About cam_common_code_mod.f90

!>@todo FIXME ask kevin if this threshold value is small enough
! to distinguish cam from waccm configurations?
real(r8), parameter :: high_top_threshold = 0.3_r8  ! pascals

All of the CAMs up to CAM6 have a p_top ~ 340 Pa.
Future CAMs will have higher tops, but 0.3 Pa should be plenty high.