Skip to content
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

fix: potential temp to temp conversion for MOM6 model_interpolate #782

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

fix: potential temp to temp conversion for MOM6 model_interpolate #782

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a612e81
feat: model_interpolate for QTY_TEMPERATURE
hkershaw-brown Dec 9, 2024
ea32617
fix: _r8 was missing from numbers
hkershaw-brown Dec 11, 2024
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
178 changes: 148 additions & 30 deletions models/MOM6/model_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -220,25 +221,27 @@ 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)

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)

Expand All @@ -247,6 +250,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
Expand Down Expand Up @@ -338,6 +351,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_EVALUTATE_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_EVALUTATE_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 PSU (model) to MSU (obersvation)
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_EVALUTATE_FAILED
return
endif
expected_obs = expected_obs/1000.0_r8

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_EVALUTATE_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
Expand Down Expand Up @@ -371,27 +444,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
Expand Down Expand Up @@ -627,28 +688,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)

Expand Down Expand Up @@ -897,6 +958,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
Expand Down
Loading