Detailed Description

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

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, tv, 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, US, m_to_BLD_units)
	Copies KPP surface boundary layer depth into BLD, in units of [Z ~> m] unless other units are specified. 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.

integer, parameter, private	nlt_shape_linear = 1
	Linear, \( G(\sigma) = 1-\sigma \).

integer, parameter, private	nlt_shape_parabolic = 2
	Parabolic, \( G(\sigma) = (1-\sigma)^2 \).

integer, parameter, private	nlt_shape_cubic = 3
	Cubic, \( G(\sigma) = 1 + (2\sigma-3) \sigma^2\).

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

integer, parameter, private	sw_method_all_sw = 0
	Use all shortwave radiation.

integer, parameter, private	sw_method_mxl_sw = 1
	Use shortwave radiation absorbed in mixing layer.

integer, parameter, private	sw_method_lv1_sw = 2
	Use shortwave radiation absorbed in layer 1.

integer, parameter, private	lt_k_constant = 1
	Constant enhance K through column.

integer, parameter, private	lt_k_scaled = 2
	Enhance K scales with G(sigma)

integer, parameter, private	lt_k_mode_constant = 1
	Prescribed enhancement for K.

integer, parameter, private	lt_k_mode_vr12 = 2
	Enhancement for K based on.

integer, parameter, private	lt_k_mode_rw16 = 3
	Enhancement for K based on.

integer, parameter, private	lt_vt2_mode_constant = 1
	Prescribed enhancement for Vt2.

integer, parameter, private	lt_vt2_mode_vr12 = 2
	Enhancement for Vt2 based on.

integer, parameter, private	lt_vt2_mode_rw16 = 3
	Enhancement for Vt2 based on.

integer, parameter, private	lt_vt2_mode_lf17 = 4
	Enhancement for Vt2 based on.


integer	id_clock_kpp_calc
	CPU time clocks.

integer	id_clock_kpp_compute_bld
	CPU time clocks.

integer	id_clock_kpp_smoothing
	CPU time clocks.

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

	cs	Control structure
[in]	g	Ocean grid
[in]	gv	Ocean vertical grid
[in]	us	A dimensional unit scaling type
	waves	Wave CS
[in]	h	Layer/level thicknesses [H ~> m or kg m-2]
[in]	ustar	Surface friction velocity [Z T-1 ~> m s-1]
[in]	buoyflux	Surface 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]	nonlocaltransheat	Temp non-local transport [m s-1]
[in,out]	nonlocaltransscalar	scalar non-local transport [m s-1]

