Feature request: MOM6 interface

[hkershaw-brown]

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:

• [x] MOM6 model_mod
• [x] run MOM6 in cesm
• [x] test nuopc cesm with DART
• [x] 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]

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]

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?

On Wed, Feb 15, 2023 at 2:27 PM hkershaw-brown @.***> wrote:

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" ;

— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/451#issuecomment-1432056976, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANDHUISFAHFPIQXLGJFCEYTWXVC53ANCNFSM6AAAAAAUSZZY5A . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[amrhein]

The kids are home today so I'm mostly away, but Helen let's put this on the agenda for our meeting with Alper tomorrow afternoon.

Dan

On Wed, Feb 15, 2023 at 2:54 PM jlaucar @.***> wrote:

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?

On Wed, Feb 15, 2023 at 2:27 PM hkershaw-brown @.***> wrote:

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" ;

— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/451#issuecomment-1432056976, or unsubscribe < https://github.com/notifications/unsubscribe-auth/ANDHUISFAHFPIQXLGJFCEYTWXVC53ANCNFSM6AAAAAAUSZZY5A

. You are receiving this because you are subscribed to this thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/451#issuecomment-1432090088, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADCQXR2V46XSTYUR6ZSL2YLWXVGCHANCNFSM6AAAAAAUSZZY5A . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[hkershaw-brown]

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]

My ignorant assumption is that yes, they could be in different 'layers' for different ensemble members. I think I recall coding this capability when writing the first draft of the forward operators for CAM when we were moving to Manhattan? Can we just implement the same capability for MOM?

On Fri, Feb 17, 2023 at 11:09 AM hkershaw-brown @.***> wrote:

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.

— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/451#issuecomment-1435053867, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANDHUIUJ6YE23TT6TFXBJ2TWX65GVANCNFSM6AAAAAAUSZZY5A . You are receiving this because you commented.Message ID: @.***>

[hkershaw-brown]

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

[jlaucar]

Maybe Helen wrote that code and I'm just remembering discussions about the issue with different layers.

On Fri, Feb 17, 2023 at 11:30 AM Jeffrey Anderson @.***> wrote:

My ignorant assumption is that yes, they could be in different 'layers' for different ensemble members. I think I recall coding this capability when writing the first draft of the forward operators for CAM when we were moving to Manhattan? Can we just implement the same capability for MOM?

On Fri, Feb 17, 2023 at 11:09 AM hkershaw-brown @.***> wrote:

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.

— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/451#issuecomment-1435053867, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANDHUIUJ6YE23TT6TFXBJ2TWX65GVANCNFSM6AAAAAAUSZZY5A . You are receiving this because you commented.Message ID: @.***>

[hkershaw-brown]

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. https://github.com/NCAR/DART/blob/224698d11846c83ed8141156b94338232dff7271/assimilation_code/modules/utilities/distributed_state_mod.f90#L35-L39

[nancycollins]

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]

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]

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]

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

[hkershaw-brown]

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]

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)

[hkershaw-brown]

boom :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]

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]

Should the thickness be part of the state and get updated by assimilation?

No, I don't think so. I talked with Frank about this. The ALE coordinate is not tied to any conservation laws as would be an isopycnal model, where e.g. if you changed temperature might have to alter layer thickness for consistency. MOM6 will do its own remapping based on the restart file it gets based on that previous set of coordinates.

If I'm wrong, hopefully something horrific will happen and we'll know quickly ;-)

On Thu, Mar 9, 2023 at 8:57 AM hkershaw-brown @.***> wrote:

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?

— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/451#issuecomment-1462309788, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADCQXRZNWAJCRBDCBWC3NODW3H4V7ANCNFSM6AAAAAAUSZZY5A . You are receiving this because you commented.Message ID: @.***>

[hkershaw-brown]

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]

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]

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]

:bug: 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]

 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 https://github.com/NCAR/DART/pull/492

[TomNicholas]

Was this closed by #467?

[hkershaw-brown]

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]

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