Skip to content

feat: Add clamping variables to MOM6, set vertical localization coord, vert interp #1033

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

Merged
hkershaw-brown merged 4 commits into from
Mar 6, 2026
+73 −34
Merged

feat: Add clamping variables to MOM6, set vertical localization coord, vert interp #1033

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
eaee042
feat: add clamping variables to MOM6
hkershaw-brown Jan 6, 2026
59dc3d1
set vertical localization coordinate
hkershaw-brown Feb 19, 2026
7a31f84
fix: vertical interp fixes
hkershaw-brown Feb 20, 2026
6cb5d7f
bump version and changelog for release
hkershaw-brown Mar 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 37 additions & 22 deletions models/MOM6/model_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,8 @@ module model_mod
loc_get_close_state => get_close_state, &
set_location, set_location_missing, &
get_location, query_location, VERTISLEVEL, &
VERTISHEIGHT, set_vertical, get_dist
VERTISHEIGHT, set_vertical, get_dist, &
set_vertical_localization_coord

use utilities_mod, only : error_handler, E_ERR, E_MSG, &
nmlfileunit, do_output, do_nml_file, do_nml_term, &
Expand Down Expand Up @@ -60,8 +61,8 @@ module model_mod
init_time => fail_init_time, &
init_conditions => fail_init_conditions, &
convert_vertical_obs, adv_1step, &
parse_variables, &
MAX_STATE_VARIABLE_FIELDS
parse_variables_clamp, &
MAX_STATE_VARIABLE_FIELDS_CLAMP

implicit none
private
Expand Down Expand Up @@ -127,7 +128,7 @@ module model_mod
character(len=256) :: ocean_geometry = 'ocean_geometry.nc'
integer :: assimilation_period_days = 1
integer :: assimilation_period_seconds = 0
character(len=vtablenamelength) :: model_state_variables(MAX_STATE_VARIABLE_FIELDS) = ' '
character(len=vtablenamelength) :: model_state_variables(MAX_STATE_VARIABLE_FIELDS_CLAMP) = ' '
character(len=NF90_MAX_NAME) :: layer_name = 'Layer'
logical :: use_pseudo_depth = .false. ! use pseudo depth instead of sum(layer thickness) for vertical location

Expand Down Expand Up @@ -165,6 +166,7 @@ subroutine static_init_model()
if (do_nml_term()) write( * , nml=model_nml)

call set_calendar_type('gregorian')
call set_vertical_localization_coord(VERTISHEIGHT)

! This time is both the minimum time you can ask the model to advance
! (for models that can be advanced by filter) and it sets the assimilation
Expand All @@ -178,7 +180,7 @@ subroutine static_init_model()
! parse_variables converts the character table that was read in from
! model_nml:model_state_variables to a state_var_type that can be passed
! to add_domain
dom_id = add_domain(template_file, parse_variables(model_state_variables))
dom_id = add_domain(template_file, parse_variables_clamp(model_state_variables))

model_size = get_domain_size(dom_id)

