Merge pull request #782 from NCAR/mom6-pot-temp

fix: potential temp to temp conversion for MOM6 model_interpolate

Author: hkershaw-brown

CHANGELOG.rst

@@ -22,7 +22,12 @@ individual files.
 
 The changes are now listed with the most recent at the top.
 
-**January 9 2024 :: Bug-fix 1D obs_diag. Tag v11.8.7**
+**January 14 2025 :: Bug-fix MOM6 potential temperature. Tag v11.8.8**
+
+- MOM6 model_interpolate for potential temperature
+- Update lorenz workshop input.nmls to v11
+
+**January 9 2025 :: Bug-fix 1D obs_diag. Tag v11.8.7**
 
 - Added a dummy dimension so 1D obs_diag output can be used with 
   MATLAB diagnostic tools

conf.py

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

models/MOM6/model_mod.f90

@@ -49,7 +49,8 @@ module model_mod
 
 use obs_kind_mod, only : get_index_for_quantity, QTY_U_CURRENT_COMPONENT, &
                          QTY_V_CURRENT_COMPONENT, QTY_LAYER_THICKNESS, &
-                         QTY_DRY_LAND, QTY_SALINITY
+                         QTY_DRY_LAND, QTY_SALINITY, QTY_TEMPERATURE, &
+                         QTY_POTENTIAL_TEMPERATURE
 
 use ensemble_manager_mod, only : ensemble_type
 
@@ -113,8 +114,8 @@ module model_mod
 integer, parameter :: NOT_IN_STATE = 12
 integer, parameter :: THICKNESS_NOT_IN_STATE = 13
 integer, parameter :: QUAD_LOCATE_FAILED = 14
-integer, parameter :: THICKNESS_QUAD_EVALUTATE_FAILED = 15
-integer, parameter :: QUAD_EVALUTATE_FAILED = 16
+integer, parameter :: THICKNESS_QUAD_EVALUATE_FAILED = 15
+integer, parameter :: QUAD_EVALUATE_FAILED = 16
 integer, parameter :: QUAD_ON_LAND = 17
 integer, parameter :: QUAD_ON_BASIN_EDGE = 18
 integer, parameter :: OBS_ABOVE_SURFACE = 20
@@ -220,25 +221,29 @@ end function get_model_size
 ! 0 unless there is some problem in computing the interpolation in
 ! which case a positive istatus should be returned.
 
-subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs, istatus)
+subroutine model_interpolate(state_handle, ens_size, location, qty_in, expected_obs, istatus)
 
 
 type(ensemble_type), intent(in) :: state_handle
 integer,             intent(in) :: ens_size
 type(location_type), intent(in) :: location
-integer,             intent(in) :: qty
+integer,             intent(in) :: qty_in
 real(r8),           intent(out) :: expected_obs(ens_size) !< array of interpolated values
 integer,            intent(out) :: istatus(ens_size)
 
+real(r8), parameter             :: CONCENTRATION_TO_PPT = 1000.0_r8
+
+integer  :: qty ! local qty
 integer  :: which_vert, four_ilons(4), four_ilats(4), lev(ens_size,2)
 integer  :: locate_status, quad_status
 real(r8) :: lev_fract(ens_size)
 real(r8) :: lon_lat_vert(3)
 real(r8) :: quad_vals(4, ens_size)
 real(r8) :: expected(ens_size, 2) ! level below and above obs
+real(r8) :: expected_pot_temp(ens_size), expected_salinity(ens_size), pressure_dbars(ens_size)
 type(quad_interp_handle) :: interp
 integer :: varid, i, e, thick_id
-integer(i8) :: th_indx, indx(ens_size)
+integer(i8) :: th_indx
 real(r8) :: depth_at_x(ens_size), thick_at_x(ens_size) ! depth, layer thickness at obs lat lon
 logical :: found(ens_size)
 
@@ -247,6 +252,16 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs
 expected_obs(:) = MISSING_R8
 istatus(:) = 1
 
