MOM6
Data Types | Variables
mom_cvmix_kpp Module Reference

Detailed Description

Provides the K-Profile Parameterization (KPP) of Large et al., 1994, via CVMix.

The K-Profile Parameterization

The K-Profile Parameterization (KPP) of Large et al., 1994, (http://dx.doi.org/10.1029/94RG01872) is implemented via the Community Vertical Mixing package, CVMix, which is called directly by this module.

The formulation and implementation of KPP is described in great detail in the CVMix manual (written by our own Steve Griffies).

KPP in a nutshell

Large et al., 1994, decompose the parameterized boundary layer turbulent flux of a scalar, \( s \), as

\[ \overline{w^\prime s^\prime} = -K \partial_z s + K \gamma_s(\sigma), \]

where \( \sigma = -z/h \) is a non-dimensional coordinate within the boundary layer of depth \( h \). \( K \) is the eddy diffusivity and is a function of position within the boundary layer as well as a function of the surface forcing:

\[ K = h w_s(\sigma) G(\sigma) . \]

Here, \( w_s \) is the vertical velocity scale of the boundary layer turbulence and \( G(\sigma) \) is a "shape function" which is described later. The last term is the "non-local transport" which involves a function \( \gamma_s(\sigma) \) that is matched to the forcing but is not actually needed in the final implementation. Instead, the entire non-local transport term can be equivalently written

\[ K \gamma_s(\sigma) = C_s G(\sigma) Q_s \]

where \( Q_s \) is the surface flux of \( s \) and \( C_s \) is a constant. The vertical structure of the redistribution (non-local) term is solely due to the shape function, \( G(\sigma) \). In our implementation of KPP, we allow the shape functions used for \( K \) and for the non-local transport to be chosen independently.

The particular shape function most widely used in the atmospheric community is

\[ G(\sigma) = \sigma (1-\sigma)^2 \]

which satisfies the boundary conditions \( G(0) = 0 \), \( G(1) = 0 \), \( G^\prime(0) = 1 \), and \( G^\prime(1) = 0 \). Large et al, 1994, alter the function so as to match interior diffusivities but we have found that this leads to inconsistencies within the formulation (see google groups thread Extreme values of non-local transport). Instead, we use either the above form, or even simpler forms that use alternative upper boundary conditions.

The KPP boundary layer depth is a function of the bulk Richardson number, Rib. But to compute Rib, we need the boundary layer depth. To address this circular logic, we compute Rib for each vertical cell in a column, assuming the BL depth equals to the depth of the given grid cell. Once we have a vertical array of Rib(k), we then call the OBLdepth routine from CVMix to compute the actual OBLdepth. We optionally then "correct" the OBLdepth by cycling through once more, this time knowing the OBLdepth from the first pass. This "correction" step is not used by NCAR. It has been found in idealized MOM6 tests to not be necessary.

See also
kpp_calculate(), kpp_applynonlocaltransport()

Data Types

type  kpp_cs
 Control structure for containing KPP parameters/data. More...
 

Functions/Subroutines

logical function, public kpp_init (paramFile, G, GV, US, diag, Time, CS, passive, Waves)
 Initialize the CVMix KPP module and set up diagnostics Returns True if KPP is to be used, False otherwise. More...
 
subroutine, public kpp_calculate (CS, G, GV, US, h, uStar, buoyFlux, Kt, Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves)
 KPP vertical diffusivity/viscosity and non-local tracer transport. More...
 
subroutine, public kpp_compute_bld (CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Waves)
 Compute OBL depth. More...
 
subroutine kpp_smooth_bld (CS, G, GV, h)
 Apply a 1-1-4-1-1 Laplacian filter one time on BLD to reduce any horizontal two-grid-point noise. More...
 
subroutine, public kpp_get_bld (CS, BLD, G)
 Copies KPP surface boundary layer depth into BLD. More...
 
subroutine, public kpp_nonlocaltransport_temp (CS, G, GV, h, nonLocalTrans, surfFlux, dt, scalar, C_p)
 Apply KPP non-local transport of surface fluxes for temperature. More...
 
subroutine, public kpp_nonlocaltransport_saln (CS, G, GV, h, nonLocalTrans, surfFlux, dt, scalar)
 Apply KPP non-local transport of surface fluxes for salinity. This routine is a useful prototype for other material tracers. More...
 
subroutine, public kpp_end (CS)
 Clear pointers, deallocate memory. More...
 

Variables

integer, parameter, private nlt_shape_cvmix = 0
 Use the CVMix profile. More...
 
integer, parameter, private nlt_shape_linear = 1
 Linear, \( G(\sigma) = 1-\sigma \). More...
 
integer, parameter, private nlt_shape_parabolic = 2
 Parabolic, \( G(\sigma) = (1-\sigma)^2 \). More...
 
integer, parameter, private nlt_shape_cubic = 3
 Cubic, \( G(\sigma) = 1 + (2\sigma-3) \sigma^2\). More...
 
integer, parameter, private nlt_shape_cubic_lmd = 4
 Original shape, \( G(\sigma) = \frac{27}{4} \sigma (1-\sigma)^2 \). More...
 
integer, parameter, private sw_method_all_sw = 0
 Use all shortwave radiation. More...
 
integer, parameter, private sw_method_mxl_sw = 1
 Use shortwave radiation absorbed in mixing layer. More...
 
integer, parameter, private sw_method_lv1_sw = 2
 Use shortwave radiation absorbed in layer 1. More...
 
integer, parameter, private lt_k_constant = 1
 Constant enhance K through column. More...
 
integer, parameter, private lt_k_scaled = 2
 Enhance K scales with G(sigma) More...
 
integer, parameter, private lt_k_mode_constant = 1
 Prescribed enhancement for K. More...
 
integer, parameter, private lt_k_mode_vr12 = 2
 Enhancement for K based on. More...
 
integer, parameter, private lt_k_mode_rw16 = 3
 Enhancement for K based on. More...
 
integer, parameter, private lt_vt2_mode_constant = 1
 Prescribed enhancement for Vt2. More...
 
integer, parameter, private lt_vt2_mode_vr12 = 2
 Enhancement for Vt2 based on. More...
 
integer, parameter, private lt_vt2_mode_rw16 = 3
 Enhancement for Vt2 based on. More...
 
integer, parameter, private lt_vt2_mode_lf17 = 4
 Enhancement for Vt2 based on. More...
 
integer ncall_smooth = 0
 

Function/Subroutine Documentation

◆ kpp_calculate()

subroutine, public mom_cvmix_kpp::kpp_calculate ( type(kpp_cs), pointer  CS,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g)), intent(in)  uStar,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(in)  buoyFlux,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout)  Kt,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout)  Ks,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout)  Kv,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout)  nonLocalTransHeat,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout)  nonLocalTransScalar,
type(wave_parameters_cs), optional, pointer  waves 
)

KPP vertical diffusivity/viscosity and non-local tracer transport.

Parameters
csControl structure
[in]gOcean grid
[in]gvOcean vertical grid
[in]usA dimensional unit scaling type
wavesWave CS
[in]hLayer/level thicknesses [H ~> m or kg m-2]
[in]ustarSurface friction velocity [Z T-1 ~> m s-1]
[in]buoyfluxSurface buoyancy flux [L2 T-3 ~> m2 s-3]
[in,out]kt(in) Vertical diffusivity of heat w/o KPP (out) Vertical diffusivity including KPP [Z2 T-1 ~> m2 s-1]
[in,out]ks(in) Vertical diffusivity of salt w/o KPP (out) Vertical diffusivity including KPP [Z2 T-1 ~> m2 s-1]
[in,out]kv(in) Vertical viscosity w/o KPP (out) Vertical viscosity including KPP [Z2 T-1 ~> m2 s-1]
[in,out]nonlocaltransheatTemp non-local transport [m s-1]
[in,out]nonlocaltransscalarscalar non-local transport [m s-1]

Definition at line 590 of file MOM_CVMix_KPP.F90.