Definition at line 604 of file MOM_CVMix_KPP.F90.

 
   ! Arguments
   type(KPP_CS),                               pointer       :: CS    !< Control structure
   type(ocean_grid_type),                      intent(in)    :: G     !< Ocean grid
   type(verticalGrid_type),                    intent(in)    :: GV    !< Ocean vertical grid
   type(unit_scale_type),                      intent(in)    :: US    !< A dimensional unit scaling type
   type(wave_parameters_CS),         optional, pointer       :: Waves !< Wave CS
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)    :: h     !< Layer/level thicknesses [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G)),           intent(in)    :: uStar !< Surface friction velocity [Z T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kt   !< (in)  Vertical diffusivity of heat w/o KPP
                                                                     !! (out) Vertical diffusivity including KPP
                                                                     !!       [Z2 T-1 ~> m2 s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Ks   !< (in)  Vertical diffusivity of salt w/o KPP
                                                                     !! (out) Vertical diffusivity including KPP
                                                                     !!       [Z2 T-1 ~> m2 s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kv   !< (in)  Vertical viscosity w/o KPP
                                                                     !! (out) Vertical viscosity including KPP
                                                                     !!       [Z2 T-1 ~> m2 s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransHeat   !< Temp non-local transport [m s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local transport [m s-1]
 
 ! Local variables
   integer :: i, j, k                             ! Loop indices
   real, dimension( G%ke )     :: cellHeight      ! Cell center heights referenced to surface [m] (negative in ocean)
   real, dimension( G%ke+1 )   :: iFaceHeight     ! Interface heights referenced to surface [m] (negative in ocean)
   real, dimension( G%ke+1, 2) :: Kdiffusivity    ! Vertical diffusivity at interfaces [m2 s-1]
   real, dimension( G%ke+1 )   :: Kviscosity      ! Vertical viscosity at interfaces [m2 s-1]
   real, dimension( G%ke+1, 2) :: nonLocalTrans   ! Non-local transport for heat/salt at interfaces [nondim]
 
   real :: surfFricVel, surfBuoyFlux
   real :: sigma, sigmaRatio
   real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim]
   real :: dh    ! The local thickness used for calculating interface positions [m]
   real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
 
   ! For Langmuir Calculations
   real :: LangEnhK     ! Langmuir enhancement for mixing coefficient
 
 
   if (cs%debug) then
     call hchksum(h, "KPP in: h",g%HI,haloshift=0, scale=gv%H_to_m)
     call hchksum(ustar, "KPP in: uStar",g%HI,haloshift=0, scale=us%Z_to_m*us%s_to_T)
     call hchksum(buoyflux, "KPP in: buoyFlux",g%HI,haloshift=0)
     call hchksum(kt, "KPP in: Kt",g%HI,haloshift=0, scale=us%Z2_T_to_m2_s)
     call hchksum(ks, "KPP in: Ks",g%HI,haloshift=0, scale=us%Z2_T_to_m2_s)
   endif
 
   nonlocaltrans(:,:) = 0.0
 
   if (cs%id_Kd_in > 0) call post_data(cs%id_Kd_in, kt, cs%diag)
 
   call cpu_clock_begin(id_clock_kpp_calc)
   buoy_scale = us%L_to_m**2*us%s_to_T**3
 
   !$OMP parallel do default(none) firstprivate(nonLocalTrans)                               &
   !$OMP                           private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight,  &
   !$OMP                           surfBuoyFlux, Kdiffusivity, Kviscosity, LangEnhK, sigma,  &
   !$OMP                           sigmaRatio)                                               &
   !$OMP                           shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, Kt, &
   !$OMP                           Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves)
   ! loop over horizontal points on processor
   do j = g%jsc, g%jec
     do i = g%isc, g%iec
 
       ! skip calling KPP for land points
       if (g%mask2dT(i,j)==0.) cycle
 
       ! things independent of position within the column
       surffricvel = us%Z_to_m*us%s_to_T * ustar(i,j)
 
       ifaceheight(1) = 0.0 ! BBL is all relative to the surface
       hcorr = 0.
       do k=1,g%ke
 
         ! cell center and cell bottom in meters (negative values in the ocean)
         dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
         dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
         hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
         dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
         cellheight(k)    = ifaceheight(k) - 0.5 * dh
         ifaceheight(k+1) = ifaceheight(k) - dh
 
       enddo ! k-loop finishes
 
       surfbuoyflux = buoy_scale*buoyflux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
                                      ! h to Monin-Obukov (default is false, ie. not used)
 
       ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports
 
       ! Unlike LMD94, we do not match to interior diffusivities. If using the original
       ! LMD94 shape function, not matching is equivalent to matching to a zero diffusivity.
 
       !BGR/ Add option for use of surface buoyancy flux with total sw flux.
       if (cs%SW_METHOD == sw_method_all_sw) then
          surfbuoyflux = buoy_scale * buoyflux(i,j,1)
       elseif (cs%SW_METHOD == sw_method_mxl_sw) then
          ! We know the actual buoyancy flux into the OBL
          surfbuoyflux  = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,int(cs%kOBL(i,j))+1))
       elseif (cs%SW_METHOD == sw_method_lv1_sw) then
          surfbuoyflux  = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,2))
       endif
 
       ! If option "MatchBoth" is selected in CVMix, MOM should be capable of matching.
       if (.not. (cs%MatchTechnique == 'MatchBoth')) then
          kdiffusivity(:,:) = 0. ! Diffusivities for heat and salt [m2 s-1]
          kviscosity(:)     = 0. ! Viscosity [m2 s-1]
       else
          kdiffusivity(:,1) = us%Z2_T_to_m2_s * kt(i,j,:)
          kdiffusivity(:,2) = us%Z2_T_to_m2_s * ks(i,j,:)
          kviscosity(:) = us%Z2_T_to_m2_s * kv(i,j,:)
       endif
 
       call cvmix_coeffs_kpp(kviscosity(:),        & ! (inout) Total viscosity [m2 s-1]
                             kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1]
                             kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1]
                             ifaceheight,       & ! (in) Height of interfaces [m]
                             cellheight,        & ! (in) Height of level centers [m]
                             kviscosity(:),        & ! (in) Original viscosity [m2 s-1]
                             kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1]
                             kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1]
                             cs%OBLdepth(i,j),  & ! (in) OBL depth [m]
                             cs%kOBL(i,j),      & ! (in) level (+fraction) of OBL extent
                             nonlocaltrans(:,1),& ! (out) Non-local heat transport [nondim]
                             nonlocaltrans(:,2),& ! (out) Non-local salt transport [nondim]
                             surffricvel,       & ! (in) Turbulent friction velocity at surface [m s-1]
                             surfbuoyflux,      & ! (in) Buoyancy flux at surface [m2 s-3]
                             g%ke,              & ! (in) Number of levels to compute coeffs for
                             g%ke,              & ! (in) Number of levels in array shape
                             cvmix_kpp_params_user=cs%KPP_params )
 
       ! safety check, Kviscosity and Kdiffusivity must be >= 0
       do k=1, g%ke+1
         if (kviscosity(k) < 0. .or. kdiffusivity(k,1) < 0.) then
           call mom_error(fatal,"KPP_calculate, after CVMix_coeffs_kpp: "// &
                    "Negative vertical viscosity or diffusivity has been detected. " // &
                    "This is likely related to the choice of MATCH_TECHNIQUE and INTERP_TYPE2." //&
                    "You might consider using the default options for these parameters." )
         endif
       enddo
 
       IF (cs%LT_K_ENHANCEMENT) then
         if (cs%LT_K_METHOD==lt_k_mode_constant) then
            langenhk = cs%KPP_K_ENH_FAC
         elseif (cs%LT_K_METHOD==lt_k_mode_vr12) then
            ! Added minimum value for La_SL, so removed maximum value for LangEnhK.
            langenhk = sqrt(1.+(1.5*cs%La_SL(i,j))**(-2) + &
                 (5.4*cs%La_SL(i,j))**(-4))
         elseif (cs%LT_K_METHOD==lt_k_mode_rw16) then
           !This maximum value is proposed in Reichl et al., 2016 JPO formula
           langenhk = min(2.25, 1. + 1./cs%La_SL(i,j))
         else
            !This shouldn't be reached.
            !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT")
            langenhk = 1.0
         endif
         do k=1,g%ke
           if (cs%LT_K_SHAPE== lt_k_constant) then
             if (cs%id_EnhK > 0) cs%EnhK(i,j,:) = langenhk
             kdiffusivity(k,1) = kdiffusivity(k,1) * langenhk
             kdiffusivity(k,2) = kdiffusivity(k,2) * langenhk
             kviscosity(k)     = kviscosity(k)   * langenhk
           elseif (cs%LT_K_SHAPE == lt_k_scaled) then
             sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
             sigmaratio = sigma * (1. - sigma)**2 / 0.148148037
             if (cs%id_EnhK > 0) cs%EnhK(i,j,k) = (1.0 + (langenhk - 1.)*sigmaratio)
             kdiffusivity(k,1) = kdiffusivity(k,1) * ( 1. + &
                                 ( langenhk - 1.)*sigmaratio)
             kdiffusivity(k,2) = kdiffusivity(k,2) * ( 1. + &
                                 ( langenhk - 1.)*sigmaratio)
             kviscosity(k) = kviscosity(k) * ( 1. + &
                                 ( langenhk - 1.)*sigmaratio)
           endif
         enddo
       endif
 
       ! Over-write CVMix NLT shape function with one of the following choices.
       ! The CVMix code has yet to update for thse options, so we compute in MOM6.
       ! Note that nonLocalTrans = Cs * G(sigma) (LMD94 notation), with
       ! Cs = 6.32739901508.
       ! Start do-loop at k=2, since k=1 is ocean surface (sigma=0)
       ! and we do not wish to double-count the surface forcing.
       ! Only compute nonlocal transport for 0 <= sigma <= 1.
       ! MOM6 recommended shape is the parabolic; it gives deeper boundary layer
       ! and no spurious extrema.
       if (surfbuoyflux < 0.0) then
         if (cs%NLT_shape == nlt_shape_cubic) then
           do k = 2, g%ke
             sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
             nonlocaltrans(k,1) = (1.0 - sigma)**2 * (1.0 + 2.0*sigma) !*
             nonlocaltrans(k,2) = nonlocaltrans(k,1)
           enddo
         elseif (cs%NLT_shape == nlt_shape_parabolic) then
           do k = 2, g%ke
             sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
             nonlocaltrans(k,1) = (1.0 - sigma)**2 !*CS%CS2
             nonlocaltrans(k,2) = nonlocaltrans(k,1)
           enddo
         elseif (cs%NLT_shape == nlt_shape_linear) then
           do k = 2, g%ke
             sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
             nonlocaltrans(k,1) = (1.0 - sigma)!*CS%CS2
             nonlocaltrans(k,2) = nonlocaltrans(k,1)
           enddo
         elseif (cs%NLT_shape == nlt_shape_cubic_lmd) then
           ! Sanity check (should agree with CVMix result using simple matching)
           do k = 2, g%ke
             sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
             nonlocaltrans(k,1) = cs%CS2 * sigma*(1.0 -sigma)**2
             nonlocaltrans(k,2) = nonlocaltrans(k,1)
           enddo
         endif
       endif
 
       ! we apply nonLocalTrans in subroutines
       ! KPP_NonLocalTransport_temp and KPP_NonLocalTransport_saln
       nonlocaltransheat(i,j,:)   = nonlocaltrans(:,1) ! temp
       nonlocaltransscalar(i,j,:) = nonlocaltrans(:,2) ! saln
 
       ! set the KPP diffusivity and viscosity to zero for testing purposes
       if (cs%KPPzeroDiffusivity) then
          kdiffusivity(:,1) = 0.0
          kdiffusivity(:,2) = 0.0
          kviscosity(:)     = 0.0
       endif
 
 
       ! compute unresolved squared velocity for diagnostics
       if (cs%id_Vt2 > 0) then
 !BGR Now computing VT2 above so can modify for LT
 !    therefore, don't repeat this operation here
 !        CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( &
 !                    cellHeight(1:G%ke),                 & ! Depth of cell center [m]
 !                    ws_cntr=Ws_1d,                      & ! Turbulent velocity scale profile, at centers [m s-1]
 !                    N_iface=CS%N(i,j,:),                & ! Buoyancy frequency at interface [s-1]
 !                    CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters
       endif
 
       ! Copy 1d data into 3d diagnostic arrays
       !/ grabbing obldepth_0d for next time step.
       cs%OBLdepthprev(i,j)=cs%OBLdepth(i,j)
       if (cs%id_sigma > 0) then
         cs%sigma(i,j,:)  = 0.
         if (cs%OBLdepth(i,j)>0.)   cs%sigma(i,j,:)  = -ifaceheight/cs%OBLdepth(i,j)
       endif
       if (cs%id_Kt_KPP > 0)   cs%Kt_KPP(i,j,:) = kdiffusivity(:,1)
       if (cs%id_Ks_KPP > 0)   cs%Ks_KPP(i,j,:) = kdiffusivity(:,2)
       if (cs%id_Kv_KPP > 0)   cs%Kv_KPP(i,j,:) = kviscosity(:)
 
       ! Update output of routine
       if (.not. cs%passiveMode) then
         if (cs%KPPisAdditive) then
           do k=1, g%ke+1
             kt(i,j,k) = kt(i,j,k) + us%m2_s_to_Z2_T * kdiffusivity(k,1)
             ks(i,j,k) = ks(i,j,k) + us%m2_s_to_Z2_T * kdiffusivity(k,2)
             kv(i,j,k) = kv(i,j,k) + us%m2_s_to_Z2_T * kviscosity(k)
             if (cs%Stokes_Mixing) waves%KvS(i,j,k) = kv(i,j,k)
           enddo
         else ! KPP replaces prior diffusivity when former is non-zero
           do k=1, g%ke+1
             if (kdiffusivity(k,1) /= 0.) kt(i,j,k) = us%m2_s_to_Z2_T * kdiffusivity(k,1)
             if (kdiffusivity(k,2) /= 0.) ks(i,j,k) = us%m2_s_to_Z2_T * kdiffusivity(k,2)
             if (kviscosity(k) /= 0.) kv(i,j,k) = us%m2_s_to_Z2_T * kviscosity(k)
             if (cs%Stokes_Mixing) waves%KvS(i,j,k) = kv(i,j,k)
           enddo
         endif
       endif
 
 
     ! end of the horizontal do-loops over the vertical columns
     enddo ! i
   enddo ! j
 
   call cpu_clock_end(id_clock_kpp_calc)
 
   if (cs%debug) then
     call hchksum(kt, "KPP out: Kt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
     call hchksum(ks, "KPP out: Ks", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
   endif
 
   ! send diagnostics to post_data
   if (cs%id_OBLdepth > 0) call post_data(cs%id_OBLdepth, cs%OBLdepth,        cs%diag)
   if (cs%id_OBLdepth_original > 0) call post_data(cs%id_OBLdepth_original,cs%OBLdepth_original,cs%diag)
   if (cs%id_sigma    > 0) call post_data(cs%id_sigma,    cs%sigma,           cs%diag)
   if (cs%id_Ws       > 0) call post_data(cs%id_Ws,       cs%Ws,              cs%diag)
   if (cs%id_Vt2      > 0) call post_data(cs%id_Vt2,      cs%Vt2,             cs%diag)
   if (cs%id_uStar    > 0) call post_data(cs%id_uStar,    ustar,              cs%diag)
   if (cs%id_buoyFlux > 0) call post_data(cs%id_buoyFlux, buoyflux,           cs%diag)
   if (cs%id_Kt_KPP   > 0) call post_data(cs%id_Kt_KPP,   cs%Kt_KPP,          cs%diag)
   if (cs%id_Ks_KPP   > 0) call post_data(cs%id_Ks_KPP,   cs%Ks_KPP,          cs%diag)
   if (cs%id_Kv_KPP   > 0) call post_data(cs%id_Kv_KPP,   cs%Kv_KPP,          cs%diag)
   if (cs%id_NLTt     > 0) call post_data(cs%id_NLTt,     nonlocaltransheat,  cs%diag)
   if (cs%id_NLTs     > 0) call post_data(cs%id_NLTs,     nonlocaltransscalar,cs%diag)
 
 

◆ 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(thermo_var_ptrs), intent(in)	tv,
		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

	cs	Control structure
[in,out]	g	Ocean grid
[in]	gv	Ocean vertical grid
[in]	us	A dimensional unit scaling type
[in]	h	Layer/level thicknesses [H ~> m or kg m-2]
[in]	temp	potential/cons temp [degC]
[in]	salt	Salinity [ppt]
[in]	u	Velocity i-component [L T-1 ~> m s-1]
[in]	v	Velocity j-component [L T-1 ~> m s-1]
[in]	tv	Thermodynamics structure.
[in]	ustar	Surface friction velocity [Z T-1 ~> m s-1]
[in]	buoyflux	Surface buoyancy flux [L2 T-3 ~> m2 s-3]
	waves	Wave CS

Definition at line 904 of file MOM_CVMix_KPP.F90.

 
   ! Arguments
   type(KPP_CS),                               pointer       :: CS    !< Control structure
   type(ocean_grid_type),                      intent(inout) :: G     !< Ocean grid
   type(verticalGrid_type),                    intent(in)    :: GV    !< Ocean vertical grid
   type(unit_scale_type),                      intent(in)    :: US    !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)    :: h     !< Layer/level thicknesses [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)    :: Temp  !< potential/cons temp [degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)    :: Salt  !< Salinity [ppt]
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)),  intent(in)    :: u     !< Velocity i-component [L T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)),  intent(in)    :: v     !< Velocity j-component [L T-1 ~> m s-1]
   type(thermo_var_ptrs),                      intent(in)    :: tv    !< Thermodynamics structure.
   real, dimension(SZI_(G),SZJ_(G)),           intent(in)    :: uStar !< Surface friction velocity [Z T-1 ~> m s-1]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
   type(wave_parameters_CS),         optional, pointer       :: Waves !< Wave CS
 
   ! Local variables
   integer :: i, j, k, km1                        ! Loop indices
   real, dimension( G%ke )     :: cellHeight      ! Cell center heights referenced to surface [m] (negative in ocean)
   real, dimension( G%ke+1 )   :: iFaceHeight     ! Interface heights referenced to surface [m] (negative in ocean)
   real, dimension( G%ke+1 )   :: N2_1d           ! Brunt-Vaisala frequency squared, at interfaces [s-2]
   real, dimension( G%ke )     :: Ws_1d           ! Profile of vertical velocity scale for scalars [m s-1]
   real, dimension( G%ke )     :: deltaRho        ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3]
   real, dimension( G%ke )     :: deltaU2         ! square of delta U (shear) in denominator of Bulk Ri [m2 s-2]
   real, dimension( G%ke )     :: surfBuoyFlux2
   real, dimension( G%ke )     :: BulkRi_1d       ! Bulk Richardson number for each layer
 
   ! for EOS calculation
   real, dimension( 3*G%ke )   :: rho_1D   ! A column of densities [R ~> kg m-3]
   real, dimension( 3*G%ke )   :: pres_1D  ! A column of pressures [R L2 T-2 ~> Pa]
   real, dimension( 3*G%ke )   :: Temp_1D
   real, dimension( 3*G%ke )   :: Salt_1D
 
   real :: surfFricVel, surfBuoyFlux, Coriolis
   real :: GoRho  ! Gravitational acceleration divided by density in MKS units [m R-1 s-2 ~> m4 kg-1 s-2]
   real :: pRef   ! The interface pressure [R L2 T-2 ~> Pa]
   real :: rho1, rhoK, Uk, Vk, sigma, sigmaRatio
 
   real :: zBottomMinusOffset   ! Height of bottom plus a little bit [m]
   real :: SLdepth_0d           ! Surface layer depth = surf_layer_ext*OBLdepth.
   real :: hTot                 ! Running sum of thickness used in the surface layer average [m]
   real :: buoy_scale           ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim]
   real :: delH                 ! Thickness of a layer [m]
   real :: surfHtemp, surfTemp  ! Integral and average of temp over the surface layer
   real :: surfHsalt, surfSalt  ! Integral and average of saln over the surface layer
   real :: surfHu, surfU        ! Integral and average of u over the surface layer
   real :: surfHv, surfV        ! Integral and average of v over the surface layer
   real :: dh    ! The local thickness used for calculating interface positions [m]
   real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
   integer :: kk, ksfc, ktmp
 
   ! For Langmuir Calculations
   real :: LangEnhW     ! Langmuir enhancement for turbulent velocity scale
   real, dimension(G%ke) :: LangEnhVt2   ! Langmuir enhancement for unresolved shear
   real, dimension(G%ke) :: U_H, V_H
   real :: MLD_GUESS, LA
   real :: surfHuS, surfHvS, surfUs, surfVs, wavedir, currentdir
   real :: VarUp, VarDn, M, VarLo, VarAvg
   real :: H10pct, H20pct,CMNFACT, USx20pct, USy20pct, enhvt2
   integer :: B
   real :: WST
 
 
   if (cs%debug) then
     call hchksum(salt, "KPP in: S",g%HI,haloshift=0)
     call hchksum(temp, "KPP in: T",g%HI,haloshift=0)
     call hchksum(u, "KPP in: u",g%HI,haloshift=0)
     call hchksum(v, "KPP in: v",g%HI,haloshift=0)
   endif
 
   call cpu_clock_begin(id_clock_kpp_compute_bld)
 
   ! some constants
   gorho = us%L_T_to_m_s**2*us%m_to_Z * gv%g_Earth / gv%Rho0
   buoy_scale = us%L_to_m**2*us%s_to_T**3
 
   ! loop over horizontal points on processor
   !$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight,  &
   !$OMP                           surfBuoyFlux, U_H, V_H, Coriolis, pRef, SLdepth_0d,       &
   !$OMP                           ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS,      &
   !$OMP                           surfHvS, hTot, delH, surftemp, surfsalt, surfu, surfv,    &
   !$OMP                           surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D,        &
   !$OMP                           Temp_1D, salt_1D, surfBuoyFlux2, MLD_GUESS, LA, rho_1D,   &
   !$OMP                           deltarho, N2_1d, ws_1d, LangEnhVT2, enhvt2, wst,          &
   !$OMP                           BulkRi_1d, zBottomMinusOffset) &
   !$OMP                           shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux,     &
   !$OMP                           Temp, Salt, waves, tv, GoRho, u, v)
   do j = g%jsc, g%jec
     do i = g%isc, g%iec
 
       ! skip calling KPP for land points
       if (g%mask2dT(i,j)==0.) cycle
 
       do k=1,g%ke
         u_h(k) = 0.5 * us%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k))
         v_h(k) = 0.5 * us%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k))
       enddo
 
       ! things independent of position within the column
       coriolis = 0.25*us%s_to_T*( (g%CoriolisBu(i,j)   + g%CoriolisBu(i-1,j-1)) + &
                                   (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1)) )
       surffricvel = us%Z_to_m*us%s_to_T * ustar(i,j)
 
       ! Bullk Richardson number computed for each cell in a column,
       ! assuming OBLdepth = grid cell depth. After Rib(k) is
       ! known for the column, then CVMix interpolates to find
       ! the actual OBLdepth. This approach avoids need to iterate
       ! on the OBLdepth calculation. It follows that used in MOM5
       ! and POP.
       ifaceheight(1) = 0.0 ! BBL is all relative to the surface
       pref = 0. ; if (associated(tv%p_surf)) pref = tv%p_surf(i,j)
       hcorr = 0.
       do k=1,g%ke
 
         ! cell center and cell bottom in meters (negative values in the ocean)
         dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
         dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
         hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
         dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
         cellheight(k)    = ifaceheight(k) - 0.5 * dh
         ifaceheight(k+1) = ifaceheight(k) - dh
 
         ! find ksfc for cell where "surface layer" sits
         sldepth_0d = cs%surf_layer_ext*max( max(-cellheight(k),-ifaceheight(2) ), cs%minOBLdepth )
         ksfc = k
         do ktmp = 1,k
           if (-1.0*ifaceheight(ktmp+1) >= sldepth_0d) then
             ksfc = ktmp
             exit
           endif
         enddo
 
         ! average temp, saln, u, v over surface layer
         ! use C-grid average to get u,v on T-points.
         surfhtemp=0.0
         surfhsalt=0.0
         surfhu   =0.0
         surfhv   =0.0
         surfhus  =0.0
         surfhvs  =0.0
         htot     =0.0
         do ktmp = 1,ksfc
 
           ! SLdepth_0d can be between cell interfaces
           delh = min( max(0.0, sldepth_0d - htot), h(i,j,ktmp)*gv%H_to_m )
 
           ! surface layer thickness
           htot = htot + delh
 
           ! surface averaged fields
           surfhtemp = surfhtemp + temp(i,j,ktmp) * delh
           surfhsalt = surfhsalt + salt(i,j,ktmp) * delh
           surfhu    = surfhu + 0.5*us%L_T_to_m_s*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delh
           surfhv    = surfhv + 0.5*us%L_T_to_m_s*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delh
           if (cs%Stokes_Mixing) then
             surfhus = surfhus + 0.5*(waves%US_x(i,j,ktmp)+waves%US_x(i-1,j,ktmp)) * delh
             surfhvs = surfhvs + 0.5*(waves%US_y(i,j,ktmp)+waves%US_y(i,j-1,ktmp)) * delh
           endif
 
         enddo
         surftemp = surfhtemp / htot
         surfsalt = surfhsalt / htot
         surfu    = surfhu    / htot
         surfv    = surfhv    / htot
         surfus   = surfhus   / htot
         surfvs   = surfhvs   / htot
 
         ! vertical shear between present layer and
         ! surface layer averaged surfU,surfV.
         ! C-grid average to get Uk and Vk on T-points.
         uk         = 0.5*us%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) - surfu
         vk         = 0.5*us%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) - surfv
 
         if (cs%Stokes_Mixing) then
           ! If momentum is mixed down the Stokes drift gradient, then
           !  the Stokes drift must be included in the bulk Richardson number
           !  calculation.
           uk =  uk + (0.5*(waves%Us_x(i,j,k)+waves%US_x(i-1,j,k)) -surfus )
           vk =  vk + (0.5*(waves%Us_y(i,j,k)+waves%Us_y(i,j-1,k)) -surfvs )
         endif
 
         deltau2(k) = uk**2 + vk**2
 
         ! pressure, temp, and saln for EOS
         ! kk+1 = surface fields
         ! kk+2 = k fields
         ! kk+3 = km1 fields
         km1  = max(1, k-1)
         kk   = 3*(k-1)
         pres_1d(kk+1) = pref
         pres_1d(kk+2) = pref
         pres_1d(kk+3) = pref
         temp_1d(kk+1) = surftemp
         temp_1d(kk+2) = temp(i,j,k)
         temp_1d(kk+3) = temp(i,j,km1)
         salt_1d(kk+1) = surfsalt
         salt_1d(kk+2) = salt(i,j,k)
         salt_1d(kk+3) = salt(i,j,km1)
 
         ! pRef is pressure at interface between k and km1 [R L2 T-2 ~> Pa].
         ! iterate pRef for next pass through k-loop.
         pref = pref + (gv%g_Earth * gv%H_to_RZ) * h(i,j,k)
 
         ! this difference accounts for penetrating SW
         surfbuoyflux2(k) = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,k+1))
 
       enddo ! k-loop finishes
 
       if (cs%LT_K_ENHANCEMENT .or. cs%LT_VT2_ENHANCEMENT) then
         mld_guess = max( 1.*us%m_to_Z, abs(us%m_to_Z*cs%OBLdepthprev(i,j) ) )
         call get_langmuir_number(la, g, gv, us, mld_guess, ustar(i,j), i, j, &
                                  h=h(i,j,:), u_h=u_h, v_h=v_h, waves=waves)
         cs%La_SL(i,j)=la
       endif
 
 
       ! compute in-situ density
       call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, tv%eqn_of_state)
 
       ! N2 (can be negative) and N (non-negative) on interfaces.
       ! deltaRho is non-local rho difference used for bulk Richardson number.
       ! CS%N is local N (with floor) used for unresolved shear calculation.
       do k = 1, g%ke
         km1 = max(1, k-1)
         kk = 3*(k-1)
         deltarho(k) = rho_1d(kk+2) - rho_1d(kk+1)
         n2_1d(k)    = (gorho * (rho_1d(kk+2) - rho_1d(kk+3)) ) / &
                       ((0.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_m)
         cs%N(i,j,k)     = sqrt( max( n2_1d(k), 0.) )
       enddo
       n2_1d(g%ke+1 ) = 0.0
       cs%N(i,j,g%ke+1 )  = 0.0
 
       ! turbulent velocity scales w_s and w_m computed at the cell centers.
       ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales
       ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass
       ! sigma=CS%surf_layer_ext for this calculation.
       call cvmix_kpp_compute_turbulent_scales( &
         cs%surf_layer_ext, & ! (in)  Normalized surface layer depth; sigma = CS%surf_layer_ext
         -cellheight,       & ! (in)  Assume here that OBL depth [m] = -cellHeight(k)
         surfbuoyflux2,     & ! (in)  Buoyancy flux at surface [m2 s-3]
         surffricvel,       & ! (in)  Turbulent friction velocity at surface [m s-1]
         w_s=ws_1d,         & ! (out) Turbulent velocity scale profile [m s-1]
         cvmix_kpp_params_user=cs%KPP_params )
 
       !Compute CVMix VT2
       cs%Vt2(i,j,:) = cvmix_kpp_compute_unresolved_shear( &
                       zt_cntr=cellheight(1:g%ke),         & ! Depth of cell center [m]
                       ws_cntr=ws_1d,                      & ! Turbulent velocity scale profile, at centers [m s-1]
                       n_iface=cs%N(i,j,:),                & ! Buoyancy frequency at interface [s-1]
                     cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
 
       !Modify CVMix VT2
       IF (cs%LT_VT2_ENHANCEMENT) then
         IF (cs%LT_VT2_METHOD==lt_vt2_mode_constant) then
           do k=1,g%ke
             langenhvt2(k) = cs%KPP_VT2_ENH_FAC
           enddo
         elseif (cs%LT_VT2_METHOD==lt_vt2_mode_vr12) then
           !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
           enhvt2 = sqrt(1.+(1.5*cs%La_SL(i,j))**(-2) + &
                    (5.4*cs%La_SL(i,j))**(-4))
           do k=1,g%ke
              langenhvt2(k) = enhvt2
           enddo
         elseif (cs%LT_VT2_METHOD==lt_vt2_mode_rw16) then
           !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
           enhvt2 = 1. + 2.3*cs%La_SL(i,j)**(-0.5)
           do k=1,g%ke
             langenhvt2(k) = enhvt2
           enddo
         elseif (cs%LT_VT2_METHOD==lt_vt2_mode_lf17) then
           cs%CS=cvmix_get_kpp_real('c_s',cs%KPP_params)
           do k=1,g%ke
             wst = (max(0.,-buoy_scale*buoyflux(i,j,1))*(-cellheight(k)))**(1./3.)
             langenhvt2(k) = sqrt((0.15*wst**3. + 0.17*surffricvel**3.* &
                  (1.+0.49*cs%La_SL(i,j)**(-2.)))  / &
                  (0.2*ws_1d(k)**3/(cs%cs*cs%surf_layer_ext*cs%vonKarman**4.)))
           enddo
         else
            !This shouldn't be reached.
            !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in Vt2")
            langenhvt2(:) = 1.0
         endif
       else
         langenhvt2(:) = 1.0
       endif
 
       do k=1,g%ke
         cs%Vt2(i,j,k)=cs%Vt2(i,j,k)*langenhvt2(k)
         if (cs%id_EnhVt2 > 0) cs%EnhVt2(i,j,k)=langenhvt2(k)
       enddo
 
       ! Calculate Bulk Richardson number from eq (21) of LMD94
       bulkri_1d = cvmix_kpp_compute_bulk_richardson( &
                   zt_cntr = cellheight(1:g%ke),      & ! Depth of cell center [m]
                   delta_buoy_cntr=gorho*deltarho,    & ! Bulk buoyancy difference, Br-B(z) [s-1]
                   delta_vsqr_cntr=deltau2,           & ! Square of resolved velocity difference [m2 s-2]
                   vt_sqr_cntr=cs%Vt2(i,j,:),         &
                   ws_cntr=ws_1d,                     & ! Turbulent velocity scale profile [m s-1]
                   n_iface=cs%N(i,j,:))               ! Buoyancy frequency [s-1]
 
 
       surfbuoyflux = buoy_scale * buoyflux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
                                      ! h to Monin-Obukov (default is false, ie. not used)
 
       call cvmix_kpp_compute_obl_depth( &
         bulkri_1d,              & ! (in) Bulk Richardson number
         ifaceheight,            & ! (in) Height of interfaces [m]
         cs%OBLdepth(i,j),       & ! (out) OBL depth [m]
         cs%kOBL(i,j),           & ! (out) level (+fraction) of OBL extent
         zt_cntr=cellheight,     & ! (in) Height of cell centers [m]
         surf_fric=surffricvel,  & ! (in) Turbulent friction velocity at surface [m s-1]
         surf_buoy=surfbuoyflux, & ! (in) Buoyancy flux at surface [m2 s-3]
         coriolis=coriolis,      & ! (in) Coriolis parameter [s-1]
         cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
 
       ! A hack to avoid KPP reaching the bottom. It was needed during development
       ! because KPP was unable to handle vanishingly small layers near the bottom.
       if (cs%deepOBLoffset>0.) then
         zbottomminusoffset = ifaceheight(g%ke+1) + min(cs%deepOBLoffset,-0.1*ifaceheight(g%ke+1))
         cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -zbottomminusoffset )
       endif
 
       ! apply some constraints on OBLdepth
       if(cs%fixedOBLdepth)  cs%OBLdepth(i,j) = cs%fixedOBLdepth_value
       cs%OBLdepth(i,j) = max( cs%OBLdepth(i,j), -ifaceheight(2) )      ! no shallower than top layer
       cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -ifaceheight(g%ke+1) ) ! no deeper than bottom
       cs%kOBL(i,j)     = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, cs%OBLdepth(i,j) )
 
 
       ! recompute wscale for diagnostics, now that we in fact know boundary layer depth
       !BGR consider if LTEnhancement is wanted for diagnostics
       if (cs%id_Ws > 0) then
           call cvmix_kpp_compute_turbulent_scales( &
             -cellheight/cs%OBLdepth(i,j),          & ! (in)  Normalized boundary layer coordinate
             cs%OBLdepth(i,j),                      & ! (in)  OBL depth [m]
             surfbuoyflux,                          & ! (in)  Buoyancy flux at surface [m2 s-3]
             surffricvel,                           & ! (in)  Turbulent friction velocity at surface [m s-1]
             w_s=ws_1d,                             & ! (out) Turbulent velocity scale profile [m s-1]
             cvmix_kpp_params_user=cs%KPP_params)     !       KPP parameters
           cs%Ws(i,j,:) = ws_1d(:)
       endif
 
       ! Diagnostics
       if (cs%id_N2     > 0)   cs%N2(i,j,:)     = n2_1d(:)
       if (cs%id_BulkDrho > 0) cs%dRho(i,j,:)   = deltarho(:)
       if (cs%id_BulkRi > 0)   cs%BulkRi(i,j,:) = bulkri_1d(:)
       if (cs%id_BulkUz2 > 0)  cs%Uz2(i,j,:)    = deltau2(:)
       if (cs%id_Tsurf  > 0)   cs%Tsurf(i,j)    = surftemp
       if (cs%id_Ssurf  > 0)   cs%Ssurf(i,j)    = surfsalt
       if (cs%id_Usurf  > 0)   cs%Usurf(i,j)    = surfu
       if (cs%id_Vsurf  > 0)   cs%Vsurf(i,j)    = surfv
 
     enddo
   enddo
 
   call cpu_clock_end(id_clock_kpp_compute_bld)
 
   ! send diagnostics to post_data
   if (cs%id_BulkRi   > 0) call post_data(cs%id_BulkRi,   cs%BulkRi,          cs%diag)
   if (cs%id_N        > 0) call post_data(cs%id_N,        cs%N,               cs%diag)
   if (cs%id_N2       > 0) call post_data(cs%id_N2,       cs%N2,              cs%diag)
   if (cs%id_Tsurf    > 0) call post_data(cs%id_Tsurf,    cs%Tsurf,           cs%diag)
   if (cs%id_Ssurf    > 0) call post_data(cs%id_Ssurf,    cs%Ssurf,           cs%diag)
   if (cs%id_Usurf    > 0) call post_data(cs%id_Usurf,    cs%Usurf,           cs%diag)
   if (cs%id_Vsurf    > 0) call post_data(cs%id_Vsurf,    cs%Vsurf,           cs%diag)
   if (cs%id_BulkDrho > 0) call post_data(cs%id_BulkDrho, cs%dRho,            cs%diag)
   if (cs%id_BulkUz2  > 0) call post_data(cs%id_BulkUz2,  cs%Uz2,             cs%diag)
   if (cs%id_EnhK     > 0) call post_data(cs%id_EnhK,     cs%EnhK,            cs%diag)
   if (cs%id_EnhVt2   > 0) call post_data(cs%id_EnhVt2,   cs%EnhVt2,          cs%diag)
   if (cs%id_La_SL    > 0) call post_data(cs%id_La_SL,    cs%La_SL,           cs%diag)
 
   ! BLD smoothing:
   if (cs%n_smooth > 0) call kpp_smooth_bld(cs,g,gv,h)
 

