diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 698243a7f6..963547f3c0 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -251,7 +251,7 @@ module MOM_diabatic_driver integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap -integer :: id_clock_kpp, id_clock_CVMix_ddiff +integer :: id_clock_kpp contains @@ -279,6 +279,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & type(diabatic_CS), pointer :: CS !< module control structure type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves + ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & ea, & ! amount of fluid entrained from the layer above within ! one time step (m for Bouss, kg/m^2 for non-Bouss) @@ -391,6 +392,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("diabatic(), MOM_diabatic_driver.F90") + if (.not. (CS%useALEalgorithm)) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "The ALE algorithm must be enabled when using MOM_diabatic_driver.") ! Offer diagnostics of various state varables at the start of diabatic ! these are mostly for debugging purposes. @@ -452,7 +455,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%id_frazil_h > 0) call post_data(CS%id_frazil_h, h, CS%diag) endif call disable_averaging(CS%diag) - endif + endif !associated(tv%T) .AND. associated(tv%frazil) + ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep call enable_averaging(dt, Time_end, CS%diag) if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, G) @@ -482,58 +486,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (associated(CS%optics)) & call set_opacity(CS%optics, fluxes, G, GV, CS%opacity_CSp) - if (CS%bulkmixedlayer) then - if (CS%debug) then - call MOM_forcing_chksum("Before mixedlayer", fluxes, G, haloshift=0) - endif - - if (CS%ML_mix_first > 0.0) then - ! This subroutine: - ! (1) Cools the mixed layer. - ! (2) Performs convective adjustment by mixed layer entrainment. - ! (3) Heats the mixed layer and causes it to detrain to - ! Monin-Obukhov depth or minimum mixed layer depth. - ! (4) Uses any remaining TKE to drive mixed layer entrainment. - ! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) - call find_uv_at_h(u, v, h, u_h, v_h, G, GV) - - call cpu_clock_begin(id_clock_mixedlayer) - if (CS%ML_mix_first < 1.0) then - ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*CS%ML_mix_first, & - eaml,ebml, G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.false.) - if (CS%salt_reject_below_ML) & - call insert_brine(h, tv, G, GV, fluxes, nkmb, CS%diabatic_aux_CSp, & - dt*CS%ML_mix_first, CS%id_brine_lay) - else - ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, & - G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.true.) - endif - - ! Keep salinity from falling below a small but positive threshold. - ! This constraint is needed for SIS1 ice model, which can extract - ! more salt than is present in the ocean. SIS2 does not suffer - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - call cpu_clock_end(id_clock_mixedlayer) - if (CS%debug) then - call MOM_state_chksum("After mixedlayer ", u, v, h, G, GV, haloshift=0) - call MOM_forcing_chksum("After mixedlayer", fluxes, G, haloshift=0) - endif - if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)") - if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G) - endif - endif ! end CS%bulkmixedlayer - - if (CS%debug) then + if (CS%debug) & call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, haloshift=0) - endif if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then @@ -609,7 +563,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call hchksum(Kd_Int, "after set_diffusivity Kd_Int",G%HI,haloshift=0) endif - if (CS%useKPP) then call cpu_clock_begin(id_clock_kpp) ! KPP needs the surface buoyancy flux but does not update state variables. @@ -731,11 +684,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! If using matching within the KPP scheme, then this step needs to provide ! a diffusivity and happen before KPP. But generally in MOM, we do not match ! KPP boundary layer to interior, so this diffusivity can be computed when convenient. - if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T)) then - call cpu_clock_begin(id_clock_CVMix_ddiff) + if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T) .and. .not. & + CS%use_CVMix_ddiff) then + call cpu_clock_begin(id_clock_differential_diff) call differential_diffuse_T_S(h, tv, visc, dt, G, GV) - call cpu_clock_end(id_clock_CVMix_ddiff) + call cpu_clock_end(id_clock_differential_diff) + if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G) @@ -752,47 +707,22 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! This block sets ea, eb from Kd or Kd_int. - ! If using ALE algorithm, set ea=eb=Kd_int on interfaces for - ! use in the tri-diagonal solver. - ! Otherwise, call entrainment_diffusive() which sets ea and eb - ! based on KD and target densities (ie. does remapping as well). - if (CS%useALEalgorithm) then + ! set ea=eb=Kd_int on interfaces for use in the tri-diagonal solver. - do j=js,je ; do i=is,ie - ea(i,j,1) = 0. - enddo ; enddo + do j=js,je ; do i=is,ie + ea(i,j,1) = 0. + enddo ; enddo !$OMP parallel do default(none) shared(is,ie,js,je,nz,h_neglect,h,ea,GV,dt,Kd_int,eb) & !$OMP private(hval) - do k=2,nz ; do j=js,je ; do i=is,ie - hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) - ea(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_int(i,j,k) - eb(i,j,k-1) = ea(i,j,k) - enddo ; enddo ; enddo - do j=js,je ; do i=is,ie - eb(i,j,nz) = 0. - enddo ; enddo - if (showCallTree) call callTree_waypoint("done setting ea,eb from Kd_int (diabatic)") - - else ! .not. CS%useALEalgorithm - ! When not using ALE, calculate layer entrainments/detrainments from - ! diffusivities and differences between layer and target densities - call cpu_clock_begin(id_clock_entrain) - ! Calculate appropriately limited diapycnal mass fluxes to account - ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb - call Entrainment_diffusive(u, v, h, tv, fluxes, dt, G, GV, CS%entrain_diffusive_CSp, & - ea, eb, kb, Kd_Lay=Kd, Kd_int=Kd_int) - call cpu_clock_end(id_clock_entrain) - if (showCallTree) call callTree_waypoint("done with Entrainment_diffusive (diabatic)") - - endif ! endif for (CS%useALEalgorithm) - - if (CS%debug) then - call MOM_forcing_chksum("after calc_entrain ", fluxes, G, haloshift=0) - call MOM_thermovar_chksum("after calc_entrain ", tv, G) - call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, haloshift=0) - call hchksum(ea, "after calc_entrain ea", G%HI, haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after calc_entrain eb", G%HI, haloshift=0, scale=GV%H_to_m) - endif + do k=2,nz ; do j=js,je ; do i=is,ie + hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) + ea(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_int(i,j,k) + eb(i,j,k-1) = ea(i,j,k) + enddo ; enddo ; enddo + do j=js,je ; do i=is,ie + eb(i,j,nz) = 0. + enddo ; enddo + if (showCallTree) call callTree_waypoint("done setting ea,eb from Kd_int (diabatic)") ! Save fields before boundary forcing is applied for tendency diagnostics if (CS%boundary_forcing_tendency_diag) then @@ -803,96 +733,90 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif - ! Apply forcing when using the ALE algorithm - if (CS%useALEalgorithm) then - call cpu_clock_begin(id_clock_remap) - - ! Changes made to following fields: h, tv%T and tv%S. - - do k=1,nz ; do j=js,je ; do i=is,ie - h_prebound(i,j,k) = h(i,j,k) - enddo ; enddo ; enddo - if (CS%use_energetic_PBL) then - - skinbuoyflux(:,:) = 0.0 - call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & - h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) - - if (CS%debug) then - call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after applyBoundaryFluxes eb",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(cTKE, "after applyBoundaryFluxes cTKE",G%HI,haloshift=0) - call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0) - call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0) - endif + ! Apply forcing + call cpu_clock_begin(id_clock_remap) - call find_uv_at_h(u, v, h, u_h, v_h, G, GV) - call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, & - CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) - - ! If visc%MLD exists, copy the ePBL's MLD into it - if (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G) - call pass_var(visc%MLD, G%domain, halo=1) - Hml(:,:) = visc%MLD(:,:) - endif + ! Changes made to following fields: h, tv%T and tv%S. + do k=1,nz ; do j=js,je ; do i=is,ie + h_prebound(i,j,k) = h(i,j,k) + enddo ; enddo ; enddo + if (CS%use_energetic_PBL) then - ! Augment the diffusivities due to those diagnosed in energetic_PBL. - do K=2,nz ; do j=js,je ; do i=is,ie + skinbuoyflux(:,:) = 0.0 + call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & + h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) - if (CS%ePBL_is_additive) then - Kd_add_here = Kd_ePBL(i,j,K) - visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + Kd_ePBL(i,j,K) - else - Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_shear(i,j,K), 0.0) - visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), Kd_ePBL(i,j,K)) - endif - Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / & - (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) - eb(i,j,k-1) = eb(i,j,k-1) + Ent_int - ea(i,j,k) = ea(i,j,k) + Ent_int - Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here + if (CS%debug) then + call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "after applyBoundaryFluxes eb",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(cTKE, "after applyBoundaryFluxes cTKE",G%HI,haloshift=0) + call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0) + call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0) + endif - ! for diagnostics - Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K) - Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_int(i,j,K) + call find_uv_at_h(u, v, h, u_h, v_h, G, GV) + call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, & + CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) - enddo ; enddo ; enddo + ! If visc%MLD exists, copy the ePBL's MLD into it + if (associated(visc%MLD)) then + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G) + call pass_var(visc%MLD, G%domain, halo=1) + Hml(:,:) = visc%MLD(:,:) + endif - if (CS%debug) then - call hchksum(ea, "after ePBL ea",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after ePBL eb",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(Kd_ePBL, "after ePBL Kd_ePBL",G%HI,haloshift=0) + ! Augment the diffusivities due to those diagnosed in energetic_PBL. + do K=2,nz ; do j=js,je ; do i=is,ie + if (CS%ePBL_is_additive) then + Kd_add_here = Kd_ePBL(i,j,K) + visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + Kd_ePBL(i,j,K) + else + Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_shear(i,j,K), 0.0) + visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), Kd_ePBL(i,j,K)) endif + Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / & + (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) + eb(i,j,k-1) = eb(i,j,k-1) + Ent_int + ea(i,j,k) = ea(i,j,k) + Ent_int + Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here - else - call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & - h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth) - - endif ! endif for CS%use_energetic_PBL - - ! diagnose the tendencies due to boundary forcing - ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme - ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards - if (CS%boundary_forcing_tendency_diag) then - call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, G, GV, CS) - if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h = h_diag) - endif - ! Boundary fluxes may have changed T, S, and h - call diag_update_remap_grids(CS%diag) + ! for diagnostics + Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K) + Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_int(i,j,K) + + enddo ; enddo ; enddo - call cpu_clock_end(id_clock_remap) if (CS%debug) then - call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, haloshift=0) - call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G) - call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, haloshift=0) + call hchksum(ea, "after ePBL ea",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "after ePBL eb",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(Kd_ePBL, "after ePBL Kd_ePBL",G%HI,haloshift=0) endif - if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)") - if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G) - endif ! endif for (CS%useALEalgorithm) + else + call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & + h, tv, CS%aggregate_FW_forcing, & + CS%evap_CFL_limit, CS%minimum_forcing_depth) + + endif ! endif for CS%use_energetic_PBL + + ! diagnose the tendencies due to boundary forcing + ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme + ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards + if (CS%boundary_forcing_tendency_diag) then + call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, G, GV, CS) + if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h = h_diag) + endif + ! Boundary fluxes may have changed T, S, and h + call diag_update_remap_grids(CS%diag) + call cpu_clock_end(id_clock_remap) + if (CS%debug) then + call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G) + call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, haloshift=0) + endif + if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)") + if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G) ! Update h according to divergence of the difference between ! ea and eb. We keep a record of the original h in hold. @@ -902,6 +826,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! enough iterations are permitted in Calculate_Entrainment. ! Even if too few iterations are allowed, it is still guarded ! against. In other words the checks are probably unnecessary. + + ! GMM, should the code below be deleted? eb(i,j,k-1) = ea(i,j,k), + ! see above, so h should not change. + !$OMP parallel do default(shared) do j=js,je do i=is,ie @@ -937,204 +865,55 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G) - ! Here, T and S are updated according to ea and eb. - ! If using the bulk mixed layer, T and S are also updated - ! by surface fluxes (in fluxes%*). - ! This is a very long block. - if (CS%bulkmixedlayer) then + ! calculate change in temperature & salinity due to dia-coordinate surface diffusion + if (associated(tv%T)) then - if (associated(tv%T)) then - call cpu_clock_begin(id_clock_tridiag) - ! Temperature and salinity (as state variables) are treated - ! differently from other tracers to insure massless layers that - ! are lighter than the mixed layer have temperatures and salinities - ! that correspond to their prescribed densities. - if (CS%massless_match_targets) then - !$OMP parallel do default (shared) private(h_tr,b1,d1,c1,b_denom_1) - do j=js,je - do i=is,ie - h_tr = hold(i,j,1) + h_neglect - b1(i) = 1.0 / (h_tr + eb(i,j,1)) - d1(i) = h_tr * b1(i) - tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1)) - tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1)) - enddo - do k=2,nkmb ; do i=is,ie - c1(i,k) = eb(i,j,k-1) * b1(i) - h_tr = hold(i,j,k) + h_neglect - b_denom_1 = h_tr + d1(i)*ea(i,j,k) - b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) - if (k kb(i,j)) then - c1(i,k) = eb(i,j,k-1) * b1(i) - h_tr = hold(i,j,k) + h_neglect - b_denom_1 = h_tr + d1(i)*ea(i,j,k) - b1(i) = 1.0 / (b_denom_1 + eb(i,j,k)) - d1(i) = b_denom_1 * b1(i) - tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1)) - tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1)) - elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j)) - ! The bottommost buffer layer might entrain all the mass from some - ! of the interior layers that are thin and lighter in the coordinate - ! density than that buffer layer. The T and S of these newly - ! massless interior layers are unchanged. - tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k) - tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k) - endif - enddo ; enddo - - do k=nz-1,nkmb,-1 ; do i=is,ie - if (k >= kb(i,j)) then - tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) - tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) - endif - enddo ; enddo - do i=is,ie ; if (kb(i,j) <= nz) then - tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j)) - tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j)) - endif ; enddo - do k=nkmb-1,1,-1 ; do i=is,ie - tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) - tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) - enddo ; enddo - enddo ! end of j loop - else ! .not. massless_match_targets - ! This simpler form allows T & S to be too dense for the layers - ! between the buffer layers and the interior. - ! Changes: T, S - if (CS%tracer_tridiag) then - call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV) - call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV) - else - call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) - endif - endif ! massless_match_targets - call cpu_clock_end(id_clock_tridiag) + if (CS%debug) then + call hchksum(ea, "before triDiagTS ea ",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb, "before triDiagTS eb ",G%HI,haloshift=0, scale=GV%H_to_m) + endif + call cpu_clock_begin(id_clock_tridiag) - endif ! endif for associated(T) - if (CS%debugConservation) call MOM_state_stats('BML tridiag', u, v, h, tv%T, tv%S, G) + ! Keep salinity from falling below a small but positive threshold. + ! This constraint is needed for SIS1 ice model, which can extract + ! more salt than is present in the ocean. SIS2 does not suffer + ! from this limitation, in which case we can let salinity=0 and still + ! have salt conserved with SIS2 ice. So for SIS2, we can run with + ! BOUND_SALINITY=False in MOM.F90. + if (associated(tv%S) .and. associated(tv%salt_deficit)) & + call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then - ! The mixed layer code has already been called, but there is some needed - ! bookkeeping. - !$OMP parallel do default(shared) + if (CS%diabatic_diff_tendency_diag) then do k=1,nz ; do j=js,je ; do i=is,ie - hold(i,j,k) = h_orig(i,j,k) - ea(i,j,k) = ea(i,j,k) + eaml(i,j,k) - eb(i,j,k) = eb(i,j,k) + ebml(i,j,k) + temp_diag(i,j,k) = tv%T(i,j,k) + saln_diag(i,j,k) = tv%S(i,j,k) enddo ; enddo ; enddo - if (CS%debug) then - call hchksum(ea, "after ea = ea + eaml",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after eb = eb + ebml",G%HI,haloshift=0, scale=GV%H_to_m) - endif endif - if (CS%ML_mix_first < 1.0) then - ! Call the mixed layer code now, perhaps for a second time. - ! This subroutine (1) Cools the mixed layer. - ! (2) Performs convective adjustment by mixed layer entrainment. - ! (3) Heats the mixed layer and causes it to detrain to - ! Monin-Obukhov depth or minimum mixed layer depth. - ! (4) Uses any remaining TKE to drive mixed layer entrainment. - ! (5) Possibly splits the buffer layer into two isopycnal layers. - - call find_uv_at_h(u, v, hold, u_h, v_h, G, GV, ea, eb) - if (CS%debug) call MOM_state_chksum("find_uv_at_h1 ", u, v, h, G, GV, haloshift=0) - - dt_mix = min(dt,dt*(1.0 - CS%ML_mix_first)) - call cpu_clock_begin(id_clock_mixedlayer) - ! Changes: h, tv%T, tv%S, ea and eb (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, & - G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.true.) - - if (CS%salt_reject_below_ML) & - call insert_brine(h, tv, G, GV, fluxes, nkmb, CS%diabatic_aux_CSp, dt_mix, & - CS%id_brine_lay) - - ! Keep salinity from falling below a small but positive threshold. - ! This constraint is needed for SIS1 ice model, which can extract - ! more salt than is present in the ocean. SIS2 does not suffer - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - - call cpu_clock_end(id_clock_mixedlayer) - if (showCallTree) call callTree_waypoint("done with 2nd bulkmixedlayer (diabatic)") - if (CS%debugConservation) call MOM_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, G) - endif - - else ! following block for when NOT using BULKMIXEDLAYER - - - ! calculate change in temperature & salinity due to dia-coordinate surface diffusion - if (associated(tv%T)) then - - if (CS%debug) then - call hchksum(ea, "before triDiagTS ea ",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "before triDiagTS eb ",G%HI,haloshift=0, scale=GV%H_to_m) - endif - call cpu_clock_begin(id_clock_tridiag) - - ! Keep salinity from falling below a small but positive threshold. - ! This constraint is needed for SIS1 ice model, which can extract - ! more salt than is present in the ocean. SIS2 does not suffer - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - - if (CS%diabatic_diff_tendency_diag) then - do k=1,nz ; do j=js,je ; do i=is,ie - temp_diag(i,j,k) = tv%T(i,j,k) - saln_diag(i,j,k) = tv%S(i,j,k) - enddo ; enddo ; enddo - endif - - ! Changes T and S via the tridiagonal solver; no change to h - if (CS%tracer_tridiag) then - call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV) - call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV) - else - call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) - endif - - ! diagnose temperature, salinity, heat, and salt tendencies - ! Note: hold here refers to the thicknesses from before the dual-entraintment when using - ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will have changed - ! In either case, tendencies should be posted on hold - if (CS%diabatic_diff_tendency_diag) then - call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, G, GV, CS) - if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, hold, CS%diag, alt_h = hold) - endif + ! Changes T and S via the tridiagonal solver; no change to h + if (CS%tracer_tridiag) then + call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV) + call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV) + else - call cpu_clock_end(id_clock_tridiag) - if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)") + call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) + endif - endif ! endif corresponding to if (associated(tv%T)) - if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G) + ! diagnose temperature, salinity, heat, and salt tendencies + ! Note: hold here refers to the thicknesses from before the dual-entraintment when using + ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will have changed + ! In either case, tendencies should be posted on hold + if (CS%diabatic_diff_tendency_diag) then + call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, G, GV, CS) + if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, hold, CS%diag, alt_h = hold) + endif + call cpu_clock_end(id_clock_tridiag) + if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)") - endif ! endif for the BULKMIXEDLAYER block + endif ! endif corresponding to if (associated(tv%T)) + if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G) if (CS%debug) then call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, haloshift=0) @@ -1143,14 +922,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call hchksum(eb, "after mixed layer eb", G%HI, scale=GV%H_to_m) endif - if (.not. CS%useALEalgorithm) then - call cpu_clock_begin(id_clock_remap) - call regularize_layers(h, tv, dt, ea, eb, G, GV, CS%regularize_layers_CSp) - call cpu_clock_end(id_clock_remap) - if (showCallTree) call callTree_waypoint("done with regularize_layers (diabatic)") - if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G) - endif - ! Whenever thickness changes let the diag manager know, as the ! target grids for vertical remapping may need to be regenerated. call diag_update_remap_grids(CS%diag) @@ -1187,6 +958,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! mixing of passive tracers from massless boundary layers to interior call cpu_clock_begin(id_clock_tracers) + if (CS%mix_boundary_tracers) then Tr_ea_BBL = sqrt(dt*CS%Kd_BBL_tr) !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) @@ -1235,17 +1007,12 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo - if (CS%useALEalgorithm) then ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied ! so hold should be h_orig - call call_tracer_column_fns(h_prebound, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug, & - evap_CFL_limit = CS%evap_CFL_limit, & - minimum_forcing_depth = CS%minimum_forcing_depth) - else - call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug, & + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) elseif (associated(visc%Kd_extra_S)) then ! extra diffusivity for passive tracers @@ -1265,56 +1032,31 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & eatr(i,j,k) = ea(i,j,k) + add_ent enddo ; enddo ; enddo - if (CS%useALEalgorithm) then ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied - call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug,& - evap_CFL_limit = CS%evap_CFL_limit, & - minimum_forcing_depth = CS%minimum_forcing_depth) - else - call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug,& + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) else - if (CS%useALEalgorithm) then ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied - call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug, & - evap_CFL_limit = CS%evap_CFL_limit, & - minimum_forcing_depth = CS%minimum_forcing_depth) - else - call call_tracer_column_fns(hold, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & + CS%optics, CS%tracer_flow_CSp, CS%debug, & + evap_CFL_limit = CS%evap_CFL_limit, & + minimum_forcing_depth = CS%minimum_forcing_depth) endif ! (CS%mix_boundary_tracers) - - call cpu_clock_end(id_clock_tracers) - ! sponges if (CS%use_sponge) then call cpu_clock_begin(id_clock_sponge) if (associated(CS%ALE_sponge_CSp)) then ! ALE sponge call apply_ALE_sponge(h, dt, G, CS%ALE_sponge_CSp, CS%Time) - else - ! Layer mode sponge - if (CS%bulkmixedlayer .and. associated(tv%eqn_of_state)) then - do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo - !$OMP parallel do default(shared) - do j=js,je - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), & - is, ie-is+1, tv%eqn_of_state) - enddo - call apply_sponge(h, dt, G, GV, ea, eb, CS%sponge_CSp, Rcv_ml) - else - call apply_sponge(h, dt, G, GV, ea, eb, CS%sponge_CSp) - endif endif + call cpu_clock_end(id_clock_sponge) if (CS%debug) then call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, haloshift=0) @@ -1337,29 +1079,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo endif -! For momentum, it is only the net flux that homogenizes within -! the mixed layer. Vertical viscosity that is proportional to the -! mixed layer turbulence is applied elsewhere. - if (CS%bulkmixedlayer) then - if (CS%debug) then - call hchksum(ea, "before net flux rearrangement ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "before net flux rearrangement eb",G%HI, scale=GV%H_to_m) - endif - !$OMP parallel do default(shared) private(net_ent) - do j=js,je - do K=2,GV%nkml ; do i=is,ie - net_ent = ea(i,j,k) - eb(i,j,k-1) - ea(i,j,k) = max(net_ent, 0.0) - eb(i,j,k-1) = max(-net_ent, 0.0) - enddo ; enddo - enddo - if (CS%debug) then - call hchksum(ea, "after net flux rearrangement ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "after net flux rearrangement eb",G%HI, scale=GV%H_to_m) - endif - endif - -! Initialize halo regions of ea, eb, and hold to default values. + ! Initialize halo regions of ea, eb, and hold to default values. !$OMP parallel do default(shared) do k=1,nz do i=is-1,ie+1 @@ -1387,84 +1107,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call cpu_clock_end(id_clock_pass) - if (.not. CS%useALEalgorithm) then - ! Use a tridiagonal solver to determine effect of the diapycnal - ! advection on velocity field. It is assumed that water leaves - ! or enters the ocean with the surface velocity. - if (CS%debug) then - call MOM_state_chksum("before u/v tridiag ", u, v, h, G, GV, haloshift=0) - call hchksum(ea, "before u/v tridiag ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "before u/v tridiag eb",G%HI, scale=GV%H_to_m) - call hchksum(hold, "before u/v tridiag hold",G%HI, scale=GV%H_to_m) - endif - call cpu_clock_begin(id_clock_tridiag) - !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval) - do j=js,je - do I=Isq,Ieq - if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,1) = u(I,j,1) - hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect - b1(I) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1))) - d1(I) = hval * b1(I) - u(I,j,1) = b1(I) * (hval * u(I,j,1)) - enddo - do k=2,nz ; do I=Isq,Ieq - if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,k) = u(I,j,k) - c1(I,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(I) - eaval = ea(i,j,k) + ea(i+1,j,k) - hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect - b1(I) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(I)*eaval)) - d1(I) = (hval + d1(I)*eaval) * b1(I) - u(I,j,k) = (hval*u(I,j,k) + eaval*u(I,j,k-1))*b1(I) - enddo ; enddo - do k=nz-1,1,-1 ; do I=Isq,Ieq - u(I,j,k) = u(I,j,k) + c1(I,k+1)*u(I,j,k+1) - if (associated(ADp%du_dt_dia)) & - ADp%du_dt_dia(I,j,k) = (u(I,j,k) - ADp%du_dt_dia(I,j,k)) * Idt - enddo ; enddo - if (associated(ADp%du_dt_dia)) then - do I=Isq,Ieq - ADp%du_dt_dia(I,j,nz) = (u(I,j,nz)-ADp%du_dt_dia(I,j,nz)) * Idt - enddo - endif - enddo - if (CS%debug) then - call MOM_state_chksum("aft 1st loop tridiag ", u, v, h, G, GV, haloshift=0) - endif - !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval) - do J=Jsq,Jeq - do i=is,ie - if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,1) = v(i,J,1) - hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect - b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1))) - d1(I) = hval * b1(I) - v(i,J,1) = b1(i) * (hval * v(i,J,1)) - enddo - do k=2,nz ; do i=is,ie - if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,k) = v(i,J,k) - c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i) - eaval = ea(i,j,k) + ea(i,j+1,k) - hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect - b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval)) - d1(i) = (hval + d1(i)*eaval) * b1(i) - v(i,J,k) = (hval*v(i,J,k) + eaval*v(i,J,k-1))*b1(i) - enddo ; enddo - do k=nz-1,1,-1 ; do i=is,ie - v(i,J,k) = v(i,J,k) + c1(i,k+1)*v(i,J,k+1) - if (associated(ADp%dv_dt_dia)) & - ADp%dv_dt_dia(i,J,k) = (v(i,J,k) - ADp%dv_dt_dia(i,J,k)) * Idt - enddo ; enddo - if (associated(ADp%dv_dt_dia)) then - do i=is,ie - ADp%dv_dt_dia(i,J,nz) = (v(i,J,nz)-ADp%dv_dt_dia(i,J,nz)) * Idt - enddo - endif - enddo - call cpu_clock_end(id_clock_tridiag) - if (CS%debug) then - call MOM_state_chksum("after u/v tridiag ", u, v, h, G, GV, haloshift=0) - endif - endif ! useALEalgorithm - call disable_averaging(CS%diag) ! Frazil formation keeps temperature above the freezing point. ! make_frazil is deliberately called at both the beginning and at @@ -1889,7 +1531,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, real :: Kd integer :: num_mode - logical :: use_temperature + logical :: use_temperature, differentialDiffusion type(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1940,7 +1582,18 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, the diffusivity from ePBL is added to all\n"//& "other diffusivities. Otherwise, the larger of kappa-\n"//& "shear and ePBL diffusivities are used.", default=.true.) + call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, & + "If true, apply parameterization of double-diffusion.", & + default=.false. ) + CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) + + if (CS%use_CVMix_ddiff .and. differentialDiffusion) then + call MOM_error(FATAL, 'diabatic_driver_init: '// & + 'Multiple double-diffusion options selected (DOUBLE_DIFFUSION and'//& + 'USE_CVMIX_DDIFF), please disable all but one option to proceed.') + endif + CS%use_kappa_shear = kappa_shear_is_used(param_file) CS%use_CVMix_shear = CVMix_shear_is_used(param_file) @@ -2402,8 +2055,8 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=CLOCK_MODULE) id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=CLOCK_ROUTINE) id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE) - id_clock_CVMix_ddiff = -1 ; if (CS%use_CVMix_ddiff) & - id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion)', grain=CLOCK_ROUTINE) + id_clock_differential_diff = -1 ; if (differentialDiffusion) & + id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=CLOCK_ROUTINE) ! initialize the auxiliary diabatic driver module call diabatic_aux_init(Time, G, GV, param_file, diag, CS%diabatic_aux_CSp, & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 903868795a..93018c9dac 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -130,9 +130,13 @@ module MOM_set_diffusivity !! shear-driven diapycnal diffusivity. logical :: use_CVMix_shear !< If true, use one of the CVMix modules to find !! shear-driven diapycnal diffusivity. + logical :: double_diffusion !< If true, enable double-diffusive mixing using an old method. logical :: use_CVMix_ddiff !< If true, enable double-diffusive mixing via CVMix. logical :: simple_TKE_to_Kd !< If true, uses a simple estimate of Kd/TKE that !! does not rely on a layer-formulation. + real :: Max_Rrho_salt_fingers !< max density ratio for salt fingering + real :: Max_salt_diff_salt_fingers !< max salt diffusivity for salt fingers (m2/s) + real :: Kv_molecular !< molecular visc for double diff convect (m2/s) character(len=200) :: inputdir type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() @@ -156,6 +160,10 @@ module MOM_set_diffusivity integer :: id_N2 = -1 integer :: id_N2_z = -1 + integer :: id_KT_extra = -1 + integer :: id_KS_extra = -1 + integer :: id_KT_extra_z = -1 + integer :: id_KS_extra_z = -1 end type set_diffusivity_CS @@ -166,13 +174,16 @@ module MOM_set_diffusivity Kd_BBL => NULL(),& ! BBL diffusivity at interfaces (m2/s) Kd_work => NULL(),& ! layer integrated work by diapycnal mixing (W/m2) maxTKE => NULL(),& ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd => NULL() ! conversion rate (~1.0 / (G_Earth + dRho_lay)) + TKE_to_Kd => NULL(),& ! conversion rate (~1.0 / (G_Earth + dRho_lay)) ! between TKE dissipated within a layer and Kd ! in that layer, in m2 s-1 / m3 s-3 = s2 m-1 + KT_extra => NULL(),& ! double diffusion diffusivity for temp (m2/s) + KS_extra => NULL() ! double diffusion diffusivity for saln (m2/s) + end type diffusivity_diags ! Clocks -integer :: id_clock_kappaShear +integer :: id_clock_kappaShear, id_clock_CVMix_ddiff contains @@ -236,7 +247,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & real, dimension(SZI_(G),SZK_(G)+1) :: & N2_int, & !< squared buoyancy frequency associated at interfaces (1/s2) - dRho_int !< locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? + dRho_int, & !< locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? + KT_extra, & !< double difusion diffusivity on temperature (m2/sec) + KS_extra ! double difusion diffusivity on salinity (m2/sec) real :: I_Rho0 ! inverse of Boussinesq density (m3/kg) real :: dissip ! local variable for dissipation calculations (W/m3) @@ -271,10 +284,10 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & use_EOS = associated(tv%eqn_of_state) - if ((CS%use_CVMix_ddiff) .and. & + if ((CS%use_CVMix_ddiff) .or. CS%double_diffusion .and. & .not.(associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S)) ) & call MOM_error(FATAL, "set_diffusivity: visc%Kd_extra_T and "//& - "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF is true.") + "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF or DOUBLE_DIFFUSION are true.") ! Set Kd, Kd_int and Kv_slow to constant values. ! If nothing else is specified, this will be the value used. @@ -299,6 +312,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%id_TKE_to_Kd > 0) then allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0 endif + if ((CS%id_KT_extra > 0) .or. (CS%id_KT_extra_z > 0)) then + allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0 + endif + if ((CS%id_KS_extra > 0) .or. (CS%id_KS_extra_z > 0)) then + allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0 + endif if ((CS%id_Kd_BBL > 0) .or. (CS%id_Kd_BBL_z > 0)) then allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0 endif @@ -375,10 +394,40 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! Add background mixing call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc%Kv_slow, j, G, GV, CS%bkgnd_mixing_csp) - ! Apply double diffusion + ! Double-diffusion (old method) + if (CS%double_diffusion) then + call double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, KT_extra, KS_extra) + do K=2,nz ; do i=is,ie + if (KS_extra(i,K) > KT_extra(i,K)) then ! salt fingering + Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KT_extra(i,K) + Kd(i,j,k) = Kd(i,j,k) + 0.5*KT_extra(i,K) + visc%Kd_extra_S(i,j,k) = KS_extra(i,K)-KT_extra(i,K) + visc%Kd_extra_T(i,j,k) = 0.0 + elseif (KT_extra(i,K) > 0.0) then ! double-diffusive convection + Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KS_extra(i,K) + Kd(i,j,k) = Kd(i,j,k) + 0.5*KS_extra(i,K) + visc%Kd_extra_T(i,j,k) = KT_extra(i,K)-KS_extra(i,K) + visc%Kd_extra_S(i,j,k) = 0.0 + else ! There is no double diffusion at this interface. + visc%Kd_extra_T(i,j,k) = 0.0 + visc%Kd_extra_S(i,j,k) = 0.0 + endif + enddo ; enddo + if (associated(dd%KT_extra)) then ; do K=1,nz+1 ; do i=is,ie + dd%KT_extra(i,j,K) = KT_extra(i,K) + enddo ; enddo ; endif + + if (associated(dd%KS_extra)) then ; do K=1,nz+1 ; do i=is,ie + dd%KS_extra(i,j,K) = KS_extra(i,K) + enddo ; enddo ; endif + endif + + ! Apply double diffusion via CVMix ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available. if (CS%use_CVMix_ddiff) then + call cpu_clock_begin(id_clock_CVMix_ddiff) call compute_ddiff_coeffs(h, tv, G, GV, j, visc%Kd_extra_T, visc%Kd_extra_S, CS%CVMix_ddiff_csp) + call cpu_clock_end(id_clock_CVMix_ddiff) endif ! Add the input turbulent diffusivity. @@ -562,11 +611,26 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif + if (CS%id_KT_extra > 0) call post_data(CS%id_KT_extra, dd%KT_extra, CS%diag) + if (CS%id_KS_extra > 0) call post_data(CS%id_KS_extra, dd%KS_extra, CS%diag) if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, dd%Kd_BBL, CS%diag) + if (CS%id_KT_extra_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_KT_extra_z + z_ptrs(num_z_diags)%p => dd%KT_extra + endif + + if (CS%id_KS_extra_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_KS_extra_z + z_ptrs(num_z_diags)%p => dd%KS_extra + endif + if (CS%id_Kd_BBL_z > 0) then num_z_diags = num_z_diags + 1 z_ids(num_z_diags) = CS%id_Kd_BBL_z + z_ptrs(num_z_diags)%p => dd%KS_extra endif if (num_z_diags > 0) & @@ -577,6 +641,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(dd%Kd_user)) deallocate(dd%Kd_user) if (associated(dd%maxTKE)) deallocate(dd%maxTKE) if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd) + if (associated(dd%KT_extra)) deallocate(dd%KT_extra) + if (associated(dd%KS_extra)) deallocate(dd%KS_extra) if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL) if (showCallTree) call callTree_leave("set_diffusivity()") @@ -929,6 +995,97 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & end subroutine find_N2 +!> This subroutine sets the additional diffusivities of temperature and +!! salinity due to double diffusion, using the same functional form as is +!! used in MOM4.1, and taken from an NCAR technical note (REF?) that updates +!! what was in Large et al. (1994). All the coefficients here should probably +!! be made run-time variables rather than hard-coded constants. +!! +!! \todo Find reference for NCAR tech note above. +subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available + !! thermodynamic fields; absent fields have NULL + !! ptrs. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: T_f !< layer temp in C with the values in massless layers + !! filled vertically by diffusion. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: S_f !< Layer salinities in PPT with values in massless + !! layers filled vertically by diffusion. + integer, intent(in) :: j !< Meridional index upon which to work. + type(set_diffusivity_CS), pointer :: CS !< Module control structure. + real, dimension(SZI_(G),SZK_(G)+1), & + intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal + !! diffusivity for temp (m2/sec). + real, dimension(SZI_(G),SZK_(G)+1), & + intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal + !! diffusivity for saln (m2/sec). + + real, dimension(SZI_(G)) :: & + dRho_dT, & ! partial derivatives of density wrt temp (kg m-3 degC-1) + dRho_dS, & ! partial derivatives of density wrt saln (kg m-3 PPT-1) + pres, & ! pressure at each interface (Pa) + Temp_int, & ! temp and saln at interfaces + Salin_int + + real :: alpha_dT ! density difference between layers due to temp diffs (kg/m3) + real :: beta_dS ! density difference between layers due to saln diffs (kg/m3) + + real :: Rrho ! vertical density ratio + real :: diff_dd ! factor for double-diffusion + real :: prandtl ! flux ratio for diffusive convection regime + + real, parameter :: Rrho0 = 1.9 ! limit for double-diffusive density ratio + real, parameter :: dsfmax = 1.e-4 ! max diffusivity in case of salt fingering + real, parameter :: Kv_molecular = 1.5e-6 ! molecular viscosity (m2/sec) + + integer :: i, k, is, ie, nz + is = G%isc ; ie = G%iec ; nz = G%ke + + if (associated(tv%eqn_of_state)) then + do i=is,ie + pres(i) = 0.0 ; Kd_T_dd(i,1) = 0.0 ; Kd_S_dd(i,1) = 0.0 + Kd_T_dd(i,nz+1) = 0.0 ; Kd_S_dd(i,nz+1) = 0.0 + enddo + do K=2,nz + do i=is,ie + pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1) + Temp_Int(i) = 0.5 * (T_f(i,j,k-1) + T_f(i,j,k)) + Salin_Int(i) = 0.5 * (S_f(i,j,k-1) + S_f(i,j,k)) + enddo + call calculate_density_derivs(Temp_int, Salin_int, pres, & + dRho_dT(:), dRho_dS(:), is, ie-is+1, tv%eqn_of_state) + + do i=is,ie + alpha_dT = -1.0*dRho_dT(i) * (T_f(i,j,k-1) - T_f(i,j,k)) + beta_dS = dRho_dS(i) * (S_f(i,j,k-1) - S_f(i,j,k)) + + if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case + Rrho = min(alpha_dT/beta_dS,Rrho0) + diff_dd = 1.0 - ((RRho-1.0)/(RRho0-1.0)) + diff_dd = dsfmax*diff_dd*diff_dd*diff_dd + Kd_T_dd(i,K) = 0.7*diff_dd + Kd_S_dd(i,K) = diff_dd + elseif ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection + Rrho = alpha_dT/beta_dS + diff_dd = Kv_molecular*0.909*exp(4.6*exp(-0.54*(1/Rrho-1))) + prandtl = 0.15*Rrho + if (Rrho > 0.5) prandtl = (1.85-0.85/Rrho)*Rrho + Kd_T_dd(i,K) = diff_dd + Kd_S_dd(i,K) = prandtl*diff_dd + else + Kd_T_dd(i,K) = 0.0 ; Kd_S_dd(i,K) = 0.0 + endif + enddo + enddo + endif + +end subroutine double_diffusion + !> This routine adds diffusion sustained by flow energy extracted by bottom drag. subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & maxTKE, kb, G, GV, CS, Kd, Kd_int, Kd_BBL) @@ -1945,6 +2102,43 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp endif endif + call get_param(param_file, mdl, "DOUBLE_DIFFUSION", CS%double_diffusion, & + "If true, increase diffusivitives for temperature or salt \n"//& + "based on double-diffusive paramaterization from MOM4/KPP.", & + default=.false.) + if (CS%double_diffusion) then + call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", CS%Max_Rrho_salt_fingers, & + "Maximum density ratio for salt fingering regime.", & + default=2.55, units="nondim") + call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", CS%Max_salt_diff_salt_fingers, & + "Maximum salt diffusivity for salt fingering regime.", & + default=1.e-4, units="m2 s-1") + call get_param(param_file, mdl, "KV_MOLECULAR", CS%Kv_molecular, & + "Molecular viscosity for calculation of fluxes under \n"//& + "double-diffusive convection.", default=1.5e-6, units="m2 s-1") + ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults. + + CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for temperature', 'm2 s-1') + + CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for salinity', 'm2 s-1') + + if (associated(diag_to_Z_CSp)) then + vd = var_desc("KT_extra", "m2 s-1", & + "Double-Diffusive Temperature Diffusivity, interpolated to z", & + z_grid='z') + CS%id_KT_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + vd = var_desc("KS_extra", "m2 s-1", & + "Double-Diffusive Salinity Diffusivity, interpolated to z",& + z_grid='z') + CS%id_KS_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + vd = var_desc("Kd_BBL", "m2 s-1", & + "Bottom Boundary Layer Diffusivity", z_grid='z') + CS%id_Kd_BBL_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + endif + endif + if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif @@ -1962,6 +2156,8 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! CVMix double diffusion mixing CS%use_CVMix_ddiff = CVMix_ddiff_init(Time, G, GV, param_file, CS%diag, CS%CVMix_ddiff_csp) + if (CS%use_CVMix_ddiff) & + id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_MODULE) end subroutine set_diffusivity_init