590 
591  ! Arguments
592  type(KPP_CS), pointer :: CS !< Control structure
593  type(ocean_grid_type), intent(in) :: G !< Ocean grid
594  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
595  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
596  type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS
597  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]
598  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity [Z T-1 ~> m s-1]
599  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
600  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kt !< (in) Vertical diffusivity of heat w/o KPP
601  !! (out) Vertical diffusivity including KPP
602  !! [Z2 T-1 ~> m2 s-1]
603  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Ks !< (in) Vertical diffusivity of salt w/o KPP
604  !! (out) Vertical diffusivity including KPP
605  !! [Z2 T-1 ~> m2 s-1]
606  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kv !< (in) Vertical viscosity w/o KPP
607  !! (out) Vertical viscosity including KPP
608  !! [Z2 T-1 ~> m2 s-1]
609  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransHeat !< Temp non-local transport [m s-1]
610  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local transport [m s-1]
611 
612 ! Local variables
613  integer :: i, j, k ! Loop indices
614  real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface [m] (negative in ocean)
615  real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m] (negative in ocean)
616  real, dimension( G%ke+1, 2) :: Kdiffusivity ! Vertical diffusivity at interfaces [m2 s-1]
617  real, dimension( G%ke+1 ) :: Kviscosity ! Vertical viscosity at interfaces [m2 s-1]
618  real, dimension( G%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces [nondim]
619 
620  real :: surfFricVel, surfBuoyFlux
621  real :: sigma, sigmaRatio
622  real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim]
623  real :: dh ! The local thickness used for calculating interface positions [m]
624  real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
625 
626  ! For Langmuir Calculations
627  real :: LangEnhK ! Langmuir enhancement for mixing coefficient
628 
629 
630 #ifdef __DO_SAFETY_CHECKS__
631  if (cs%debug) then
632  call hchksum(h, "KPP in: h",g%HI,haloshift=0, scale=gv%H_to_m)
633  call hchksum(ustar, "KPP in: uStar",g%HI,haloshift=0, scale=us%Z_to_m*us%s_to_T)
634  call hchksum(buoyflux, "KPP in: buoyFlux",g%HI,haloshift=0)
635  call hchksum(kt, "KPP in: Kt",g%HI,haloshift=0, scale=us%Z2_T_to_m2_s)
636  call hchksum(ks, "KPP in: Ks",g%HI,haloshift=0, scale=us%Z2_T_to_m2_s)
637  endif
638 #endif
639 
640  nonlocaltrans(:,:) = 0.0
641 
642  if (cs%id_Kd_in > 0) call post_data(cs%id_Kd_in, kt, cs%diag)
643 
644  buoy_scale = us%L_to_m**2*us%s_to_T**3
645 
646  !$OMP parallel do default(none) firstprivate(nonLocalTrans) &
647  !$OMP private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, &
648  !$OMP surfBuoyFlux, Kdiffusivity, Kviscosity, LangEnhK, sigma, &
649  !$OMP sigmaRatio) &
650  !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, Kt, &
651  !$OMP Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves)
652  ! loop over horizontal points on processor
653  do j = g%jsc, g%jec
654  do i = g%isc, g%iec
655 
656  ! skip calling KPP for land points
657  if (g%mask2dT(i,j)==0.) cycle
658 
659  ! things independent of position within the column
660  surffricvel = us%Z_to_m*us%s_to_T * ustar(i,j)
661 
662  ifaceheight(1) = 0.0 ! BBL is all relative to the surface
663  hcorr = 0.
664  do k=1,g%ke
665 
666  ! cell center and cell bottom in meters (negative values in the ocean)
667  dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
668  dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
669  hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
670  dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
671  cellheight(k) = ifaceheight(k) - 0.5 * dh
672  ifaceheight(k+1) = ifaceheight(k) - dh
673 
674  enddo ! k-loop finishes
675 
676  surfbuoyflux = buoy_scale*buoyflux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
677  ! h to Monin-Obukov (default is false, ie. not used)
678 
679  ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports
680 
681  ! Unlike LMD94, we do not match to interior diffusivities. If using the original
682  ! LMD94 shape function, not matching is equivalent to matching to a zero diffusivity.
683 
684  !BGR/ Add option for use of surface buoyancy flux with total sw flux.
685  if (cs%SW_METHOD == sw_method_all_sw) then
686  surfbuoyflux = buoy_scale * buoyflux(i,j,1)
687  elseif (cs%SW_METHOD == sw_method_mxl_sw) then
688  ! We know the actual buoyancy flux into the OBL
689  surfbuoyflux = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,int(cs%kOBL(i,j))+1))
690  elseif (cs%SW_METHOD == sw_method_lv1_sw) then
691  surfbuoyflux = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,2))
692  endif
693 
694  ! If option "MatchBoth" is selected in CVMix, MOM should be capable of matching.
695  if (.not. (cs%MatchTechnique == 'MatchBoth')) then
696  kdiffusivity(:,:) = 0. ! Diffusivities for heat and salt [m2 s-1]
697  kviscosity(:) = 0. ! Viscosity [m2 s-1]
698  else
699  kdiffusivity(:,1) = us%Z2_T_to_m2_s * kt(i,j,:)
700  kdiffusivity(:,2) = us%Z2_T_to_m2_s * ks(i,j,:)
701  kviscosity(:) = us%Z2_T_to_m2_s * kv(i,j,:)
702  endif
703 
704  call cvmix_coeffs_kpp(kviscosity(:), & ! (inout) Total viscosity [m2 s-1]
705  kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1]
706  kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1]
707  ifaceheight, & ! (in) Height of interfaces [m]
708  cellheight, & ! (in) Height of level centers [m]
709  kviscosity(:), & ! (in) Original viscosity [m2 s-1]
710  kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1]
711  kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1]
712  cs%OBLdepth(i,j), & ! (in) OBL depth [m]
713  cs%kOBL(i,j), & ! (in) level (+fraction) of OBL extent
714  nonlocaltrans(:,1),& ! (out) Non-local heat transport [nondim]
715  nonlocaltrans(:,2),& ! (out) Non-local salt transport [nondim]
716  surffricvel, & ! (in) Turbulent friction velocity at surface [m s-1]
717  surfbuoyflux, & ! (in) Buoyancy flux at surface [m2 s-3]
718  g%ke, & ! (in) Number of levels to compute coeffs for
719  g%ke, & ! (in) Number of levels in array shape
720  cvmix_kpp_params_user=cs%KPP_params )
721 
722  ! safety check, Kviscosity and Kdiffusivity must be >= 0
723  do k=1, g%ke+1
724  if (kviscosity(k) < 0. .or. kdiffusivity(k,1) < 0.) then
725  call mom_error(fatal,"KPP_calculate, after CVMix_coeffs_kpp: "// &
726  "Negative vertical viscosity or diffusivity has been detected. " // &
727  "This is likely related to the choice of MATCH_TECHNIQUE and INTERP_TYPE2." //&
728  "You might consider using the default options for these parameters." )
729  endif
730  enddo
731 
732  IF (cs%LT_K_ENHANCEMENT) then
733  if (cs%LT_K_METHOD==lt_k_mode_constant) then
734  langenhk = cs%KPP_K_ENH_FAC
735  elseif (cs%LT_K_METHOD==lt_k_mode_vr12) then
736  ! Added minimum value for La_SL, so removed maximum value for LangEnhK.
737  langenhk = sqrt(1.+(1.5*cs%La_SL(i,j))**(-2) + &
738  (5.4*cs%La_SL(i,j))**(-4))
739  elseif (cs%LT_K_METHOD==lt_k_mode_rw16) then
740  !This maximum value is proposed in Reichl et al., 2016 JPO formula
741  langenhk = min(2.25, 1. + 1./cs%La_SL(i,j))
742  else
743  !This shouldn't be reached.
744  !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT")
745  langenhk = 1.0
746  endif
747  do k=1,g%ke
748  if (cs%LT_K_SHAPE== lt_k_constant) then
749  if (cs%id_EnhK > 0) cs%EnhK(i,j,:) = langenhk
750  kdiffusivity(k,1) = kdiffusivity(k,1) * langenhk
751  kdiffusivity(k,2) = kdiffusivity(k,2) * langenhk
752  kviscosity(k) = kviscosity(k) * langenhk
753  elseif (cs%LT_K_SHAPE == lt_k_scaled) then
754  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
755  sigmaratio = sigma * (1. - sigma)**2 / 0.148148037
756  if (cs%id_EnhK > 0) cs%EnhK(i,j,k) = (1.0 + (langenhk - 1.)*sigmaratio)
757  kdiffusivity(k,1) = kdiffusivity(k,1) * ( 1. + &
758  ( langenhk - 1.)*sigmaratio)
759  kdiffusivity(k,2) = kdiffusivity(k,2) * ( 1. + &
760  ( langenhk - 1.)*sigmaratio)
761  kviscosity(k) = kviscosity(k) * ( 1. + &
762  ( langenhk - 1.)*sigmaratio)
763  endif
764  enddo
765  endif
766 
767  ! Over-write CVMix NLT shape function with one of the following choices.
768  ! The CVMix code has yet to update for thse options, so we compute in MOM6.
769  ! Note that nonLocalTrans = Cs * G(sigma) (LMD94 notation), with
770  ! Cs = 6.32739901508.
771  ! Start do-loop at k=2, since k=1 is ocean surface (sigma=0)
772  ! and we do not wish to double-count the surface forcing.
773  ! Only compute nonlocal transport for 0 <= sigma <= 1.
774  ! MOM6 recommended shape is the parabolic; it gives deeper boundary layer
775  ! and no spurious extrema.
776  if (surfbuoyflux < 0.0) then
777  if (cs%NLT_shape == nlt_shape_cubic) then
778  do k = 2, g%ke
779  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
780  nonlocaltrans(k,1) = (1.0 - sigma)**2 * (1.0 + 2.0*sigma) !*
781  nonlocaltrans(k,2) = nonlocaltrans(k,1)
782  enddo
783  elseif (cs%NLT_shape == nlt_shape_parabolic) then
784  do k = 2, g%ke
785  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
786  nonlocaltrans(k,1) = (1.0 - sigma)**2 !*CS%CS2
787  nonlocaltrans(k,2) = nonlocaltrans(k,1)
788  enddo
789  elseif (cs%NLT_shape == nlt_shape_linear) then
790  do k = 2, g%ke
791  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
792  nonlocaltrans(k,1) = (1.0 - sigma)!*CS%CS2
793  nonlocaltrans(k,2) = nonlocaltrans(k,1)
794  enddo
795  elseif (cs%NLT_shape == nlt_shape_cubic_lmd) then
796  ! Sanity check (should agree with CVMix result using simple matching)
797  do k = 2, g%ke
798  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
799  nonlocaltrans(k,1) = cs%CS2 * sigma*(1.0 -sigma)**2
800  nonlocaltrans(k,2) = nonlocaltrans(k,1)
801  enddo
802  endif
803  endif
804 
805  ! we apply nonLocalTrans in subroutines
806  ! KPP_NonLocalTransport_temp and KPP_NonLocalTransport_saln
807  nonlocaltransheat(i,j,:) = nonlocaltrans(:,1) ! temp
808  nonlocaltransscalar(i,j,:) = nonlocaltrans(:,2) ! saln
809 
810  ! set the KPP diffusivity and viscosity to zero for testing purposes
811  if (cs%KPPzeroDiffusivity) then
812  kdiffusivity(:,1) = 0.0
813  kdiffusivity(:,2) = 0.0
814  kviscosity(:) = 0.0
815  endif
816 
817 
818  ! compute unresolved squared velocity for diagnostics
819  if (cs%id_Vt2 > 0) then
820 !BGR Now computing VT2 above so can modify for LT
821 ! therefore, don't repeat this operation here
822 ! CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( &
823 ! cellHeight(1:G%ke), & ! Depth of cell center [m]
824 ! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1]
825 ! N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1]
826 ! CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters
827  endif
828 
829  ! Copy 1d data into 3d diagnostic arrays
830  !/ grabbing obldepth_0d for next time step.
831  cs%OBLdepthprev(i,j)=cs%OBLdepth(i,j)
832  if (cs%id_sigma > 0) then
833  cs%sigma(i,j,:) = 0.
834  if (cs%OBLdepth(i,j)>0.) cs%sigma(i,j,:) = -ifaceheight/cs%OBLdepth(i,j)
835  endif
836  if (cs%id_Kt_KPP > 0) cs%Kt_KPP(i,j,:) = kdiffusivity(:,1)
837  if (cs%id_Ks_KPP > 0) cs%Ks_KPP(i,j,:) = kdiffusivity(:,2)
838  if (cs%id_Kv_KPP > 0) cs%Kv_KPP(i,j,:) = kviscosity(:)
839 
840  ! Update output of routine
841  if (.not. cs%passiveMode) then
842  if (cs%KPPisAdditive) then
843  do k=1, g%ke+1
844  kt(i,j,k) = kt(i,j,k) + us%m2_s_to_Z2_T * kdiffusivity(k,1)
845  ks(i,j,k) = ks(i,j,k) + us%m2_s_to_Z2_T * kdiffusivity(k,2)
846  kv(i,j,k) = kv(i,j,k) + us%m2_s_to_Z2_T * kviscosity(k)
847  if (cs%Stokes_Mixing) waves%KvS(i,j,k) = kv(i,j,k)
848  enddo
849  else ! KPP replaces prior diffusivity when former is non-zero
850  do k=1, g%ke+1
851  if (kdiffusivity(k,1) /= 0.) kt(i,j,k) = us%m2_s_to_Z2_T * kdiffusivity(k,1)
852  if (kdiffusivity(k,2) /= 0.) ks(i,j,k) = us%m2_s_to_Z2_T * kdiffusivity(k,2)
853  if (kviscosity(k) /= 0.) kv(i,j,k) = us%m2_s_to_Z2_T * kviscosity(k)
854  if (cs%Stokes_Mixing) waves%KvS(i,j,k) = kv(i,j,k)
855  enddo
856  endif
857  endif
858 
859 
860  ! end of the horizontal do-loops over the vertical columns
861  enddo ! i
862  enddo ! j
863 
864 
865 #ifdef __DO_SAFETY_CHECKS__
866  if (cs%debug) then
867  call hchksum(kt, "KPP out: Kt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
868  call hchksum(ks, "KPP out: Ks", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
869  endif
870 #endif
871 
872  ! send diagnostics to post_data
873  if (cs%id_OBLdepth > 0) call post_data(cs%id_OBLdepth, cs%OBLdepth, cs%diag)
874  if (cs%id_OBLdepth_original > 0) call post_data(cs%id_OBLdepth_original,cs%OBLdepth_original,cs%diag)
875  if (cs%id_sigma > 0) call post_data(cs%id_sigma, cs%sigma, cs%diag)
876  if (cs%id_Ws > 0) call post_data(cs%id_Ws, cs%Ws, cs%diag)
877  if (cs%id_Vt2 > 0) call post_data(cs%id_Vt2, cs%Vt2, cs%diag)
878  if (cs%id_uStar > 0) call post_data(cs%id_uStar, ustar, cs%diag)
879  if (cs%id_buoyFlux > 0) call post_data(cs%id_buoyFlux, buoyflux, cs%diag)
880  if (cs%id_Kt_KPP > 0) call post_data(cs%id_Kt_KPP, cs%Kt_KPP, cs%diag)
881  if (cs%id_Ks_KPP > 0) call post_data(cs%id_Ks_KPP, cs%Ks_KPP, cs%diag)
882  if (cs%id_Kv_KPP > 0) call post_data(cs%id_Kv_KPP, cs%Kv_KPP, cs%diag)
883  if (cs%id_NLTt > 0) call post_data(cs%id_NLTt, nonlocaltransheat, cs%diag)
884  if (cs%id_NLTs > 0) call post_data(cs%id_NLTs, nonlocaltransscalar,cs%diag)
885 
886 

References lt_k_constant, lt_k_mode_constant, lt_k_mode_rw16, lt_k_mode_vr12, lt_k_scaled, mom_error_handler::mom_error(), nlt_shape_cubic, nlt_shape_cubic_lmd, nlt_shape_linear, nlt_shape_parabolic, sw_method_all_sw, sw_method_lv1_sw, and sw_method_mxl_sw.

Here is the call graph for this function:

◆ kpp_compute_bld()

subroutine, public mom_cvmix_kpp::kpp_compute_bld ( type(kpp_cs), pointer  CS,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  Temp,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  Salt,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
type(eos_type), pointer  EOS,
real, dimension(szi_(g),szj_(g)), intent(in)  uStar,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(in)  buoyFlux,
type(wave_parameters_cs), optional, pointer  Waves 
)

Compute OBL depth.

Parameters
csControl structure
[in,out]gOcean grid
[in]gvOcean vertical grid
[in]usA dimensional unit scaling type
[in]hLayer/level thicknesses [H ~> m or kg m-2]
[in]temppotential/cons temp [degC]
[in]saltSalinity [ppt]
[in]uVelocity i-component [L T-1 ~> m s-1]
[in]vVelocity j-component [L T-1 ~> m s-1]
eosEquation of state
[in]ustarSurface friction velocity [Z T-1 ~> m s-1]
[in]buoyfluxSurface buoyancy flux [L2 T-3 ~> m2 s-3]
wavesWave CS

Definition at line 892 of file MOM_CVMix_KPP.F90.

892 
893  ! Arguments
894  type(KPP_CS), pointer :: CS !< Control structure
895  type(ocean_grid_type), intent(inout) :: G !< Ocean grid
896  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
897  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
898  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]
899  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Temp !< potential/cons temp [degC]
900  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Salt !< Salinity [ppt]
901  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Velocity i-component [L T-1 ~> m s-1]
902  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Velocity j-component [L T-1 ~> m s-1]
903  type(EOS_type), pointer :: EOS !< Equation of state
904  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity [Z T-1 ~> m s-1]
905  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
906  type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS
907 
908  ! Local variables
909  integer :: i, j, k, km1 ! Loop indices
910  real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface [m] (negative in ocean)
911  real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m] (negative in ocean)
912  real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces [s-2]
913  real, dimension( G%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars [m s-1]
914  real, dimension( G%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number
915  real, dimension( G%ke ) :: deltaU2 ! square of delta U (shear) in denominator of Bulk Ri [m2 s-2]
916  real, dimension( G%ke ) :: surfBuoyFlux2
917  real, dimension( G%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer
918 
919  ! for EOS calculation
920  real, dimension( 3*G%ke ) :: rho_1D
921  real, dimension( 3*G%ke ) :: pres_1D
922  real, dimension( 3*G%ke ) :: Temp_1D
923  real, dimension( 3*G%ke ) :: Salt_1D
924 
925  real :: surfFricVel, surfBuoyFlux, Coriolis
926  real :: GoRho ! Gravitational acceleration divided by density in MKS units [m4 s-2]
927  real :: pRef, rho1, rhoK, Uk, Vk, sigma, sigmaRatio
928 
929  real :: zBottomMinusOffset ! Height of bottom plus a little bit [m]
930  real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth.
931  real :: hTot ! Running sum of thickness used in the surface layer average [m]
932  real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim]
933  real :: delH ! Thickness of a layer [m]
934  real :: surfHtemp, surfTemp ! Integral and average of temp over the surface layer
935  real :: surfHsalt, surfSalt ! Integral and average of saln over the surface layer
936  real :: surfHu, surfU ! Integral and average of u over the surface layer
937  real :: surfHv, surfV ! Integral and average of v over the surface layer
938  real :: dh ! The local thickness used for calculating interface positions [m]
939  real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
940  integer :: kk, ksfc, ktmp
941 
942  ! For Langmuir Calculations
943  real :: LangEnhW ! Langmuir enhancement for turbulent velocity scale
944  real, dimension(G%ke) :: LangEnhVt2 ! Langmuir enhancement for unresolved shear
945  real, dimension(G%ke) :: U_H, V_H
946  real :: MLD_GUESS, LA
947  real :: surfHuS, surfHvS, surfUs, surfVs, wavedir, currentdir
948  real :: VarUp, VarDn, M, VarLo, VarAvg
949  real :: H10pct, H20pct,CMNFACT, USx20pct, USy20pct, enhvt2
950  integer :: B
951  real :: WST
952 
953 
954 #ifdef __DO_SAFETY_CHECKS__
955  if (cs%debug) then
956  call hchksum(salt, "KPP in: S",g%HI,haloshift=0)
957  call hchksum(temp, "KPP in: T",g%HI,haloshift=0)
958  call hchksum(u, "KPP in: u",g%HI,haloshift=0)
959  call hchksum(v, "KPP in: v",g%HI,haloshift=0)
960  endif
961 #endif
962 
963  ! some constants
964  gorho = gv%mks_g_Earth / (us%R_to_kg_m3*gv%Rho0)
965  buoy_scale = us%L_to_m**2*us%s_to_T**3
966 
967  ! loop over horizontal points on processor
968  !$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, &
969  !$OMP surfBuoyFlux, U_H, V_H, u, v, Coriolis, pRef, SLdepth_0d, &
970  !$OMP ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS, &
971  !$OMP surfHvS, hTot, delH, surftemp, surfsalt, surfu, surfv, &
972  !$OMP surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D, &
973  !$OMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_GUESS, LA, rho_1D, &
974  !$OMP deltarho, N2_1d, ws_1d, LangEnhVT2, enhvt2, wst, &
975  !$OMP BulkRi_1d, zBottomMinusOffset) &
976  !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, &
977  !$OMP Temp, Salt, waves, EOS, GoRho)
978  do j = g%jsc, g%jec
979  do i = g%isc, g%iec
980 
981  ! skip calling KPP for land points
982  if (g%mask2dT(i,j)==0.) cycle
983 
984  do k=1,g%ke
985  u_h(k) = 0.5 * us%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k))
986  v_h(k) = 0.5 * us%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k))
987  enddo
988 
989  ! things independent of position within the column
990  coriolis = 0.25*us%s_to_T*( (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
991  (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1)) )
992  surffricvel = us%Z_to_m*us%s_to_T * ustar(i,j)
993 
994  ! Bullk Richardson number computed for each cell in a column,
995  ! assuming OBLdepth = grid cell depth. After Rib(k) is
996  ! known for the column, then CVMix interpolates to find
997  ! the actual OBLdepth. This approach avoids need to iterate
998  ! on the OBLdepth calculation. It follows that used in MOM5
999  ! and POP.
1000  ifaceheight(1) = 0.0 ! BBL is all relative to the surface
1001  pref = 0.
1002  hcorr = 0.
1003  do k=1,g%ke
1004 
1005  ! cell center and cell bottom in meters (negative values in the ocean)
1006  dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
1007  dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
1008  hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
1009  dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
1010  cellheight(k) = ifaceheight(k) - 0.5 * dh
1011  ifaceheight(k+1) = ifaceheight(k) - dh
1012 
1013  ! find ksfc for cell where "surface layer" sits
1014  sldepth_0d = cs%surf_layer_ext*max( max(-cellheight(k),-ifaceheight(2) ), cs%minOBLdepth )
1015  ksfc = k
1016  do ktmp = 1,k
1017  if (-1.0*ifaceheight(ktmp+1) >= sldepth_0d) then
1018  ksfc = ktmp
1019  exit
1020  endif
1021  enddo
1022 
1023  ! average temp, saln, u, v over surface layer
1024  ! use C-grid average to get u,v on T-points.
1025  surfhtemp=0.0
1026  surfhsalt=0.0
1027  surfhu =0.0
1028  surfhv =0.0
1029  surfhus =0.0
1030  surfhvs =0.0
1031  htot =0.0
1032  do ktmp = 1,ksfc
1033 
1034  ! SLdepth_0d can be between cell interfaces
1035  delh = min( max(0.0, sldepth_0d - htot), h(i,j,ktmp)*gv%H_to_m )
1036 
1037  ! surface layer thickness
1038  htot = htot + delh
1039 
1040  ! surface averaged fields
1041  surfhtemp = surfhtemp + temp(i,j,ktmp) * delh
1042  surfhsalt = surfhsalt + salt(i,j,ktmp) * delh
1043  surfhu = surfhu + 0.5*us%L_T_to_m_s*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delh
1044  surfhv = surfhv + 0.5*us%L_T_to_m_s*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delh
1045  if (cs%Stokes_Mixing) then
1046  surfhus = surfhus + 0.5*(waves%US_x(i,j,ktmp)+waves%US_x(i-1,j,ktmp)) * delh
1047  surfhvs = surfhvs + 0.5*(waves%US_y(i,j,ktmp)+waves%US_y(i,j-1,ktmp)) * delh
1048  endif
1049 
1050  enddo
1051  surftemp = surfhtemp / htot
1052  surfsalt = surfhsalt / htot
1053  surfu = surfhu / htot
1054  surfv = surfhv / htot
1055  surfus = surfhus / htot
1056  surfvs = surfhvs / htot
1057 
1058  ! vertical shear between present layer and
1059  ! surface layer averaged surfU,surfV.
1060  ! C-grid average to get Uk and Vk on T-points.
1061  uk = 0.5*us%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) - surfu
1062  vk = 0.5*us%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) - surfv
1063 
1064  if (cs%Stokes_Mixing) then
1065  ! If momentum is mixed down the Stokes drift gradient, then
1066  ! the Stokes drift must be included in the bulk Richardson number
1067  ! calculation.
1068  uk = uk + (0.5*(waves%Us_x(i,j,k)+waves%US_x(i-1,j,k)) -surfus )
1069  vk = vk + (0.5*(waves%Us_y(i,j,k)+waves%Us_y(i,j-1,k)) -surfvs )
1070  endif
1071 
1072  deltau2(k) = uk**2 + vk**2
1073 
1074  ! pressure, temp, and saln for EOS
1075  ! kk+1 = surface fields
1076  ! kk+2 = k fields
1077  ! kk+3 = km1 fields
1078  km1 = max(1, k-1)
1079  kk = 3*(k-1)
1080  pres_1d(kk+1) = pref
1081  pres_1d(kk+2) = pref
1082  pres_1d(kk+3) = pref
1083  temp_1d(kk+1) = surftemp
1084  temp_1d(kk+2) = temp(i,j,k)
1085  temp_1d(kk+3) = temp(i,j,km1)
1086  salt_1d(kk+1) = surfsalt
1087  salt_1d(kk+2) = salt(i,j,k)
1088  salt_1d(kk+3) = salt(i,j,km1)
1089 
1090  ! pRef is pressure at interface between k and km1.
1091  ! iterate pRef for next pass through k-loop.
1092  pref = pref + gv%H_to_Pa * h(i,j,k)
1093 
1094  ! this difference accounts for penetrating SW
1095  surfbuoyflux2(k) = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,k+1))
1096 
1097  enddo ! k-loop finishes
1098 
1099  if (cs%LT_K_ENHANCEMENT .or. cs%LT_VT2_ENHANCEMENT) then
1100  mld_guess = max( 1.*us%m_to_Z, abs(us%m_to_Z*cs%OBLdepthprev(i,j) ) )
1101  call get_langmuir_number(la, g, gv, us, mld_guess, ustar(i,j), i, j, &
1102  h=h(i,j,:), u_h=u_h, v_h=v_h, waves=waves)
1103  cs%La_SL(i,j)=la
1104  endif
1105 
1106 
1107  ! compute in-situ density
1108  call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, 1, 3*g%ke, eos)
1109 
1110  ! N2 (can be negative) and N (non-negative) on interfaces.
1111  ! deltaRho is non-local rho difference used for bulk Richardson number.
1112  ! CS%N is local N (with floor) used for unresolved shear calculation.
1113  do k = 1, g%ke
1114  km1 = max(1, k-1)
1115  kk = 3*(k-1)
1116  deltarho(k) = rho_1d(kk+2) - rho_1d(kk+1)
1117  n2_1d(k) = (gorho * (rho_1d(kk+2) - rho_1d(kk+3)) ) / &
1118  ((0.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_m)
1119  cs%N(i,j,k) = sqrt( max( n2_1d(k), 0.) )
1120  enddo
1121  n2_1d(g%ke+1 ) = 0.0
1122  cs%N(i,j,g%ke+1 ) = 0.0
1123 
1124  ! turbulent velocity scales w_s and w_m computed at the cell centers.
1125  ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales
1126  ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass
1127  ! sigma=CS%surf_layer_ext for this calculation.
1128  call cvmix_kpp_compute_turbulent_scales( &
1129  cs%surf_layer_ext, & ! (in) Normalized surface layer depth; sigma = CS%surf_layer_ext
1130  -cellheight, & ! (in) Assume here that OBL depth [m] = -cellHeight(k)
1131  surfbuoyflux2, & ! (in) Buoyancy flux at surface [m2 s-3]
1132  surffricvel, & ! (in) Turbulent friction velocity at surface [m s-1]
1133  w_s=ws_1d, & ! (out) Turbulent velocity scale profile [m s-1]
1134  cvmix_kpp_params_user=cs%KPP_params )
1135 
1136  !Compute CVMix VT2
1137  cs%Vt2(i,j,:) = cvmix_kpp_compute_unresolved_shear( &
1138  zt_cntr=cellheight(1:g%ke), & ! Depth of cell center [m]
1139  ws_cntr=ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1]
1140  n_iface=cs%N(i,j,:), & ! Buoyancy frequency at interface [s-1]
1141  cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
1142 
1143  !Modify CVMix VT2
1144  IF (cs%LT_VT2_ENHANCEMENT) then
1145  IF (cs%LT_VT2_METHOD==lt_vt2_mode_constant) then
1146  do k=1,g%ke
1147  langenhvt2(k) = cs%KPP_VT2_ENH_FAC
1148  enddo
1149  elseif (cs%LT_VT2_METHOD==lt_vt2_mode_vr12) then
1150  !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
1151  enhvt2 = sqrt(1.+(1.5*cs%La_SL(i,j))**(-2) + &
1152  (5.4*cs%La_SL(i,j))**(-4))
1153  do k=1,g%ke
1154  langenhvt2(k) = enhvt2
1155  enddo
1156  elseif (cs%LT_VT2_METHOD==lt_vt2_mode_rw16) then
1157  !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
1158  enhvt2 = 1. + 2.3*cs%La_SL(i,j)**(-0.5)
1159  do k=1,g%ke
1160  langenhvt2(k) = enhvt2
1161  enddo
1162  elseif (cs%LT_VT2_METHOD==lt_vt2_mode_lf17) then
1163  cs%CS=cvmix_get_kpp_real('c_s',cs%KPP_params)
1164  do k=1,g%ke
1165  wst = (max(0.,-buoy_scale*buoyflux(i,j,1))*(-cellheight(k)))**(1./3.)
1166  langenhvt2(k) = sqrt((0.15*wst**3. + 0.17*surffricvel**3.* &
1167  (1.+0.49*cs%La_SL(i,j)**(-2.))) / &
1168  (0.2*ws_1d(k)**3/(cs%cs*cs%surf_layer_ext*cs%vonKarman**4.)))
1169  enddo
1170  else
1171  !This shouldn't be reached.
1172  !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in Vt2")
1173  langenhvt2(:) = 1.0
1174  endif
1175  else
1176  langenhvt2(:) = 1.0
1177  endif
1178 
1179  do k=1,g%ke
1180  cs%Vt2(i,j,k)=cs%Vt2(i,j,k)*langenhvt2(k)
1181  if (cs%id_EnhVt2 > 0) cs%EnhVt2(i,j,k)=langenhvt2(k)
1182  enddo
1183 
1184  ! Calculate Bulk Richardson number from eq (21) of LMD94
1185  bulkri_1d = cvmix_kpp_compute_bulk_richardson( &
1186  zt_cntr = cellheight(1:g%ke), & ! Depth of cell center [m]
1187  delta_buoy_cntr=gorho*deltarho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
1188  delta_vsqr_cntr=deltau2, & ! Square of resolved velocity difference [m2 s-2]
1189  vt_sqr_cntr=cs%Vt2(i,j,:), &
1190  ws_cntr=ws_1d, & ! Turbulent velocity scale profile [m s-1]
1191  n_iface=cs%N(i,j,:)) ! Buoyancy frequency [s-1]
1192 
1193 
1194  surfbuoyflux = buoy_scale * buoyflux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
1195  ! h to Monin-Obukov (default is false, ie. not used)
1196 
1197  call cvmix_kpp_compute_obl_depth( &
1198  bulkri_1d, & ! (in) Bulk Richardson number
1199  ifaceheight, & ! (in) Height of interfaces [m]
1200  cs%OBLdepth(i,j), & ! (out) OBL depth [m]
1201  cs%kOBL(i,j), & ! (out) level (+fraction) of OBL extent
1202  zt_cntr=cellheight, & ! (in) Height of cell centers [m]
1203  surf_fric=surffricvel, & ! (in) Turbulent friction velocity at surface [m s-1]
1204  surf_buoy=surfbuoyflux, & ! (in) Buoyancy flux at surface [m2 s-3]
1205  coriolis=coriolis, & ! (in) Coriolis parameter [s-1]
1206  cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
1207 
1208  ! A hack to avoid KPP reaching the bottom. It was needed during development
1209  ! because KPP was unable to handle vanishingly small layers near the bottom.
1210  if (cs%deepOBLoffset>0.) then
1211  zbottomminusoffset = ifaceheight(g%ke+1) + min(cs%deepOBLoffset,-0.1*ifaceheight(g%ke+1))
1212  cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -zbottomminusoffset )
1213  endif
1214 
1215  ! apply some constraints on OBLdepth
1216  if(cs%fixedOBLdepth) cs%OBLdepth(i,j) = cs%fixedOBLdepth_value
1217  cs%OBLdepth(i,j) = max( cs%OBLdepth(i,j), -ifaceheight(2) ) ! no shallower than top layer
1218  cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -ifaceheight(g%ke+1) ) ! no deeper than bottom
1219  cs%kOBL(i,j) = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, cs%OBLdepth(i,j) )
1220 
1221 !*************************************************************************
1222 ! smg: remove code below
1223 
1224 ! Following "correction" step has been found to be unnecessary.
1225 ! Code should be removed after further testing.
1226 ! BGR: 03/15/2018-> Restructured code (Vt2 changed to compute from call in MOM_CVMix_KPP now)
1227 ! I have not taken this restructuring into account here.
1228 ! Do we ever run with correctSurfLayerAvg?
1229 ! smg's suggested testing and removal is advised, in the meantime
1230 ! I have added warning if correctSurfLayerAvg is attempted.
1231  ! if (CS%correctSurfLayerAvg) then
1232 
1233  ! SLdepth_0d = CS%surf_layer_ext * CS%OBLdepth(i,j)
1234  ! hTot = h(i,j,1)
1235  ! surfTemp = Temp(i,j,1) ; surfHtemp = surfTemp * hTot
1236  ! surfSalt = Salt(i,j,1) ; surfHsalt = surfSalt * hTot
1237  ! surfU = 0.5*US%L_T_to_m_s*(u(i,j,1)+u(i-1,j,1)) ; surfHu = surfU * hTot
1238  ! surfV = 0.5*US%L_T_to_m_s*(v(i,j,1)+v(i,j-1,1)) ; surfHv = surfV * hTot
1239  ! pRef = 0.0
1240 
1241  ! do k = 2, G%ke
1242 
1243  ! ! Recalculate differences with surface layer
1244  ! Uk = 0.5*US%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) - surfU
1245  ! Vk = 0.5*US%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) - surfV
1246  ! deltaU2(k) = Uk**2 + Vk**2
1247  ! pRef = pRef + GV%H_to_Pa * h(i,j,k)
1248  ! call calculate_density(surfTemp, surfSalt, pRef, rho1, EOS)
1249  ! call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS)
1250  ! deltaRho(k) = rhoK - rho1
1251 
1252  ! ! Surface layer averaging (needed for next k+1 iteration of this loop)
1253  ! if (hTot < SLdepth_0d) then
1254  ! delH = min( max(0., SLdepth_0d - hTot), h(i,j,k)*GV%H_to_m )
1255  ! hTot = hTot + delH
1256  ! surfHtemp = surfHtemp + Temp(i,j,k) * delH ; surfTemp = surfHtemp / hTot
1257  ! surfHsalt = surfHsalt + Salt(i,j,k) * delH ; surfSalt = surfHsalt / hTot
1258  ! surfHu = surfHu + 0.5*US%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) * delH ; surfU = surfHu / hTot
1259  ! surfHv = surfHv + 0.5*US%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) * delH ; surfV = surfHv / hTot
1260  ! endif
1261 
1262  ! enddo
1263 
1264  ! BulkRi_1d = CVMix_kpp_compute_bulk_Richardson( &
1265  ! cellHeight(1:G%ke), & ! Depth of cell center [m]
1266  ! GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
1267  ! deltaU2, & ! Square of resolved velocity difference [m2 s-2]
1268  ! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1]
1269  ! N_iface=CS%N ) ! Buoyancy frequency [s-1]
1270 
1271  ! surfBuoyFlux = buoy_scale*buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
1272  ! ! h to Monin-Obukov (default is false, ie. not used)
1273 
1274  ! call CVMix_kpp_compute_OBL_depth( &
1275  ! BulkRi_1d, & ! (in) Bulk Richardson number
1276  ! iFaceHeight, & ! (in) Height of interfaces [m]
1277  ! CS%OBLdepth(i,j), & ! (out) OBL depth [m]
1278  ! CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent
1279  ! zt_cntr=cellHeight, & ! (in) Height of cell centers [m]
1280  ! surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
1281  ! surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3]
1282  ! Coriolis=Coriolis, & ! (in) Coriolis parameter [s-1]
1283  ! CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters
1284 
1285  ! if (CS%deepOBLoffset>0.) then
1286  ! zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1))
1287  ! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset )
1288  ! CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
1289  ! endif
1290 
1291  ! ! apply some constraints on OBLdepth
1292  ! if(CS%fixedOBLdepth) CS%OBLdepth(i,j) = CS%fixedOBLdepth_value
1293  ! CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer
1294  ! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deep than bottom
1295  ! CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
1296 
1297  ! endif ! endif for "correction" step
1298 
1299 ! smg: remove code above
1300 ! **********************************************************************
1301 
1302  ! recompute wscale for diagnostics, now that we in fact know boundary layer depth
1303  !BGR consider if LTEnhancement is wanted for diagnostics
1304  if (cs%id_Ws > 0) then
1305  call cvmix_kpp_compute_turbulent_scales( &
1306  -cellheight/cs%OBLdepth(i,j), & ! (in) Normalized boundary layer coordinate
1307  cs%OBLdepth(i,j), & ! (in) OBL depth [m]
1308  surfbuoyflux, & ! (in) Buoyancy flux at surface [m2 s-3]
1309  surffricvel, & ! (in) Turbulent friction velocity at surface [m s-1]
1310  w_s=ws_1d, & ! (out) Turbulent velocity scale profile [m s-1]
1311  cvmix_kpp_params_user=cs%KPP_params) ! KPP parameters
1312  cs%Ws(i,j,:) = ws_1d(:)
1313  endif
1314 
1315  ! Diagnostics
1316  if (cs%id_N2 > 0) cs%N2(i,j,:) = n2_1d(:)
1317  if (cs%id_BulkDrho > 0) cs%dRho(i,j,:) = deltarho(:)
1318  if (cs%id_BulkRi > 0) cs%BulkRi(i,j,:) = bulkri_1d(:)
1319  if (cs%id_BulkUz2 > 0) cs%Uz2(i,j,:) = deltau2(:)
1320  if (cs%id_Tsurf > 0) cs%Tsurf(i,j) = surftemp
1321  if (cs%id_Ssurf > 0) cs%Ssurf(i,j) = surfsalt
1322  if (cs%id_Usurf > 0) cs%Usurf(i,j) = surfu
1323  if (cs%id_Vsurf > 0) cs%Vsurf(i,j) = surfv
1324 
1325  enddo
1326  enddo
1327 
1328  ! send diagnostics to post_data
1329  if (cs%id_BulkRi > 0) call post_data(cs%id_BulkRi, cs%BulkRi, cs%diag)
1330  if (cs%id_N > 0) call post_data(cs%id_N, cs%N, cs%diag)
1331  if (cs%id_N2 > 0) call post_data(cs%id_N2, cs%N2, cs%diag)
1332  if (cs%id_Tsurf > 0) call post_data(cs%id_Tsurf, cs%Tsurf, cs%diag)
1333  if (cs%id_Ssurf > 0) call post_data(cs%id_Ssurf, cs%Ssurf, cs%diag)
1334  if (cs%id_Usurf > 0) call post_data(cs%id_Usurf, cs%Usurf, cs%diag)
1335  if (cs%id_Vsurf > 0) call post_data(cs%id_Vsurf, cs%Vsurf, cs%diag)
1336  if (cs%id_BulkDrho > 0) call post_data(cs%id_BulkDrho, cs%dRho, cs%diag)
1337  if (cs%id_BulkUz2 > 0) call post_data(cs%id_BulkUz2, cs%Uz2, cs%diag)
1338  if (cs%id_EnhK > 0) call post_data(cs%id_EnhK, cs%EnhK, cs%diag)
1339  if (cs%id_EnhVt2 > 0) call post_data(cs%id_EnhVt2, cs%EnhVt2, cs%diag)
1340  if (cs%id_La_SL > 0) call post_data(cs%id_La_SL, cs%La_SL, cs%diag)
1341 
1342  ! BLD smoothing:
1343  if (cs%n_smooth > 0) call kpp_smooth_bld(cs,g,gv,h)
1344 

