Skip to content

Commit

Permalink
+Add 11 runtime params for determine_temperature
Browse files Browse the repository at this point in the history 
  Added 11 new runtime parameters (DETERMINE_TEMP_ADJUST_T_AND_S,
DETERMINE_TEMP_T_MIN, DETERMINE_TEMP_T_MAX, DETERMINE_TEMP_S_MIN,
DETERMINE_TEMP_S_MAX, DETERMINE_TEMP_T_TOLERANCE, DETERMINE_TEMP_S_TOLERANCE,
DETERMINE_TEMP_RHO_TOLERANCE, DETERMINE_TEMP_DT_DS_WEIGHT,
DETERMINE_TEMP_T_ADJ_RANGE, and DETERMINE_TEMP_S_ADJ_RANGE) to replace hard
coded dimensional parameters used in the routine determine_temperature.  This
change also requires that a param_file_type argument and the logical argument
just_read are passed to determine temperature.  By default, all answers are
bitwise identical but there are up to 10 new entries in the MOM_parameter_doc
files for some layer-mode configurations with INIT_LAYERS_FROM_Z_FILE and
FIT_TO_TARGET_DENSITY_IC set to true.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Dec 20, 2022
1 parent e9300da commit dda6d45
Show file tree
Hide file tree
Showing 2 changed files with 70 additions and 22 deletions.
8 changes: 6 additions & 2 deletions src/initialization/MOM_state_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2737,6 +2737,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
units="1e-3", default=1.0e-3, scale=US%ppt_to_S, do_not_log=just_read)

if (just_read) then
if ((.not.useALEremapping) .and. adjust_temperature) &
! This call is just here to read and log the determine_temperature parameters
call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), eos, tv%P_Ref, 0, &
h, 0, G, GV, US, PF, just_read=.true.)
call cpu_clock_end(id_clock_routine)
return ! All run-time parameters have been read, so return.
endif
Expand Down Expand Up @@ -2957,8 +2961,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
if (adjust_temperature) then
! Finally adjust to target density
ks = 1 ; if (separate_mixed_layer) ks = GV%nk_rho_varies + 1
call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), tv%P_Ref, niter, &
h, ks, G, GV, US, eos)
call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), eos, tv%P_Ref, niter, &
h, ks, G, GV, US, PF, just_read)
endif

endif ! useALEremapping
Expand Down
84 changes: 64 additions & 20 deletions src/tracer/MOM_tracer_Z_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ module MOM_tracer_Z_init
! This file is part of MOM6. See LICENSE.md for the license.

use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe
! use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_grid, only : ocean_grid_type
use MOM_io, only : MOM_read_data, get_var_sizes, read_attribute, read_variable
use MOM_io, only : open_file_to_read, close_file_to_read
Expand Down Expand Up @@ -556,23 +556,27 @@ end function find_limited_slope

!> This subroutine determines the potential temperature and salinity that
!! is consistent with the target density using provided initial guess
subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G, GV, US, &
EOS, h_massless)
subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_start, G, GV, US, &
PF, just_read, h_massless)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: temp !< potential temperature [C ~> degC]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: salt !< salinity [S ~> ppt]
real, dimension(SZK_(GV)), intent(in) :: R_tgt !< desired potential density [R ~> kg m-3].
type(EOS_type), intent(in) :: EOS !< seawater equation of state control structure
real, intent(in) :: p_ref !< reference pressure [R L2 T-2 ~> Pa].
integer, intent(in) :: niter !< maximum number of iterations
integer, intent(in) :: k_start !< starting index (i.e. below the buffer layer)
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< layer thickness, used only to avoid working on
!! massless layers [H ~> m or kg m-2]
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(EOS_type), intent(in) :: EOS !< seawater equation of state control structure
type(param_file_type), intent(in) :: PF !< A structure indicating the open file
!! to parse for model parameter values.
logical, intent(in) :: just_read !< If true, this call will only read
!! parameters without changing T or S.
real, optional, intent(in) :: h_massless !< A threshold below which a layer is
!! determined to be massless [H ~> m or kg m-2]

Expand Down Expand Up @@ -600,29 +604,69 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G,
! when old_fit is true [C ~> degC]
real :: max_s_adj ! The largest permitted salinity changes with each iteration
! when old_fit is true [S ~> ppt]
logical :: adjust_salt, old_fit
! This include declares and sets the variable "version".
# include "version_variable.h"
character(len=40) :: mdl = "determine_temperature" ! This subroutine's name.
logical :: adjust_salt, fit_together
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
integer :: i, j, k, is, ie, js, je, nz, itt

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

! These hard coded parameters need to be set properly.
S_min = 0.5*US%ppt_to_S ; S_max = 65.0*US%ppt_to_S
T_max = 31.0*US%degC_to_C ; T_min = -2.0*US%degC_to_C
max_t_adj = 1.0*US%degC_to_C
max_s_adj = 0.5*US%ppt_to_S
tol_T = 1.0e-4*US%degC_to_C
tol_S = 1.0e-4*US%ppt_to_S
tol_rho = 1.0e-4*US%kg_m3_to_R
old_fit = .true. ! reproduces siena behavior
! ### The algorithms of determine_temperature subroutine needs to be reexamined.

