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.

Sign up for GitHub

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

Jump to bottom

Feature request: MOM6 interface #451

Closed
4 tasks done
hkershaw-brown opened this issue Feb 6, 2023 · 25 comments
Closed
4 tasks done

Feature request: MOM6 interface #451

hkershaw-brown opened this issue Feb 6, 2023 · 25 comments
Assignees
@hkershaw-brown
Labels
mom6 Modular Ocean Model

Comments

@hkershaw-brown
Copy link
Member

hkershaw-brown commented Feb 6, 2023 •
edited
Loading

Use case

Data assimilation with MOM6 and DART

Is your feature request related to a problem?

Currently no interface for MOM6 for dart.

Describe your preferred solution

To do:

  • MOM6 model_mod
  • run MOM6 in cesm
  • test nuopc cesm with DART
  • scripting for DART-MOM6 run in CESM

Describe any alternatives you have considered

There is some work here https://github.com/amrhein/DART/tree/mom6 where Dan and I looked at starting with the POP model_mod. POP has the KMU KMT files for bathymetry and uses its own structure for finding the enclosing grid point (vs. the new quad_interp_mod). MOM6 has variables on T,U,V.

Here is the branch for the latest work on MOM6 model_mod.
https://github.com/NCAR/DART/tree/mom6-cesm released in v10.7

Note restart files we have, have 1d lon, lat vs. 2d lon,lat. Watch out for this, we may need to have two versions of the quad interp calls. The 1d lon,lat are nominal. The 2d lon,lats are in the geolon variables.
@hkershaw-brown hkershaw-brown self-assigned this Feb 6, 2023
@hkershaw-brown
Copy link
Member Author

Q. For vertical location should we use Layer pseudo-depth?

 double Layer(Layer) ;
                Layer:long_name = "Layer pseudo-depth, -z*" ;
                Layer:units = "meter" ;
                Layer:axis = "Z" ;
                Layer:positive = "up" ;

Q. Is Layer pseudo-depth per ensemble member?
Q. Layer pseudo-depth is only 1D (Layer), h Layer Thickness is 3D (Layer, lath, lonh). What is going on here?

double h(Time, Layer, lath, lonh) ;
                h:long_name = "Layer Thickness" ;
                h:units = "m" ;
                h:checksum = "B012521BB74AB9DA" ;

@jlaucar
Copy link
Contributor

jlaucar commented Feb 15, 2023 via email

@amrhein
Copy link
Collaborator

amrhein commented Feb 15, 2023 via email

hkershaw-brown added a commit that referenced this issue Feb 17, 2023
@hkershaw-brown
find level bounds for pseudo-depth using 1D interf
71eabd1 
ignoring per ensemble member pseudo-depths for now, see #451
for discussion on depth
@hkershaw-brown
Copy link
Member Author

I think this is similar to what goes on with eta (or hybrid eta/pressure) coordinates in the atmosphere. The actual depth in meters of a given grid point can vary with other aspects of the state. In atmospheric models, this is normally the surface pressure. We probably need to talk with the modelers to understand this better. As I recall, our response in atmospheric models has been to use the ensemble mean to do the computations of vertical location?

just for completeness, after chatting with Dan and Alper,
Yes inside filter_assim, get_close vertical calculations are with the ensemble mean. The subtlety is for the forward operator, e.g. can an observation at depth X be in level 3 for one ensemble member and level 4 for another.

@jlaucar
Copy link
Contributor

jlaucar commented Feb 17, 2023 via email

@hkershaw-brown
Copy link
Member Author

yup you are correct.
For now, I'm got the 1D pseudo-depth being used (so all the same).

@jlaucar
Copy link
Contributor

jlaucar commented Feb 17, 2023 via email

@hkershaw-brown
Copy link
Member Author

yeah it is kind of pain with the vectorized forward operators.
There is get_state_array so you can grab a different state index across the ensemble with one call.

DART/assimilation_code/modules/utilities/distributed_state_mod.f90

Lines 35 to 39 in 224698d

!---------------------------------------------------------
!> Gets all copies of an element of the state vector from the process who owns it
!> Assumes ensemble complete. This differes from get_state as it now works on an
!> array of state indices rather than a single index.
subroutine get_state_array(x, my_index, state_ens_handle)

@nancycollins
Copy link
Collaborator

re: forwards operators and different vertical levels in different ensemble members

the wrf model_mod has this code already (written by helen originally) and mpas has it, too.

@hkershaw-brown
Copy link
Member Author

For masking out land.

ocean_geometry.nc (example in /glade/scratch/hkershaw/c.e22.CMOM.T62_t061.nuopc.001)

