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