+if (qty_in == QTY_TEMPERATURE) then
+   qty = QTY_POTENTIAL_TEMPERATURE  ! model has potential temperature
+   if (get_varid_from_kind(dom_id, QTY_SALINITY) < 0) then ! Require salinity to convert to temperature
+      istatus = NOT_IN_STATE
+      return
+   end if
+else
+   qty = qty_in
+endif
+
 varid = get_varid_from_kind(dom_id, qty)
 if (varid < 0) then ! not in state
    istatus = NOT_IN_STATE
@@ -310,7 +325,7 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs
                               thick_at_x, &
                               quad_status)
    if (quad_status /= 0) then
-      istatus(:) = THICKNESS_QUAD_EVALUTATE_FAILED
+      istatus(:) = THICKNESS_QUAD_EVALUATE_FAILED
       return
    endif
 
@@ -338,6 +353,66 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs
    return
 endif
 
+
+select case (qty_in)
+   case (QTY_TEMPERATURE)
+      ! convert from potential temperature to temperature
+
+      call state_on_quad(four_ilons, four_ilats, lon_lat_vert, ens_size, lev, lev_fract, interp, state_handle, varid, expected_pot_temp, quad_status)
+      if (quad_status /= 0) then
+         istatus(:) = QUAD_EVALUATE_FAILED
+         return
+      endif
+      call state_on_quad(four_ilons, four_ilats, lon_lat_vert, ens_size, lev, lev_fract, interp, state_handle, get_varid_from_kind(dom_id, QTY_SALINITY), expected_salinity, quad_status)
+      if (quad_status /= 0) then
+         istatus(:) = QUAD_EVALUATE_FAILED
+         return
+      endif
+
+      pressure_dbars =  0.059808_r8*(exp(-0.025_r8*depth_at_x) - 1.0_r8)  &
+                        + 0.100766_r8*depth_at_x + 2.28405e-7_r8*lon_lat_vert(3)**2
+      expected_obs = sensible_temp(expected_pot_temp, expected_salinity, pressure_dbars)
+
+   case (QTY_SALINITY) ! convert from g of salt per kg of seawater (model) to kg of salt per kg of seawater (observation)
+      call state_on_quad(four_ilons, four_ilats, lon_lat_vert, ens_size, lev, lev_fract, interp, state_handle, varid, expected_obs, quad_status)
+      if (quad_status /= 0) then 
+         istatus(:) = QUAD_EVALUATE_FAILED
+         return
+      endif
+      expected_obs = expected_obs/CONCENTRATION_TO_PPT
+ 
+   case default
+      call state_on_quad(four_ilons, four_ilats, lon_lat_vert, ens_size, lev, lev_fract, interp, state_handle, varid, expected_obs, quad_status)
+      if (quad_status /= 0) then
+         istatus(:) = QUAD_EVALUATE_FAILED
+         return
+      endif
+end select
+
+istatus(:) = 0
+
+end subroutine model_interpolate
+
+!------------------------------------------------------------------
+! Interpolate on the quad, between two levels
+subroutine state_on_quad(four_ilons, four_ilats, lon_lat_vert, ens_size, lev, lev_fract, interp, state_handle, varid, expected_obs, quad_status)
+
+integer,  intent(in) :: four_ilons(4), four_ilats(4) ! indices into lon, lat
+real(r8), intent(in) :: lon_lat_vert(3) ! lon, lat, vert of obs
+integer,  intent(in) :: ens_size
+integer,  intent(in) :: lev(ens_size,2) ! levels below and above obs
+real(r8), intent(in) :: lev_fract(ens_size) 
+type(quad_interp_handle), intent(in) :: interp
+type(ensemble_type), intent(in) :: state_handle
+integer,  intent(in) :: varid ! which state variable
+real(r8), intent(out) :: expected_obs(ens_size)
+integer,  intent(out) :: quad_status
+
+integer :: i, e
+integer(i8) :: indx(ens_size)
+real(r8) :: quad_vals(4, ens_size)
+real(r8) :: expected(ens_size, 2) ! state value at level below and above obs
+
 do i = 1, 2 
    !HK which corner of the quad is which?
    ! corner1
@@ -371,27 +446,15 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs
                               quad_vals, & ! 4 corners x ens_size
                               expected(:,i), &
                               quad_status)