double wet(lath, lonh) ;
               wet:long_name = "land or ocean?" ;
               wet:units = "nondim" ;

double D(lath, lonh) ;
                D:long_name = "Basin Depth" ;
                D:units = "meter" ;

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Feb 23, 2023 •
edited
Loading

MOM6 does a check checksums for the variables in restart files.

Example, here is me filling 'Temp' with junk:

ncap2 -s 'Temp(:,:,:)=-999.9' c.e22.CMOM.T62_t061.nuopc.001.mom6.r.0001-01-11-00000.nc temp.nc
mv temp.nc c.e22.CMOM.T62_t061.nuopc.001.mom6.r.0001-01-11-00000.nc

This is the error you get ./case.submit

cesm.log.8706019.chadmin1.ib0.cheyenne.ucar.edu.230223-070401

72:FATAL from PE     0: MOM_restart(restore_state): Checksum of input field Temp 10311FFFFFCEF0C8 does not match value 9DFB15BB7F911B95 store
d in ./c.e22.CMOM.T62_t061.nuopc.001.mom6.r.0001-01-11-00000.nc

DART is going to change the variables, so the checksums aren't going to match. What is the best thing to do here? Switch off the checking checksums? Update the checksum attribute so the restart file is consistent? Remove the checksum attribute from the restart file?

I think I should be able to switch of the checksum check by namelist in user_nml_mom
but I'm missing something.

checksum_required = .false. Edit: the option is:

RESTART_CHECKSUMS_REQUIRED = False

@hkershaw-brown
Copy link
Member Author

Hi @alperaltuntas
I'm using this issue to keep a note of dart-mom6 questions as we hit them
Cheers,
Helen

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Mar 6, 2023 •
edited
Loading

Advice from Alper on the restart files appearing to be variable(lat,lon):

Even though it seems like the outputs are on a lat-lon grid, they actually are on the native grid.

To get the native coordinates, you should read in geolat and geolon (or geolatb or geolonb depending on the staggering of a field). These are in the ocean_geometry.nc output file in rundir.

ocean_geometry.nc

        double geolatb(latq, lonq) ;
                geolatb:long_name = "latitude at corner (Bu) points" ;
                geolatb:units = "degree" ;
        double geolonb(latq, lonq) ;
                geolonb:long_name = "longitude at corner (Bu) points" ;
                geolonb:units = "degree" ;
        double geolat(lath, lonh) ;
                geolat:long_name = "latitude at tracer (T) points" ;
                geolat:units = "degree" ;
        double geolon(lath, lonh) ;
                geolon:long_name = "longitude at tracer (T) points" ;
                geolon:units = "degree" ;

At some point hopefully the MOM6 grid docs get fleshed out: https://mom6.readthedocs.io/en/main/grids.html

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Mar 6, 2023 •
edited
Loading

this would be cooler if it had the locations of the u and v grid points

geolatu(lath, lonq) ?
geolatv(latq, lonh) ?

Edit: just for kicks, plotting (geolon, geolat)
Screen Shot 2023-03-06 at 5 52 32 PM

@hkershaw-brown
Copy link
Member Author

boom 💥 c.e22.CMOM.T62_t061.nuopc.001.mom6.static.nc is what you want.

example: /glade/scratch/hkershaw/c.e22.CMOM.T62_t061.nuopc.001/run/c.e22.CMOM.T62_t061.nuopc.001.mom6.static.nc

        float geolon(yh, xh) ;
                geolon:long_name = "Longitude of tracer (T) points" ;
                geolon:units = "degrees_east" ;
                geolon:missing_value = 1.e+20f ;
                geolon:_FillValue = 1.e+20f ;
                geolon:cell_methods = "time: point" ;
        float geolat(yh, xh) ;
                geolat:long_name = "Latitude of tracer (T) points" ;
                geolat:units = "degrees_north" ;
                geolat:missing_value = 1.e+20f ;
                geolat:_FillValue = 1.e+20f ;
                geolat:cell_methods = "time: point" ;