References mom_wave_interface::get_langmuir_number(), kpp_smooth_bld(), lt_vt2_mode_constant, lt_vt2_mode_lf17, lt_vt2_mode_rw16, and lt_vt2_mode_vr12.

Here is the call graph for this function:

◆ kpp_end()

subroutine, public mom_cvmix_kpp::kpp_end ( type(kpp_cs), pointer  CS)

Clear pointers, deallocate memory.

Parameters
csControl structure

Definition at line 1603 of file MOM_CVMix_KPP.F90.

1603  type(KPP_CS), pointer :: CS !< Control structure
1604 
1605  if (.not.associated(cs)) return
1606 
1607  deallocate(cs)
1608 

◆ kpp_get_bld()

subroutine, public mom_cvmix_kpp::kpp_get_bld ( type(kpp_cs), pointer  CS,
real, dimension(szi_(g),szj_(g)), intent(inout)  BLD,
type(ocean_grid_type), intent(in)  G 
)

Copies KPP surface boundary layer depth into BLD.

Parameters
csControl structure for this module
[in]gGrid structure
[in,out]bldbnd. layer depth [m]

Definition at line 1469 of file MOM_CVMix_KPP.F90.

1469  type(KPP_CS), pointer :: CS !< Control structure for
1470  !! this module
1471  type(ocean_grid_type), intent(in) :: G !< Grid structure
1472  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD!< bnd. layer depth [m]
1473  ! Local variables
1474  integer :: i,j
1475  do j = g%jsc, g%jec ; do i = g%isc, g%iec
1476  bld(i,j) = cs%OBLdepth(i,j)
1477  enddo ; enddo