Expand Down Expand Up @@ -272,6 +274,7 @@ subroutine model_interpolate(state_handle, ens_size, location, qty_in, expected_
! unpack the location type into lon, lat, vert, vert_type
lon_lat_vert = get_location(location)
which_vert = nint(query_location(location))
if (which_vert /= VERTISHEIGHT) call error_handler(E_ERR, 'model_interpolate', 'only supports VERTISHEIGHT')

! get the indices for the 4 corners of the quad in the horizontal
call quad_lon_lat_locate(interp, lon_lat_vert(1), lon_lat_vert(2), &
Expand All @@ -290,7 +293,7 @@ subroutine model_interpolate(state_handle, ens_size, location, qty_in, expected_
if (use_pseudo_depth) then

! Get the bounding vertical levels and the fraction between bottom and top
call find_level_bounds(lon_lat_vert(3), which_vert, levz, levz_fract, locate_status)
call find_level_bounds(lon_lat_vert(3), levz, levz_fract, locate_status)
if (locate_status /= 0) then
istatus(:) = locate_status
return
Expand All @@ -306,7 +309,7 @@ subroutine model_interpolate(state_handle, ens_size, location, qty_in, expected_
! HK @todo Do you need to use t_grid interp for thickness four_ilons, four_ilats?
found(:) = .false.
depth_at_x(:) = 0
FIND_LAYER: do i = 2, nz
FIND_LAYER: do i = 1, nz

do corner = 1, 4
th_indx = get_dart_vector_index(four_ilons(corner), four_ilats(corner), i, dom_id, thick_id)
Expand All @@ -330,10 +333,17 @@ subroutine model_interpolate(state_handle, ens_size, location, qty_in, expected_

do e = 1, ens_size
if (lon_lat_vert(3) < depth_at_x(e)) then
lev(e,1) = i ! layer_below
lev(e,2) = i-1 ! layer_above
lev_fract(e) = (depth_at_x(e) - lon_lat_vert(3)) / thick_at_x(e)
found(e) = .true.
if (i ==1) then
lev(e,1) = 1 ! in surface layer
lev(e,2) = 1 ! in surface layer
lev_fract(e) = 0.0_r8
found(e) = .true.
else
lev(e,1) = i ! layer_below
lev(e,2) = i-1 ! layer_above
lev_fract(e) = (depth_at_x(e) - lon_lat_vert(3)) / thick_at_x(e)
found(e) = .true.
endif
if (all(found)) exit FIND_LAYER
endif
enddo
Expand Down Expand Up @@ -519,6 +529,7 @@ subroutine convert_vertical_state(state_handle, num, locs, loc_qtys, loc_indx, &
real(r8) :: depth(1)

! assert(which_vert == VERTISHEIGHT)
if (which_vert /= VERTISHEIGHT) call error_handler(E_ERR,'model_mod convert_vertical_state', 'only supports VERTISHEIGHT')

if (use_pseudo_depth) then
do ii = 1, num
Expand Down Expand Up @@ -1081,26 +1092,30 @@ end function read_model_time
! --- 1 i=3
!

subroutine find_level_bounds(vert_loc, which_vert, lev, lev_fract, istatus)
subroutine find_level_bounds(vert_loc, lev, lev_fract, istatus)

real(r8), intent(in) :: vert_loc ! observation location
integer, intent(in) :: which_vert ! obs vertical coordinate
integer, intent(out) :: lev(2) ! bottom, top
real(r8), intent(out) :: lev_fract
integer, intent(out) :: istatus

integer :: i

! HK assert(which_vert == height meters) ?

do i = 2, nz
do i = 1, nz
if(vert_loc < zstar(i)) then
lev(1) = i ! bottom
lev(2) = i-1 ! top

lev_fract = (zstar(lev(1)) - vert_loc) / (zstar(lev(2)) - zstar(lev(1)))
istatus = 0
return
if (i == 1) then
lev(1) = 1 ! obs is in surface layer
lev(2) = 1 ! obs is in surface layer
lev_fract = 0.0_r8
istatus = 0
return
else
lev(1) = i ! bottom
lev(2) = i-1 ! top
lev_fract = (zstar(lev(1)) - vert_loc) / (zstar(lev(1)) - zstar(lev(2)))
istatus = 0
return
endif
endif
enddo

Expand Down
21 changes: 15 additions & 6 deletions models/MOM6/readme.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,12 +74,12 @@ The namelist options for DART-MOM6 are as follows:
&model_nml
template_file = 'mom6.r.nc',
ocean_geometry = 'ocean_geometry.nc',
static_file = 'c.e22.GMOM.T62_g16.nuopc.001.mom6.static.nc',
model_state_variables = 'Salt ', 'QTY_SALINITY ', 'UPDATE',
'Temp ', 'QTY_POTENTIAL_TEMPERATURE', 'UPDATE',
'u ', 'QTY_U_CURRENT_COMPONENT ', 'UPDATE',
'v ', 'QTY_V_CURRENT_COMPONENT ', 'UPDATE',
'h ', 'QTY_LAYER_THICKNESS ', 'UPDATE',
static_file = 'static.nc',
model_state_variables = 'Salt ', 'QTY_SALINITY ', 'NA', 'NA', 'UPDATE',
'Temp ', 'QTY_POTENTIAL_TEMPERATURE', 'NA', 'NA', 'UPDATE',
'u ', 'QTY_U_CURRENT_COMPONENT ', 'NA', 'NA', 'UPDATE',
'v ', 'QTY_V_CURRENT_COMPONENT ', 'NA', 'NA', 'UPDATE',
'h ', 'QTY_LAYER_THICKNESS ', 'NA', 'NA', 'UPDATE',
assimilation_period_days = 1
assimilation_period_seconds = 0
use_pseudo_depth = .false. ! use pseudo depth instead of sum(layer thickness) for vertical location
Expand All @@ -104,6 +104,15 @@ The namelist options for DART-MOM6 are as follows:
geolat_v(:,:) Latitude of meridional velocity (Cv) point


* ``model_state_variables`` defines the list of model variables from the MOM6 restart file that will be included in the DART state. Each row in the table should have the following fields:

- **NetCDF variable name**: Name of the variable in the MOM6 restart file (e.g., 'Salt').
- **DART Quantity**: The DART quantity associated with the variable (e.g., 'QTY_SALINITY').
- **Clamping lower bound**: Minimum allowed value for the variable when writing out restarts (use 'NA' for no bound).
- **Clamping upper bound**: Maximum allowed value for the variable when writing out restarts (use 'NA' for no bound).
- **UPDATE or NO_COPY_BACK**: Use 'UPDATE' to allow DART to update this variable during assimilation, or 'NO_COPY_BACK' to prevent updates (variable will be read but not written back).


Vertical Coordinate
-------------------

Expand Down
10 changes: 5 additions & 5 deletions models/MOM6/work/input.nml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,11 +142,11 @@
template_file = 'mom6.r.nc',
static_file = 'mom6.static.nc',
ocean_geometry = 'ocean_geometry.nc',
model_state_variables = 'Salt ', 'QTY_SALINITY ', 'UPDATE',
'Temp ', 'QTY_POTENTIAL_TEMPERATURE', 'UPDATE',
'u ', 'QTY_U_CURRENT_COMPONENT ', 'UPDATE',
'v ', 'QTY_V_CURRENT_COMPONENT ', 'UPDATE',
'h ', 'QTY_LAYER_THICKNESS ', 'UPDATE',
model_state_variables = 'Salt ', 'QTY_SALINITY ', 'NA', 'NA', 'UPDATE',
'Temp ', 'QTY_POTENTIAL_TEMPERATURE', 'NA', 'NA', 'UPDATE',
'u ', 'QTY_U_CURRENT_COMPONENT ', 'NA', 'NA', 'UPDATE',
'v ', 'QTY_V_CURRENT_COMPONENT ', 'NA', 'NA', 'UPDATE',
'h ', 'QTY_LAYER_THICKNESS ', 'NA', 'NA', 'UPDATE',
/

&utilities_nml
Expand Down