...
        float geolon_u(yh, xq) ;
                geolon_u:long_name = "Longitude of zonal velocity (Cu) points" ;
                geolon_u:units = "degrees_east" ;
                geolon_u:missing_value = 1.e+20f ;
                geolon_u:_FillValue = 1.e+20f ;
                geolon_u:cell_methods = "time: point" ;
                geolon_u:interp_method = "none" ;
        float geolat_u(yh, xq) ;
                geolat_u:long_name = "Latitude of zonal velocity (Cu) points" ;
                geolat_u:units = "degrees_north" ;
                geolat_u:missing_value = 1.e+20f ;
                geolat_u:_FillValue = 1.e+20f ;
                geolat_u:cell_methods = "time: point" ;
                geolat_u:interp_method = "none" ;
        float geolon_v(yq, xh) ;
                geolon_v:long_name = "Longitude of meridional velocity (Cv) points" ;
                geolon_v:units = "degrees_east" ;
                geolon_v:missing_value = 1.e+20f ;
                geolon_v:_FillValue = 1.e+20f ;
                geolon_v:cell_methods = "time: point" ;
                geolon_v:interp_method = "none" ;
        float geolat_v(yq, xh) ;
                geolat_v:long_name = "Latitude of meridional velocity (Cv) points" ;
                geolat_v:units = "degrees_north" ;
                geolat_v:missing_value = 1.e+20f ;
                geolat_v:_FillValue = 1.e+20f ;
                geolat_v:cell_methods = "time: point" ;
                geolat_v:interp_method = "none" ;
...  wet wet wet
        float wet(yh, xh) ;
                wet:long_name = "0 if land, 1 if ocean at tracer points" ;
                wet:missing_value = 1.e+20f ;
                wet:_FillValue = 1.e+20f ;
                wet:cell_methods = "time: point" ;
                wet:cell_measures = "area: areacello" ;
        float wet_u(yh, xq) ;
                wet_u:long_name = "0 if land, 1 if ocean at zonal velocity (Cu) points" ;
                wet_u:missing_value = 1.e+20f ;
                wet_u:_FillValue = 1.e+20f ;
                wet_u:cell_methods = "time: point" ;
                wet_u:interp_method = "none" ;
        float wet_v(yq, xh) ;
                wet_v:long_name = "0 if land, 1 if ocean at meridional velocity (Cv) points" ;
                wet_v:missing_value = 1.e+20f ;
                wet_v:_FillValue = 1.e+20f ;
                wet_v:cell_methods = "time: point" ;
                wet_v:interp_method = "none" ;

hkershaw-brown added a commit that referenced this issue Mar 8, 2023
@hkershaw-brown
using full 2D geolon, geolat for T,U,V quad interpolation
d93a56f 
see #451 for discusion on which files the grids are stored in. Be
aware MOM6 is in development.

Note wet (ocean or land) is available for T,U,V. Still just using T
in this commit.

Also fixes on_t_grid if statement logic
@hkershaw-brown
Copy link
Member Author

Note on layer thickness which varies across the ensemble.

For this example run:
/glade/scratch/hkershaw/c.e22.GMOM.T62_g16.nuopc.001/run

restart files have 60 layers (61 interfaces)

netcdf c.e22.GMOM.T62_g16.nuopc.001.mom6.r.0001-01-06-00000 {
dimensions:
        lath = 384 ;
        lonh = 320 ;
        latq = 384 ;
        lonq = 320 ;
        Layer = 60 ;
        Interface = 61 ;
        Time = UNLIMITED ; // (1 currently)

and h as layer thickness

        double h(Time, Layer, lath, lonh) ;
                h:long_name = "Layer Thickness" ;
                h:units = "m" ;

history files have fewer layers?

netcdf c.e22.GMOM.T62_g16.nuopc.001.mom6.h_0001_09 {
dimensions:
        xq = 320 ;
        yh = 384 ;
        z_l = 34 ;
        z_i = 35 ;

where h is layer thickness

        float h(time, z_l, yh, xh) ;
                h:long_name = "Layer Thickness" ;
                h:units = "m" ;

Note, layer thickness may be called e in some MOM6 versions, I think it is h here.

I'm going with restart file layers since we are using the restart file data as the state.
Should the thickness be part of the state and get updated by assimilation?

@amrhein
Copy link
Collaborator

amrhein commented Mar 9, 2023 via email

@hkershaw-brown hkershaw-brown mentioned this issue Mar 29, 2023
Merged
15 tasks
@hkershaw-brown
Copy link
Member Author

Notes on the scripting for MOM6

Looking at the cam-fv various scripts I'm not sure what the goal is of having a no_assimilate.csh and an assimilate.csh
vs. DATA_ASSIMILATION=FALSE/TRUE

e.g setup_advanced

4)  The default initial configuration is to assimilate.
    Verify the ${caseroot}/input.nml contents.
    Assimilation can be turned off by
    ./xmlchange DATA_ASSIMILATION_SCRIPT=${caseroot}/no_assimilate.csh
    DART can be turned off by
    ./xmlchange DATA_ASSIMILATION=FALSE

DART.config.template

2) If you want to turn DART on or off, edit the 
   env_run.xml: DATA_ASSIMILATION_* to specify which DART script 
   to use after the model forecast;
   assimilate.csh, no_assimilate.csh, or perfect_model.csh.