Referenced by mom_lateral_boundary_diffusion::lateral_boundary_diffusion(), and mom_neutral_diffusion::neutral_diffusion_calc_coeffs().

Here is the caller graph for this function:

◆ kpp_init()

logical function, public mom_cvmix_kpp::kpp_init ( type(param_file_type), intent(in)  paramFile,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(diag_ctrl), intent(in), target  diag,
type(time_type), intent(in)  Time,
type(kpp_cs), pointer  CS,
logical, intent(out), optional  passive,
type(wave_parameters_cs), optional, pointer  Waves 
)

Initialize the CVMix KPP module and set up diagnostics Returns True if KPP is to be used, False otherwise.

Parameters
[in]paramfileFile parser
[in]gOcean grid
[in]gvVertical grid structure.
[in]usA dimensional unit scaling type
[in]diagDiagnostics
[in]timeModel time
csControl structure
[out]passiveCopy of passiveMode
wavesWave CS

Definition at line 181 of file MOM_CVMix_KPP.F90.

181 
182  ! Arguments
183  type(param_file_type), intent(in) :: paramFile !< File parser
184  type(ocean_grid_type), intent(in) :: G !< Ocean grid
185  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
186  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
187  type(diag_ctrl), target, intent(in) :: diag !< Diagnostics
188  type(time_type), intent(in) :: Time !< Model time
189  type(KPP_CS), pointer :: CS !< Control structure
190  logical, optional, intent(out) :: passive !< Copy of %passiveMode
191  type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS
192 
193  ! Local variables
194 #include "version_variable.h"
195  character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module
196  character(len=20) :: string !< local temporary string
197  logical :: CS_IS_ONE=.false. !< Logical for setting Cs based on Non-local
198  logical :: lnoDGat1=.false. !< True => G'(1) = 0 (shape function)
199  !! False => compute G'(1) as in LMD94
200  if (associated(cs)) call mom_error(fatal, 'MOM_CVMix_KPP, KPP_init: '// &
201  'Control structure has already been initialized')
202 
203  ! Read parameters
204  call log_version(paramfile, mdl, version, 'This is the MOM wrapper to CVMix:KPP\n' // &
205  'See http://cvmix.github.io/')
206  call get_param(paramfile, mdl, "USE_KPP", kpp_init, &
207  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "// &
208  "to calculate diffusivities and non-local transport in the OBL.", &
209  default=.false.)
210  ! Forego remainder of initialization if not using this scheme
211  if (.not. kpp_init) return
212  allocate(cs)
213 
214  call openparameterblock(paramfile,'KPP')
215  call get_param(paramfile, mdl, 'PASSIVE', cs%passiveMode, &
216  'If True, puts KPP into a passive-diagnostic mode.', &
217  default=.false.)
218  !BGR: Note using PASSIVE for KPP creates warning for PASSIVE from Convection
219  ! should we create a separate flag?
220  if (present(passive)) passive=cs%passiveMode ! This is passed back to the caller so
221  ! the caller knows to not use KPP output
222  call get_param(paramfile, mdl, 'APPLY_NONLOCAL_TRANSPORT', cs%applyNonLocalTrans, &
223  'If True, applies the non-local transport to heat and scalars. '// &
224  'If False, calculates the non-local transport and tendencies but '//&
225  'purely for diagnostic purposes.', &
226  default=.not. cs%passiveMode)
227  call get_param(paramfile, mdl, 'N_SMOOTH', cs%n_smooth, &
228  'The number of times the 1-1-4-1-1 Laplacian filter is applied on '// &
229  'OBL depth.', &
230  default=0)
231  if (cs%n_smooth > 0) then
232  call get_param(paramfile, mdl, 'DEEPEN_ONLY_VIA_SMOOTHING', cs%deepen_only, &
233  'If true, apply OBLdepth smoothing at a cell only if the OBLdepth '// &
234  'gets deeper via smoothing.', &
235  default=.false.)
236  endif
237  call get_param(paramfile, mdl, 'RI_CRIT', cs%Ri_crit, &
238  'Critical bulk Richardson number used to define depth of the '// &
239  'surface Ocean Boundary Layer (OBL).', &
240  units='nondim', default=0.3)
241  call get_param(paramfile, mdl, 'VON_KARMAN', cs%vonKarman, &
242  'von Karman constant.', &
243  units='nondim', default=0.40)
244  call get_param(paramfile, mdl, 'ENHANCE_DIFFUSION', cs%enhance_diffusion, &
245  'If True, adds enhanced diffusion at the based of the boundary layer.', &
246  default=.true.)
247  call get_param(paramfile, mdl, 'INTERP_TYPE', cs%interpType, &
248  'Type of interpolation to determine the OBL depth.\n'// &
249  'Allowed types are: linear, quadratic, cubic.', &
250  default='quadratic')
251  call get_param(paramfile, mdl, 'INTERP_TYPE2', cs%interpType2, &
252  'Type of interpolation to compute diff and visc at OBL_depth.\n'// &
253  'Allowed types are: linear, quadratic, cubic or LMD94.', &
254  default='LMD94')
255  call get_param(paramfile, mdl, 'COMPUTE_EKMAN', cs%computeEkman, &
256  'If True, limit OBL depth to be no deeper than Ekman depth.', &
257  default=.false.)
258  call get_param(paramfile, mdl, 'COMPUTE_MONIN_OBUKHOV', cs%computeMoninObukhov, &
259  'If True, limit the OBL depth to be no deeper than '// &
260  'Monin-Obukhov depth.', &
261  default=.false.)
262  call get_param(paramfile, mdl, 'CS', cs%cs, &
263  'Parameter for computing velocity scale function.', &
264  units='nondim', default=98.96)
265  call get_param(paramfile, mdl, 'CS2', cs%cs2, &
266  'Parameter for computing non-local term.', &
267  units='nondim', default=6.32739901508)
268  call get_param(paramfile, mdl, 'DEEP_OBL_OFFSET', cs%deepOBLoffset, &
269  'If non-zero, the distance above the bottom to which the OBL is clipped '// &
270  'if it would otherwise reach the bottom. The smaller of this and 0.1D is used.', &
271  units='m',default=0.)
272  call get_param(paramfile, mdl, 'FIXED_OBLDEPTH', cs%fixedOBLdepth, &
273  'If True, fix the OBL depth to FIXED_OBLDEPTH_VALUE '// &
274  'rather than using the OBL depth from CVMix. '// &
275  'This option is just for testing purposes.', &
276  default=.false.)
277  call get_param(paramfile, mdl, 'FIXED_OBLDEPTH_VALUE', cs%fixedOBLdepth_value, &
278  'Value for the fixed OBL depth when fixedOBLdepth==True. '// &
279  'This parameter is for just for testing purposes. '// &
280  'It will over-ride the OBLdepth computed from CVMix.', &
281  units='m',default=30.0)
282  call get_param(paramfile, mdl, 'SURF_LAYER_EXTENT', cs%surf_layer_ext, &
283  'Fraction of OBL depth considered in the surface layer.', &
284  units='nondim',default=0.10)
285  call get_param(paramfile, mdl, 'MINIMUM_OBL_DEPTH', cs%minOBLdepth, &
286  'If non-zero, a minimum depth to use for KPP OBL depth. Independent of '// &
287  'this parameter, the OBL depth is always at least as deep as the first layer.', &
288  units='m',default=0.)
289  call get_param(paramfile, mdl, 'MINIMUM_VT2', cs%minVtsqr, &
290  'Min of the unresolved velocity Vt2 used in Rib CVMix calculation.\n'// &
291  'Scaling: MINIMUM_VT2 = const1*d*N*ws, with d=1m, N=1e-5/s, ws=1e-6 m/s.', &
292  units='m2/s2',default=1e-10)
293 
294 ! smg: for removal below
295  call get_param(paramfile, mdl, 'CORRECT_SURFACE_LAYER_AVERAGE', cs%correctSurfLayerAvg, &
296  'If true, applies a correction step to the averaging of surface layer '// &
297  'properties. This option is obsolete.', default=.false.)
298  if (cs%correctSurfLayerAvg) &
299  call mom_error(fatal,'Correct surface layer average disabled in code. To recover \n'// &
300  ' feature will require code intervention.')
301  call get_param(paramfile, mdl, 'FIRST_GUESS_SURFACE_LAYER_DEPTH', cs%surfLayerDepth, &
302  'The first guess at the depth of the surface layer used for averaging '// &
303  'the surface layer properties. If =0, the top model level properties '// &
304  'will be used for the surface layer. If CORRECT_SURFACE_LAYER_AVERAGE=True, a '// &
305  'subsequent correction is applied. This parameter is obsolete', units='m', default=0.)
306 ! smg: for removal above
307 
308  call get_param(paramfile, mdl, 'NLT_SHAPE', string, &
309  'MOM6 method to set nonlocal transport profile. '// &
310  'Over-rides the result from CVMix. Allowed values are: \n'// &
311  '\t CVMix - Uses the profiles from CVMix specified by MATCH_TECHNIQUE\n'//&
312  '\t LINEAR - A linear profile, 1-sigma\n'// &
313  '\t PARABOLIC - A parablic profile, (1-sigma)^2\n'// &
314  '\t CUBIC - A cubic profile, (1-sigma)^2(1+2*sigma)\n'// &
315  '\t CUBIC_LMD - The original KPP profile', &
316  default='CVMix')
317  select case ( trim(string) )
318  case ("CVMix") ; cs%NLT_shape = nlt_shape_cvmix
319  case ("LINEAR") ; cs%NLT_shape = nlt_shape_linear
320  case ("PARABOLIC") ; cs%NLT_shape = nlt_shape_parabolic
321  case ("CUBIC") ; cs%NLT_shape = nlt_shape_cubic
322  case ("CUBIC_LMD") ; cs%NLT_shape = nlt_shape_cubic_lmd
323  case default ; call mom_error(fatal,"KPP_init: "// &
324  "Unrecognized NLT_SHAPE option"//trim(string))
325  end select
326  call get_param(paramfile, mdl, 'MATCH_TECHNIQUE', cs%MatchTechnique, &
327  'CVMix method to set profile function for diffusivity and NLT, '// &
328  'as well as matching across OBL base. Allowed values are: \n'// &
329  '\t SimpleShapes = sigma*(1-sigma)^2 for both diffusivity and NLT\n'// &
330  '\t MatchGradient = sigma*(1-sigma)^2 for NLT; diffusivity profile from matching\n'//&
331  '\t MatchBoth = match gradient for both diffusivity and NLT\n'// &
332  '\t ParabolicNonLocal = sigma*(1-sigma)^2 for diffusivity; (1-sigma)^2 for NLT', &
333  default='SimpleShapes')
334  if (cs%MatchTechnique == 'ParabolicNonLocal') then
335  ! This forces Cs2 (Cs in non-local computation) to equal 1 for parabolic non-local option.
336  ! May be used during CVMix initialization.
337  cs_is_one=.true.
338  endif
339  if (cs%MatchTechnique == 'ParabolicNonLocal' .or. cs%MatchTechnique == 'SimpleShapes') then
340  ! if gradient won't be matched, lnoDGat1=.true.
341  lnodgat1=.true.
342  endif
343 
344  ! safety check to avoid negative diff/visc
345  if (cs%MatchTechnique == 'MatchBoth' .and. (cs%interpType2 == 'cubic' .or. &
346  cs%interpType2 == 'quadratic')) then
347  call mom_error(fatal,"If MATCH_TECHNIQUE=MatchBoth, INTERP_TYPE2 must be set to \n"//&
348  "linear or LMD94 (recommended) to avoid negative viscosity and diffusivity.\n"//&
349  "Please select one of these valid options." )
350  endif
351 
352  call get_param(paramfile, mdl, 'KPP_ZERO_DIFFUSIVITY', cs%KPPzeroDiffusivity, &
353  'If True, zeroes the KPP diffusivity and viscosity; for testing purpose.',&
354  default=.false.)
355  call get_param(paramfile, mdl, 'KPP_IS_ADDITIVE', cs%KPPisAdditive, &
356  'If true, adds KPP diffusivity to diffusivity from other schemes.\n'//&
357  'If false, KPP is the only diffusivity wherever KPP is non-zero.', &
358  default=.true.)
359  call get_param(paramfile, mdl, 'KPP_SHORTWAVE_METHOD',string, &
360  'Determines contribution of shortwave radiation to KPP surface '// &
361  'buoyancy flux. Options include:\n'// &
362  ' ALL_SW: use total shortwave radiation\n'// &
363  ' MXL_SW: use shortwave radiation absorbed by mixing layer\n'// &
364  ' LV1_SW: use shortwave radiation absorbed by top model layer', &
365  default='MXL_SW')
366  select case ( trim(string) )
367  case ("ALL_SW") ; cs%SW_METHOD = sw_method_all_sw
368  case ("MXL_SW") ; cs%SW_METHOD = sw_method_mxl_sw
369  case ("LV1_SW") ; cs%SW_METHOD = sw_method_lv1_sw
370  case default ; call mom_error(fatal,"KPP_init: "// &
371  "Unrecognized KPP_SHORTWAVE_METHOD option"//trim(string))
372  end select
373  call get_param(paramfile, mdl, 'CVMix_ZERO_H_WORK_AROUND', cs%min_thickness, &
374  'A minimum thickness used to avoid division by small numbers in the vicinity '// &
375  'of vanished layers. This is independent of MIN_THICKNESS used in other parts of MOM.', &
376  units='m', default=0.)
377 
378 !/BGR: New options for including Langmuir effects
379 !/ 1. Options related to enhancing the mixing coefficient
380  call get_param(paramfile, mdl, "USE_KPP_LT_K", cs%LT_K_Enhancement, &
381  'Flag for Langmuir turbulence enhancement of turbulent'//&
382  'mixing coefficient.', units="", default=.false.)
383  call get_param(paramfile, mdl, "STOKES_MIXING", cs%STOKES_MIXING, &
384  'Flag for Langmuir turbulence enhancement of turbulent'//&
385  'mixing coefficient.', units="", default=.false.)
386  if (cs%LT_K_Enhancement) then
387  call get_param(paramfile, mdl, 'KPP_LT_K_SHAPE', string, &
388  'Vertical dependence of LT enhancement of mixing. '// &
389  'Valid options are: \n'// &
390  '\t CONSTANT = Constant value for full OBL\n'// &
391  '\t SCALED = Varies based on normalized shape function.', &
392  default='CONSTANT')
393  select case ( trim(string))
394  case ("CONSTANT") ; cs%LT_K_SHAPE = lt_k_constant
395  case ("SCALED") ; cs%LT_K_SHAPE = lt_k_scaled
396  case default ; call mom_error(fatal,"KPP_init: "//&
397  "Unrecognized KPP_LT_K_SHAPE option: "//trim(string))
398  end select
399  call get_param(paramfile, mdl, "KPP_LT_K_METHOD", string , &
400  'Method to enhance mixing coefficient in KPP. '// &
401  'Valid options are: \n'// &
402  '\t CONSTANT = Constant value (KPP_K_ENH_FAC) \n'// &
403  '\t VR12 = Function of Langmuir number based on VR12\n'// &
404  '\t RW16 = Function of Langmuir number based on RW16', &
405  default='CONSTANT')
406  select case ( trim(string))
407  case ("CONSTANT") ; cs%LT_K_METHOD = lt_k_mode_constant
408  case ("VR12") ; cs%LT_K_METHOD = lt_k_mode_vr12
409  case ("RW16") ; cs%LT_K_METHOD = lt_k_mode_rw16
410  case default ; call mom_error(fatal,"KPP_init: "//&
411  "Unrecognized KPP_LT_K_METHOD option: "//trim(string))
412  end select
413  if (cs%LT_K_METHOD==lt_k_mode_constant) then
414  call get_param(paramfile, mdl, "KPP_K_ENH_FAC",cs%KPP_K_ENH_FAC , &
415  'Constant value to enhance mixing coefficient in KPP.', &
416  default=1.0)
417  endif
418  endif
419 !/ 2. Options related to enhancing the unresolved Vt2/entrainment in Rib
420  call get_param(paramfile, mdl, "USE_KPP_LT_VT2", cs%LT_Vt2_Enhancement, &
421  'Flag for Langmuir turbulence enhancement of Vt2'//&
422  'in Bulk Richardson Number.', units="", default=.false.)
423  if (cs%LT_Vt2_Enhancement) then
424  call get_param(paramfile, mdl, "KPP_LT_VT2_METHOD",string , &
425  'Method to enhance Vt2 in KPP. '// &
426  'Valid options are: \n'// &
427  '\t CONSTANT = Constant value (KPP_VT2_ENH_FAC) \n'// &
428  '\t VR12 = Function of Langmuir number based on VR12\n'// &
429  '\t RW16 = Function of Langmuir number based on RW16\n'// &
430  '\t LF17 = Function of Langmuir number based on LF17', &
431  default='CONSTANT')
432  select case ( trim(string))
433  case ("CONSTANT") ; cs%LT_VT2_METHOD = lt_vt2_mode_constant
434  case ("VR12") ; cs%LT_VT2_METHOD = lt_vt2_mode_vr12
435  case ("RW16") ; cs%LT_VT2_METHOD = lt_vt2_mode_rw16
436  case ("LF17") ; cs%LT_VT2_METHOD = lt_vt2_mode_lf17
437  case default ; call mom_error(fatal,"KPP_init: "//&
438  "Unrecognized KPP_LT_VT2_METHOD option: "//trim(string))
439  end select
440  if (cs%LT_VT2_METHOD==lt_vt2_mode_constant) then
441  call get_param(paramfile, mdl, "KPP_VT2_ENH_FAC",cs%KPP_VT2_ENH_FAC , &
442  'Constant value to enhance VT2 in KPP.', &
443  default=1.0)
444  endif
445  endif
446 
447  call closeparameterblock(paramfile)
448  call get_param(paramfile, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)
449 
450  call cvmix_init_kpp( ri_crit=cs%Ri_crit, &
451  minobldepth=cs%minOBLdepth, &
452  minvtsqr=cs%minVtsqr, &
453  vonkarman=cs%vonKarman, &
454  surf_layer_ext=cs%surf_layer_ext, &
455  interp_type=cs%interpType, &
456  interp_type2=cs%interpType2, &
457  lekman=cs%computeEkman, &
458  lmonob=cs%computeMoninObukhov, &
459  matchtechnique=cs%MatchTechnique, &
460  lenhanced_diff=cs%enhance_diffusion,&
461  lnonzero_surf_nonlocal=cs_is_one ,&
462  lnodgat1=lnodgat1 ,&
463  cvmix_kpp_params_user=cs%KPP_params )
464 
465  ! Register diagnostics
466  cs%diag => diag
467  cs%id_OBLdepth = register_diag_field('ocean_model', 'KPP_OBLdepth', diag%axesT1, time, &
468  'Thickness of the surface Ocean Boundary Layer calculated by [CVMix] KPP', 'meter', &
469  cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
470  cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
471  ! CMOR names are placeholders; must be modified by time period
472  ! for CMOR compliance. Diag manager will be used for omlmax and
473  ! omldamax.
474  if (cs%n_smooth > 0) then
475  cs%id_OBLdepth_original = register_diag_field('ocean_model', 'KPP_OBLdepth_original', diag%axesT1, time, &
476  'Thickness of the surface Ocean Boundary Layer without smoothing calculated by [CVMix] KPP', 'meter', &
477  cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
478  cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
479  endif
480  cs%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, time, &
481  'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', 'kg/m3')
482  cs%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, time, &
483  'Square of bulk difference in resolved velocity used in Bulk Richardson number via [CVMix] KPP', 'm2/s2')
484  cs%id_BulkRi = register_diag_field('ocean_model', 'KPP_BulkRi', diag%axesTL, time, &
485  'Bulk Richardson number used to find the OBL depth used by [CVMix] KPP', 'nondim')
486  cs%id_Sigma = register_diag_field('ocean_model', 'KPP_sigma', diag%axesTi, time, &
487  'Sigma coordinate used by [CVMix] KPP', 'nondim')
488  cs%id_Ws = register_diag_field('ocean_model', 'KPP_Ws', diag%axesTL, time, &
489  'Turbulent vertical velocity scale for scalars used by [CVMix] KPP', 'm/s')
490  cs%id_N = register_diag_field('ocean_model', 'KPP_N', diag%axesTi, time, &
491  '(Adjusted) Brunt-Vaisala frequency used by [CVMix] KPP', '1/s')
492  cs%id_N2 = register_diag_field('ocean_model', 'KPP_N2', diag%axesTi, time, &
493  'Square of Brunt-Vaisala frequency used by [CVMix] KPP', '1/s2')
494  cs%id_Vt2 = register_diag_field('ocean_model', 'KPP_Vt2', diag%axesTL, time, &
495  'Unresolved shear turbulence used by [CVMix] KPP', 'm2/s2')
496  cs%id_uStar = register_diag_field('ocean_model', 'KPP_uStar', diag%axesT1, time, &
497  'Friction velocity, u*, as used by [CVMix] KPP', 'm/s', conversion=us%Z_to_m*us%s_to_T)
498  cs%id_buoyFlux = register_diag_field('ocean_model', 'KPP_buoyFlux', diag%axesTi, time, &
499  'Surface (and penetrating) buoyancy flux, as used by [CVMix] KPP', 'm2/s3', conversion=us%L_to_m**2*us%s_to_T**3)
500  cs%id_QminusSW = register_diag_field('ocean_model', 'KPP_QminusSW', diag%axesT1, time, &
501  'Net temperature flux ignoring short-wave, as used by [CVMix] KPP', 'K m/s')
502  cs%id_netS = register_diag_field('ocean_model', 'KPP_netSalt', diag%axesT1, time, &
503  'Effective net surface salt flux, as used by [CVMix] KPP', 'ppt m/s')
504  cs%id_Kt_KPP = register_diag_field('ocean_model', 'KPP_Kheat', diag%axesTi, time, &
505  'Heat diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
506  cs%id_Kd_in = register_diag_field('ocean_model', 'KPP_Kd_in', diag%axesTi, time, &
507  'Diffusivity passed to KPP', 'm2/s', conversion=us%Z2_T_to_m2_s)
508  cs%id_Ks_KPP = register_diag_field('ocean_model', 'KPP_Ksalt', diag%axesTi, time, &
509  'Salt diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
510  cs%id_Kv_KPP = register_diag_field('ocean_model', 'KPP_Kv', diag%axesTi, time, &
511  'Vertical viscosity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
512  cs%id_NLTt = register_diag_field('ocean_model', 'KPP_NLtransport_heat', diag%axesTi, time, &
513  'Non-local transport (Cs*G(sigma)) for heat, as calculated by [CVMix] KPP', 'nondim')
514  cs%id_NLTs = register_diag_field('ocean_model', 'KPP_NLtransport_salt', diag%axesTi, time, &
515  'Non-local tranpsort (Cs*G(sigma)) for scalars, as calculated by [CVMix] KPP', 'nondim')
516  cs%id_NLT_dTdt = register_diag_field('ocean_model', 'KPP_NLT_dTdt', diag%axesTL, time, &
517  'Temperature tendency due to non-local transport of heat, as calculated by [CVMix] KPP', 'K/s')
518  cs%id_NLT_dSdt = register_diag_field('ocean_model', 'KPP_NLT_dSdt', diag%axesTL, time, &
519  'Salinity tendency due to non-local transport of salt, as calculated by [CVMix] KPP', 'ppt/s')
520  cs%id_NLT_temp_budget = register_diag_field('ocean_model', 'KPP_NLT_temp_budget', diag%axesTL, time, &
521  'Heat content change due to non-local transport, as calculated by [CVMix] KPP', 'W/m^2')
522  cs%id_NLT_saln_budget = register_diag_field('ocean_model', 'KPP_NLT_saln_budget', diag%axesTL, time, &
523  'Salt content change due to non-local transport, as calculated by [CVMix] KPP', 'kg/(sec*m^2)')
524  cs%id_Tsurf = register_diag_field('ocean_model', 'KPP_Tsurf', diag%axesT1, time, &
525  'Temperature of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'C')
526  cs%id_Ssurf = register_diag_field('ocean_model', 'KPP_Ssurf', diag%axesT1, time, &
527  'Salinity of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'ppt')
528  cs%id_Usurf = register_diag_field('ocean_model', 'KPP_Usurf', diag%axesCu1, time, &
529  'i-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'm/s')
530  cs%id_Vsurf = register_diag_field('ocean_model', 'KPP_Vsurf', diag%axesCv1, time, &
531  'j-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'm/s')
532  cs%id_EnhK = register_diag_field('ocean_model', 'EnhK', diag%axesTI, time, &
533  'Langmuir number enhancement to K as used by [CVMix] KPP','nondim')
534  cs%id_EnhVt2 = register_diag_field('ocean_model', 'EnhVt2', diag%axesTL, time, &
535  'Langmuir number enhancement to Vt2 as used by [CVMix] KPP','nondim')
536  cs%id_La_SL = register_diag_field('ocean_model', 'KPP_La_SL', diag%axesT1, time, &
537  'Surface-layer Langmuir number computed in [CVMix] KPP','nondim')
538 
539  allocate( cs%N( szi_(g), szj_(g), szk_(g)+1 ) )
540  cs%N(:,:,:) = 0.
541  allocate( cs%OBLdepth( szi_(g), szj_(g) ) )
542  cs%OBLdepth(:,:) = 0.
543  allocate( cs%kOBL( szi_(g), szj_(g) ) )
544  cs%kOBL(:,:) = 0.
545  allocate( cs%La_SL( szi_(g), szj_(g) ) )
546  cs%La_SL(:,:) = 0.
547  allocate( cs%Vt2( szi_(g), szj_(g), szk_(g) ) )
548  cs%Vt2(:,:,:) = 0.
549  if (cs%id_OBLdepth_original > 0) allocate( cs%OBLdepth_original( szi_(g), szj_(g) ) )
550 
551  allocate( cs%OBLdepthprev( szi_(g), szj_(g) ) ) ; cs%OBLdepthprev(:,:) = 0.0
552  if (cs%id_BulkDrho > 0) allocate( cs%dRho( szi_(g), szj_(g), szk_(g) ) )
553  if (cs%id_BulkDrho > 0) cs%dRho(:,:,:) = 0.
554  if (cs%id_BulkUz2 > 0) allocate( cs%Uz2( szi_(g), szj_(g), szk_(g) ) )
555  if (cs%id_BulkUz2 > 0) cs%Uz2(:,:,:) = 0.
556  if (cs%id_BulkRi > 0) allocate( cs%BulkRi( szi_(g), szj_(g), szk_(g) ) )
557  if (cs%id_BulkRi > 0) cs%BulkRi(:,:,:) = 0.
558  if (cs%id_Sigma > 0) allocate( cs%sigma( szi_(g), szj_(g), szk_(g)+1 ) )
559  if (cs%id_Sigma > 0) cs%sigma(:,:,:) = 0.
560  if (cs%id_Ws > 0) allocate( cs%Ws( szi_(g), szj_(g), szk_(g) ) )
561  if (cs%id_Ws > 0) cs%Ws(:,:,:) = 0.
562  if (cs%id_N2 > 0) allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) )
563  if (cs%id_N2 > 0) cs%N2(:,:,:) = 0.
564  if (cs%id_Kt_KPP > 0) allocate( cs%Kt_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
565  if (cs%id_Kt_KPP > 0) cs%Kt_KPP(:,:,:) = 0.
566  if (cs%id_Ks_KPP > 0) allocate( cs%Ks_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
567  if (cs%id_Ks_KPP > 0) cs%Ks_KPP(:,:,:) = 0.
568  if (cs%id_Kv_KPP > 0) allocate( cs%Kv_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
569  if (cs%id_Kv_KPP > 0) cs%Kv_KPP(:,:,:) = 0.
570  if (cs%id_Tsurf > 0) allocate( cs%Tsurf( szi_(g), szj_(g)) )
571  if (cs%id_Tsurf > 0) cs%Tsurf(:,:) = 0.
572  if (cs%id_Ssurf > 0) allocate( cs%Ssurf( szi_(g), szj_(g)) )
573  if (cs%id_Ssurf > 0) cs%Ssurf(:,:) = 0.
574  if (cs%id_Usurf > 0) allocate( cs%Usurf( szib_(g), szj_(g)) )
575  if (cs%id_Usurf > 0) cs%Usurf(:,:) = 0.
576  if (cs%id_Vsurf > 0) allocate( cs%Vsurf( szi_(g), szjb_(g)) )
577  if (cs%id_Vsurf > 0) cs%Vsurf(:,:) = 0.
578  if (cs%id_EnhVt2 > 0) allocate( cs%EnhVt2( szi_(g), szj_(g), szk_(g)) )
579  if (cs%id_EnhVt2 > 0) cs%EnhVt2(:,:,:) = 0.
580  if (cs%id_EnhK > 0) allocate( cs%EnhK( szi_(g), szj_(g), szk_(g)+1 ) )
581  if (cs%id_EnhK > 0) cs%EnhK(:,:,:) = 0.
582 
583 

References mom_file_parser::closeparameterblock(), lt_k_constant, lt_k_mode_constant, lt_k_mode_rw16, lt_k_mode_vr12, lt_k_scaled, lt_vt2_mode_constant, lt_vt2_mode_lf17, lt_vt2_mode_rw16, lt_vt2_mode_vr12, mom_error_handler::mom_error(), nlt_shape_cubic, nlt_shape_cubic_lmd, nlt_shape_cvmix, nlt_shape_linear, nlt_shape_parabolic, mom_file_parser::openparameterblock(), mom_diag_mediator::register_diag_field(), sw_method_all_sw, sw_method_lv1_sw, and sw_method_mxl_sw.

Here is the call graph for this function:

◆ kpp_nonlocaltransport_saln()

subroutine, public mom_cvmix_kpp::kpp_nonlocaltransport_saln ( type(kpp_cs), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(in)  nonLocalTrans,
real, dimension(szi_(g),szj_(g)), intent(in)  surfFlux,
real, intent(in)  dt,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  scalar 
)

Apply KPP non-local transport of surface fluxes for salinity. This routine is a useful prototype for other material tracers.

Parameters
[in]csControl structure
[in]gOcean grid
[in]gvOcean vertical grid
[in]hLayer/level thickness [H ~> m or kg m-2]
[in]nonlocaltransNon-local transport [nondim]
[in]surffluxSurface flux of scalar [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
[in]dtTime-step [s]
[in,out]scalarScalar (scalar units [conc])

Definition at line 1543 of file MOM_CVMix_KPP.F90.

1543 
1544  type(KPP_CS), intent(in) :: CS !< Control structure
1545  type(ocean_grid_type), intent(in) :: G !< Ocean grid
1546  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1547  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
1548  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim]
1549  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of scalar
1550  !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
1551  real, intent(in) :: dt !< Time-step [s]
1552  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: scalar !< Scalar (scalar units [conc])
1553 
1554  integer :: i, j, k
1555  real, dimension( SZI_(G), SZJ_(G), SZK_(G) ) :: dtracer
1556 
1557 
1558  dtracer(:,:,:) = 0.0
1559  !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux)
1560  do k = 1, g%ke
1561  do j = g%jsc, g%jec
1562  do i = g%isc, g%iec
1563  dtracer(i,j,k) = ( nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1) ) / &
1564  ( h(i,j,k) + gv%H_subroundoff ) * surfflux(i,j)
1565  enddo
1566  enddo
1567  enddo
1568 
1569  ! Update tracer due to non-local redistribution of surface flux
1570  if (cs%applyNonLocalTrans) then
1571  do k = 1, g%ke
1572  do j = g%jsc, g%jec
1573  do i = g%isc, g%iec
1574  scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
1575  enddo
1576  enddo
1577  enddo
1578  endif
1579 
1580  ! Diagnostics
1581  if (cs%id_netS > 0) call post_data(cs%id_netS, surfflux, cs%diag)
1582  if (cs%id_NLT_dSdt > 0) call post_data(cs%id_NLT_dSdt, dtracer, cs%diag)
1583  if (cs%id_NLT_saln_budget > 0) then
1584  dtracer(:,:,:) = 0.0
1585  do k = 1, g%ke
1586  do j = g%jsc, g%jec
1587  do i = g%isc, g%iec
1588  dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * &
1589  surfflux(i,j) * gv%H_to_kg_m2
1590  enddo
1591  enddo
1592  enddo
1593  call post_data(cs%id_NLT_saln_budget, dtracer, cs%diag)
1594  endif
1595 

Referenced by mom_diabatic_driver::diabatic_ale(), mom_diabatic_driver::diabatic_ale_legacy(), and mom_diabatic_driver::layered_diabatic().

Here is the caller graph for this function:

◆ kpp_nonlocaltransport_temp()

subroutine, public mom_cvmix_kpp::kpp_nonlocaltransport_temp ( type(kpp_cs), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(in)  nonLocalTrans,
real, dimension(szi_(g),szj_(g)), intent(in)  surfFlux,
real, intent(in)  dt,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  scalar,
real, intent(in)  C_p 
)

Apply KPP non-local transport of surface fluxes for temperature.

Parameters
[in]csControl structure
[in]gOcean grid
[in]gvOcean vertical grid
[in]hLayer/level thickness [H ~> m or kg m-2]
[in]nonlocaltransNon-local transport [nondim]
[in]surffluxSurface flux of scalar [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
[in]dtTime-step [s]
[in,out]scalartemperature
[in]c_pSeawater specific heat capacity [J kg-1 degC-1]

Definition at line 1483 of file MOM_CVMix_KPP.F90.

1483 
1484  type(KPP_CS), intent(in) :: CS !< Control structure
1485  type(ocean_grid_type), intent(in) :: G !< Ocean grid
1486  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1487  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
1488  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim]
1489  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of scalar
1490  !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
1491  real, intent(in) :: dt !< Time-step [s]
1492  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: scalar !< temperature
1493  real, intent(in) :: C_p !< Seawater specific heat capacity [J kg-1 degC-1]
1494 
1495  integer :: i, j, k
1496  real, dimension( SZI_(G), SZJ_(G), SZK_(G) ) :: dtracer
1497 
1498 
1499  dtracer(:,:,:) = 0.0
1500  !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux)
1501  do k = 1, g%ke
1502  do j = g%jsc, g%jec
1503  do i = g%isc, g%iec
1504  dtracer(i,j,k) = ( nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1) ) / &
1505  ( h(i,j,k) + gv%H_subroundoff ) * surfflux(i,j)
1506  enddo
1507  enddo
1508  enddo
1509 
1510  ! Update tracer due to non-local redistribution of surface flux
1511  if (cs%applyNonLocalTrans) then
1512  do k = 1, g%ke
1513  do j = g%jsc, g%jec
1514  do i = g%isc, g%iec
1515  scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
1516  enddo
1517  enddo
1518  enddo
1519  endif
1520 
1521  ! Diagnostics
1522  if (cs%id_QminusSW > 0) call post_data(cs%id_QminusSW, surfflux, cs%diag)
1523  if (cs%id_NLT_dTdt > 0) call post_data(cs%id_NLT_dTdt, dtracer, cs%diag)
1524  if (cs%id_NLT_temp_budget > 0) then
1525  dtracer(:,:,:) = 0.0
1526  do k = 1, g%ke
1527  do j = g%jsc, g%jec
1528  do i = g%isc, g%iec
1529  dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * &
1530  surfflux(i,j) * c_p * gv%H_to_kg_m2
1531  enddo
1532  enddo
1533  enddo
1534  call post_data(cs%id_NLT_temp_budget, dtracer, cs%diag)
1535  endif
1536 

Referenced by mom_diabatic_driver::diabatic_ale(), mom_diabatic_driver::diabatic_ale_legacy(), and mom_diabatic_driver::layered_diabatic().

Here is the caller graph for this function:

◆ kpp_smooth_bld()

subroutine mom_cvmix_kpp::kpp_smooth_bld ( type(kpp_cs), pointer  CS,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h 
)
private

Apply a 1-1-4-1-1 Laplacian filter one time on BLD to reduce any horizontal two-grid-point noise.

Parameters
csControl structure
[in,out]gOcean grid
[in]gvOcean vertical grid
[in]hLayer/level thicknesses [H ~> m or kg m-2]

Definition at line 1350 of file MOM_CVMix_KPP.F90.

1350  ! Arguments
1351  type(KPP_CS), pointer :: CS !< Control structure
1352  type(ocean_grid_type), intent(inout) :: G !< Ocean grid
1353  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1354  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]
1355 
1356  ! local
1357  real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface [m]
1358  ! (negative in the ocean)
1359  real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m]
1360  ! (negative in the ocean)
1361  real :: wc, ww, we, wn, ws ! averaging weights for smoothing
1362  real :: dh ! The local thickness used for calculating interface positions [m]
1363  real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
1364  real :: pref
1365  integer :: i, j, k, s
1366 
1367  ncall_smooth = ncall_smooth+1
1368 
1369  if (cs%id_OBLdepth_original > 0) cs%OBLdepth_original = cs%OBLdepth
1370 
1371  cs%n_smooth = 2
1372  do s=1,cs%n_smooth
1373 
1374  ! Update halos
1375  !call pass_var(CS%OBLdepth, G%Domain, halo=4)
1376  if (s==1 .or. ncall_smooth<3) then
1377  call pass_var(cs%OBLdepth, g%Domain, halo=4)
1378  endif
1379 
1380  cs%OBLdepthprev = cs%OBLdepth
1381 
1382  ! apply smoothing on OBL depth
1383  do j = g%jsc-2, g%jec+2
1384  do i = g%isc-2, g%iec+2
1385 
1386  ! skip land points
1387  if (g%mask2dT(i,j)==0.) cycle
1388 
1389  ! compute weights
1390  !! fail
1391  ww = 0.125 * g%mask2dT(i-1,j)
1392  we = 0.125 * g%mask2dT(i+1,j)
1393  ws = 0.125 * g%mask2dT(i,j-1)
1394  wn = 0.125 * g%mask2dT(i,j+1)
1395 
1396  !! pass
1397  !ww = 0.25 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i-1,j)
1398  !we = 0.25 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i+1,j)
1399  !ws = 0.0 ! 0.125 * G%mask2dT(i,j-1)
1400  !wn = 0.0 ! 0.125 * G%mask2dT(i,j+1)
1401 
1402  !! fail
1403  !!ww = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i-1,j)
1404  !!we = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i+1,j)
1405  !!ws = 0.25 ! 0.125 * G%mask2dT(i,j-1)
1406  !!wn = 0.25 ! 0.125 * G%mask2dT(i,j+1)
1407 
1408  !! pass
1409  !!ww = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i-1,j)
1410  !!we = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i+1,j)
1411  !!ws = 0.25 ! 0.125 * G%mask2dT(i,j-1)
1412  !!wn = 0.0 ! 0.125 * G%mask2dT(i,j+1)
1413 
1414  !! fail
1415  !!ww = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i-1,j)
1416  !!we = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i+1,j)
1417  !!ws = 0.0 ! 0.125 * G%mask2dT(i,j-1)
1418  !!wn = 0.25 ! 0.125 * G%mask2dT(i,j+1)
1419 
1420  wc = 1.0 - (ww+we+wn+ws)
1421 
1422  cs%OBLdepth(i,j) = wc * cs%OBLdepthprev(i,j) &
1423  + ww * cs%OBLdepthprev(i-1,j) &
1424  + we * cs%OBLdepthprev(i+1,j) &
1425  + ws * cs%OBLdepthprev(i,j-1) &
1426  + wn * cs%OBLdepthprev(i,j+1)
1427 
1428  enddo
1429  enddo
1430 
1431  enddo ! s-loop
1432 
1433  ! Update kOBL for smoothed OBL depths
1434  do j = g%jsc, g%jec
1435  do i = g%isc, g%iec
1436 
1437  ! skip land points
1438  if (g%mask2dT(i,j)==0.) cycle
1439 
1440  ifaceheight(1) = 0.0 ! BBL is all relative to the surface
1441  hcorr = 0.
1442  do k=1,g%ke
1443 
1444  ! cell center and cell bottom in meters (negative values in the ocean)
1445  dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
1446  dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
1447  hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
1448  dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
1449  cellheight(k) = ifaceheight(k) - 0.5 * dh
1450  ifaceheight(k+1) = ifaceheight(k) - dh
1451  enddo
1452 
1453  ! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing.
1454  if (cs%deepen_only) cs%OBLdepth(i,j) = max(cs%OBLdepth(i,j),cs%OBLdepth_original(i,j))
1455 
1456  ! prevent OBL depths deeper than the bathymetric depth
1457  cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -ifaceheight(g%ke+1) ) ! no deeper than bottom
1458 
1459  cs%kOBL(i,j) = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, cs%OBLdepth(i,j) )
1460  enddo
1461  enddo
1462 