◆ kpp_end()

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

Clear pointers, deallocate memory.

Parameters

cs	Control structure

Definition at line 1514 of file MOM_CVMix_KPP.F90.

   type(KPP_CS), pointer :: CS !< Control structure
 
   if (.not.associated(cs)) return
 
   deallocate(cs)

◆ kpp_get_bld()

subroutine, public mom_cvmix_kpp::kpp_get_bld	(	type(kpp_cs), pointer	CS,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)	BLD,
		type(ocean_grid_type), intent(in)	G,
		type(unit_scale_type), intent(in)	US,
		real, intent(in), optional	m_to_BLD_units
	)

Copies KPP surface boundary layer depth into BLD, in units of [Z ~> m] unless other units are specified.

Parameters

	cs	Control structure for this module
[in]	g	Grid structure
[in]	us	A dimensional unit scaling type
[in,out]	bld	Boundary layer depth [Z ~> m] or other units
[in]	m_to_bld_units	A conversion factor from meters to the desired units for BLD

Definition at line 1367 of file MOM_CVMix_KPP.F90.

   type(KPP_CS),                     pointer     :: CS  !< Control structure for
                                                        !! this module
   type(ocean_grid_type),            intent(in)  :: G   !< Grid structure
   type(unit_scale_type),            intent(in)  :: US  !< A dimensional unit scaling type
   real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Boundary layer depth [Z ~> m] or other units
   real,                   optional, intent(in)  :: m_to_BLD_units !< A conversion factor from meters
                                                        !! to the desired units for BLD
   ! Local variables
   real :: scale  ! A dimensional rescaling factor
   integer :: i,j
 
   scale = us%m_to_Z ; if (present(m_to_bld_units)) scale = m_to_bld_units
 
   !$OMP parallel do default(none) shared(BLD, CS, G, scale)
   do j = g%jsc, g%jec ; do i = g%isc, g%iec
     bld(i,j) = scale * cs%OBLdepth(i,j)
   enddo ; enddo
 