-   if (quad_status /= 0) then
-      istatus(:) = QUAD_EVALUTATE_FAILED
-      return
-   else
-      istatus = 0
-   endif
+   if (quad_status /= 0) return
 
 enddo
 
 ! Interpolate between levels
 ! expected_obs = bot_val + lev_fract * (top_val - bot_val)
 expected_obs = expected(:,1) + lev_fract(:) * (expected(:,2) - expected(:,1))
 
-if (qty == QTY_SALINITY) then  ! convert from PSU (model) to MSU (obersvation)
-   expected_obs = expected_obs/1000.0_r8
-endif
-
-
-end subroutine model_interpolate
-
-
+end subroutine state_on_quad
 
 !------------------------------------------------------------------
 ! Returns the smallest increment in time that the model is capable 
@@ -627,28 +690,28 @@ subroutine read_horizontal_grid()
 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
+where (geolon == fillval) geolon = 72.51_r8
+where (geolat == fillval) geolat = 42.56_r8
 
 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
+where (geolon_u == fillval) geolon_u = 72.51_r8
+where (geolat_u == fillval) geolat_u = 42.56_r8
 
 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
+where (geolon_v == fillval) geolon_v = 72.51_r8
+where (geolat_v == fillval) geolat_v = 42.56_r8
 
 ! 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)
+geolon = mod(geolon, 360.0_r8)
+geolon_u = mod(geolon_u, 360.0_r8)
+geolon_v = mod(geolon_v, 360.0_r8)
 
-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
+where (geolon < 0.0) geolon = geolon + 360.0_r8
+where (geolon_u < 0.0) geolon_u = geolon_u + 360.0_r8
+where (geolon_v < 0.0) geolon_v = geolon_v + 360.0_r8
 
 call nc_close_file(ncid)
 
@@ -897,6 +960,63 @@ function get_interp_handle(qty)
 
 end function
 
+!------------------------------------------------------------
+! calculate sensible (in-situ) temperature from 
+! local pressure, salinity, and potential temperature
+elemental function sensible_temp(potemp, s, lpres)
+
+real(r8), intent(in)  :: potemp ! potential temperature in C
+real(r8), intent(in)  :: s      ! salinity Practical Salinity Scale 1978 (PSS-78)
+real(r8), intent(in)  :: lpres  ! pressure in decibars
+real(r8) :: sensible_temp ! in-situ (sensible) temperature (C)
+
+integer  :: i,j,n
+real(r8) :: dp,p,q,r1,r2,r3,r4,r5,s1,t,x
+
+s1 = s - 35.0_r8
+p  = 0.0_r8
+t  = potemp
+
+dp = lpres - p
+n  = int (abs(dp)/1000.0_r8) + 1
+dp = dp/n
+
+do i=1,n
+   do j=1,4
+
+      r1 = ((-2.1687e-16_r8 * t + 1.8676e-14_r8) * t - 4.6206e-13_r8) * p
+      r2 = (2.7759e-12_r8*t - 1.1351e-10_r8) * s1
+      r3 = ((-5.4481e-14_r8 * t + 8.733e-12_r8) * t - 6.7795e-10_r8) * t
+      r4 = (r1 + (r2 + r3 + 1.8741e-8_r8)) * p + (-4.2393e-8_r8 * t+1.8932e-6_r8) * s1
+      r5 = r4 + ((6.6228e-10_r8 * t-6.836e-8_r8) * t + 8.5258e-6_r8) * t + 3.5803e-5_r8
+
+      x  = dp*r5
+
+      if (j == 1) then
+         t = t + 0.5_r8 * x
+         q = x
+         p = p + 0.5_r8 * dp
+      
+      else if (j == 2) then
+         t = t + 0.29298322_r8 * (x-q)
+         q = 0.58578644_r8 * x + 0.121320344_r8 * q
+
+      else if (j == 3) then
+         t = t + 1.707106781_r8 * (x-q)
+         q = 3.414213562_r8*x - 4.121320344_r8*q
+         p = p + 0.5_r8*dp
+
+      else ! j must == 4
+         t = t + (x - 2.0_r8 * q) / 6.0_r8
+
+      endif
+
+   enddo ! j loop
+enddo ! i loop
+
+sensible_temp = t
+
+end function sensible_temp
 
 !------------------------------------------------------------------
 ! Verify that the namelist was filled in correctly, and check