References ncall_smooth.

Referenced by kpp_compute_bld().

Here is the caller graph for this function:

Variable Documentation

◆ lt_k_constant

integer, parameter, private mom_cvmix_kpp::lt_k_constant = 1
private

Constant enhance K through column.

Definition at line 53 of file MOM_CVMix_KPP.F90.

53 integer, private, parameter :: LT_K_CONSTANT = 1, & !< Constant enhance K through column
54  lt_k_scaled = 2, & !< Enhance K scales with G(sigma)
55  lt_k_mode_constant = 1, & !< Prescribed enhancement for K
56  lt_k_mode_vr12 = 2, & !< Enhancement for K based on
57  !! Van Roekel et al., 2012
58  lt_k_mode_rw16 = 3, & !< Enhancement for K based on
59  !! Reichl et al., 2016
60  lt_vt2_mode_constant = 1, & !< Prescribed enhancement for Vt2
61  lt_vt2_mode_vr12 = 2, & !< Enhancement for Vt2 based on
62  !! Van Roekel et al., 2012
63  lt_vt2_mode_rw16 = 3, & !< Enhancement for Vt2 based on
64  !! Reichl et al., 2016
65  lt_vt2_mode_lf17 = 4 !< Enhancement for Vt2 based on

Referenced by kpp_calculate(), and kpp_init().