@hkershaw-brown
Copy link
Member Author

Looking at where the DATA_ASSIMILATION_SCRIPT is called in case_run.py:
If can we pass the case (would need a cime change) to _do_assimilate
https://github.com/ESMCI/cime/blob/42bc2916478a060f6daca7871299966dced52f6c/CIME/case/case_run.py#L411

and use the case module to access al the xml variables (DATA_ASSIMILATION_SCRIPT in python)?
https://github.com/ESMCI/cime/blob/42bc2916478a060f6daca7871299966dced52f6c/CIME/utils.py#L667-L674

@hkershaw-brown
Copy link
Member Author

There is a data assimilation true/false for each component, but only one data_assimilation_script

((cesm2_3_beta11))$ ./xmlquery DATA_ASSIMILATION_SCRIPT
	DATA_ASSIMILATION_SCRIPT: assimilate.sh
((cesm2_3_beta11))$ ./xmlquery DATA_ASSIMILATION
	DATA_ASSIMILATION: ['CPL:TRUE', 'ATM:TRUE', 'LND:TRUE', 'ICE:TRUE', 'OCN:TRUE', 'ROF:TRUE', 'GLC:TRUE', 'WAV:TRUE']
((cesm2_3_beta11))$ ./xmlquery DATA_ASSIMILATION_LND
	DATA_ASSIMILATION_LND: TRUE

data_assimilation logical a function in case_run.py
https://github.com/ESMCI/cime/blob/42bc2916478a060f6daca7871299966dced52f6c/CIME/case/case_run.py#L463-L470

hkershaw-brown added a commit that referenced this issue May 12, 2023
@hkershaw-brown
toy assimilate script for MOM6
a1115c5 
See #451 for discussion
@hkershaw-brown
Copy link
Member Author

🐛 MOM6 model_mod read_time is currently just reading the time from the MOM6 restart file. This is years from 1 (or maybe 0?). The obs_seq files have the dart time (1601 + X).
For read_model_time, should we convert the model_time to dart time?

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Jun 6, 2023 •
edited
Loading 

 PE 0:  filter trace: Before computing prior observation values
mlx5: r5i3n3: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
000000e9 00000000 00000000 00000000
00000000 00008813 100036c4 03421ad2
mlx5: r5i3n3: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
000000f9 00000000 00000000 00000000
00000000 00008813 1000390d 00038ed2
MPT ERROR: 07:16:51: Rank 265: Send completion to rank 233, status 10: remote memory region protection error
	265->233 uses r5i3n3:mlx5_0:1:0x80fe -> r5i0n16
MPT ERROR: 07:16:51: Rank 267: Send completion to rank 249, status 10: remote memory region protection error
	267->249 uses r5i3n3:mlx5_0:1:0x80fe -> r5i0n16
MPT ERROR: Rank 267(g:267) is aborting with error code 0.
	Process ID: 25429, Host: r5i3n3, Program: /glade/scratch/hkershaw/c.T62_g16.Alper.ens3/bld/filter
	MPT Version: HPE MPT 2.25  08/14/21 03:05:20

Not sure if this is a bug on our end or Cheyenne network problems. June 6 2023
Edit: or both.

Edit: this was a bug, not catching observations that are too deep: fixed in #492

@TomNicholas
Copy link

Was this closed by #467?

@hkershaw-brown
Copy link
Member Author

Was this closed by #467?

Hi @TomNicholas the scripting for running a MOM6-DART with CESM is not done yet.
We have a CESM 2.3 (and 3.0) target with DART integrated into cime.

Let us know if you are interested in running MOM6-DART or if you have scripts/code that you want to contribute.

@hkershaw-brown hkershaw-brown added the mom6 Modular Ocean Model label Jun 20, 2023
@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Nov 10, 2023 •
edited
Loading

Derecho: CESM tag to use cesm2_3_alpha16g
mom6-scripting https://github.com/NCAR/DART/tree/mom6-scripting

@hkershaw-brown hkershaw-brown mentioned this issue Dec 13, 2023
Closed
15 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@hkershaw-brown hkershaw-brown

Labels
mom6 Modular Ocean Model
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

MOM6 assimilate script for CESM NCAR/DART
5 participants
@amrhein @nancycollins @hkershaw-brown @TomNicholas @jlaucar