diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index b450127156..89ed3bada9 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1258,9 +1258,6 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, call hchksum(Kd_heat, "after set_diffusivity Kd_heat", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) endif - ! Store the diagnosed typical diffusivity at interfaces. - if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_heat, CS%diag) - ! Set diffusivities for heat and salt separately, and possibly change the meaning of Kd_heat. if (CS%double_diffuse) then ! Add contributions from double diffusion @@ -1520,6 +1517,14 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) + if (CS%id_Kd_int > 0) then + if (CS%double_diffuse .or. CS%useKPP) then + do K=1,nz ; do j=js,je ; do i=is,ie + Kd_heat(i,j,k) = min(Kd_heat(i,j,k), Kd_salt(i,j,k)) + enddo ; enddo ; enddo + endif + call post_data(CS%id_Kd_int, Kd_heat, CS%diag) + endif if (CS%id_ea_t > 0) call post_data(CS%id_ea_t, ent_t(:,:,1:nz), CS%diag) if (CS%id_eb_t > 0) call post_data(CS%id_eb_t, ent_t(:,:,2:nz+1), CS%diag) @@ -3275,8 +3280,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di CS%id_Kd_int = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, Time, & 'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) if (CS%use_energetic_PBL) then - CS%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, Time, & - 'ePBL diapycnal diffusivity at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) + CS%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, Time, & + 'ePBL diapycnal diffusivity at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) endif CS%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, Time, & diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 2cf0ba1efe..c471b61717 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -464,11 +464,6 @@ subroutine read_Z_edges(filename, tr_name, z_edges, nz_out, has_edges, & end subroutine read_Z_edges -!### `find_overlap` and `find_limited_slope` were previously part of -! MOM_diag_to_Z.F90, and are nearly identical to `find_overlap` in -! `midas_vertmap.F90` with some slight differences. We keep it here for -! reproducibility, but the two should be merged at some point - !> Determines the layers bounded by interfaces e that overlap !! with the depth range between Z_top and Z_bot, and the fractional weights !! of each layer. It also calculates the normalized relative depths of the range @@ -620,15 +615,13 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, ! 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 + logical :: domore(SZK_(GV)) ! Records which layers need additional iterations + logical :: adjust_salt, fit_together, convergence_bug, do_any 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 - ! ### The algorithms of determine_temperature subroutine needs to be reexamined. - - call log_version(PF, mdl, version, "") ! We should switch the default to the newer method which simultaneously adjusts @@ -638,7 +631,13 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, "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_CONVERGENCE_BUG", convergence_bug, & + "If true, use layout-dependent tests on the changes in temperature and salinity "//& + "to determine when the iterations have converged when DETERMINE_TEMP_ADJUST_T_AND_S "//& + "is false. For realistic equations of state and the default values of the "//& + "various tolerances, this bug does not impact the solutions.", & + default=.true., do_not_log=just_read) !### Change the default to false. + 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) @@ -653,10 +652,12 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, units="ppt", 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) + units="degC", default=1.0e-4, scale=US%degC_to_C, & + do_not_log=just_read.or.(.not.convergence_bug)) call get_param(PF, mdl, "DETERMINE_TEMP_S_TOLERANCE", tol_S, & "The convergence tolerance for temperature in determine_temperature.", & - units="ppt", default=1.0e-4, scale=US%ppt_to_S, do_not_log=just_read) + units="ppt", default=1.0e-4, scale=US%ppt_to_S, & + do_not_log=just_read.or.(.not.convergence_bug)) 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) @@ -689,49 +690,69 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, T(:,:) = temp(:,j,:) S(:,:) = salt(:,j,:) dT(:,:) = 0.0 + domore(:) = .true. adjust_salt = .true. iter_loop: do itt = 1,niter - do k=1,nz + do k=k_start,nz ; if (domore(k)) then + domore(k) = .false. 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), & EOS, EOSdom ) - enddo - do k=k_start,nz ; do i=is,ie -! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln) then - if (abs(rho(i,k)-R_tgt(k))>tol_rho) 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 - I_denom = 1.0 / (drho_dS(i,k)**2 + dT_dS_gauge**2*drho_dT(i,k)**2) - dS(i,k) = (R_tgt(k)-rho(i,k)) * drho_dS(i,k) * I_denom - dT(i,k) = (R_tgt(k)-rho(i,k)) * dT_dS_gauge**2*drho_dT(i,k) * I_denom + do i=is,ie +! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln) then + if (abs(rho(i,k)-R_tgt(k))>tol_rho) then + domore(k) = .true. + 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 + I_denom = 1.0 / (drho_dS(i,k)**2 + dT_dS_gauge**2*drho_dT(i,k)**2) + dS(i,k) = (R_tgt(k)-rho(i,k)) * drho_dS(i,k) * I_denom + dT(i,k) = (R_tgt(k)-rho(i,k)) * dT_dS_gauge**2*drho_dT(i,k) * I_denom - T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) - S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) + T(i,k) = max(min(T(i,k)+dT(i,k), T_max), T_min) + S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) + endif endif + enddo + endif ; enddo + if (convergence_bug) then + ! If this test does anything, it is layout-dependent. + if (maxval(abs(dT)) < tol_T) then + adjust_salt = .false. + exit iter_loop endif - enddo ; enddo - if (maxval(abs(dT)) < tol_T) then - adjust_salt = .false. - exit iter_loop endif + + do_any = .false. + do k=k_start,nz ; if (domore(k)) do_any = .true. ; enddo + if (.not.do_any) exit iter_loop ! Further iterations will not change anything. enddo iter_loop if (adjust_salt .and. .not.fit_together) then ; do itt = 1,niter - do k=1,nz + do k=k_start,nz ; if (domore(k)) then + domore(k) = .false. 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), & EOS, EOSdom ) - enddo - do k=k_start,nz ; do i=is,ie -! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln ) then - if (abs(rho(i,k)-R_tgt(k)) > tol_rho) then - dS(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dS(i,k), max_s_adj), -max_s_adj) - S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) - endif - enddo ; enddo - if (maxval(abs(dS)) < tol_S) exit + do i=is,ie +! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln ) then + if (abs(rho(i,k)-R_tgt(k)) > tol_rho) then + dS(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dS(i,k), max_s_adj), -max_s_adj) + S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) + domore(k) = .true. + endif + enddo + endif ; enddo + + if (convergence_bug) then + ! If this test does anything, it is layout-dependent. + if (maxval(abs(dS)) < tol_S) exit + endif + + do_any = .false. + do k=k_start,nz ; if (domore(k)) do_any = .true. ; enddo + if (.not.do_any) exit ! Further iterations will not change anything enddo ; endif temp(:,j,:) = T(:,:)