◆ lt_k_mode_constant

integer, parameter, private mom_cvmix_kpp::lt_k_mode_constant = 1
private

Prescribed enhancement for K.

Definition at line 53 of file MOM_CVMix_KPP.F90.

Referenced by kpp_calculate(), and kpp_init().

◆ lt_k_mode_rw16

integer, parameter, private mom_cvmix_kpp::lt_k_mode_rw16 = 3
private

Enhancement for K based on.

Definition at line 53 of file MOM_CVMix_KPP.F90.

Referenced by kpp_calculate(), and kpp_init().

◆ lt_k_mode_vr12

integer, parameter, private mom_cvmix_kpp::lt_k_mode_vr12 = 2
private

Enhancement for K based on.

Definition at line 53 of file MOM_CVMix_KPP.F90.

Referenced by kpp_calculate(), and kpp_init().

◆ lt_k_scaled

integer, parameter, private mom_cvmix_kpp::lt_k_scaled = 2
private

Enhance K scales with G(sigma)

Definition at line 53 of file MOM_CVMix_KPP.F90.

Referenced by kpp_calculate(), and kpp_init().

◆ lt_vt2_mode_constant

integer, parameter, private mom_cvmix_kpp::lt_vt2_mode_constant = 1
private

Prescribed enhancement for Vt2.