◆ 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]	paramfile	File parser
[in]	g	Ocean grid
[in]	gv	Vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	diag	Diagnostics
[in]	time	Model time
	cs	Control structure
[out]	passive	Copy of passiveMode
	waves	Wave CS

Definition at line 185 of file MOM_CVMix_KPP.F90.

 
   ! Arguments
   type(param_file_type),   intent(in)    :: paramFile !< File parser
   type(ocean_grid_type),   intent(in)    :: G         !< Ocean grid
   type(verticalGrid_type), intent(in)    :: GV        !< Vertical grid structure.
   type(unit_scale_type),   intent(in)    :: US        !< A dimensional unit scaling type
   type(diag_ctrl), target, intent(in)    :: diag      !< Diagnostics
   type(time_type),         intent(in)    :: Time      !< Model time
   type(KPP_CS),            pointer       :: CS        !< Control structure
   logical,       optional, intent(out)   :: passive   !< Copy of %passiveMode
   type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS
 
   ! Local variables
 # include "version_variable.h"
   character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module
   character(len=20) :: string          !< local temporary string
   logical :: CS_IS_ONE=.false.         !< Logical for setting Cs based on Non-local
   logical :: lnoDGat1=.false.          !< True => G'(1) = 0 (shape function)
                                        !! False => compute G'(1) as in LMD94
   if (associated(cs)) call mom_error(fatal, 'MOM_CVMix_KPP, KPP_init: '// &
            'Control structure has already been initialized')
 
   ! Read parameters
   call get_param(paramfile, mdl, "USE_KPP", kpp_init, default=.false., do_not_log=.true.)
   call log_version(paramfile, mdl, version, 'This is the MOM wrapper to CVMix:KPP\n' // &
             'See http://cvmix.github.io/', all_default=.not.kpp_init)
   call get_param(paramfile, mdl, "USE_KPP", kpp_init, &
                  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "// &
                  "to calculate diffusivities and non-local transport in the OBL.",     &
                  default=.false.)
   ! Forego remainder of initialization if not using this scheme
   if (.not. kpp_init) return
   allocate(cs)
 
   call openparameterblock(paramfile,'KPP')
   call get_param(paramfile, mdl, 'PASSIVE', cs%passiveMode,           &
                  'If True, puts KPP into a passive-diagnostic mode.', &
                   default=.false.)
   !BGR: Note using PASSIVE for KPP creates warning for PASSIVE from Convection
   !     should we create a separate flag?
   if (present(passive)) passive=cs%passiveMode ! This is passed back to the caller so
                                                ! the caller knows to not use KPP output
   call get_param(paramfile, mdl, 'APPLY_NONLOCAL_TRANSPORT', cs%applyNonLocalTrans,  &
                  'If True, applies the non-local transport to heat and scalars. '//  &
                  'If False, calculates the non-local transport and tendencies but '//&
                  'purely for diagnostic purposes.',                                   &
                  default=.not. cs%passiveMode)
   call get_param(paramfile, mdl, 'N_SMOOTH', cs%n_smooth,  &
                  'The number of times the 1-1-4-1-1 Laplacian filter is applied on '//  &
                  'OBL depth.',   &
                  default=0)
   if (cs%n_smooth > g%domain%nihalo) then
     call mom_error(fatal,'KPP smoothing number (N_SMOOTH) cannot be greater than NIHALO.')
   elseif (cs%n_smooth > g%domain%njhalo) then
     call mom_error(fatal,'KPP smoothing number (N_SMOOTH) cannot be greater than NJHALO.')
   endif
   if (cs%n_smooth > 0) then
     call get_param(paramfile, mdl, 'DEEPEN_ONLY_VIA_SMOOTHING', cs%deepen_only,  &
                    'If true, apply OBLdepth smoothing at a cell only if the OBLdepth '// &
                    'gets deeper via smoothing.',   &
                    default=.false.)
     id_clock_kpp_smoothing = cpu_clock_id('(Ocean KPP BLD smoothing)', grain=clock_routine)
   endif
   call get_param(paramfile, mdl, 'RI_CRIT', cs%Ri_crit,                            &
                  'Critical bulk Richardson number used to define depth of the '// &
                  'surface Ocean Boundary Layer (OBL).',                            &
                  units='nondim', default=0.3)
   call get_param(paramfile, mdl, 'VON_KARMAN', cs%vonKarman, &
                  'von Karman constant.',                     &
                  units='nondim', default=0.40)
   call get_param(paramfile, mdl, 'ENHANCE_DIFFUSION', cs%enhance_diffusion,              &
                  'If True, adds enhanced diffusion at the based of the boundary layer.', &
                  default=.true.)
   call get_param(paramfile, mdl, 'INTERP_TYPE', cs%interpType,           &
                  'Type of interpolation to determine the OBL depth.\n'// &
                  'Allowed types are: linear, quadratic, cubic.',         &
                  default='quadratic')
   call get_param(paramfile, mdl, 'INTERP_TYPE2', cs%interpType2,           &
                  'Type of interpolation to compute diff and visc at OBL_depth.\n'// &
                  'Allowed types are: linear, quadratic, cubic or LMD94.',         &
                  default='LMD94')
   call get_param(paramfile, mdl, 'COMPUTE_EKMAN', cs%computeEkman,             &
                  'If True, limit OBL depth to be no deeper than Ekman depth.', &
                  default=.false.)
   call get_param(paramfile, mdl, 'COMPUTE_MONIN_OBUKHOV', cs%computeMoninObukhov, &
                  'If True, limit the OBL depth to be no deeper than '//          &
                  'Monin-Obukhov depth.',                                          &
                  default=.false.)
   call get_param(paramfile, mdl, 'CS', cs%cs,                        &
                  'Parameter for computing velocity scale function.', &
                  units='nondim', default=98.96)
   call get_param(paramfile, mdl, 'CS2', cs%cs2,                        &
                  'Parameter for computing non-local term.', &
                  units='nondim', default=6.32739901508)
   call get_param(paramfile, mdl, 'DEEP_OBL_OFFSET', cs%deepOBLoffset,                             &
                  'If non-zero, the distance above the bottom to which the OBL is clipped '//     &
                  'if it would otherwise reach the bottom. The smaller of this and 0.1D is used.', &
                  units='m',default=0.)
   call get_param(paramfile, mdl, 'FIXED_OBLDEPTH', cs%fixedOBLdepth,       &
                  'If True, fix the OBL depth to FIXED_OBLDEPTH_VALUE '//  &
                  'rather than using the OBL depth from CVMix. '//         &
                  'This option is just for testing purposes.',              &
                  default=.false.)
   call get_param(paramfile, mdl, 'FIXED_OBLDEPTH_VALUE', cs%fixedOBLdepth_value,  &
                  'Value for the fixed OBL depth when fixedOBLdepth==True. '//   &
                  'This parameter is for just for testing purposes. '//          &
                  'It will over-ride the OBLdepth computed from CVMix.',           &
                  units='m',default=30.0)
   call get_param(paramfile, mdl, 'SURF_LAYER_EXTENT', cs%surf_layer_ext,   &
                  'Fraction of OBL depth considered in the surface layer.', &
                  units='nondim',default=0.10)
   call get_param(paramfile, mdl, 'MINIMUM_OBL_DEPTH', cs%minOBLdepth,                            &
                  'If non-zero, a minimum depth to use for KPP OBL depth. Independent of '//     &
                  'this parameter, the OBL depth is always at least as deep as the first layer.', &
                  units='m',default=0.)
   call get_param(paramfile, mdl, 'MINIMUM_VT2', cs%minVtsqr,                                   &
                  'Min of the unresolved velocity Vt2 used in Rib CVMix calculation.\n'//  &
                  'Scaling: MINIMUM_VT2 = const1*d*N*ws, with d=1m, N=1e-5/s, ws=1e-6 m/s.',    &
                  units='m2/s2',default=1e-10)
 
 ! smg: for removal below
   call get_param(paramfile, mdl, 'CORRECT_SURFACE_LAYER_AVERAGE', cs%correctSurfLayerAvg,   &
                  'If true, applies a correction step to the averaging of surface layer '// &
                  'properties. This option is obsolete.', default=.false.)
   if (cs%correctSurfLayerAvg) &
     call mom_error(fatal,'Correct surface layer average disabled in code.  To recover \n'// &
                        ' feature will require code intervention.')
   call get_param(paramfile, mdl, 'FIRST_GUESS_SURFACE_LAYER_DEPTH', cs%surfLayerDepth,              &
                  'The first guess at the depth of the surface layer used for averaging '//         &
                  'the surface layer properties. If =0, the top model level properties '//          &
                  'will be used for the surface layer. If CORRECT_SURFACE_LAYER_AVERAGE=True, a '// &
                  'subsequent correction is applied. This parameter is obsolete', units='m', default=0.)
 ! smg: for removal above
 
   call get_param(paramfile, mdl, 'NLT_SHAPE', string, &
                  'MOM6 method to set nonlocal transport profile. '//                          &
                  'Over-rides the result from CVMix.  Allowed values are: \n'//                 &
                  '\t CVMix     - Uses the profiles from CVMix specified by MATCH_TECHNIQUE\n'//&
                  '\t LINEAR    - A linear profile, 1-sigma\n'//                                &
                  '\t PARABOLIC - A parablic profile, (1-sigma)^2\n'//                          &
                  '\t CUBIC     - A cubic profile, (1-sigma)^2(1+2*sigma)\n'//                  &
                  '\t CUBIC_LMD - The original KPP profile',                                    &
                  default='CVMix')
   select case ( trim(string) )
     case ("CVMix")     ; cs%NLT_shape = nlt_shape_cvmix
     case ("LINEAR")    ; cs%NLT_shape = nlt_shape_linear
     case ("PARABOLIC") ; cs%NLT_shape = nlt_shape_parabolic
     case ("CUBIC")     ; cs%NLT_shape = nlt_shape_cubic
     case ("CUBIC_LMD") ; cs%NLT_shape = nlt_shape_cubic_lmd
     case default ; call mom_error(fatal,"KPP_init: "// &
                    "Unrecognized NLT_SHAPE option"//trim(string))
   end select
   call get_param(paramfile, mdl, 'MATCH_TECHNIQUE', cs%MatchTechnique,                                    &
                  'CVMix method to set profile function for diffusivity and NLT, '//                      &
                  'as well as matching across OBL base. Allowed values are: \n'//                          &
                  '\t SimpleShapes      = sigma*(1-sigma)^2 for both diffusivity and NLT\n'//              &
                  '\t MatchGradient     = sigma*(1-sigma)^2 for NLT; diffusivity profile from matching\n'//&
                  '\t MatchBoth         = match gradient for both diffusivity and NLT\n'//                 &
                  '\t ParabolicNonLocal = sigma*(1-sigma)^2 for diffusivity; (1-sigma)^2 for NLT',         &
                  default='SimpleShapes')
   if (cs%MatchTechnique == 'ParabolicNonLocal') then
     ! This forces Cs2 (Cs in non-local computation) to equal 1 for parabolic non-local option.
     !  May be used during CVMix initialization.
     cs_is_one=.true.
   endif
   if (cs%MatchTechnique == 'ParabolicNonLocal' .or. cs%MatchTechnique == 'SimpleShapes') then
     ! if gradient won't be matched, lnoDGat1=.true.
     lnodgat1=.true.
   endif
 
   ! safety check to avoid negative diff/visc
   if (cs%MatchTechnique == 'MatchBoth' .and. (cs%interpType2 == 'cubic' .or. &
       cs%interpType2 == 'quadratic')) then
     call mom_error(fatal,"If MATCH_TECHNIQUE=MatchBoth, INTERP_TYPE2 must be set to \n"//&
                "linear or LMD94 (recommended) to avoid negative viscosity and diffusivity.\n"//&
                "Please select one of these valid options." )
   endif
 
   call get_param(paramfile, mdl, 'KPP_ZERO_DIFFUSIVITY', cs%KPPzeroDiffusivity,            &
                  'If True, zeroes the KPP diffusivity and viscosity; for testing purpose.',&
                  default=.false.)
   call get_param(paramfile, mdl, 'KPP_IS_ADDITIVE', cs%KPPisAdditive,                &
                  'If true, adds KPP diffusivity to diffusivity from other schemes.\n'//&
                  'If false, KPP is the only diffusivity wherever KPP is non-zero.',  &
                  default=.true.)
   call get_param(paramfile, mdl, 'KPP_SHORTWAVE_METHOD',string,                      &
                  'Determines contribution of shortwave radiation to KPP surface '// &
                  'buoyancy flux.  Options include:\n'//                             &
                  '  ALL_SW: use total shortwave radiation\n'//                      &
                  '  MXL_SW: use shortwave radiation absorbed by mixing layer\n'//  &
                  '  LV1_SW: use shortwave radiation absorbed by top model layer',  &
                  default='MXL_SW')
   select case ( trim(string) )
     case ("ALL_SW") ; cs%SW_METHOD = sw_method_all_sw
     case ("MXL_SW") ; cs%SW_METHOD = sw_method_mxl_sw
     case ("LV1_SW") ; cs%SW_METHOD = sw_method_lv1_sw
     case default ; call mom_error(fatal,"KPP_init: "// &
                    "Unrecognized KPP_SHORTWAVE_METHOD option"//trim(string))
   end select
   call get_param(paramfile, mdl, 'CVMix_ZERO_H_WORK_AROUND', cs%min_thickness,                           &
                  'A minimum thickness used to avoid division by small numbers in the vicinity '//       &
                  'of vanished layers. This is independent of MIN_THICKNESS used in other parts of MOM.', &
                  units='m', default=0.)
 
 !/BGR: New options for including Langmuir effects
 !/ 1. Options related to enhancing the mixing coefficient
   call get_param(paramfile, mdl, "USE_KPP_LT_K", cs%LT_K_Enhancement, &
        'Flag for Langmuir turbulence enhancement of turbulent'//&
        'mixing coefficient.', units="", default=.false.)
   call get_param(paramfile, mdl, "STOKES_MIXING", cs%STOKES_MIXING, &
        'Flag for Langmuir turbulence enhancement of turbulent'//&
        'mixing coefficient.', units="", default=.false.)
   if (cs%LT_K_Enhancement) then
     call get_param(paramfile, mdl, 'KPP_LT_K_SHAPE', string,                 &
                  'Vertical dependence of LT enhancement of mixing. '//     &
                  'Valid options are: \n'//                                   &
                  '\t CONSTANT = Constant value for full OBL\n'//             &
                  '\t SCALED   = Varies based on normalized shape function.', &
                  default='CONSTANT')
     select case ( trim(string))
       case ("CONSTANT") ; cs%LT_K_SHAPE = lt_k_constant
       case ("SCALED")   ; cs%LT_K_SHAPE = lt_k_scaled
       case default ; call mom_error(fatal,"KPP_init: "//&
                     "Unrecognized KPP_LT_K_SHAPE option: "//trim(string))
     end select
     call get_param(paramfile, mdl, "KPP_LT_K_METHOD", string ,                   &
                    'Method to enhance mixing coefficient in KPP. '//           &
                    'Valid options are: \n'//                                     &
                    '\t CONSTANT = Constant value (KPP_K_ENH_FAC) \n'//           &
                    '\t VR12     = Function of Langmuir number based on VR12\n'// &
                    '\t RW16     = Function of Langmuir number based on RW16',    &
                    default='CONSTANT')
     select case ( trim(string))
       case ("CONSTANT") ; cs%LT_K_METHOD = lt_k_mode_constant
       case ("VR12")     ; cs%LT_K_METHOD = lt_k_mode_vr12
       case ("RW16")     ; cs%LT_K_METHOD = lt_k_mode_rw16
       case default      ; call mom_error(fatal,"KPP_init: "//&
                     "Unrecognized KPP_LT_K_METHOD option: "//trim(string))
     end select
     if (cs%LT_K_METHOD==lt_k_mode_constant) then
       call get_param(paramfile, mdl, "KPP_K_ENH_FAC",cs%KPP_K_ENH_FAC ,     &
                    'Constant value to enhance mixing coefficient in KPP.',  &
                    default=1.0)
     endif
   endif
 !/ 2. Options related to enhancing the unresolved Vt2/entrainment in Rib
   call get_param(paramfile, mdl, "USE_KPP_LT_VT2", cs%LT_Vt2_Enhancement, &
        'Flag for Langmuir turbulence enhancement of Vt2'//&
        'in Bulk Richardson Number.', units="", default=.false.)
   if (cs%LT_Vt2_Enhancement) then
     call get_param(paramfile, mdl, "KPP_LT_VT2_METHOD",string ,                  &
                    'Method to enhance Vt2 in KPP. '//                          &
                    'Valid options are: \n'//                                     &
                    '\t CONSTANT = Constant value (KPP_VT2_ENH_FAC) \n'//         &
                    '\t VR12     = Function of Langmuir number based on VR12\n'// &
                    '\t RW16     = Function of Langmuir number based on RW16\n'// &
                    '\t LF17     = Function of Langmuir number based on LF17',    &
                    default='CONSTANT')
     select case ( trim(string))
       case ("CONSTANT") ; cs%LT_VT2_METHOD = lt_vt2_mode_constant
       case ("VR12")     ; cs%LT_VT2_METHOD = lt_vt2_mode_vr12
       case ("RW16")     ; cs%LT_VT2_METHOD = lt_vt2_mode_rw16
       case ("LF17")     ; cs%LT_VT2_METHOD = lt_vt2_mode_lf17
       case default      ; call mom_error(fatal,"KPP_init: "//&
                     "Unrecognized KPP_LT_VT2_METHOD option: "//trim(string))
     end select
     if (cs%LT_VT2_METHOD==lt_vt2_mode_constant) then
       call get_param(paramfile, mdl, "KPP_VT2_ENH_FAC",cs%KPP_VT2_ENH_FAC ,     &
                    'Constant value to enhance VT2 in KPP.',  &
                    default=1.0)
     endif
   endif
 
   call closeparameterblock(paramfile)
   call get_param(paramfile, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)
 
   call cvmix_init_kpp( ri_crit=cs%Ri_crit,                 &
                        minobldepth=cs%minOBLdepth,         &
                        minvtsqr=cs%minVtsqr,               &
                        vonkarman=cs%vonKarman,             &
                        surf_layer_ext=cs%surf_layer_ext,   &
                        interp_type=cs%interpType,          &
                        interp_type2=cs%interpType2,        &
                        lekman=cs%computeEkman,             &
                        lmonob=cs%computeMoninObukhov,      &
                        matchtechnique=cs%MatchTechnique,   &
                        lenhanced_diff=cs%enhance_diffusion,&
                        lnonzero_surf_nonlocal=cs_is_one   ,&
                        lnodgat1=lnodgat1                  ,&
                        cvmix_kpp_params_user=cs%KPP_params )
 
   ! Register diagnostics
   cs%diag => diag
   cs%id_OBLdepth = register_diag_field('ocean_model', 'KPP_OBLdepth', diag%axesT1, time, &
       'Thickness of the surface Ocean Boundary Layer calculated by [CVMix] KPP', 'meter', &
       cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
       cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
       ! CMOR names are placeholders; must be modified by time period
       ! for CMOR compliance. Diag manager will be used for omlmax and
       ! omldamax.
   if (cs%n_smooth > 0) then
     cs%id_OBLdepth_original = register_diag_field('ocean_model', 'KPP_OBLdepth_original', diag%axesT1, time, &
         'Thickness of the surface Ocean Boundary Layer without smoothing calculated by [CVMix] KPP', 'meter', &
         cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
         cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
   endif
   cs%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, time, &
       'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', &
       'kg/m3', conversion=us%R_to_kg_m3)
   cs%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, time, &
       'Square of bulk difference in resolved velocity used in Bulk Richardson number via [CVMix] KPP', 'm2/s2')
   cs%id_BulkRi = register_diag_field('ocean_model', 'KPP_BulkRi', diag%axesTL, time, &
       'Bulk Richardson number used to find the OBL depth used by [CVMix] KPP', 'nondim')
   cs%id_Sigma = register_diag_field('ocean_model', 'KPP_sigma', diag%axesTi, time, &
       'Sigma coordinate used by [CVMix] KPP', 'nondim')
   cs%id_Ws = register_diag_field('ocean_model', 'KPP_Ws', diag%axesTL, time, &
       'Turbulent vertical velocity scale for scalars used by [CVMix] KPP', 'm/s')
   cs%id_N = register_diag_field('ocean_model', 'KPP_N', diag%axesTi, time, &
       '(Adjusted) Brunt-Vaisala frequency used by [CVMix] KPP', '1/s')
   cs%id_N2 = register_diag_field('ocean_model', 'KPP_N2', diag%axesTi, time, &
       'Square of Brunt-Vaisala frequency used by [CVMix] KPP', '1/s2')
   cs%id_Vt2 = register_diag_field('ocean_model', 'KPP_Vt2', diag%axesTL, time, &
       'Unresolved shear turbulence used by [CVMix] KPP', 'm2/s2')
   cs%id_uStar = register_diag_field('ocean_model', 'KPP_uStar', diag%axesT1, time, &
       'Friction velocity, u*, as used by [CVMix] KPP', 'm/s', conversion=us%Z_to_m*us%s_to_T)
   cs%id_buoyFlux = register_diag_field('ocean_model', 'KPP_buoyFlux', diag%axesTi, time, &
       'Surface (and penetrating) buoyancy flux, as used by [CVMix] KPP', 'm2/s3', conversion=us%L_to_m**2*us%s_to_T**3)
   cs%id_QminusSW = register_diag_field('ocean_model', 'KPP_QminusSW', diag%axesT1, time, &
       'Net temperature flux ignoring short-wave, as used by [CVMix] KPP', 'K m/s')
   cs%id_netS = register_diag_field('ocean_model', 'KPP_netSalt', diag%axesT1, time, &
       'Effective net surface salt flux, as used by [CVMix] KPP', 'ppt m/s')
   cs%id_Kt_KPP = register_diag_field('ocean_model', 'KPP_Kheat', diag%axesTi, time, &
       'Heat diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
   cs%id_Kd_in = register_diag_field('ocean_model', 'KPP_Kd_in', diag%axesTi, time, &
       'Diffusivity passed to KPP', 'm2/s', conversion=us%Z2_T_to_m2_s)
   cs%id_Ks_KPP = register_diag_field('ocean_model', 'KPP_Ksalt', diag%axesTi, time, &
       'Salt diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
   cs%id_Kv_KPP = register_diag_field('ocean_model', 'KPP_Kv', diag%axesTi, time, &
       'Vertical viscosity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
   cs%id_NLTt = register_diag_field('ocean_model', 'KPP_NLtransport_heat', diag%axesTi, time, &
       'Non-local transport (Cs*G(sigma)) for heat, as calculated by [CVMix] KPP', 'nondim')
   cs%id_NLTs = register_diag_field('ocean_model', 'KPP_NLtransport_salt', diag%axesTi, time, &
       'Non-local tranpsort (Cs*G(sigma)) for scalars, as calculated by [CVMix] KPP', 'nondim')
   cs%id_NLT_dTdt = register_diag_field('ocean_model', 'KPP_NLT_dTdt', diag%axesTL, time, &
       'Temperature tendency due to non-local transport of heat, as calculated by [CVMix] KPP', 'K/s')
   cs%id_NLT_dSdt = register_diag_field('ocean_model', 'KPP_NLT_dSdt', diag%axesTL, time, &
       'Salinity tendency due to non-local transport of salt, as calculated by [CVMix] KPP', 'ppt/s')
   cs%id_NLT_temp_budget = register_diag_field('ocean_model', 'KPP_NLT_temp_budget', diag%axesTL, time, &
       'Heat content change due to non-local transport, as calculated by [CVMix] KPP', 'W/m^2')
   cs%id_NLT_saln_budget = register_diag_field('ocean_model', 'KPP_NLT_saln_budget', diag%axesTL, time, &
       'Salt content change due to non-local transport, as calculated by [CVMix] KPP', 'kg/(sec*m^2)')
   cs%id_Tsurf = register_diag_field('ocean_model', 'KPP_Tsurf', diag%axesT1, time, &
       'Temperature of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'C')
   cs%id_Ssurf = register_diag_field('ocean_model', 'KPP_Ssurf', diag%axesT1, time, &
       'Salinity of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'ppt')
   cs%id_Usurf = register_diag_field('ocean_model', 'KPP_Usurf', diag%axesCu1, time, &
       'i-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'm/s')
   cs%id_Vsurf = register_diag_field('ocean_model', 'KPP_Vsurf', diag%axesCv1, time, &
       'j-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'm/s')
   cs%id_EnhK = register_diag_field('ocean_model', 'EnhK', diag%axesTI, time, &
       'Langmuir number enhancement to K as used by [CVMix] KPP','nondim')
   cs%id_EnhVt2 = register_diag_field('ocean_model', 'EnhVt2', diag%axesTL, time, &
       'Langmuir number enhancement to Vt2 as used by [CVMix] KPP','nondim')
   cs%id_La_SL = register_diag_field('ocean_model', 'KPP_La_SL', diag%axesT1, time, &
       'Surface-layer Langmuir number computed in [CVMix] KPP','nondim')
 
   allocate( cs%N( szi_(g), szj_(g), szk_(g)+1 ) )
   cs%N(:,:,:) = 0.
   allocate( cs%OBLdepth( szi_(g), szj_(g) ) )
   cs%OBLdepth(:,:) = 0.
   allocate( cs%kOBL( szi_(g), szj_(g) ) )
   cs%kOBL(:,:) = 0.
   allocate( cs%La_SL( szi_(g), szj_(g) ) )
   cs%La_SL(:,:) = 0.
   allocate( cs%Vt2( szi_(g), szj_(g), szk_(g) ) )
   cs%Vt2(:,:,:) = 0.
   if (cs%id_OBLdepth_original > 0) allocate( cs%OBLdepth_original( szi_(g), szj_(g) ) )
 
   allocate( cs%OBLdepthprev( szi_(g), szj_(g) ) ) ; cs%OBLdepthprev(:,:) = 0.0
   if (cs%id_BulkDrho > 0) allocate( cs%dRho( szi_(g), szj_(g), szk_(g) ) )
   if (cs%id_BulkDrho > 0) cs%dRho(:,:,:) = 0.
   if (cs%id_BulkUz2 > 0)  allocate( cs%Uz2( szi_(g), szj_(g), szk_(g) ) )
   if (cs%id_BulkUz2 > 0)  cs%Uz2(:,:,:) = 0.
   if (cs%id_BulkRi > 0)   allocate( cs%BulkRi( szi_(g), szj_(g), szk_(g) ) )
   if (cs%id_BulkRi > 0)   cs%BulkRi(:,:,:) = 0.
   if (cs%id_Sigma > 0)    allocate( cs%sigma( szi_(g), szj_(g), szk_(g)+1 ) )
   if (cs%id_Sigma > 0)    cs%sigma(:,:,:) = 0.
   if (cs%id_Ws > 0)       allocate( cs%Ws( szi_(g), szj_(g), szk_(g) ) )
   if (cs%id_Ws > 0)       cs%Ws(:,:,:) = 0.
   if (cs%id_N2 > 0)       allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) )
   if (cs%id_N2 > 0)       cs%N2(:,:,:) = 0.
   if (cs%id_Kt_KPP > 0)   allocate( cs%Kt_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
   if (cs%id_Kt_KPP > 0)   cs%Kt_KPP(:,:,:) = 0.
   if (cs%id_Ks_KPP > 0)   allocate( cs%Ks_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
   if (cs%id_Ks_KPP > 0)   cs%Ks_KPP(:,:,:) = 0.
   if (cs%id_Kv_KPP > 0)   allocate( cs%Kv_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
   if (cs%id_Kv_KPP > 0)   cs%Kv_KPP(:,:,:) = 0.
   if (cs%id_Tsurf > 0)    allocate( cs%Tsurf( szi_(g), szj_(g)) )
   if (cs%id_Tsurf > 0)    cs%Tsurf(:,:) = 0.
   if (cs%id_Ssurf > 0)    allocate( cs%Ssurf( szi_(g), szj_(g)) )
   if (cs%id_Ssurf > 0)    cs%Ssurf(:,:) = 0.
   if (cs%id_Usurf > 0)    allocate( cs%Usurf( szib_(g), szj_(g)) )
   if (cs%id_Usurf > 0)    cs%Usurf(:,:) = 0.
   if (cs%id_Vsurf > 0)    allocate( cs%Vsurf( szi_(g), szjb_(g)) )
   if (cs%id_Vsurf > 0)    cs%Vsurf(:,:) = 0.
   if (cs%id_EnhVt2 > 0)    allocate( cs%EnhVt2( szi_(g), szj_(g), szk_(g)) )
   if (cs%id_EnhVt2 > 0)    cs%EnhVt2(:,:,:) = 0.
   if (cs%id_EnhK > 0)    allocate( cs%EnhK( szi_(g), szj_(g), szk_(g)+1 ) )
   if (cs%id_EnhK > 0)    cs%EnhK(:,:,:) = 0.
 
   id_clock_kpp_calc = cpu_clock_id('Ocean KPP calculate)', grain=clock_module)
   id_clock_kpp_compute_bld = cpu_clock_id('(Ocean KPP comp BLD)', grain=clock_routine)
 

◆ 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( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(in)	nonLocalTrans,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)	surfFlux,
		real, intent(in)	dt,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), 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]	cs	Control structure
[in]	g	Ocean grid
[in]	gv	Ocean vertical grid
[in]	h	Layer/level thickness [H ~> m or kg m-2]
[in]	nonlocaltrans	Non-local transport [nondim]
[in]	surfflux	Surface flux of scalar [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
[in]	dt	Time-step [s]
[in,out]	scalar	Scalar (scalar units [conc])

Definition at line 1452 of file MOM_CVMix_KPP.F90.

 
   type(KPP_CS),                               intent(in)    :: CS            !< Control structure
   type(ocean_grid_type),                      intent(in)    :: G             !< Ocean grid
   type(verticalGrid_type),                    intent(in)    :: GV            !< Ocean vertical grid
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)    :: h             !< Layer/level thickness [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: nonLocalTrans !< Non-local transport [nondim]
   real, dimension(SZI_(G),SZJ_(G)),           intent(in)    :: surfFlux      !< Surface flux of scalar
                                                                         !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
   real,                                       intent(in)    :: dt            !< Time-step [s]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: scalar        !< Scalar (scalar units [conc])
 
   integer :: i, j, k
   real, dimension( SZI_(G), SZJ_(G), SZK_(G) ) :: dtracer
 
 
   dtracer(:,:,:) = 0.0
   !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux)
   do k = 1, g%ke
     do j = g%jsc, g%jec
       do i = g%isc, g%iec
         dtracer(i,j,k) = ( nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1) ) / &
                          ( h(i,j,k) + gv%H_subroundoff ) * surfflux(i,j)
       enddo
     enddo
   enddo
 
   !  Update tracer due to non-local redistribution of surface flux
   if (cs%applyNonLocalTrans) then
     !$OMP parallel do default(none) shared(G, dt, scalar, dtracer)
     do k = 1, g%ke
       do j = g%jsc, g%jec
         do i = g%isc, g%iec
           scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
         enddo
       enddo
     enddo
   endif
 
   ! Diagnostics
   if (cs%id_netS            > 0) call post_data(cs%id_netS,     surfflux, cs%diag)
   if (cs%id_NLT_dSdt        > 0) call post_data(cs%id_NLT_dSdt, dtracer,  cs%diag)
   if (cs%id_NLT_saln_budget > 0) then
     dtracer(:,:,:) = 0.0
     !$OMP parallel do default(none) shared(G, GV, dtracer, nonLocalTrans, surfFlux)
     do k = 1, g%ke
       do j = g%jsc, g%jec
         do i = g%isc, g%iec
           dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * &
                            surfflux(i,j) * gv%H_to_kg_m2
         enddo
       enddo
     enddo
     call post_data(cs%id_NLT_saln_budget, dtracer, cs%diag)
   endif
 

◆ 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( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(in)	nonLocalTrans,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)	surfFlux,
		real, intent(in)	dt,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)	scalar,
		real, intent(in)	C_p
	)

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

Parameters

[in]	cs	Control structure
[in]	g	Ocean grid
[in]	gv	Ocean vertical grid
[in]	h	Layer/level thickness [H ~> m or kg m-2]
[in]	nonlocaltrans	Non-local transport [nondim]
[in]	surfflux	Surface flux of scalar [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
[in]	dt	Time-step [s]
[in,out]	scalar	temperature
[in]	c_p	Seawater specific heat capacity [J kg-1 degC-1]

Definition at line 1390 of file MOM_CVMix_KPP.F90.

 
   type(KPP_CS),                               intent(in)    :: CS     !< Control structure
   type(ocean_grid_type),                      intent(in)    :: G      !< Ocean grid
   type(verticalGrid_type),                    intent(in)    :: GV     !< Ocean vertical grid
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)    :: h      !< Layer/level thickness [H ~> m or kg m-2]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: nonLocalTrans !< Non-local transport [nondim]
   real, dimension(SZI_(G),SZJ_(G)),           intent(in)    :: surfFlux  !< Surface flux of scalar
                                                                       !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
   real,                                       intent(in)    :: dt     !< Time-step [s]
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: scalar !< temperature
   real,                                       intent(in)    :: C_p    !< Seawater specific heat capacity [J kg-1 degC-1]
 
   integer :: i, j, k
   real, dimension( SZI_(G), SZJ_(G), SZK_(G) ) :: dtracer
 
 
   dtracer(:,:,:) = 0.0
   !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux)
   do k = 1, g%ke
     do j = g%jsc, g%jec
       do i = g%isc, g%iec
         dtracer(i,j,k) = ( nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1) ) / &
                          ( h(i,j,k) + gv%H_subroundoff ) * surfflux(i,j)
       enddo
     enddo
   enddo
 
   !  Update tracer due to non-local redistribution of surface flux
   if (cs%applyNonLocalTrans) then
     !$OMP parallel do default(none) shared(dt, scalar, dtracer, G)
     do k = 1, g%ke
       do j = g%jsc, g%jec
         do i = g%isc, g%iec
           scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
         enddo
       enddo
     enddo
   endif
 
   ! Diagnostics
   if (cs%id_QminusSW        > 0) call post_data(cs%id_QminusSW, surfflux, cs%diag)
   if (cs%id_NLT_dTdt        > 0) call post_data(cs%id_NLT_dTdt, dtracer,  cs%diag)
   if (cs%id_NLT_temp_budget > 0) then
     dtracer(:,:,:) = 0.0
     !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, surfFlux, C_p, G, GV)
     do k = 1, g%ke
       do j = g%jsc, g%jec
         do i = g%isc, g%iec
           dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * &
                            surfflux(i,j) * c_p * gv%H_to_kg_m2
         enddo
       enddo
     enddo
     call post_data(cs%id_NLT_temp_budget, dtracer, cs%diag)
   endif
 

◆ 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( g %isd: g %ied, g %jsd: g %jed, g %ke), 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

	cs	Control structure
[in,out]	g	Ocean grid
[in]	gv	Ocean vertical grid
[in]	h	Layer/level thicknesses [H ~> m or kg m-2]

Definition at line 1285 of file MOM_CVMix_KPP.F90.

   ! Arguments
   type(KPP_CS),                           pointer       :: CS   !< Control structure
   type(ocean_grid_type),                  intent(inout) :: G    !< Ocean grid
   type(verticalGrid_type),                intent(in)    :: GV   !< Ocean vertical grid
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h    !< Layer/level thicknesses [H ~> m or kg m-2]
 
   ! local
   real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_prev     ! OBLdepth before s.th smoothing iteration
   real, dimension( G%ke )          :: cellHeight        ! Cell center heights referenced to surface [m]
                                                         ! (negative in the ocean)
   real, dimension( G%ke+1 )        :: iFaceHeight       ! Interface heights referenced to surface [m]
                                                         ! (negative in the ocean)
   real :: wc, ww, we, wn, ws ! averaging weights for smoothing
   real :: dh                 ! The local thickness used for calculating interface positions [m]
   real :: hcorr              ! A cumulative correction arising from inflation of vanished layers [m]
   integer :: i, j, k, s
 
   call cpu_clock_begin(id_clock_kpp_smoothing)
 
   ! Update halos
   call pass_var(cs%OBLdepth, g%Domain, halo=cs%n_smooth)
 
   if (cs%id_OBLdepth_original > 0) cs%OBLdepth_original = cs%OBLdepth
 
   do s=1,cs%n_smooth
 
     obldepth_prev = cs%OBLdepth
 
     ! apply smoothing on OBL depth
     !$OMP parallel do default(none) shared(G, GV, CS, h, OBLdepth_prev) &
     !$OMP                           private(wc, ww, we, wn, ws, dh, hcorr, cellHeight, iFaceHeight)
     do j = g%jsc, g%jec
       do i = g%isc, g%iec
 
          ! skip land points
         if (g%mask2dT(i,j)==0.) cycle
 
         ifaceheight(1) = 0.0 ! BBL is all relative to the surface
         hcorr = 0.
         do k=1,g%ke
 
           ! cell center and cell bottom in meters (negative values in the ocean)
           dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
           dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
           hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
           dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
           cellheight(k)    = ifaceheight(k) - 0.5 * dh
           ifaceheight(k+1) = ifaceheight(k) - dh
         enddo
 
         ! compute weights
         ww = 0.125 * g%mask2dT(i-1,j)
         we = 0.125 * g%mask2dT(i+1,j)
         ws = 0.125 * g%mask2dT(i,j-1)
         wn = 0.125 * g%mask2dT(i,j+1)
         wc = 1.0 - (ww+we+wn+ws)
 
         cs%OBLdepth(i,j) =  wc * obldepth_prev(i,j)   &
                           + ww * obldepth_prev(i-1,j) &
                           + we * obldepth_prev(i+1,j) &
                           + ws * obldepth_prev(i,j-1) &
                           + wn * obldepth_prev(i,j+1)
 
         ! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing.
         if (cs%deepen_only) cs%OBLdepth(i,j) = max(cs%OBLdepth(i,j), obldepth_prev(i,j))
 
         ! prevent OBL depths deeper than the bathymetric depth
         cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -ifaceheight(g%ke+1) ) ! no deeper than bottom
         cs%kOBL(i,j)     = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, cs%OBLdepth(i,j) )
       enddo
     enddo
 
   enddo ! s-loop
 
   call cpu_clock_end(id_clock_kpp_smoothing)
 

Detailed Description

Data Types

Functions/Subroutines

Variables

Function/Subroutine Documentation

◆ kpp_calculate()

◆ kpp_compute_bld()

◆ kpp_end()

◆ kpp_get_bld()

◆ kpp_init()

◆ kpp_nonlocaltransport_saln()

◆ kpp_nonlocaltransport_temp()

◆ kpp_smooth_bld()