Merge pull request #727 from NCAR/mom6-missing-grid

MOM6 static file mask missing values

Author: hkershaw-brown

CHANGELOG.rst

@@ -22,7 +22,12 @@ individual files.
 
 The changes are now listed with the most recent at the top.
 
-**September 10 2024:: MARBL_column. Tag 11.8.0**
+**September 27 2024 :: MOM6 mask bug-fix. Tag 11.8.1**
+
+- Fix for MOM6 CESM3 workhorse 2/3 degree grid TL319_t232 to 
+  mask missing geolon|lat|u|v|t values  
+
+**September 10 2024 :: MARBL_column. Tag 11.8.0**
 
 - Interface for MARBL_column for DART: 
 

conf.py

@@ -21,7 +21,7 @@
 author = 'Data Assimilation Research Section'
 
 # The full version, including alpha/beta/rc tags
-release = '11.8.0'
+release = '11.8.1'
 root_doc = 'index'
 
 # -- General configuration ---------------------------------------------------

models/MOM6/model_mod.f90

@@ -32,7 +32,7 @@ module model_mod
                                  nc_begin_define_mode, nc_end_define_mode, &
                                  nc_open_file_readonly, nc_close_file, &
                                  nc_get_variable, nc_get_variable_size, &
-                                 NF90_MAX_NAME
+                                 NF90_MAX_NAME, nc_get_attribute_from_variable
 
 use        quad_utils_mod,  only : quad_interp_handle, init_quad_interp, &
                                    set_quad_coords, quad_lon_lat_locate, &
@@ -97,6 +97,7 @@ module model_mod
 real(r8), allocatable :: geolon(:,:), geolat(:,:),     & ! T
                          geolon_u(:,:), geolat_u(:,:), & ! U
                          geolon_v(:,:), geolat_v(:,:)    ! V
+logical, allocatable  :: mask(:,:), mask_u(:,:), mask_v(:,:) ! geolat/lon/u/v has missing values
 type(quad_interp_handle) :: interp_t_grid,  &
                             interp_u_grid,  &
                             interp_v_grid
@@ -595,6 +596,7 @@ subroutine read_horizontal_grid()
 
 integer :: ncid
 integer :: nxy(2) ! (nx,ny)
+real(r8) :: fillval
 
 character(len=*), parameter :: routine = 'read_horizontal_grid'
 
@@ -606,22 +608,48 @@ subroutine read_horizontal_grid()
 allocate(geolon(nx,ny), geolat(nx,ny))      ! T grid
 allocate(geolon_u(nx,ny), geolat_u(nx,ny))  ! U grid
 allocate(geolon_v(nx,ny), geolat_v(nx,ny))  ! V grid
+allocate(mask(nx,ny))  ! missing values
+allocate(mask_u(nx,ny))  ! missing values
+allocate(mask_v(nx,ny))  ! missing values
 
 call nc_get_variable(ncid, 'geolon', geolon, routine)
 call nc_get_variable(ncid, 'geolon_u', geolon_u, routine)
 call nc_get_variable(ncid, 'geolon_v', geolon_v, routine)
 
-
-! mom6 example file has longitude > 360
-! DART uses [0,360]
-where(geolon > 360.0_r8 )   geolon   = geolon   - 360.0_r8
-where(geolon_u > 360.0_r8 ) geolon_u = geolon_u - 360.0_r8
-where(geolon_v > 360.0_r8 ) geolon_v = geolon_v - 360.0_r8
-
 call nc_get_variable(ncid, 'geolat', geolat, routine)
 call nc_get_variable(ncid, 'geolat_u', geolat_u, routine)
 call nc_get_variable(ncid, 'geolat_v', geolat_v, routine)
 
+! mom6 has missing values in the grid
+! and set missing value to a land point to prevent set_location erroring
+mask(:,:) = .false.
+mask_u(:,:) = .false.
+mask_v(:,:) = .false.
+call nc_get_attribute_from_variable(ncid, 'geolon', '_FillValue', fillval)
+where (geolon == fillval) mask = .true.  
+where (geolon == fillval) geolon = 72.51
+where (geolat == fillval) geolat = 42.56
+
+call nc_get_attribute_from_variable(ncid, 'geolon_u', '_FillValue', fillval)
+where (geolon_u == fillval) mask_u = .true.  
+where (geolon_u == fillval) geolon_u = 72.51
+where (geolat_u == fillval) geolat_u = 42.56
+
+call nc_get_attribute_from_variable(ncid, 'geolon_v', '_FillValue', fillval)
+where (geolon_v == fillval) mask_v = .true.  
+where (geolon_v == fillval) geolon_v = 72.51
+where (geolat_v == fillval) geolat_v = 42.56
+
+! mom6 example files have longitude > 360 and longitudes < 0
+! DART uses [0,360]
+geolon = mod(geolon, 360.0)
+geolon_u = mod(geolon_u, 360.0)
+geolon_v = mod(geolon_v, 360.0)
+
+where (geolon < 0.0) geolon = geolon + 360
+where (geolon_u < 0.0) geolon_u = geolon_u + 360
+where (geolon_v < 0.0) geolon_v = geolon_v + 360
+
 call nc_close_file(ncid)
 
 end subroutine read_horizontal_grid
@@ -831,15 +859,15 @@ subroutine setup_interpolation()
                       global=.true., spans_lon_zero=.true., pole_wrap=.true., &
                       interp_handle=interp_t_grid)
 
-call set_quad_coords(interp_t_grid, geolon, geolat)
+call set_quad_coords(interp_t_grid, geolon, geolat, mask)
 
 ! U
 call init_quad_interp(GRID_QUAD_FULLY_IRREGULAR, nx, ny, &
                       QUAD_LOCATED_CELL_CENTERS, &
                       global=.true., spans_lon_zero=.true., pole_wrap=.true., &
                       interp_handle=interp_u_grid)
 
-call set_quad_coords(interp_u_grid, geolon_u, geolat_u)
+call set_quad_coords(interp_u_grid, geolon_u, geolat_u, mask_u)
 
 
 ! V
@@ -848,7 +876,7 @@ subroutine setup_interpolation()
                       global=.true., spans_lon_zero=.true., pole_wrap=.true., &
                       interp_handle=interp_v_grid)
 
-call set_quad_coords(interp_v_grid, geolon_v, geolat_v)
+call set_quad_coords(interp_v_grid, geolon_v, geolat_v, mask_v)
 
 end subroutine setup_interpolation
 

models/MOM6/work/input.nml

@@ -20,8 +20,8 @@
 
    obs_seq_in_file_name       = "obs_seq.in",
    obs_seq_out_file_name      = "obs_seq.out",
-   init_time_days             = 0,
-   init_time_seconds          = 0,
+   init_time_days             = -1,
+   init_time_seconds          = -1,
    first_obs_days             = -1,
    first_obs_seconds          = -1,
    last_obs_days              = -1,
@@ -93,7 +93,7 @@
    /
 
 &assim_tools_nml
-   cutoff                          = 1000000.0
+   cutoff                          = 0.03
    sort_obs_inc                    = .false.,
    spread_restoration              = .false.,
    sampling_error_correction       = .false.,