Definition at line 53 of file MOM_CVMix_KPP.F90.

Referenced by kpp_compute_bld(), and kpp_init().

◆ lt_vt2_mode_lf17

integer, parameter, private mom_cvmix_kpp::lt_vt2_mode_lf17 = 4
private

Enhancement for Vt2 based on.

Definition at line 53 of file MOM_CVMix_KPP.F90.

Referenced by kpp_compute_bld(), and kpp_init().

◆ lt_vt2_mode_rw16

integer, parameter, private mom_cvmix_kpp::lt_vt2_mode_rw16 = 3
private

Enhancement for Vt2 based on.

Definition at line 53 of file MOM_CVMix_KPP.F90.

Referenced by kpp_compute_bld(), and kpp_init().

◆ lt_vt2_mode_vr12

integer, parameter, private mom_cvmix_kpp::lt_vt2_mode_vr12 = 2
private

Enhancement for Vt2 based on.

Definition at line 53 of file MOM_CVMix_KPP.F90.

Referenced by kpp_compute_bld(), and kpp_init().

◆ ncall_smooth

integer mom_cvmix_kpp::ncall_smooth = 0
private

Definition at line 68 of file MOM_CVMix_KPP.F90.

68 integer :: ncall_smooth = 0

Referenced by kpp_smooth_bld().

