diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 97425b3ea..fe8172245 100644
--- a/CHANGELOG.rst
+++ b/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
diff --git a/conf.py b/conf.py
index 59df22338..7f26293ab 100644
--- a/conf.py
+++ b/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 ---------------------------------------------------
diff --git a/models/MOM6/model_mod.f90 b/models/MOM6/model_mod.f90
index e6ebcd03f..e0b88ae28 100644
--- a/models/MOM6/model_mod.f90
+++ b/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,12 +446,7 @@ 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
@@ -384,14 +454,7 @@ subroutine model_interpolate(state_handle, ens_size, location, qty, expected_obs
! 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