diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index e686261fdf..19f3d87429 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -83,7 +83,8 @@ subroutine find_obsolete_params(param_file) call obsolete_real(param_file, "ETA_TOLERANCE_AUX", only_warn=.true.) call obsolete_real(param_file, "BT_MASS_SOURCE_LIMIT", 0.0) - + call obsolete_real(param_file, "FIRST_GUESS_SURFACE_LAYER_DEPTH") + call obsolete_logical(param_file, "CORRECT_SURFACE_LAYER_AVERAGE") call obsolete_int(param_file, "SEAMOUNT_LENGTH_SCALE", hint="Use SEAMOUNT_X_LENGTH_SCALE instead.") call obsolete_logical(param_file, "MSTAR_FIXED", hint="Instead use MSTAR_MODE.") diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 1e068509b1..d73bba1551 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -73,10 +73,10 @@ module MOM_CVMix_KPP type, public :: KPP_CS ; private ! Parameters - real :: Ri_crit !< Critical bulk Richardson number (defines OBL depth) - real :: vonKarman !< von Karman constant (dimensionless) - real :: cs !< Parameter for computing velocity scale function (dimensionless) - real :: cs2 !< Parameter for multiplying by non-local term + real :: Ri_crit !< Critical bulk Richardson number (defines OBL depth) [nondim] + real :: vonKarman !< von Karman constant (dimensionless) [nondim] + real :: cs !< Parameter for computing velocity scale function (dimensionless) [nondim] + real :: cs2 !< Parameter for multiplying by non-local term [nondim] ! This is active for NLT_SHAPE_CUBIC_LMD only logical :: enhance_diffusion !< If True, add enhanced diffusivity at base of boundary layer. character(len=32) :: interpType !< Type of interpolation to compute bulk Richardson number @@ -85,12 +85,13 @@ module MOM_CVMix_KPP logical :: computeMoninObukhov !< If True, compute Monin-Obukhov limit for OBLdepth logical :: passiveMode !< If True, makes KPP passive meaning it does NOT alter the diffusivity real :: deepOBLoffset !< If non-zero, is a distance from the bottom that the OBL can not - !! penetrate through [m] - real :: minOBLdepth !< If non-zero, is a minimum depth for the OBL [m] + !! penetrate through [Z ~> m] + real :: minOBLdepth !< If non-zero, is a minimum depth for the OBL [Z ~> m] real :: surf_layer_ext !< Fraction of OBL depth considered in the surface layer [nondim] - real :: minVtsqr !< Min for the squared unresolved velocity used in Rib CVMix calculation [m2 s-2] + real :: minVtsqr !< Min for the squared unresolved velocity used in Rib CVMix + !! calculation [L2 T-2 ~> m2 s-2] logical :: fixedOBLdepth !< If True, will fix the OBL depth at fixedOBLdepth_value - real :: fixedOBLdepth_value !< value for the fixed OBL depth when fixedOBLdepth==True. + real :: fixedOBLdepth_value !< value for the fixed OBL depth when fixedOBLdepth==True [Z ~> m] logical :: debug !< If True, calculate checksums and write debugging information character(len=30) :: MatchTechnique !< Method used in CVMix for setting diffusivity and NLT profile functions integer :: NLT_shape !< MOM6 over-ride of CVMix NLT shape function @@ -103,21 +104,17 @@ module MOM_CVMix_KPP !! If False, will replace initial diffusivity wherever KPP diffusivity !! is non-zero. real :: min_thickness !< A minimum thickness used to avoid division by small numbers - !! in the vicinity of vanished layers. - ! smg: obsolete below - logical :: correctSurfLayerAvg !< If true, applies a correction to the averaging of surface layer properties - real :: surfLayerDepth !< A guess at the depth of the surface layer (which should 0.1 of OBLdepth) [m] - ! smg: obsolete above + !! in the vicinity of vanished layers [Z ~> m] integer :: SW_METHOD !< Sets method for using shortwave radiation in surface buoyancy flux logical :: LT_K_Enhancement !< Flags if enhancing mixing coefficients due to LT integer :: LT_K_Shape !< Integer for constant or shape function enhancement integer :: LT_K_Method !< Integer for mixing coefficients LT method - real :: KPP_K_ENH_FAC !< Factor to multiply by K if Method is CONSTANT + real :: KPP_K_ENH_FAC !< Factor to multiply by K if Method is CONSTANT [nondim] logical :: LT_Vt2_Enhancement !< Flags if enhancing Vt2 due to LT integer :: LT_VT2_METHOD !< Integer for Vt2 LT method - real :: KPP_VT2_ENH_FAC !< Factor to multiply by VT2 if Method is CONSTANT + real :: KPP_VT2_ENH_FAC !< Factor to multiply by VT2 if Method is CONSTANT [nondim] logical :: STOKES_MIXING !< Flag if model is mixing down Stokes gradient - !! This is relavent for which current to use in RiB + !! This is relevant for which current to use in RiB !> CVMix parameters type(CVMix_kpp_params_type), pointer :: KPP_params => NULL() @@ -143,15 +140,15 @@ module MOM_CVMix_KPP !>@} ! Diagnostics arrays - real, allocatable, dimension(:,:) :: OBLdepth !< Depth (positive) of OBL [m] + real, allocatable, dimension(:,:) :: OBLdepth !< Depth (positive) of ocean boundary layer (OBL) [m] real, allocatable, dimension(:,:) :: OBLdepth_original !< Depth (positive) of OBL [m] without smoothing - real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent + real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent [nondim] real, allocatable, dimension(:,:) :: OBLdepthprev !< previous Depth (positive) of OBL [m] - real, allocatable, dimension(:,:) :: La_SL !< Langmuir number used in KPP + real, allocatable, dimension(:,:) :: La_SL !< Langmuir number used in KPP [nondim] real, allocatable, dimension(:,:,:) :: dRho !< Bulk difference in density [R ~> kg m-3] real, allocatable, dimension(:,:,:) :: Uz2 !< Square of bulk difference in resolved velocity [m2 s-2] - real, allocatable, dimension(:,:,:) :: BulkRi !< Bulk Richardson number for each layer (dimensionless) - real, allocatable, dimension(:,:,:) :: sigma !< Sigma coordinate (dimensionless) + real, allocatable, dimension(:,:,:) :: BulkRi !< Bulk Richardson number for each layer [nondim] + real, allocatable, dimension(:,:,:) :: sigma !< Sigma coordinate (dimensionless) [nondim] real, allocatable, dimension(:,:,:) :: Ws !< Turbulent velocity scale for scalars [m s-1] real, allocatable, dimension(:,:,:) :: N !< Brunt-Vaisala frequency [s-1] real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency [s-2] @@ -161,10 +158,10 @@ module MOM_CVMix_KPP real, allocatable, dimension(:,:,:) :: Kv_KPP !< Viscosity due to KPP [m2 s-1] real, allocatable, dimension(:,:) :: Tsurf !< Temperature of surface layer [C ~> degC] real, allocatable, dimension(:,:) :: Ssurf !< Salinity of surface layer [S ~> ppt] - real, allocatable, dimension(:,:) :: Usurf !< i-velocity of surface layer [m s-1] - real, allocatable, dimension(:,:) :: Vsurf !< j-velocity of surface layer [m s-1] - real, allocatable, dimension(:,:,:) :: EnhK !< Enhancement for mixing coefficient - real, allocatable, dimension(:,:,:) :: EnhVt2 !< Enhancement for Vt2 + real, allocatable, dimension(:,:) :: Usurf !< i-velocity of surface layer [L T-1 ~> m s-1] + real, allocatable, dimension(:,:) :: Vsurf !< j-velocity of surface layer [L T-1 ~> m s-1] + real, allocatable, dimension(:,:,:) :: EnhK !< Enhancement for mixing coefficient [nondim] + real, allocatable, dimension(:,:,:) :: EnhVt2 !< Enhancement for Vt2 [nondim] end type KPP_CS @@ -194,8 +191,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) # include "version_variable.h" character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module character(len=20) :: string !< local temporary string - character(len=20) :: langmuir_mixing_opt = 'NONE' !< langmuir mixing opt to be passed to CVMix, e.g., LWF16 - character(len=20) :: langmuir_entrainment_opt = 'NONE' !< langmuir entrainment opt to be passed to CVMix, e.g., LWF16 + character(len=20) :: langmuir_mixing_opt = 'NONE' !< Langmuir mixing option to be passed to CVMix, e.g., LWF16 + character(len=20) :: langmuir_entrainment_opt = 'NONE' !< Langmuir entrainment option to be + !! passed to CVMix, e.g., LWF16 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 @@ -228,8 +226,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) '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.', & + '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.') @@ -277,7 +274,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) 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.) + units='m', default=0., scale=US%m_to_Z) 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. '// & @@ -287,32 +284,18 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) '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) + units='m', default=30.0, scale=US%m_to_Z) 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) + 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.) + units='m', default=0., scale=US%m_to_Z) 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 + units='m2/s2', default=1e-10, scale=US%m_s_to_L_T**2) call get_param(paramFile, mdl, 'NLT_SHAPE', string, & 'MOM6 method to set nonlocal transport profile. '// & @@ -382,7 +365,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) 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.) + units='m', default=0., scale=US%m_to_Z) !/BGR: New options for including Langmuir effects !/ 1. Options related to enhancing the mixing coefficient @@ -430,9 +413,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) "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) + call get_param(paramFile, mdl, "KPP_K_ENH_FAC", CS%KPP_K_ENH_FAC, & + 'Constant value to enhance mixing coefficient in KPP.', & + units="nondim", default=1.0) endif endif !/ 2. Options related to enhancing the unresolved Vt2/entrainment in Rib @@ -470,9 +453,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) "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 , & + call get_param(paramFile, mdl, "KPP_VT2_ENH_FAC", CS%KPP_VT2_ENH_FAC, & 'Constant value to enhance VT2 in KPP.', & - default=1.0) + units="nondim", default=1.0) endif endif @@ -481,8 +464,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) 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, & + minOBLdepth=US%Z_to_m*CS%minOBLdepth, & + minVtsqr=US%L_T_to_m_s**2*CS%minVtsqr, & vonKarman=CS%vonKarman, & surf_layer_ext=CS%surf_layer_ext, & interp_type=CS%interpType, & @@ -547,13 +530,17 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) 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_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', conversion=US%C_to_degC) + 'Temperature of surface layer (10% of OBL depth) as passed to [CVMix] KPP', & + 'C', conversion=US%C_to_degC) 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', conversion=US%S_to_ppt) + 'Salinity of surface layer (10% of OBL depth) as passed to [CVMix] KPP', & + 'ppt', conversion=US%S_to_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') + 'i-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', & + 'm/s', conversion=US%L_T_to_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') + 'j-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', & + 'm/s', conversion=US%L_T_to_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, & @@ -616,7 +603,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, buoyFlux, Kt, Ks, Kv, & type(wave_parameters_CS), pointer :: Waves !< Wave CS for Langmuir turbulence real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult !< Langmuir enhancement multiplier -! Local variables + ! Local variables integer :: i, j, k ! Loop indices real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [m] (negative in ocean) real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m] (negative in ocean) @@ -624,14 +611,16 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, buoyFlux, Kt, Ks, Kv, & real, dimension( GV%ke+1 ) :: Kviscosity ! Vertical viscosity at interfaces [m2 s-1] real, dimension( GV%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces [nondim] - real :: surfFricVel, surfBuoyFlux - real :: sigma, sigmaRatio + real :: surfFricVel ! Surface friction velocity in MKS units [m s-1] + real :: surfBuoyFlux ! Surface buoyancy flux in MKS units [m2 s-3] + real :: sigma ! Fractional vertical position within the boundary layer [nondim] + real :: sigmaRatio ! A cubic function of sigma [nondim] real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1] - real :: dh ! The local thickness used for calculating interface positions [m] - real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m] + real :: dh ! The local thickness used for calculating interface positions [Z ~> m] + real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m] ! For Langmuir Calculations - real :: LangEnhK ! Langmuir enhancement for mixing coefficient + real :: LangEnhK ! Langmuir enhancement for mixing coefficient [nondim] if (CS%Stokes_Mixing .and. .not.associated(Waves)) call MOM_error(FATAL, & "KPP_calculate: The Waves control structure must be associated if STOKES_MIXING is True.") @@ -672,17 +661,17 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, buoyFlux, Kt, Ks, Kv, & do k=1,GV%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 = h(i,j,k) * GV%H_to_Z ! 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 + cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh + iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*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) + ! h to Monin-Obukhov (default is false, ie. not used) ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports @@ -820,8 +809,8 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, buoyFlux, Kt, Ks, Kv, & ! we apply nonLocalTrans in subroutines ! KPP_NonLocalTransport_temp and KPP_NonLocalTransport_saln - nonLocalTransHeat(i,j,:) = nonLocalTrans(:,1) ! temp - nonLocalTransScalar(i,j,:) = nonLocalTrans(:,2) ! saln + nonLocalTransHeat(i,j,:) = nonLocalTrans(:,1) ! temperature + nonLocalTransScalar(i,j,:) = nonLocalTrans(:,2) ! salinity ! set the KPP diffusivity and viscosity to zero for testing purposes if (CS%KPPzeroDiffusivity) then @@ -906,48 +895,54 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl 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_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3] type(wave_parameters_CS), pointer :: Waves !< Wave CS for Langmuir turbulence - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult!< Langmuir enhancement factor + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult !< Langmuir enhancement factor ! Local variables - integer :: i, j, k, km1 ! Loop indices + ! Variables in MKS units for passing to CVMix routines real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [m] (negative in ocean) real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m] (negative in ocean) real, dimension( GV%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces [s-2] real, dimension( GV%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars [m s-1] real, dimension( GV%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3] real, dimension( GV%ke ) :: deltaU2 ! square of delta U (shear) in denominator of Bulk Ri [m2 s-2] - real, dimension( GV%ke ) :: surfBuoyFlux2 + real, dimension( GV%ke ) :: surfBuoyFlux2 ! Surface buoyancy flux in MKS units [m2 s-3] real, dimension( GV%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer [nondim] + real :: surfFricVel ! Surface friction velocity in MKS units [m s-1] + real :: surfBuoyFlux ! Surface buoyancy flux in MKS units [m2 s-3] + real :: Coriolis ! Coriolis parameter at tracer points [s-1] + real :: zBottomMinusOffset ! Height of bottom plus a little bit [m] - ! for EOS calculation + ! Variables for EOS calculations real, dimension( 3*GV%ke ) :: rho_1D ! A column of densities [R ~> kg m-3] real, dimension( 3*GV%ke ) :: pres_1D ! A column of pressures [R L2 T-2 ~> Pa] real, dimension( 3*GV%ke ) :: Temp_1D ! A column of temperatures [C ~> degC] real, dimension( 3*GV%ke ) :: Salt_1D ! A column of salinities [S ~> ppt] - 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 :: Uk, Vk - - 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 ~> 1] - real :: delH ! Thickness of a layer [m] - real :: surfHtemp, surfTemp ! Integral and average of temp over the surface layer [C ~> degC] - real :: surfHsalt, surfSalt ! Integral and average of saln over the surface layer [S ~> ppt] - 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 + real :: GoRho ! Gravitational acceleration in MKS units divided by density [m s-2 R-1 ~> m4 kg-1 s-2] + real :: pRef ! The interface pressure [R L2 T-2 ~> Pa] + real :: Uk, Vk ! Layer velocities relative to their averages in the surface layer [L T-1 ~> m s-1] + real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth [Z ~> m] + real :: hTot ! Running sum of thickness used in the surface layer average [Z ~> m] + real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1] + real :: delH ! Thickness of a layer [Z ~> m] + real :: surfTemp ! Average of temperature over the surface layer [C ~> degC] + real :: surfHtemp ! Integral of temperature over the surface layer [Z C ~> m degC] + real :: surfSalt ! Average of salinity over the surface layer [S ~> ppt] + real :: surfHsalt ! Integral of salinity over the surface layer [Z S ~> m ppt] + real :: surfHu, surfHv ! Integral of u and v over the surface layer [Z L T-1 ~> m2 s-1] + real :: surfU, surfV ! Average of u and v over the surface layer [Z T-1 ~> m s-1] + real :: dh ! The local thickness used for calculating interface positions [Z ~> m] + real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m] ! For Langmuir Calculations - real :: LangEnhVt2 ! Langmuir enhancement for unresolved shear - real, dimension(GV%ke) :: U_H, V_H - real :: MLD_GUESS, LA - real :: surfHuS, surfHvS, surfUs, surfVs + real :: LangEnhVt2 ! Langmuir enhancement for unresolved shear [nondim] + real, dimension(GV%ke) :: U_H, V_H ! Velocities at tracer points [L T-1 ~> m s-1] + real :: MLD_guess ! A guess at the mixed layer depth for calculating the Langmuir number [Z ~> m] + real :: LA ! The local Langmuir number [nondim] + real :: surfHuS, surfHvS ! Stokes drift velocities integrated over the boundary layer [Z L T-1 ~> m2 s-1] + real :: surfUs, surfVs ! Stokes drift velocities averaged over the boundary layer [Z T-1 ~> m s-1] + + integer :: i, j, k, km1, kk, ksfc, ktmp ! Loop indices if (CS%Stokes_Mixing .and. .not.associated(Waves)) call MOM_error(FATAL, & "KPP_compute_BLD: The Waves control structure must be associated if STOKES_MIXING is True.") @@ -971,7 +966,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl !$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 Temp_1D, salt_1D, surfBuoyFlux2, MLD_guess, LA, rho_1D, & !$OMP deltarho, N2_1d, ws_1d, LangEnhVT2, & !$OMP BulkRi_1d, zBottomMinusOffset) & !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, & @@ -983,8 +978,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if (G%mask2dT(i,j)==0.) cycle do k=1,GV%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)) + U_H(k) = 0.5 * (u(i,j,k)+u(i-1,j,k)) + V_H(k) = 0.5 * (v(i,j,k)+v(i,j-1,k)) enddo ! things independent of position within the column @@ -1004,36 +999,36 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl do k=1,GV%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 = h(i,j,k) * GV%H_to_Z ! 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 + cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh + iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*dh ! find ksfc for cell where "surface layer" sits - SLdepth_0d = CS%surf_layer_ext*max( max(-cellHeight(k),-iFaceHeight(2) ), CS%minOBLdepth ) + SLdepth_0d = CS%surf_layer_ext*max( US%m_to_Z*max(-cellHeight(k),-iFaceHeight(2) ), CS%minOBLdepth ) ksfc = k do ktmp = 1,k - if (-1.0*iFaceHeight(ktmp+1) >= SLdepth_0d) then + if (-1.0*iFaceHeight(ktmp+1) >= US%Z_to_m*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 + ! average temperature, salinity, u and v over surface layer + ! use C-grid average to get u and 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 ) + delH = min( max(0.0, SLdepth_0d - hTot), h(i,j,ktmp)*GV%H_to_Z ) ! surface layer thickness hTot = hTot + delH @@ -1041,11 +1036,11 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! 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 + surfHu = surfHu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delH + surfHv = surfHv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delH if (CS%Stokes_Mixing) then - surfHus = surfHus + 0.5*US%L_T_to_m_s*(Waves%US_x(i,j,ktmp)+Waves%US_x(i-1,j,ktmp)) * delH - surfHvs = surfHvs + 0.5*US%L_T_to_m_s*(Waves%US_y(i,j,ktmp)+Waves%US_y(i,j-1,ktmp)) * delH + 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 @@ -1056,23 +1051,22 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl surfUs = surfHus / hTot surfVs = surfHvs / hTot - ! vertical shear between present layer and - ! surface layer averaged surfU,surfV. + ! vertical shear between present layer and surface layer averaged surfU and 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 + Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU + Vk = 0.5*(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*US%L_T_to_m_s*(Waves%Us_x(i,j,k)+Waves%US_x(i-1,j,k)) - surfUs ) - Vk = Vk + (0.5*US%L_T_to_m_s*(Waves%Us_y(i,j,k)+Waves%Us_y(i,j-1,k)) - surfVs ) + 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 + deltaU2(k) = US%L_T_to_m_s**2 * (Uk**2 + Vk**2) - ! pressure, temp, and saln for EOS + ! pressure, temperature, and salinity for calling the equation of state ! kk+1 = surface fields ! kk+2 = k fields ! kk+3 = km1 fields @@ -1098,7 +1092,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl enddo ! k-loop finishes if ( (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) .and. .not. present(lamult)) then - MLD_GUESS = max( 1.*US%m_to_Z, abs(US%m_to_Z*CS%OBLdepthprev(i,j) ) ) + 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 @@ -1127,12 +1121,12 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! 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 ) + 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 ) ! Determine the enhancement factor for unresolved shear IF (CS%LT_VT2_ENHANCEMENT) then @@ -1172,25 +1166,25 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl 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 + 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(GV%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(GV%ke+1)) + zBottomMinusOffset = iFaceHeight(GV%ke+1) + min(US%Z_to_m*CS%deepOBLoffset, -0.1*iFaceHeight(GV%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 + if (CS%fixedOBLdepth) CS%OBLdepth(i,j) = US%Z_to_m*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(GV%ke+1) ) ! no deeper than bottom CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) @@ -1211,14 +1205,14 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl ! 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( & + 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(:) + CS%Ws(i,j,:) = Ws_1d(:) endif ! Diagnostics @@ -1229,7 +1223,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl 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 + if (CS%id_Vsurf > 0) CS%Vsurf(i,j) = surfV enddo enddo @@ -1252,17 +1246,18 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if (CS%id_Vt2 > 0) call post_data(CS%id_Vt2, CS%Vt2, CS%diag) ! BLD smoothing: - if (CS%n_smooth > 0) call KPP_smooth_BLD(CS,G,GV,h) + if (CS%n_smooth > 0) call KPP_smooth_BLD(CS, G, GV, US, h) end subroutine KPP_compute_BLD !> Apply a 1-1-4-1-1 Laplacian filter one time on BLD to reduce any horizontal two-grid-point noise -subroutine KPP_smooth_BLD(CS,G,GV,h) +subroutine KPP_smooth_BLD(CS, G, GV, US, h) ! 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_(GV)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2] ! local @@ -1272,8 +1267,8 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m] ! (negative in the ocean) real :: wc, ww, we, wn, ws ! averaging weights for smoothing [nondim] - real :: dh ! The local thickness used for calculating interface positions [m] - real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m] + real :: dh ! The local thickness used for calculating interface positions [Z ~> m] + real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m] integer :: i, j, k, s call cpu_clock_begin(id_clock_KPP_smoothing) @@ -1288,7 +1283,7 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) OBLdepth_prev = CS%OBLdepth ! apply smoothing on OBL depth - !$OMP parallel do default(none) shared(G, GV, CS, h, OBLdepth_prev) & + !$OMP parallel do default(none) shared(G, GV, US, 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 @@ -1301,12 +1296,12 @@ subroutine KPP_smooth_BLD(CS,G,GV,h) do k=1,GV%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 = h(i,j,k) * GV%H_to_Z ! 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 + cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh + iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*dh enddo ! compute weights @@ -1347,9 +1342,9 @@ subroutine KPP_get_BLD(CS, BLD, G, US, m_to_BLD_units) 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 + !! to the desired units for BLD [various] ! Local variables - real :: scale ! A dimensional rescaling factor + real :: scale ! A dimensional rescaling factor in [Z m-1 ~> 1] or other units. integer :: i,j scale = US%m_to_Z ; if (present(m_to_BLD_units)) scale = m_to_BLD_units @@ -1376,11 +1371,12 @@ subroutine KPP_NonLocalTransport(CS, G, GV, h, nonLocalTrans, surfFlux, & type(tracer_type), pointer, intent(in) :: tr_ptr !< tracer_type has diagnostic ids on it real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< Scalar (scalar units [conc]) real, optional, intent(in) :: flux_scale !< Scale factor to get surfFlux - !! into proper units + !! into proper units [various] integer :: i, j, k real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dtracer ! Rate of tracer change [conc T-1 ~> conc s-1] - real, dimension(SZI_(G),SZJ_(G)) :: surfFlux_loc + real, dimension(SZI_(G),SZJ_(G)) :: surfFlux_loc ! An optionally rescaled surface flux of the scalar + ! in [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1] or other units ! term used to scale if (present(flux_scale)) then diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 58d6e3417a..a0b24dee70 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -27,14 +27,14 @@ module MOM_CVMix_conv type, public :: CVMix_conv_cs ; private ! Parameters - real :: kd_conv_const !< diffusivity constant used in convective regime [m2 s-1] - real :: kv_conv_const !< viscosity constant used in convective regime [m2 s-1] + real :: kd_conv_const !< diffusivity constant used in convective regime [Z2 T-1 ~> m2 s-1] + real :: kv_conv_const !< viscosity constant used in convective regime [Z2 T-1 ~> m2 s-1] real :: bv_sqr_conv !< Threshold for squared buoyancy frequency - !! needed to trigger Brunt-Vaisala parameterization [s-2] - real :: min_thickness !< Minimum thickness allowed [m] + !! needed to trigger Brunt-Vaisala parameterization [T-2 ~> s-2] + real :: min_thickness !< Minimum thickness allowed [Z ~> m] logical :: debug !< If true, turn on debugging - ! Daignostic handles and pointers + ! Diagnostic handles and pointers type(diag_ctrl), pointer :: diag => NULL() !< Pointer to diagnostics control structure !>@{ Diagnostics handles integer :: id_N2 = -1, id_kd_conv = -1, id_kv_conv = -1 @@ -55,13 +55,13 @@ logical function CVMix_conv_init(Time, G, GV, US, param_file, diag, CS) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. - type(CVMix_conv_cs), intent(inout) :: CS !< CVMix convetction control struct + type(CVMix_conv_cs), intent(inout) :: CS !< CVMix convection control structure real :: prandtl_conv !< Turbulent Prandtl number used in convective instabilities. logical :: useEPBL !< If True, use the ePBL boundary layer scheme. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" ! Read parameters call get_param(param_file, mdl, "USE_CVMix_CONVECTION", CVMix_conv_init, default=.false., do_not_log=.true.) @@ -90,7 +90,8 @@ logical function CVMix_conv_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) - call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, do_not_log=.True.) + call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, & + units="m", scale=US%m_to_Z, default=0.001, do_not_log=.True.) call openParameterBlock(param_file,'CVMix_CONVECTION') @@ -102,12 +103,12 @@ logical function CVMix_conv_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv_const, & "Diffusivity used in convective regime. Corresponding viscosity "//& "(KV_CONV) will be set to KD_CONV * PRANDTL_CONV.", & - units='m2/s', default=1.00) + units='m2/s', default=1.00, scale=US%m2_s_to_Z2_T) call get_param(param_file, mdl, 'BV_SQR_CONV', CS%bv_sqr_conv, & "Threshold for squared buoyancy frequency needed to trigger "//& "Brunt-Vaisala parameterization.", & - units='1/s^2', default=0.0) + units='1/s^2', default=0.0, scale=US%T_to_s**2) call closeParameterBlock(param_file) @@ -123,10 +124,10 @@ logical function CVMix_conv_init(Time, G, GV, US, param_file, diag, CS) CS%id_kv_conv = register_diag_field('ocean_model', 'kv_conv', diag%axesTi, Time, & 'Additional viscosity added by MOM_CVMix_conv module', 'm2/s', conversion=US%Z2_T_to_m2_s) - call CVMix_init_conv(convect_diff=CS%kd_conv_const, & - convect_visc=CS%kv_conv_const, & + call CVMix_init_conv(convect_diff=US%Z2_T_to_m2_s*CS%kd_conv_const, & + convect_visc=US%Z2_T_to_m2_s*CS%kv_conv_const, & lBruntVaisala=.true., & - BVsqr_convect=CS%bv_sqr_conv) + BVsqr_convect=US%s_to_T**2*CS%bv_sqr_conv) end function CVMix_conv_init @@ -139,7 +140,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - type(CVMix_conv_cs), intent(in) :: CS !< CVMix convection control struct + type(CVMix_conv_cs), intent(in) :: CS !< CVMix convection control structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hbl !< Depth of ocean boundary layer [Z ~> m] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & intent(inout) :: Kd !< Diapycnal diffusivity at each interface that @@ -167,14 +168,14 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) kd_conv, & !< Diffusivity added by convection for diagnostics [Z2 T-1 ~> m2 s-1] kv_conv, & !< Viscosity added by convection for diagnostics [Z2 T-1 ~> m2 s-1] N2_3d !< Squared buoyancy frequency for diagnostics [T-2 ~> s-2] - integer :: kOBL !< level of OBL extent - real :: g_o_rho0 ! Gravitational acceleration divided by density times unit convserion factors + integer :: kOBL !< level of ocean boundary layer extent + real :: g_o_rho0 ! Gravitational acceleration divided by density times unit conversion factors ! [Z s-2 R-1 ~> m4 s-2 kg-1] real :: pref ! Interface pressures [R L2 T-2 ~> Pa] real :: rhok, rhokm1 ! In situ densities of the layers above and below at the interface pressure [R ~> kg m-3] real :: hbl_KPP ! The depth of the ocean boundary as used by KPP [m] real :: dz ! A thickness [Z ~> m] - real :: dh, hcorr ! Two thicknesses [m] + real :: dh, hcorr ! Limited thicknesses and a cumulative correction [Z ~> m] integer :: i, j, k g_o_rho0 = US%L_to_Z**2*US%s_to_T**2 * GV%g_Earth / GV%Rho0 @@ -213,12 +214,12 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) hcorr = 0.0 ! compute heights at cell center and interfaces do k=1,GV%ke - dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment, in the units used by CVMix. + dh = h(i,j,k) * GV%H_to_Z ! Nominal thickness to use for increment, in the units of heights 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 + dh = max(dh, CS%min_thickness) ! Limited increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh + iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*dh enddo ! gets index of the level and interface above hbl diff --git a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 index 413b87f631..6e2c76ba8d 100644 --- a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 @@ -28,14 +28,14 @@ module MOM_CVMix_ddiff ! Parameters real :: strat_param_max !< maximum value for the stratification parameter [nondim] real :: kappa_ddiff_s !< leading coefficient in formula for salt-fingering regime - !! for salinity diffusion [m2 s-1] + !! for salinity diffusion [Z2 T-1 ~> m2 s-1] real :: ddiff_exp1 !< interior exponent in salt-fingering regime formula [nondim] real :: ddiff_exp2 !< exterior exponent in salt-fingering regime formula [nondim] - real :: mol_diff !< molecular diffusivity [m2 s-1] + real :: mol_diff !< molecular diffusivity [Z2 T-1 ~> m2 s-1] real :: kappa_ddiff_param1 !< exterior coefficient in diffusive convection regime [nondim] real :: kappa_ddiff_param2 !< middle coefficient in diffusive convection regime [nondim] real :: kappa_ddiff_param3 !< interior coefficient in diffusive convection regime [nondim] - real :: min_thickness !< Minimum thickness allowed [m] + real :: min_thickness !< Minimum thickness allowed [Z ~> m] character(len=4) :: diff_conv_type !< type of diffusive convection to use. Options are Marmorino & !! Caldwell 1976 ("MC76"; default) and Kelley 1988, 1990 ("K90") logical :: debug !< If true, turn on debugging @@ -57,8 +57,8 @@ logical function CVMix_ddiff_init(Time, G, GV, US, param_file, diag, CS) type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. type(CVMix_ddiff_cs), pointer :: CS !< This module's control structure. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" if (associated(CS)) then call MOM_error(WARNING, "CVMix_ddiff_init called with an associated "// & @@ -82,7 +82,8 @@ logical function CVMix_ddiff_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) - call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, do_not_log=.True.) + call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, & + units="m", scale=US%m_to_Z, default=0.001, do_not_log=.True.) call openParameterBlock(param_file,'CVMIX_DDIFF') @@ -91,8 +92,8 @@ logical function CVMix_ddiff_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=2.55) call get_param(param_file, mdl, "KAPPA_DDIFF_S", CS%kappa_ddiff_s, & - "Leading coefficient in formula for salt-fingering regime "//& - "for salinity diffusion.", units="m2 s-1", default=1.0e-4) + "Leading coefficient in formula for salt-fingering regime for salinity diffusion.", & + units="m2 s-1", default=1.0e-4, scale=US%m2_s_to_Z2_T) call get_param(param_file, mdl, "DDIFF_EXP1", CS%ddiff_exp1, & "Interior exponent in salt-fingering regime formula.", & @@ -116,7 +117,7 @@ logical function CVMix_ddiff_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "MOL_DIFF", CS%mol_diff, & "Molecular diffusivity used in CVMix double diffusion.", & - units="m2 s-1", default=1.5e-6) + units="m2 s-1", default=1.5e-6, scale=US%m2_s_to_Z2_T) call get_param(param_file, mdl, "DIFF_CONV_TYPE", CS%diff_conv_type, & "type of diffusive convection to use. Options are Marmorino \n" //& @@ -126,10 +127,10 @@ logical function CVMix_ddiff_init(Time, G, GV, US, param_file, diag, CS) call closeParameterBlock(param_file) call cvmix_init_ddiff(strat_param_max=CS%strat_param_max, & - kappa_ddiff_s=CS%kappa_ddiff_s, & + kappa_ddiff_s=US%Z2_T_to_m2_s*CS%kappa_ddiff_s, & ddiff_exp1=CS%ddiff_exp1, & ddiff_exp2=CS%ddiff_exp2, & - mol_diff=CS%mol_diff, & + mol_diff=US%Z2_T_to_m2_s*CS%mol_diff, & kappa_ddiff_param1=CS%kappa_ddiff_param1, & kappa_ddiff_param2=CS%kappa_ddiff_param2, & kappa_ddiff_param3=CS%kappa_ddiff_param3, & @@ -160,21 +161,21 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho) ! Local variables real, dimension(SZK_(GV)) :: & cellHeight, & !< Height of cell centers [m] - dRho_dT, & !< partial derivatives of density wrt temp [R C-1 ~> kg m-3 degC-1] - dRho_dS, & !< partial derivatives of density wrt saln [R S-1 ~> kg m-3 ppt-1] + dRho_dT, & !< partial derivatives of density with temperature [R C-1 ~> kg m-3 degC-1] + dRho_dS, & !< partial derivatives of density with salinity [R S-1 ~> kg m-3 ppt-1] pres_int, & !< pressure at each interface [R L2 T-2 ~> Pa] temp_int, & !< temp and at interfaces [C ~> degC] salt_int, & !< salt at at interfaces [S ~> ppt] alpha_dT, & !< alpha*dT across interfaces [kg m-3] beta_dS, & !< beta*dS across interfaces [kg m-3] - dT, & !< temp. difference between adjacent layers [C ~> degC] - dS !< salt difference between adjacent layers [S ~> ppt] + dT, & !< temperature difference between adjacent layers [C ~> degC] + dS !< salinity difference between adjacent layers [S ~> ppt] real, dimension(SZK_(GV)+1) :: & Kd1_T, & !< Diapycanal diffusivity of temperature [m2 s-1]. Kd1_S !< Diapycanal diffusivity of salinity [m2 s-1]. real, dimension(SZK_(GV)+1) :: iFaceHeight !< Height of interfaces [m] - real :: dh, hcorr + real :: dh, hcorr ! Limited thicknesses and a cumulative correction [Z ~> m] integer :: i, k ! initialize dummy variables @@ -184,7 +185,7 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho) ! GMM, I am leaving some code commented below. We need to pass BLD to - ! this soubroutine to avoid adding diffusivity above that. This needs + ! this subroutine to avoid adding diffusivity above that. This needs ! to be done once we re-structure the order of the calls. !if (.not. associated(hbl)) then ! allocate(hbl(SZI_(G), SZJ_(G))); @@ -234,16 +235,16 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho) hcorr = 0.0 ! compute heights at cell center and interfaces do k=1,GV%ke - dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = h(i,j,k) * GV%H_to_Z ! Nominal thickness to use for increment, in height units 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 + cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh + iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*dh enddo ! gets index of the level and interface above hbl - !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) + !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight, hbl(i,j)) Kd1_T(:) = 0.0 ; Kd1_S(:) = 0.0 call CVMix_coeffs_ddiff(Tdiff_out=Kd1_T(:), & @@ -277,7 +278,7 @@ logical function CVMix_ddiff_is_used(param_file) end function CVMix_ddiff_is_used -!> Clear pointers and dealocate memory +!> Clear pointers and deallocate memory ! NOTE: Placeholder destructor subroutine CVMix_ddiff_end(CS) type(CVMix_ddiff_cs), pointer :: CS !< Control structure for this module that diff --git a/src/parameterizations/vertical/MOM_CVMix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 index 7ec45dbe11..b69cd2daae 100644 --- a/src/parameterizations/vertical/MOM_CVMix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -28,19 +28,19 @@ module MOM_CVMix_shear ! vary with the Boussinesq approximation, the Boussinesq variant is given first. !> Control structure including parameters for CVMix interior shear schemes. -type, public :: CVMix_shear_cs ! TODO: private +type, public :: CVMix_shear_cs ; private logical :: use_LMD94 !< Flags to use the LMD94 scheme logical :: use_PP81 !< Flags to use Pacanowski and Philander (JPO 1981) logical :: smooth_ri !< If true, smooth Ri using a 1-2-1 filter - real :: Ri_zero !< LMD94 critical Richardson number - real :: Nu_zero !< LMD94 maximum interior diffusivity - real :: KPP_exp !< Exponent of unitless factor of diff. - !! for KPP internal shear mixing scheme. + real :: Ri_zero !< LMD94 critical Richardson number [nondim] + real :: Nu_zero !< LMD94 maximum interior diffusivity [Z2 T-1 ~> m2 s-1] + real :: KPP_exp !< Exponent of unitless factor of diffusivities + !! for KPP internal shear mixing scheme [nondim] real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency [T-2 ~> s-2] real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency [T-2 ~> s-2] - real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number + real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number [nondim] real, allocatable, dimension(:,:,:) :: ri_grad_smooth !< Gradient Richardson number - !! after smoothing + !! after smoothing [nondim] character(10) :: Mix_Scheme !< Mixing scheme name (string) type(diag_ctrl), pointer :: diag => NULL() !< Pointer to the diagnostics control structure @@ -137,7 +137,7 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS ) S2 = US%L_to_Z**2*(DU*DU+DV*DV)/(DZ*DZ) Ri_Grad(k) = max(0., N2) / max(S2, 1.e-10*US%T_to_s**2) - ! fill 3d arrays, if user asks for diagsnostics + ! fill 3d arrays, if user asks for diagnostics if (CS%id_N2 > 0) CS%N2(i,j,k) = N2 if (CS%id_S2 > 0) CS%S2(i,j,k) = S2 @@ -264,22 +264,22 @@ logical function CVMix_shear_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "NU_ZERO", CS%Nu_Zero, & "Leading coefficient in KPP shear mixing.", & - units="nondim", default=5.e-3) + units="m2 s-1", default=5.e-3, scale=US%m2_s_to_Z2_T) call get_param(param_file, mdl, "RI_ZERO", CS%Ri_Zero, & "Critical Richardson for KPP shear mixing, "// & "NOTE this the internal mixing and this is "// & - "not for setting the boundary layer depth." & - ,units="nondim", default=0.8) + "not for setting the boundary layer depth.", & + units="nondim", default=0.8) call get_param(param_file, mdl, "KPP_EXP", CS%KPP_exp, & "Exponent of unitless factor of diffusivities, "// & - "for KPP internal shear mixing scheme." & - ,units="nondim", default=3.0) + "for KPP internal shear mixing scheme.", & + units="nondim", default=3.0) call get_param(param_file, mdl, "SMOOTH_RI", CS%smooth_ri, & "If true, vertically smooth the Richardson "// & "number by applying a 1-2-1 filter once.", & default = .false.) call cvmix_init_shear(mix_scheme=CS%Mix_Scheme, & - KPP_nu_zero=CS%Nu_Zero, & + KPP_nu_zero=US%Z2_T_to_m2_s*CS%Nu_Zero, & KPP_Ri_zero=CS%Ri_zero, & KPP_exp=CS%KPP_exp) @@ -332,7 +332,7 @@ logical function CVMix_shear_is_used(param_file) CVMix_shear_is_used = (LMD94 .or. PP81) end function CVMix_shear_is_used -!> Clear pointers and dealocate memory +!> Clear pointers and deallocate memory subroutine CVMix_shear_end(CS) type(CVMix_shear_cs), intent(inout) :: CS !< Control structure for this module that !! will be deallocated in this subroutine diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 75798ed466..fa3dbe4b87 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -62,7 +62,7 @@ module MOM_tidal_mixing real, allocatable :: N2_meanz(:,:) !< vertically averaged buoyancy frequency [T-2 ~> s-2] real, allocatable :: Polzin_decay_scale_scaled(:,:) !< vertical scale of decay for tidal dissipation [Z ~> m] real, allocatable :: Polzin_decay_scale(:,:) !< vertical decay scale for tidal dissipation with Polzin [Z ~> m] - real, allocatable :: Simmons_coeff_2d(:,:) !< The Simmons et al mixing coefficient + real, allocatable :: Simmons_coeff_2d(:,:) !< The Simmons et al mixing coefficient [nondim] end type !> Control structure with parameters for the tidal mixing module. @@ -129,13 +129,14 @@ module MOM_tidal_mixing logical :: use_CVMix_tidal = .false. !< true if CVMix is to be used for determining !! diffusivity due to tidal mixing - real :: min_thickness !< Minimum thickness allowed [m] + real :: min_thickness !< Minimum thickness allowed [Z ~> m] ! CVMix-specific parameters integer :: CVMix_tidal_scheme = -1 !< 1 for Simmons, 2 for Schmittner type(CVMix_tidal_params_type) :: CVMix_tidal_params !< A CVMix-specific type with parameters for tidal mixing type(CVMix_global_params_type) :: CVMix_glb_params !< CVMix-specific for Prandtl number only - real :: tidal_max_coef !< CVMix-specific maximum allowable tidal diffusivity. [m^2/s] + real :: tidal_max_coef !< CVMix-specific maximum allowable tidal + !! diffusivity. [Z2 T-1 ~> m2 s-1] real :: tidal_diss_lim_tc !< CVMix-specific dissipation limit depth for !! tidal-energy-constituent data [Z ~> m]. type(remapping_CS) :: remap_CS !< The control structure for remapping @@ -443,8 +444,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di call get_param(param_file, mdl, "INT_TIDE_DECAY_SCALE", CS%Int_tide_decay_scale, & "The decay scale away from the bottom for tidal TKE with "//& "the new coding when INT_TIDE_DISSIPATION is used.", & - !units="m", default=0.0) - units="m", default=500.0, scale=US%m_to_Z) ! TODO: confirm this new default + units="m", default=500.0, scale=US%m_to_Z) call get_param(param_file, mdl, "MU_ITIDES", CS%Mu_itides, & "A dimensionless turbulent mixing efficiency used with "//& "INT_TIDE_DISSIPATION, often 0.2.", units="nondim", default=0.2) @@ -564,12 +564,10 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di call get_param(param_file, mdl, "GAMMA_NIKURASHIN",CS%Gamma_lee, & "The fraction of the lee wave energy that is dissipated "//& - "locally with LEE_WAVE_DISSIPATION.", units="nondim", & - default=0.3333) + "locally with LEE_WAVE_DISSIPATION.", units="nondim", default=0.3333) call get_param(param_file, mdl, "DECAY_SCALE_FACTOR_LEE",CS%Decay_scale_factor_lee, & "Scaling for the vertical decay scale of the local "//& - "dissipation of lee wave dissipation.", units="nondim", & - default=1.0) + "dissipation of lee wave dissipation.", units="nondim", default=1.0) else CS%Decay_scale_factor_lee = -9.e99 ! This should never be used if CS%Lee_wave_dissipation = False endif @@ -581,18 +579,17 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di !call openParameterBlock(param_file,'CVMix_TIDAL') call get_param(param_file, mdl, "TIDAL_MAX_COEF", CS%tidal_max_coef, & "largest acceptable value for tidal diffusivity", & - units="m^2/s", default=50e-4) ! the default is 50e-4 in CVMix, 100e-4 in POP. + units="m^2/s", default=50e-4, scale=US%m2_s_to_Z2_T) ! the default is 50e-4 in CVMix, 100e-4 in POP. call get_param(param_file, mdl, "TIDAL_DISS_LIM_TC", CS%tidal_diss_lim_tc, & "Min allowable depth for dissipation for tidal-energy-constituent data. "//& "No dissipation contribution is applied above TIDAL_DISS_LIM_TC.", & units="m", default=0.0, scale=US%m_to_Z) - call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, & - do_not_log=.True.) + call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, & + units="m", default=0.001, scale=US%m_to_Z, do_not_log=.True.) call get_param(param_file, mdl, "PRANDTL_TIDAL", prandtl_tidal, & "Prandtl number used by CVMix tidal mixing schemes "//& "to convert vertical diffusivities into viscosities.", & - units="nondim", default=1.0, & - do_not_log=.true.) + units="nondim", default=1.0, do_not_log=.true.) call CVMix_put(CS%CVMix_glb_params,'Prandtl',prandtl_tidal) call get_param(param_file, mdl, "TIDAL_ENERGY_TYPE",tidal_energy_type, & @@ -615,7 +612,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di mix_scheme = CVMix_tidal_scheme_str, & efficiency = CS%Mu_itides, & vertical_decay_scale = CS%int_tide_decay_scale*US%Z_to_m, & - max_coefficient = CS%tidal_max_coef, & + max_coefficient = CS%tidal_max_coef*US%Z2_T_to_m2_s, & local_mixing_frac = CS%Gamma_itides, & depth_cutoff = CS%min_zbot_itides*US%Z_to_m) @@ -777,19 +774,21 @@ subroutine calculate_CVMix_tidal(h, j, N2_int, G, GV, US, CS, Kv, Kd_lay, Kd_int ! Local variables real, dimension(SZK_(GV)+1) :: Kd_tidal ! tidal diffusivity [m2 s-1] real, dimension(SZK_(GV)+1) :: Kv_tidal ! tidal viscosity [m2 s-1] - real, dimension(SZK_(GV)+1) :: vert_dep ! vertical deposition + real, dimension(SZK_(GV)+1) :: vert_dep ! vertical deposition [nondim] real, dimension(SZK_(GV)+1) :: iFaceHeight ! Height of interfaces [m] real, dimension(SZK_(GV)+1) :: SchmittnerSocn real, dimension(SZK_(GV)) :: cellHeight ! Height of cell centers [m] real, dimension(SZK_(GV)) :: tidal_qe_md ! Tidal dissipation energy interpolated from 3d input ! to model coordinates real, dimension(SZK_(GV)+1) :: N2_int_i ! De-scaled interface buoyancy frequency [s-2] - real, dimension(SZK_(GV)) :: Schmittner_coeff + real, dimension(SZK_(GV)) :: Schmittner_coeff ! A coefficient in the Schmittner et al (2014) mixing + ! parameterization [nondim] real, dimension(SZK_(GV)) :: h_m ! Cell thickness [m] real, allocatable, dimension(:,:) :: exp_hab_zetar + real :: dh, hcorr ! Limited thicknesses and a cumulative correction [Z ~> m] + real :: Simmons_coeff ! A coefficient in the Simmons et al (2004) mixing parameterization [nondim] integer :: i, k, is, ie - real :: dh, hcorr, Simmons_coeff real, parameter :: rho_fw = 1000.0 ! fresh water density [kg m-3] ! TODO: when coupled, get this from CESM (SHR_CONST_RHOFW) @@ -803,14 +802,14 @@ subroutine calculate_CVMix_tidal(h, j, N2_int, G, GV, US, CS, Kv, Kd_lay, Kd_int iFaceHeight = 0.0 ! BBL is all relative to the surface hcorr = 0.0 + ! Compute cell center depth and cell bottom in meters (negative values in the ocean) do k=1,GV%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, rescaled to m for use by CVMix. + dh = h(i,j,k) * GV%H_to_Z ! Nominal thickness to use for increment, in the units of heights 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 + dh = max(dh, CS%min_thickness) ! Limited increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh + iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*dh enddo call CVMix_compute_Simmons_invariant( nlev = GV%ke, & @@ -889,16 +888,17 @@ subroutine calculate_CVMix_tidal(h, j, N2_int, G, GV, US, CS, Kv, Kd_lay, Kd_int if (G%mask2dT(i,j)<1) cycle - iFaceHeight = 0.0 ! BBL is all relative to the surface + iFaceHeight(:) = 0.0 ! BBL is all relative to the surface hcorr = 0.0 + ! Compute heights at cell center and interfaces, and rescale layer thicknesses do k=1,GV%ke h_m(k) = h(i,j,k)*GV%H_to_m ! Rescale thicknesses to m for use by CVmix. - ! cell center and cell bottom in meters (negative values in the ocean) - dh = h_m(k) + hcorr ! Nominal thickness less the accumulated error (could temporarily make dh<0) + dh = h(i,j,k) * GV%H_to_Z ! Nominal thickness to use for increment, in the units of heights + 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 + dh = max(dh, CS%min_thickness) ! Limited increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh + iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*dh enddo SchmittnerSocn = 0.0 ! TODO: compute this