◆ nlt_shape_cubic

integer, parameter, private mom_cvmix_kpp::nlt_shape_cubic = 3
private

Cubic, \( G(\sigma) = 1 + (2\sigma-3) \sigma^2\).

Definition at line 46 of file MOM_CVMix_KPP.F90.

46 integer, private, parameter :: NLT_SHAPE_CUBIC = 3 !< Cubic, \f$ G(\sigma) = 1 + (2\sigma-3) \sigma^2\f$

Referenced by kpp_calculate(), and kpp_init().

◆ nlt_shape_cubic_lmd

integer, parameter, private mom_cvmix_kpp::nlt_shape_cubic_lmd = 4
private

Original shape, \( G(\sigma) = \frac{27}{4} \sigma (1-\sigma)^2 \).

Definition at line 47 of file MOM_CVMix_KPP.F90.

47 integer, private, parameter :: NLT_SHAPE_CUBIC_LMD = 4 !< Original shape,

Referenced by kpp_calculate(), and kpp_init().

◆ nlt_shape_cvmix

integer, parameter, private mom_cvmix_kpp::nlt_shape_cvmix = 0
private

Use the CVMix profile.

Definition at line 43 of file MOM_CVMix_KPP.F90.

43 integer, private, parameter :: NLT_SHAPE_CVMix = 0 !< Use the CVMix profile

Referenced by kpp_init().

◆ nlt_shape_linear

integer, parameter, private mom_cvmix_kpp::nlt_shape_linear = 1
private

Linear, \( G(\sigma) = 1-\sigma \).

Definition at line 44 of file MOM_CVMix_KPP.F90.

44 integer, private, parameter :: NLT_SHAPE_LINEAR = 1 !< Linear, \f$ G(\sigma) = 1-\sigma \f$

Referenced by kpp_calculate(), and kpp_init().

◆ nlt_shape_parabolic

integer, parameter, private mom_cvmix_kpp::nlt_shape_parabolic = 2
private

Parabolic, \( G(\sigma) = (1-\sigma)^2 \).

Definition at line 45 of file MOM_CVMix_KPP.F90.

45 integer, private, parameter :: NLT_SHAPE_PARABOLIC = 2 !< Parabolic, \f$ G(\sigma) = (1-\sigma)^2 \f$

Referenced by kpp_calculate(), and kpp_init().

◆ sw_method_all_sw

integer, parameter, private mom_cvmix_kpp::sw_method_all_sw = 0
private

Use all shortwave radiation.

Definition at line 50 of file MOM_CVMix_KPP.F90.

50 integer, private, parameter :: SW_METHOD_ALL_SW = 0 !< Use all shortwave radiation

Referenced by kpp_calculate(), and kpp_init().

◆ sw_method_lv1_sw

integer, parameter, private mom_cvmix_kpp::sw_method_lv1_sw = 2
private

Use shortwave radiation absorbed in layer 1.

Definition at line 52 of file MOM_CVMix_KPP.F90.

52 integer, private, parameter :: SW_METHOD_LV1_SW = 2 !< Use shortwave radiation absorbed in layer 1

Referenced by kpp_calculate(), and kpp_init().

◆ sw_method_mxl_sw

integer, parameter, private mom_cvmix_kpp::sw_method_mxl_sw = 1
private

Use shortwave radiation absorbed in mixing layer.

Definition at line 51 of file MOM_CVMix_KPP.F90.

51 integer, private, parameter :: SW_METHOD_MXL_SW = 1 !< Use shortwave radiation absorbed in mixing layer

Referenced by kpp_calculate(), and kpp_init().