diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9e358ab17e..8d45114a39 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3299,17 +3299,18 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & CS%mixedlayer_restrat = mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, & CS%mixedlayer_restrat_CSp, restart_CSp) - if (CS%mixedlayer_restrat) then - if (GV%Boussinesq .and. associated(CS%visc%h_ML)) then - ! This is here to allow for a transition of restart files between model versions. - call get_param(param_file, "MOM", "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, & - default=.false., do_not_log=.true.) - if (MLE_use_PBL_MLD .and. .not.query_initialized(CS%visc%h_ML, "h_ML", restart_CSp) .and. & - associated(CS%visc%MLD)) then - do j=js,je ; do i=is,ie ; CS%visc%h_ML(i,j) = GV%Z_to_H * CS%visc%MLD(i,j) ; enddo ; enddo - endif + + if (GV%Boussinesq .and. associated(CS%visc%h_ML)) then + ! This is here to allow for a transition of restart files between model versions. + call get_param(param_file, "MOM", "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, & + default=.false., do_not_log=.true.) + if (MLE_use_PBL_MLD .and. .not.query_initialized(CS%visc%h_ML, "h_ML", restart_CSp) .and. & + associated(CS%visc%MLD)) then + do j=js,je ; do i=is,ie ; CS%visc%h_ML(i,j) = GV%Z_to_H * CS%visc%MLD(i,j) ; enddo ; enddo endif + endif + if (CS%mixedlayer_restrat) then if (.not.(bulkmixedlayer .or. CS%use_ALE_algorithm)) & call MOM_error(FATAL, "MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.") ! When DIABATIC_FIRST=False and using CS%visc%ML in mixedlayer_restrat we need to update after a restart @@ -4056,7 +4057,7 @@ subroutine extract_surface_state(CS, sfc_state_in) endif endif - if (CS%debug) call MOM_surface_chksum("Post extract_sfc", sfc_state, G, US, haloshift=0) + if (CS%debug) call MOM_surface_chksum("Post extract_sfc", sfc_state, G, US, haloshift=0, symmetric=.true.) ! Rotate sfc_state back onto the input grid, sfc_state_in if (CS%rotate_index) then diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index 2eb84a1c74..6b89323475 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -167,7 +167,7 @@ subroutine MOM_surface_chksum(mesg, sfc_state, G, US, haloshift, symmetric) logical :: sym sym = .false. ; if (present(symmetric)) sym = symmetric - hs = 1 ; if (present(haloshift)) hs = haloshift + hs = 0 ; if (present(haloshift)) hs = haloshift if (allocated(sfc_state%SST)) call hchksum(sfc_state%SST, mesg//" SST", G%HI, haloshift=hs, & unscale=US%C_to_degC) @@ -182,6 +182,14 @@ subroutine MOM_surface_chksum(mesg, sfc_state, G, US, haloshift, symmetric) unscale=US%L_T_to_m_s) if (allocated(sfc_state%frazil)) call hchksum(sfc_state%frazil, mesg//" frazil", G%HI, & haloshift=hs, unscale=US%Q_to_J_kg*US%RZ_to_kg_m2) + if (allocated(sfc_state%melt_potential)) call hchksum(sfc_state%melt_potential, mesg//" melt_potential", & + G%HI, haloshift=hs, unscale=US%Q_to_J_kg*US%RZ_to_kg_m2) + if (allocated(sfc_state%ocean_mass)) call hchksum(sfc_state%ocean_mass, mesg//" ocean_mass", & + G%HI, haloshift=hs, unscale=US%RZ_to_kg_m2) + if (allocated(sfc_state%ocean_heat)) call hchksum(sfc_state%ocean_heat, mesg//" ocean_heat", & + G%HI, haloshift=hs, unscale=US%C_to_degC*US%RZ_to_kg_m2) + if (allocated(sfc_state%ocean_salt)) call hchksum(sfc_state%ocean_salt, mesg//" ocean_salt", & + G%HI, haloshift=hs, unscale=US%S_to_ppt*US%RZ_to_kg_m2) end subroutine MOM_surface_chksum diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 0180cdd2c0..a66b53656c 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -1526,6 +1526,7 @@ subroutine EOS_init(param_file, EOS, US) "If true, use a bug in the calculation of the second derivatives of density "//& "with temperature and with temperature and pressure that causes some terms "//& "to be only 2/3 of what they should be.", default=.false.) + call EOS_manual_init(EOS, form_of_EOS=EOS_WRIGHT, use_Wright_2nd_deriv_bug=EOS%use_Wright_2nd_deriv_bug) endif EOS_quad_default = .not.((EOS%form_of_EOS == EOS_LINEAR) .or. & @@ -1645,6 +1646,8 @@ subroutine EOS_manual_init(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Co select type (t => EOS%type) type is (linear_EOS) call t%set_params_linear(Rho_T0_S0, dRho_dT, dRho_dS) + type is (buggy_Wright_EOS) + call t%set_params_buggy_Wright(use_Wright_2nd_deriv_bug) end select endif if (present(form_of_TFreeze)) EOS%form_of_TFreeze = form_of_TFreeze diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index e71ae1791d..f71b5c9a2c 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -12,6 +12,7 @@ module MOM_EOS_Wright public buggy_Wright_EOS public int_density_dz_wright, int_spec_vol_dp_wright public avg_spec_vol_buggy_Wright +public set_params_buggy_Wright !>@{ Parameters in the Wright equation of state using the reduced range formula, which is a fit to the UNESCO ! equation of state for the restricted range: -2 < theta < 30 [degC], 28 < S < 38 [PSU], 0 < p < 5e7 [Pa]. @@ -39,6 +40,8 @@ module MOM_EOS_Wright !> The EOS_base implementation of the Wright 1997 equation of state with some bugs type, extends (EOS_base) :: buggy_Wright_EOS + real :: three = 3.0 !< A constant that can be adjusted to recreate some bugs [nondim] + contains !> Implementation of the in-situ density as an elemental function [kg m-3] procedure :: density_elem => density_elem_buggy_Wright @@ -59,6 +62,9 @@ module MOM_EOS_Wright !> Implementation of the range query function procedure :: EOS_fit_range => EOS_fit_range_buggy_Wright + !> Instance specific function to set internal parameters + procedure :: set_params_buggy_Wright => set_params_buggy_Wright + !> Local implementation of generic calculate_density_array for efficiency procedure :: calculate_density_array => calculate_density_array_buggy_Wright !> Local implementation of generic calculate_spec_vol_array for efficiency @@ -228,13 +234,16 @@ elemental subroutine calculate_density_second_derivs_elem_buggy_Wright(this, T, real :: z11 ! A local work variable [Pa m2 s-2 PSU-1] = [kg m s-4 PSU-1] real :: z2_2 ! A local work variable [m4 s-4] real :: z2_3 ! A local work variable [m6 s-6] + real :: six ! A constant that can be adjusted from 6. to 4. to recreate a bug [nondim] ! Based on the above expression with common terms factored, there probably exists a more numerically stable ! and/or efficient expression + six = 2.0*this%three ! When recreating a bug from the original version of this routine, six = 4. + z0 = T*(b1 + b5*S + T*(b2 + b3*T)) z1 = (b0 + pressure + b4*S + z0) - z3 = (b1 + b5*S + T*(2.*b2 + 2.*b3*T)) ! BUG: This should be z3 = b1 + b5*S + T*(2.*b2 + 3.*b3*T) + z3 = (b1 + b5*S + T*(2.*b2 + this%three*b3*T)) ! When recreating a bug here this%three = 2. z4 = (c0 + c4*S + T*(c1 + c5*S + T*(c2 + c3*T))) z5 = (b1 + b5*S + T*(b2 + b3*T) + T*(b2 + 2.*b3*T)) z6 = c1 + c5*S + T*(c2 + c3*T) + T*(c2 + 2.*c3*T) @@ -249,8 +258,7 @@ elemental subroutine calculate_density_second_derivs_elem_buggy_Wright(this, T, drho_ds_ds = (z10*(c4 + c5*T) - a2*z10*z1 - z10*z7)/z2_2 - (2.*(c4 + c5*T + z9*z10 + a2*z1)*z11)/z2_3 drho_ds_dt = (z10*z6 - z1*(c5 + a2*z5) + b5*z4 - z5*z7)/z2_2 - (2.*(z6 + z9*z5 + a1*z1)*z11)/z2_3 - ! BUG: In the following line: (2.*b2 + 4.*b3*T) should be (2.*b2 + 6.*b3*T) - drho_dt_dt = (z3*z6 - z1*(2.*c2 + 6.*c3*T + a1*z5) + (2.*b2 + 4.*b3*T)*z4 - z5*z8)/z2_2 - & + drho_dt_dt = (z3*z6 - z1*(2.*c2 + 6.*c3*T + a1*z5) + (2.*b2 + six*b3*T)*z4 - z5*z8)/z2_2 - & (2.*(z6 + z9*z5 + a1*z1)*(z3*z4 - z1*z8))/z2_3 drho_ds_dp = (-c4 - c5*T - 2.*a2*z1)/z2_2 - (2.*z9*z11)/z2_3 drho_dt_dp = (-c1 - c5*S - T*(2.*c2 + 3.*c3*T) - 2.*a1*z1)/z2_2 - (2.*z9*(z3*z4 - z1*z8))/z2_3 @@ -933,6 +941,23 @@ subroutine calculate_spec_vol_array_buggy_Wright(this, T, S, pressure, specvol, end subroutine calculate_spec_vol_array_buggy_Wright +!> Set coefficients that can correct bugs un the buggy Wright equation of state. +subroutine set_params_buggy_Wright(this, use_Wright_2nd_deriv_bug) + class(buggy_Wright_EOS), intent(inout) :: this !< This EOS + logical, optional, intent(in) :: use_Wright_2nd_deriv_bug !< If true, use a buggy + !! buggy version of the calculations of the second + !! derivative of density with temperature and with temperature and + !! pressure. This bug is corrected in the default version. + + this%three = 3.0 + if (present(use_Wright_2nd_deriv_bug)) then + if (use_Wright_2nd_deriv_bug) then ; this%three = 2.0 + else ; this%three = 3.0 ; endif + endif + +end subroutine set_params_buggy_Wright + + !> \namespace mom_eos_wright !! !! \section section_EOS_Wright Wright equation of state diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 541e349d29..c28e2e5896 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -14,7 +14,6 @@ module MOM_diag_mediator use MOM_diag_manager_infra, only : send_data_infra, MOM_diag_field_add_attribute, EAST, NORTH use MOM_diag_manager_infra, only : register_diag_field_infra, register_static_field_infra use MOM_diag_manager_infra, only : get_MOM_diag_field_id, DIAG_FIELD_NOT_FOUND -use MOM_diag_manager_infra, only : diag_send_complete_infra use MOM_diag_remap, only : diag_remap_ctrl, diag_remap_update, diag_remap_calc_hmask use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap use MOM_diag_remap, only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field @@ -2079,10 +2078,8 @@ end subroutine enable_averages subroutine disable_averaging(diag_cs) type(diag_ctrl), intent(inout) :: diag_CS !< Structure used to regulate diagnostic output - call diag_send_complete_infra() diag_cs%time_int = 0.0 diag_cs%ave_enabled = .false. - end subroutine disable_averaging !> Call this subroutine to determine whether the averaging is diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 6d8afa1493..18bbc45e13 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -676,6 +676,14 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim call cpu_clock_begin(id_clock_kpp) ! total vertical viscosity in the interior is represented via visc%Kv_shear + ! NOTE: The following do not require initialization, but their checksums do + ! require initialization, and past versions were initialized to zero. + KPP_NLTheat(:,:,:) = 0. + KPP_NLTscalar(:,:,:) = 0. + KPP_buoy_flux(:,:,:) = 0. + KPP_temp_flux(:,:) = 0. + KPP_salt_flux(:,:) = 0. + ! KPP needs the surface buoyancy flux but does not update state variables. ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes. ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux @@ -1285,6 +1293,14 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k) enddo ; enddo ; enddo + ! NOTE: The following do not require initialization, but their checksums do + ! require initialization, and past versions were initialized to zero. + KPP_NLTheat(:,:,:) = 0. + KPP_NLTscalar(:,:,:) = 0. + KPP_buoy_flux(:,:,:) = 0. + KPP_temp_flux(:,:) = 0. + KPP_salt_flux(:,:) = 0. + ! KPP needs the surface buoyancy flux but does not update state variables. ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes. ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux @@ -1890,6 +1906,15 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e if (CS%useKPP) then call cpu_clock_begin(id_clock_kpp) + + ! NOTE: The following do not require initialization, but their checksums do + ! require initialization, and past versions were initialized to zero. + KPP_NLTheat(:,:,:) = 0. + KPP_NLTscalar(:,:,:) = 0. + KPP_buoy_flux(:,:,:) = 0. + KPP_temp_flux(:,:) = 0. + KPP_salt_flux(:,:) = 0. + ! KPP needs the surface buoyancy flux but does not update state variables. ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes. ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 88845a390f..ff2d178adc 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1180,7 +1180,8 @@ subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV) ! Each cell extends from x=-1/2 to 1/2, and has a topography ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. - crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + !crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + crv_3 = (Dp + Dm - (2.0*D_vel)) ; crv = 3.0*crv_3 slope = Dp - Dm ! Calculate the volume above which the entire cell is open and the volume at which the @@ -1207,10 +1208,13 @@ subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV) ! vol_below(K) = 0.5*L^2*(slope + crv/3*(3-4L)). if (a2x48_apb3*vol_below(K) < 1e-8) then ! Could be 1e-7? ! There is a very good approximation here for massless layers. - L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0) + !L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0) + L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + (ax2_3apb*L0)) else + !L(K) = apb_4a * (1.0 - & + ! 2.0 * cos(C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3)) L(K) = apb_4a * (1.0 - & - 2.0 * cos(C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3)) + 2.0 * cos(C1_3*acos((a2x48_apb3*vol_below(K)) - 1.0) - C2pi_3)) endif ! To check the answers. ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index b77949342d..e927f2f89d 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -34,6 +34,8 @@ module MOM_tracer_advect logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: usePPM !< If true, use PPM instead of PLM logical :: useHuynh !< If true, use the Huynh scheme for PPM interface values + logical :: useHuynhStencilBug = .false. !< If true, use the incorrect stencil width. + !! This is provided for compatibility with legacy simuations. type(group_pass_type) :: pass_uhr_vhr_t_hprev !< A structure used for group passes end type tracer_advect_CS @@ -94,6 +96,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first ! can be simply discarded [H L2 ~> m3 or kg]. real :: landvolfill ! An arbitrary? nonzero cell volume [H L2 ~> m3 or kg]. + logical :: use_PPM_stencil ! If true, use the correct PPM stencil width. real :: Idt ! 1/dt [T-1 ~> s-1]. logical :: domore_u(SZJ_(G),SZK_(GV)) ! domore_u and domore_v indicate whether there is more logical :: domore_v(SZJB_(G),SZK_(GV)) ! advection to be done in the corresponding row or column. @@ -112,7 +115,7 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB landvolfill = 1.0e-20 ! This is arbitrary, but must be positive. - stencil = 2 ! The scheme's stencil; 2 for PLM and PPM:H3 + stencil = 2 ! The scheme's stencil; 2 for PLM if (.not. associated(CS)) call MOM_error(FATAL, "MOM_tracer_advect: "// & "tracer_advect_init must be called before advect_tracer.") @@ -123,8 +126,8 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first x_first = (MOD(G%first_direction,2) == 0) ! increase stencil size for Colella & Woodward PPM -! if (CS%usePPM .and. .not. CS%useHuynh) stencil = 3 - if (CS%usePPM) stencil = 3 + use_PPM_stencil = CS%usePPM .and. .not. CS%useHuynhStencilBug + if (use_PPM_stencil) stencil = 3 ntr = Reg%ntr Idt = 1.0 / dt @@ -1130,6 +1133,16 @@ subroutine tracer_advect_init(Time, G, US, param_file, diag, CS) "Unknown TRACER_ADVECTION_SCHEME = "//trim(mesg)) end select + if (CS%useHuynh) then + call get_param(param_file, mdl, "USE_HUYNH_STENCIL_BUG", & + CS%useHuynhStencilBug, & + desc="If true, use a stencil width of 2 in PPM:H3 tracer advection. " & + // "This is incorrect and will produce regressions in certain " & + // "configurations, but may be required to reproduce results in " & + // "legacy simulations.", & + default=.false.) + endif + id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=CLOCK_MODULE) id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=CLOCK_ROUTINE) id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=CLOCK_ROUTINE)