From e2e52ea46a7f70d9e460b71d29014c4bf071feb2 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 24 Apr 2024 09:35:43 -0400 Subject: [PATCH] +*Use visc%h_ML in tracer diffusion routines Pass visc%h_ML to applyBoundaryFluxesInOut, which avoids the need for the thickness conversions inside of that routine. Answers are bitwise identical for Boussinesq cases, but they will change in non-Boussinesq cases because the layer specific volumes have changed between the time the boundary layer thicknesses have been calculated and when they are being used in the brine plume parameterization. Also use visc%h_ML in place of calls to KPP_get_BLD or energetic_PBL_get_MLD in step_MOM_dyn_split_RK2, hor_bnd_diffusion and neutral_diffusion_calc_coeffs. This change requires a new vertvisc_type argument to tracer_hordiff, hor_bnd_diffusion and neutral_diffusion_calc_coeffs. Boussinesq answers are bitwise identical, but non-Boussinesq answers change because they now use the actual model specific volumes rather than some prescribed constant value to to the relevant unit conversions. The logic in set_visc_register_restarts was also updated so that visc%MLD and visc%h_ML are allocated and registered for restarts in the cases where they will be used. All answers are bitwise identical in Boussinesq mode, but in non-Boussinesq mode there are changes in answers in cases that use neutral tracer diffusion that excludes the boundary layer, horizontal tracer diffusion in the boundary layer or FPMIX. There are also changes to the arguments of several publicly visible routines. --- src/core/MOM.F90 | 8 ++-- src/core/MOM_dynamics_split_RK2.F90 | 9 +---- src/core/MOM_dynamics_split_RK2b.F90 | 9 +---- .../vertical/MOM_CVMix_ddiff.F90 | 4 +- .../vertical/MOM_diabatic_aux.F90 | 39 +++++++------------ .../vertical/MOM_diabatic_driver.F90 | 8 ++-- .../vertical/MOM_set_viscosity.F90 | 29 +++++++++++--- src/tracer/MOM_hor_bnd_diffusion.F90 | 24 ++++++++---- src/tracer/MOM_neutral_diffusion.F90 | 18 ++++++--- src/tracer/MOM_tracer_hor_diff.F90 | 16 ++++---- 10 files changed, 87 insertions(+), 77 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 238e4f4f09..bdfe55424f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1432,7 +1432,7 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) call advect_tracer(h, CS%uhtr, CS%vhtr, CS%OBC, CS%t_dyn_rel_adv, G, GV, US, & CS%tracer_adv_CSp, CS%tracer_Reg, x_first_in=x_first) if (CS%debug) call MOM_tracer_chksum("Post-advect ", CS%tracer_Reg, G) - call tracer_hordiff(h, CS%t_dyn_rel_adv, CS%MEKE, CS%VarMix, G, GV, US, & + call tracer_hordiff(h, CS%t_dyn_rel_adv, CS%MEKE, CS%VarMix, CS%visc, G, GV, US, & CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) if (CS%debug) call MOM_tracer_chksum("Post-diffuse ", CS%tracer_Reg, G) if (showCallTree) call callTree_waypoint("finished tracer advection/diffusion (step_MOM)") @@ -1881,7 +1881,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif - call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, & + call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, CS%visc, G, GV, US, & CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) endif endif @@ -1908,7 +1908,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS call calc_depth_function(G, CS%VarMix) call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC) endif - call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, & + call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, CS%visc, G, GV, US, & CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) endif endif @@ -1965,7 +1965,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS CS%h, eatr, ebtr, uhtr, vhtr) ! Perform offline diffusion if requested if (.not. skip_diffusion) then - call tracer_hordiff(h_end, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, & + call tracer_hordiff(h_end, dt_offline, CS%MEKE, CS%VarMix, CS%visc, G, GV, US, & CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) endif diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 0557ec7cd5..85469811b2 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -71,9 +71,6 @@ module MOM_dynamics_split_RK2 use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units use MOM_wave_interface, only : wave_parameters_CS, Stokes_PGF -use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS -use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS -use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member implicit none ; private @@ -139,8 +136,6 @@ module MOM_dynamics_split_RK2 real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce !< pbce times eta gives the baroclinic pressure !! anomaly in each layer due to free surface height !! anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]. - type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to ge - type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure real, pointer, dimension(:,:) :: taux_bot => NULL() !< frictional x-bottom stress from the ocean !! to the seafloor [R L Z T-2 ~> Pa] @@ -717,9 +712,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f if (CS%fpmix) then hbl(:,:) = 0.0 - if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) - if (ASSOCIATED(CS%energetic_PBL_CSp)) & - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) + if (associated(visc%h_ML)) hbl(:,:) = visc%h_ML(:,:) call vertFPmix(up, vp, uold, vold, hbl, h, forces, & dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, & diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index 44a0b0bf5c..42f5590ba6 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -73,9 +73,6 @@ module MOM_dynamics_split_RK2b use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units use MOM_wave_interface, only: wave_parameters_CS, Stokes_PGF -use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS -use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS -use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member implicit none ; private @@ -147,8 +144,6 @@ module MOM_dynamics_split_RK2b real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: dv_av_inst !< The barotropic meridional velocity increment !! between filtered and instantaneous velocities !! [L T-1 ~> m s-1] - type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to ge - type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure real, pointer, dimension(:,:) :: taux_bot => NULL() !< frictional x-bottom stress from the ocean !! to the seafloor [R L Z T-2 ~> Pa] @@ -734,9 +729,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc ! if (CS%fpmix) then ! hbl(:,:) = 0.0 - ! if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) - ! if (ASSOCIATED(CS%energetic_PBL_CSp)) & - ! call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) + ! if (associated(visc%h_ML)) hbl(:,:) = visc%h_ML(:,:) ! call vertFPmix(up, vp, uold, vold, hbl, h, forces, & ! dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC) ! call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, & diff --git a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 index af17e0287f..958b2478f3 100644 --- a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 @@ -245,8 +245,8 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho) iFaceHeight(k+1) = iFaceHeight(k) - dh enddo - ! gets index of the level and interface above hbl - !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight, GV%Z_to_H*hbl(i,j)) + ! gets index of the level and interface above hbl in [H ~> m or kg m-2] + !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight, hbl(i,j)) Kd1_T(:) = 0.0 ; Kd1_S(:) = 0.0 call CVMix_coeffs_ddiff(Tdiff_out=Kd1_T(:), & diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 1ba5aef392..aa31024b24 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1064,7 +1064,7 @@ end subroutine diagnoseMLDbyEnergy subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, & aggregate_FW_forcing, evap_CFL_limit, & minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, & - SkinBuoyFlux, MLD) + SkinBuoyFlux, MLD_h) type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -1094,7 +1094,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !! salinity [R-1 S-1 ~> m3 kg-1 ppt-1]. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3]. - real, pointer, dimension(:,:), optional :: MLD !< Mixed layer depth for brine plumes [Z ~> m] + real, dimension(:,:), & + optional, pointer :: MLD_h !< Mixed layer thickness for brine plumes [H ~> m or kg m-2] ! Local variables integer, parameter :: maxGroundings = 5 @@ -1137,8 +1138,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] netMassInOut_rate, & ! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1] mixing_depth, & ! The mixing depth for brine plumes [H ~> m or kg m-2] - MLD_H, & ! The mixed layer depth for brine plumes in thickness units [H ~> m or kg m-2] - MLD_Z, & ! Running sum of distance from the surface for finding MLD_H [Z ~> m] total_h ! Total thickness of the water column [H ~> m or kg m-2] real, dimension(SZI_(G), SZK_(GV)) :: & h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2] @@ -1204,11 +1203,15 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0 endif - if (CS%do_brine_plume .and. .not. associated(MLD)) then + if (CS%do_brine_plume .and. .not.present(MLD_h)) then call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//& - "Brine plume parameterization requires a mixed-layer depth,\n"//& + "Brine plume parameterization requires a mixed-layer depth argument,\n"//& "currently coming from the energetic PBL scheme.") endif + if (CS%do_brine_plume .and. .not.associated(MLD_h)) then + call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//& + "Brine plume parameterization requires an associated mixed-layer depth.") + endif if (CS%do_brine_plume .and. .not. associated(fluxes%salt_left_behind)) then call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//& "Brine plume parameterization requires DO_BRINE_PLUME\n"//& @@ -1232,7 +1235,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, & !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho,& !$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, & - !$OMP EnthalpyConst,MLD) & + !$OMP EnthalpyConst,MLD_h) & !$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, & !$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, & !$OMP IforcingDepthScale,g_conv,dSpV_dT,dSpV_dS, & @@ -1242,7 +1245,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !$OMP drhodt,drhods,pen_sw_bnd_rate, & !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, & !$OMP mixing_depth,A_brine,fraction_left_brine, & - !$OMP plume_fraction,dK,MLD_H,MLD_Z,total_h) & + !$OMP plume_fraction,dK,total_h) & !$OMP firstprivate(SurfPressure,plume_flux) do j=js,je ! Work in vertical slices for efficiency @@ -1368,24 +1371,10 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! C/ update temp due to penetrative SW if (CS%do_brine_plume) then ! Find the plume mixing depth. - if (GV%Boussinesq .or. .not.allocated(tv%SpV_avg)) then - do i=is,ie ; MLD_H(i) = GV%Z_to_H * MLD(i,j) ; total_h(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie ; total_h(i) = total_h(i) + h(i,j,k) ; enddo ; enddo - else - do i=is,ie ; MLD_H(i) = 0.0 ; MLD_Z(i) = 0.0 ; total_h(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie - total_h(i) = total_h(i) + h(i,j,k) - if (MLD_Z(i) + GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) < MLD(i,j)) then - MLD_H(i) = MLD_H(i) + h(i,j,k) - MLD_Z(i) = MLD_Z(i) + GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) - elseif (MLD_Z(i) < MLD(i,j)) then ! This is the last layer in the mixed layer - MLD_H(i) = MLD_H(i) + GV%RZ_to_H * (MLD(i,j) - MLD_Z(i)) / tv%SpV_avg(i,j,k) - MLD_Z(i) = MLD(i,j) - endif - enddo ; enddo - endif + do i=is,ie ; total_h(i) = 0.0 ; enddo + do k=1,nz ; do i=is,ie ; total_h(i) = total_h(i) + h(i,j,k) ; enddo ; enddo do i=is,ie - mixing_depth(i) = min( max(MLD_H(i) - minimum_forcing_depth, minimum_forcing_depth), & + mixing_depth(i) = min( max(MLD_h(i,j) - minimum_forcing_depth, minimum_forcing_depth), & max(total_h(i), GV%angstrom_h) ) + GV%H_subroundoff A_brine(i) = (CS%brine_plume_n + 1) / (mixing_depth(i) ** (CS%brine_plume_n + 1)) enddo diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 4c96c0a9c4..f967c005e6 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -830,7 +830,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim skinbuoyflux(:,:) = 0.0 call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=visc%MLD) + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD_h=visc%h_ML) if (CS%debug) then call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_mks) @@ -890,7 +890,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim else call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD=visc%MLD) + CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD_h=visc%h_ML) ! Find the vertical distances across layers, which may have been modified by the net surface flux call thickness_to_dz(h, tv, dz, G, GV, US) @@ -1378,7 +1378,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, skinbuoyflux(:,:) = 0.0 call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=visc%MLD) + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD_h=visc%h_ML) if (CS%debug) then call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_MKS) @@ -1424,7 +1424,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, else call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD=visc%MLD) + CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD_h=visc%h_ML) endif ! endif for CS%use_energetic_PBL diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index dba32b24d3..b2821a4ea8 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -2680,6 +2680,7 @@ subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_C ! Local variables logical :: use_kappa_shear, KS_at_vertex logical :: adiabatic, useKPP, useEPBL + logical :: do_brine_plume, use_hor_bnd_diff, use_neutral_diffusion, use_fpmix logical :: use_CVMix_shear, MLE_use_PBL_MLD, MLE_use_Bodner, use_CVMix_conv integer :: isd, ied, jsd, jed, nz real :: hfreeze !< If hfreeze > 0 [Z ~> m], melt potential will be computed. @@ -2743,25 +2744,43 @@ subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_C call safe_alloc_ptr(visc%Kv_slow, isd, ied, jsd, jed, nz+1) endif - ! visc%MLD is used to communicate the state of the (e)PBL or KPP to the rest of the model + ! visc%MLD and visc%h_ML are used to communicate the state of the (e)PBL or KPP to the rest of the model call get_param(param_file, mdl, "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, & default=.false., do_not_log=.true.) - ! visc%MLD needs to be allocated when melt potential is computed (HFREEZE>0) + ! visc%h_ML needs to be allocated when melt potential is computed (HFREEZE>0) or one of + ! several other parameterizations are in use. call get_param(param_file, mdl, "HFREEZE", hfreeze, & units="m", default=-1.0, scale=US%m_to_Z, do_not_log=.true.) + call get_param(param_file, mdl, "DO_BRINE_PLUME", do_brine_plume, & + "If true, use a brine plume parameterization from Nguyen et al., 2009.", & + default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "USE_HORIZONTAL_BOUNDARY_DIFFUSION", use_hor_bnd_diff, & + default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", use_neutral_diffusion, & + default=.false., do_not_log=.true.) + if (use_neutral_diffusion) & + call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", use_neutral_diffusion, & + default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "FPMIX", use_fpmix, & + default=.false., do_not_log=.true.) - if (hfreeze >= 0.0 .or. MLE_use_PBL_MLD) then + if (MLE_use_PBL_MLD) then call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed) endif - if (hfreeze >= 0.0 .or. MLE_use_PBL_MLD) then + if ((hfreeze >= 0.0) .or. MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. & + use_neutral_diffusion .or. use_hor_bnd_diff) then call safe_alloc_ptr(visc%h_ML, isd, ied, jsd, jed) endif if (MLE_use_PBL_MLD) then call register_restart_field(visc%MLD, "MLD", .false., restart_CS, & "Instantaneous active mixing layer depth", units="m", conversion=US%Z_to_m) + endif + if (MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. & + use_neutral_diffusion .or. use_hor_bnd_diff) then call register_restart_field(visc%h_ML, "h_ML", .false., restart_CS, & - "Instantaneous active mixing layer thickness", units=get_thickness_units(GV), conversion=GV%H_to_mks) + "Instantaneous active mixing layer thickness", & + units=get_thickness_units(GV), conversion=GV%H_to_mks) endif ! visc%sfc_buoy_flx is used to communicate the state of the (e)PBL or KPP to the rest of the model diff --git a/src/tracer/MOM_hor_bnd_diffusion.F90 b/src/tracer/MOM_hor_bnd_diffusion.F90 index 163d8a480f..13e91e8973 100644 --- a/src/tracer/MOM_hor_bnd_diffusion.F90 +++ b/src/tracer/MOM_hor_bnd_diffusion.F90 @@ -20,6 +20,7 @@ module MOM_hor_bnd_diffusion use MOM_spatial_means, only : global_mass_integral use MOM_tracer_registry, only : tracer_registry_type, tracer_type use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : vertvisc_type use MOM_verticalGrid, only : verticalGrid_type use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS @@ -161,7 +162,7 @@ end function hor_bnd_diffusion_init !! 2) calculate diffusive tracer fluxes (F) in the HBD grid using a layer by layer approach !! 3) remap fluxes to the native grid !! 4) update tracer by adding the divergence of F -subroutine hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) +subroutine hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, visc, CS) type(ocean_grid_type), intent(inout) :: G !< Grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -171,6 +172,8 @@ subroutine hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) real, intent(in) :: dt !< Tracer time step * I_numitts !! (I_numitts in tracer_hordiff) [T ~> s] type(tracer_registry_type), pointer :: Reg !< Tracer registry + type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities, + !! boundary layer properties and related fields type(hbd_CS), pointer :: CS !< Control structure for this module ! Local variables @@ -203,10 +206,14 @@ subroutine hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) call cpu_clock_begin(id_clock_hbd) Idt = 1./dt - if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) - if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, & - m_to_MLD_units=GV%m_to_H) - call pass_var(hbl,G%Domain) + + if (associated(visc%h_ML)) then + hbl(:,:) = visc%h_ML(:,:) + else + call MOM_error(FATAL, "hor_bnd_diffusion requires that visc%h_ML is associated.") + endif + ! This halo update is probably not necessary because visc%h_ML has valid halo data. + call pass_var(hbl, G%Domain, halo=1) ! build HBD grid call hbd_grid(SURFACE, G, GV, hbl, h, CS) @@ -336,10 +343,11 @@ subroutine hbd_grid(boundary, G, GV, hbl, h, CS) integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] type(ocean_grid_type), intent(inout) :: G !< Grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G)) :: hbl !< Boundary layer depth [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: hbl !< Boundary layer depth [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Layer thickness in the native grid [H ~> m or kg m-2] - type(hbd_CS), pointer :: CS !< Horizontal diffusion control structure + intent(in) :: h !< Layer thickness in the native grid [H ~> m or kg m-2] + type(hbd_CS), pointer :: CS !< Horizontal diffusion control structure ! Local variables real, allocatable :: dz_top(:) !< temporary HBD grid given by merge_interfaces [H ~> m or kg m-2] diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 094825d031..402a008244 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -20,6 +20,7 @@ module MOM_neutral_diffusion use MOM_remapping, only : average_value_ppoly, remappingSchemesDoc, remappingDefaultScheme use MOM_tracer_registry, only : tracer_registry_type, tracer_type use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : vertvisc_type use MOM_verticalGrid, only : verticalGrid_type use polynomial_functions, only : evaluation_polynomial, first_derivative_polynomial use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation @@ -142,7 +143,7 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), target, intent(in) :: EOS !< Equation of state - type(diabatic_CS), pointer :: diabatic_CSp!< KPP control structure needed to get BLD + type(diabatic_CS), pointer :: diabatic_CSp!< diabatic control structure needed to get BLD type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure ! Local variables @@ -333,13 +334,15 @@ end function neutral_diffusion_init !> Calculate remapping factors for u/v columns used to map adjoining columns to !! a shared coordinate space. -subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) +subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, visc, CS, p_surf) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: T !< Potential temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: S !< Salinity [S ~> ppt] + type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities, + !! boundary layer properties and related fields type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: p_surf !< Surface pressure to include in pressures used !! for equation of state calculations [R L2 T-2 ~> Pa] @@ -369,10 +372,13 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) ! Check if hbl needs to be extracted if (CS%interior_only) then - if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, CS%hbl, G, US, m_to_BLD_units=GV%m_to_H) - if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, CS%hbl, G, US, & - m_to_MLD_units=GV%m_to_H) - call pass_var(CS%hbl,G%Domain) + if (associated(visc%h_ML)) then + CS%hbl(:,:) = visc%h_ML(:,:) + else + call MOM_error(FATAL, "hor_bnd_diffusion requires that visc%h_ML is associated.") + endif + call pass_var(CS%hbl, G%Domain, halo=1) + ! get k-indices and zeta do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1 if (G%mask2dT(i,j) > 0.0) then diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 732a42e44b..c3ce48b06a 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -27,7 +27,7 @@ module MOM_tracer_hor_diff use MOM_hor_bnd_diffusion, only : hor_bnd_diffusion, hor_bnd_diffusion_end use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs +use MOM_variables, only : thermo_var_ptrs, vertvisc_type use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -113,7 +113,7 @@ module MOM_tracer_hor_diff !! using the diffusivity in CS%KhTr, or using space-dependent diffusivity. !! Multiple iterations are used (if necessary) so that there is no limit !! on the acceptable time increment. -subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y) +subroutine tracer_hordiff(h, dt, MEKE, VarMix, visc, G, GV, US, CS, Reg, tv, do_online_flag, read_khdt_x, read_khdt_y) type(ocean_grid_type), intent(inout) :: G !< Grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -121,6 +121,8 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online real, intent(in) :: dt !< time step [T ~> s] type(MEKE_type), intent(in) :: MEKE !< MEKE fields type(VarMix_CS), intent(in) :: VarMix !< Variable mixing type + type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities, + !! boundary layer properties and related fields type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(tracer_hor_diff_CS), pointer :: CS !< module control structure type(tracer_registry_type), pointer :: Reg !< registered tracers @@ -442,7 +444,7 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online if (itt>1) then ! Update halos for subsequent iterations call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) endif - call hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, I_numitts*dt, Reg, & + call hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, I_numitts*dt, Reg, visc, & CS%hor_bnd_diffusion_CSp) enddo ! itt endif @@ -457,9 +459,9 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online ! would be inside the itt-loop. -AJA if (associated(tv%p_surf)) then - call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, CS%neutral_diffusion_CSp, p_surf=tv%p_surf) + call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp, p_surf=tv%p_surf) else - call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, CS%neutral_diffusion_CSp) + call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp) endif do k=1,nz+1 @@ -499,9 +501,9 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online call do_group_pass(CS%pass_t, G%Domain, clock=id_clock_pass) if (CS%recalc_neutral_surf) then if (associated(tv%p_surf)) then - call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, CS%neutral_diffusion_CSp, p_surf=tv%p_surf) + call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp, p_surf=tv%p_surf) else - call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, CS%neutral_diffusion_CSp) + call neutral_diffusion_calc_coeffs(G, GV, US, h, tv%T, tv%S, visc, CS%neutral_diffusion_CSp) endif endif endif