dT_dS_gauge = 10.0*US%degC_to_C*US%S_to_ppt ! 10 degC is weighted equivalently to 1 ppt.

! ### The whole determine_temperature subroutine needs to be reexamined, both the algorithms
! and the extensive use of hard-coded dimensional parameters.
call log_version(PF, mdl, version, "")

! We will switch to the newer method which simultaneously adjusts
! We should switch the default to the newer method which simultaneously adjusts
! temp and salt based on the ratio of the thermal and haline coefficients, once it is tested.
call get_param(PF, mdl, "DETERMINE_TEMP_ADJUST_T_AND_S", fit_together, &
"If true, simltaneously adjust the estimates of the temperature and salinity "//&
"based on the ratio of the thermal and haline coefficients. Otherwise try to "//&
"match the density by only adjusting temperatures within a maximum range before "//&
"revising estimates of the salinity.", default=.false., do_not_log=just_read)
! These hard coded parameters need to be set properly.
call get_param(PF, mdl, "DETERMINE_TEMP_T_MIN", T_min, &
"The minimum temperature that can be found by determine_temperature.", &
units="degC", default=-2.0, scale=US%degC_to_C, do_not_log=just_read)
call get_param(PF, mdl, "DETERMINE_TEMP_T_MAX", T_max, &
"The maximum temperature that can be found by determine_temperature.", &
units="degC", default=31.0, scale=US%degC_to_C, do_not_log=just_read)
call get_param(PF, mdl, "DETERMINE_TEMP_S_MIN", S_min, &
"The minimum salinity that can be found by determine_temperature.", &
units="1e-3", default=0.5, scale=US%ppt_to_S, do_not_log=just_read)
call get_param(PF, mdl, "DETERMINE_TEMP_S_MAX", S_max, &
"The maximum salinity that can be found by determine_temperature.", &
units="1e-3", default=65.0, scale=US%ppt_to_S, do_not_log=just_read)
call get_param(PF, mdl, "DETERMINE_TEMP_T_TOLERANCE", tol_T, &
"The convergence tolerance for temperature in determine_temperature.", &
units="degC", default=1.0e-4, scale=US%degC_to_C, do_not_log=just_read)
call get_param(PF, mdl, "DETERMINE_TEMP_S_TOLERANCE", tol_S, &
"The convergence tolerance for temperature in determine_temperature.", &
units="1e-3", default=1.0e-4, scale=US%ppt_to_S, do_not_log=just_read)
call get_param(PF, mdl, "DETERMINE_TEMP_RHO_TOLERANCE", tol_rho, &
"The convergence tolerance for density in determine_temperature.", &
units="kg m-3", default=1.0e-4, scale=US%kg_m3_to_R, do_not_log=just_read)
if (fit_together) then
! By default 10 degC is weighted equivalently to 1 ppt when minimizing changes.
call get_param(PF, mdl, "DETERMINE_TEMP_DT_DS_WEIGHT", dT_dS_gauge, &
"When extrapolating T & S to match the layer target densities, this "//&
"factor (in deg C / PSU) is combined with the derivatives of density "//&
"with T & S to determine what direction is orthogonal to density contours. "//&
"It could be based on a typical value of (dR/dS) / (dR/dT) in oceanic profiles.", &
units="degC PSU-1", default=10.0, scale=US%degC_to_C*US%S_to_ppt)
else
call get_param(PF, mdl, "DETERMINE_TEMP_T_ADJ_RANGE", max_t_adj, &
"The maximum amount by which the initial layer temperatures can be "//&
"modified in determine_temperature.", &
units="degC", default=1.0, scale=US%degC_to_C, do_not_log=just_read)
call get_param(PF, mdl, "DETERMINE_TEMP_S_ADJ_RANGE", max_S_adj, &
"The maximum amount by which the initial layer salinities can be "//&
"modified in determine_temperature.", &
units="1e-3", default=0.5, scale=US%ppt_to_S, do_not_log=just_read)
endif

if (just_read) return ! All run-time parameters have been read, so return.

press(:) = p_ref
EOSdom(:) = EOS_domain(G%HI)
Expand All @@ -643,7 +687,7 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G,
do k=k_start,nz ; do i=is,ie
! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln) then
if (abs(rho(i,k)-R_tgt(k))>tol_rho) then
if (old_fit) then
if (.not.fit_together) then
dT(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj)
T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min)
else
Expand All @@ -662,7 +706,7 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G,
endif
enddo iter_loop

if (adjust_salt .and. old_fit) then ; do itt = 1,niter
if (adjust_salt .and. .not.fit_together) then ; do itt = 1,niter
do k=1,nz
call calculate_density(T(:,k), S(:,k), press, rho(:,k), EOS, EOSdom )
call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), &
Expand Down

0 comments on commit dda6d45

Please sign in to comment.