portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> Interface to CVMix convection scheme. | |
| 6 | module MOM_CVMix_conv | |
| 7 | ||
| 8 | use MOM_debugging, only : hchksum | |
| 9 | use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field | |
| 10 | use MOM_diag_mediator, only : post_data | |
| 11 | use MOM_EOS, only : calculate_density | |
| 12 | use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE | |
| 13 | use MOM_file_parser, only : openParameterBlock, closeParameterBlock | |
| 14 | use MOM_file_parser, only : get_param, log_version, param_file_type | |
| 15 | use MOM_grid, only : ocean_grid_type | |
| 16 | use MOM_interface_heights, only : thickness_to_dz | |
| 17 | use MOM_unit_scaling, only : unit_scale_type | |
| 18 | use MOM_variables, only : thermo_var_ptrs | |
| 19 | use MOM_verticalGrid, only : verticalGrid_type | |
| 20 | use CVMix_convection, only : CVMix_init_conv, CVMix_coeffs_conv | |
| 21 | use CVMix_kpp, only : CVMix_kpp_compute_kOBL_depth | |
| 22 | ||
| 23 | implicit none ; private | |
| 24 | ||
| 25 | #include <MOM_memory.h> | |
| 26 | ||
| 27 | public CVMix_conv_init, calculate_CVMix_conv, CVMix_conv_is_used | |
| 28 | ||
| 29 | !> Control structure including parameters for CVMix convection. | |
| 30 | type, public :: CVMix_conv_cs ; private | |
| 31 | ||
| 32 | ! Parameters | |
| 33 | real :: kd_conv_const !< diffusivity constant used in convective regime [Z2 T-1 ~> m2 s-1] | |
| 34 | real :: kv_conv_const !< viscosity constant used in convective regime [Z2 T-1 ~> m2 s-1] | |
| 35 | real :: bv_sqr_conv !< Threshold for squared buoyancy frequency | |
| 36 | !! needed to trigger Brunt-Vaisala parameterization [T-2 ~> s-2] | |
| 37 | real :: min_thickness !< Minimum thickness allowed [Z ~> m] | |
| 38 | logical :: debug !< If true, turn on debugging | |
| 39 | ||
| 40 | ! Diagnostic handles and pointers | |
| 41 | type(diag_ctrl), pointer :: diag => NULL() !< Pointer to diagnostics control structure | |
| 42 | !>@{ Diagnostics handles | |
| 43 | integer :: id_N2 = -1, id_kd_conv = -1, id_kv_conv = -1 | |
| 44 | !>@} | |
| 45 | ||
| 46 | end type CVMix_conv_cs | |
| 47 | ||
| 48 | character(len=40) :: mdl = "MOM_CVMix_conv" !< This module's name. | |
| 49 | ||
| 50 | contains | |
| 51 | ||
| 52 | !> Initialized the CVMix convection mixing routine. | |
| 53 | 1 | logical function CVMix_conv_init(Time, G, GV, US, param_file, diag, CS) |
| 54 | ||
| 55 | type(time_type), intent(in) :: Time !< The current time. | |
| 56 | type(ocean_grid_type), intent(in) :: G !< Grid structure. | |
| 57 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. | |
| 58 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 59 | type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle | |
| 60 | type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. | |
| 61 | type(CVMix_conv_cs), intent(inout) :: CS !< CVMix convection control structure | |
| 62 | ||
| 63 | real :: prandtl_conv !< Turbulent Prandtl number used in convective instabilities [nondim] | |
| 64 | logical :: useEPBL !< If True, use the ePBL boundary layer scheme. | |
| 65 | ||
| 66 | ! This include declares and sets the variable "version". | |
| 67 | # include "version_variable.h" | |
| 68 | ||
| 69 | ! Read parameters | |
| 70 | 1 | call get_param(param_file, mdl, "USE_CVMix_CONVECTION", CVMix_conv_init, default=.false., do_not_log=.true.) |
| 71 | call log_version(param_file, mdl, version, & | |
| 72 | "Parameterization of enhanced mixing due to convection via CVMix", & | |
| 73 | 1 | all_default=.not.CVMix_conv_init) |
| 74 | call get_param(param_file, mdl, "USE_CVMix_CONVECTION", CVMix_conv_init, & | |
| 75 | "If true, turns on the enhanced mixing due to convection "//& | |
| 76 | "via CVMix. This scheme increases diapycnal diffs./viscs. "//& | |
| 77 | "at statically unstable interfaces. Relevant parameters are "//& | |
| 78 | "contained in the CVMix_CONVECTION% parameter block.", & | |
| 79 | 1 | default=.false.) |
| 80 | ||
| 81 | 1 | if (.not. CVMix_conv_init) return |
| 82 | ||
| 83 | call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useEPBL, default=.false., & | |
| 84 | 0 | do_not_log=.true.) |
| 85 | ||
| 86 | ! Warn user if EPBL is being used, since in this case mixing due to convection will | |
| 87 | ! be aplied in the boundary layer | |
| 88 | 0 | if (useEPBL) then |
| 89 | call MOM_error(WARNING, 'MOM_CVMix_conv_init: '// & | |
| 90 | 'CVMix convection may not be properly applied when ENERGETICS_SFC_PBL = True '//& | |
| 91 | 0 | 'as convective mixing might occur in the boundary layer.') |
| 92 | endif | |
| 93 | ||
| 94 | 0 | call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) |
| 95 | ||
| 96 | call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, & | |
| 97 | 0 | units="m", scale=US%m_to_Z, default=0.001, do_not_log=.True.) |
| 98 | ||
| 99 | 0 | call openParameterBlock(param_file,'CVMix_CONVECTION') |
| 100 | ||
| 101 | call get_param(param_file, mdl, "PRANDTL_CONV", prandtl_conv, & | |
| 102 | "The turbulent Prandtl number applied to convective "//& | |
| 103 | "instabilities (i.e., used to convert KD_CONV into KV_CONV)", & | |
| 104 | 0 | units="nondim", default=1.0) |
| 105 | ||
| 106 | call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv_const, & | |
| 107 | "Diffusivity used in convective regime. Corresponding viscosity "//& | |
| 108 | "(KV_CONV) will be set to KD_CONV * PRANDTL_CONV.", & | |
| 109 | 0 | units='m2/s', default=1.00, scale=US%m2_s_to_Z2_T) |
| 110 | ||
| 111 | call get_param(param_file, mdl, 'BV_SQR_CONV', CS%bv_sqr_conv, & | |
| 112 | "Threshold for squared buoyancy frequency needed to trigger "//& | |
| 113 | "Brunt-Vaisala parameterization.", & | |
| 114 | 0 | units='1/s^2', default=0.0, scale=US%T_to_s**2) |
| 115 | ||
| 116 | 0 | call closeParameterBlock(param_file) |
| 117 | ||
| 118 | ! set kv_conv_const based on kd_conv_const and prandtl_conv | |
| 119 | 0 | CS%kv_conv_const = CS%kd_conv_const * prandtl_conv |
| 120 | ||
| 121 | ! Register diagnostics | |
| 122 | 0 | CS%diag => diag |
| 123 | CS%id_N2 = register_diag_field('ocean_model', 'N2_conv', diag%axesTi, Time, & | |
| 124 | 0 | 'Square of Brunt-Vaisala frequency used by MOM_CVMix_conv module', '1/s2', conversion=US%s_to_T**2) |
| 125 | CS%id_kd_conv = register_diag_field('ocean_model', 'kd_conv', diag%axesTi, Time, & | |
| 126 | 0 | 'Additional diffusivity added by MOM_CVMix_conv module', 'm2/s', conversion=US%Z2_T_to_m2_s) |
| 127 | CS%id_kv_conv = register_diag_field('ocean_model', 'kv_conv', diag%axesTi, Time, & | |
| 128 | 0 | 'Additional viscosity added by MOM_CVMix_conv module', 'm2/s', conversion=US%Z2_T_to_m2_s) |
| 129 | ||
| 130 | call CVMix_init_conv(convect_diff=US%Z2_T_to_m2_s*CS%kd_conv_const, & | |
| 131 | convect_visc=US%Z2_T_to_m2_s*CS%kv_conv_const, & | |
| 132 | lBruntVaisala=.true., & | |
| 133 | 0 | BVsqr_convect=US%s_to_T**2*CS%bv_sqr_conv) |
| 134 | ||
| 135 | 0 | end function CVMix_conv_init |
| 136 | ||
| 137 | !> Subroutine for calculating enhanced diffusivity/viscosity | |
| 138 | !! due to convection via CVMix | |
| 139 | 0 | subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) |
| 140 | ||
| 141 | type(ocean_grid_type), intent(in) :: G !< Grid structure. | |
| 142 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. | |
| 143 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 144 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. | |
| 145 | type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. | |
| 146 | type(CVMix_conv_cs), intent(in) :: CS !< CVMix convection control structure | |
| 147 | real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hbl !< Depth of ocean boundary layer [Z ~> m] | |
| 148 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & | |
| 149 | intent(inout) :: Kd !< Diapycnal diffusivity at each interface | |
| 150 | !! that will be incremented here | |
| 151 | !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 152 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & | |
| 153 | intent(inout) :: Kv !< Viscosity at each interface that will be | |
| 154 | !! incremented here [H Z T-1 ~> m2 s-1 or Pa s] | |
| 155 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & | |
| 156 | optional, intent(inout) :: Kd_aux !< A second diapycnal diffusivity at each | |
| 157 | !! interface that will also be incremented | |
| 158 | !! here [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 159 | ||
| 160 | ! local variables | |
| 161 | 0 | real, dimension(SZK_(GV)) :: rho_lwr !< Adiabatic Water Density [kg m-3], this is a dummy |
| 162 | !! variable since here convection is always | |
| 163 | !! computed based on Brunt Vaisala. | |
| 164 | 0 | real, dimension(SZK_(GV)) :: rho_1d !< water density in a column [kg m-3], this is also |
| 165 | !! a dummy variable, same reason as above. | |
| 166 | 0 | real, dimension(SZK_(GV)+1) :: N2 !< Squared buoyancy frequency [s-2] |
| 167 | 0 | real, dimension(SZK_(GV)+1) :: kv_col !< Viscosities at interfaces in the column [m2 s-1] |
| 168 | 0 | real, dimension(SZK_(GV)+1) :: kd_col !< Diffusivities at interfaces in the column [m2 s-1] |
| 169 | 0 | real, dimension(SZK_(GV)+1) :: iFaceHeight !< Height of interfaces [Z ~> m] |
| 170 | 0 | real, dimension(SZK_(GV)) :: cellHeight !< Height of cell centers [Z ~> m] |
| 171 | 0 | real, dimension(SZI_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m] |
| 172 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & | |
| 173 | 0 | kd_conv, & !< Diffusivity added by convection for diagnostics [Z2 T-1 ~> m2 s-1] |
| 174 | 0 | kv_conv, & !< Viscosity added by convection for diagnostics [Z2 T-1 ~> m2 s-1] |
| 175 | 0 | N2_3d !< Squared buoyancy frequency for diagnostics [T-2 ~> s-2] |
| 176 | integer :: kOBL !< level of ocean boundary layer extent | |
| 177 | real :: g_o_rho0 ! Gravitational acceleration, perhaps divided by density, times unit conversion factors | |
| 178 | ! [H s-2 R-1 ~> m4 s-2 kg-1 or m s-2] | |
| 179 | real :: pref ! Interface pressures [R L2 T-2 ~> Pa] | |
| 180 | real :: rhok, rhokm1 ! In situ densities of the layers above and below at the interface pressure [R ~> kg m-3] | |
| 181 | real :: dh_int ! The distance between layer centers [H ~> m or kg m-2] | |
| 182 | real :: dh, hcorr ! Limited thicknesses and a cumulative correction [Z ~> m] | |
| 183 | integer :: i, j, k | |
| 184 | ||
| 185 | 0 | if (GV%Boussinesq) then |
| 186 | 0 | g_o_rho0 = (US%s_to_T**2*GV%Z_to_H) * GV%g_Earth_Z_T2 / GV%Rho0 |
| 187 | else | |
| 188 | 0 | g_o_rho0 = (US%s_to_T**2*GV%RZ_to_H) * GV%g_Earth_Z_T2 |
| 189 | endif | |
| 190 | ||
| 191 | ! initialize dummy variables | |
| 192 | 0 | rho_lwr(:) = 0.0 ; rho_1d(:) = 0.0 |
| 193 | ||
| 194 | ! set N2 to zero at the top- and bottom-most interfaces | |
| 195 | 0 | N2(1) = 0.0 ; N2(GV%ke+1) = 0.0 |
| 196 | ||
| 197 | 0 | if (CS%id_N2 > 0) N2_3d(:,:,:) = 0.0 |
| 198 | 0 | if (CS%id_kv_conv > 0) Kv_conv(:,:,:) = 0.0 |
| 199 | 0 | if (CS%id_kd_conv > 0) Kd_conv(:,:,:) = 0.0 |
| 200 | ||
| 201 | 0 | do j = G%jsc, G%jec |
| 202 | ||
| 203 | ! Find the vertical distances across layers. | |
| 204 | 0 | call thickness_to_dz(h, tv, dz, j, G, GV) |
| 205 | ||
| 206 | 0 | do i = G%isc, G%iec |
| 207 | ||
| 208 | ! skip calling at land points | |
| 209 | !if (G%mask2dT(i,j) == 0.) cycle | |
| 210 | ||
| 211 | 0 | pRef = 0. ; if (associated(tv%p_surf)) pRef = tv%p_surf(i,j) |
| 212 | ! Compute Brunt-Vaisala frequency (static stability) on interfaces | |
| 213 | 0 | do K=2,GV%ke |
| 214 | ||
| 215 | ! pRef is pressure at interface between k and km1 [R L2 T-2 ~> Pa]. | |
| 216 | 0 | pRef = pRef + (GV%H_to_RZ*GV%g_Earth) * h(i,j,k) |
| 217 | 0 | call calculate_density(tv%t(i,j,k), tv%s(i,j,k), pRef, rhok, tv%eqn_of_state) |
| 218 | 0 | call calculate_density(tv%t(i,j,k-1), tv%s(i,j,k-1), pRef, rhokm1, tv%eqn_of_state) |
| 219 | ||
| 220 | 0 | dh_int = 0.5*(h(i,j,k-1) + h(i,j,k)) + GV%H_subroundoff |
| 221 | 0 | N2(K) = g_o_rho0 * (rhok - rhokm1) / dh_int ! Can be negative |
| 222 | ||
| 223 | enddo | |
| 224 | ||
| 225 | 0 | iFaceHeight(1) = 0.0 ! BBL is all relative to the surface |
| 226 | 0 | hcorr = 0.0 |
| 227 | ! compute heights at cell center and interfaces | |
| 228 | 0 | do k=1,GV%ke |
| 229 | 0 | dh = dz(i,k) ! Nominal thickness to use for increment, in the units of heights |
| 230 | 0 | dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) |
| 231 | 0 | hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 |
| 232 | 0 | dh = max(dh, CS%min_thickness) ! Limited increment dh>=min_thickness |
| 233 | 0 | cellHeight(k) = iFaceHeight(k) - 0.5 * dh |
| 234 | 0 | iFaceHeight(k+1) = iFaceHeight(k) - dh |
| 235 | enddo | |
| 236 | ||
| 237 | ! gets index of the level and interface above hbl | |
| 238 | 0 | kOBL = CVMix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight, hbl(i,j)) |
| 239 | ||
| 240 | 0 | kv_col(:) = 0.0 ; kd_col(:) = 0.0 |
| 241 | call CVMix_coeffs_conv(Mdiff_out=kv_col(:), & | |
| 242 | Tdiff_out=kd_col(:), & | |
| 243 | Nsqr=N2(:), & | |
| 244 | dens=rho_1d(:), & | |
| 245 | dens_lwr=rho_lwr(:), & | |
| 246 | nlev=GV%ke, & | |
| 247 | max_nlev=GV%ke, & | |
| 248 | 0 | OBL_ind=kOBL) |
| 249 | ||
| 250 | ! Increment the diffusivity outside of the boundary layer. | |
| 251 | 0 | do K=max(1,kOBL+1),GV%ke+1 |
| 252 | 0 | Kd(i,j,K) = Kd(i,j,K) + GV%m2_s_to_HZ_T * kd_col(K) |
| 253 | enddo | |
| 254 | 0 | if (present(Kd_aux)) then |
| 255 | ! Increment the other diffusivity outside of the boundary layer. | |
| 256 | 0 | do K=max(1,kOBL+1),GV%ke+1 |
| 257 | 0 | Kd_aux(i,j,K) = Kd_aux(i,j,K) + GV%m2_s_to_HZ_T * kd_col(K) |
| 258 | enddo | |
| 259 | endif | |
| 260 | ||
| 261 | ! Increment the viscosity outside of the boundary layer. | |
| 262 | 0 | do K=max(1,kOBL+1),GV%ke+1 |
| 263 | 0 | Kv(i,j,K) = Kv(i,j,K) + GV%m2_s_to_HZ_T * kv_col(K) |
| 264 | enddo | |
| 265 | ||
| 266 | ! Store 3-d arrays for diagnostics. | |
| 267 | 0 | if (CS%id_kv_conv > 0) then |
| 268 | ! Do not apply mixing due to convection within the boundary layer | |
| 269 | 0 | do K=max(1,kOBL+1),GV%ke+1 |
| 270 | 0 | Kv_conv(i,j,K) = US%m2_s_to_Z2_T * kv_col(K) |
| 271 | enddo | |
| 272 | endif | |
| 273 | 0 | if (CS%id_kd_conv > 0) then |
| 274 | ! Do not apply mixing due to convection within the boundary layer | |
| 275 | 0 | do K=max(1,kOBL+1),GV%ke+1 |
| 276 | 0 | Kd_conv(i,j,K) = US%m2_s_to_Z2_T * kd_col(K) |
| 277 | enddo | |
| 278 | endif | |
| 279 | ||
| 280 | 0 | if (CS%id_N2 > 0) then ; do k=2,GV%ke ; N2_3d(i,j,K) = US%T_to_s**2*N2(K) ; enddo ; endif |
| 281 | ||
| 282 | enddo | |
| 283 | enddo | |
| 284 | ||
| 285 | 0 | if (CS%debug) then |
| 286 | ! if (CS%id_N2 > 0) call hchksum(N2_3d, "MOM_CVMix_conv: N2", G%HI, haloshift=0, unscale=US%s_to_T**2) | |
| 287 | ! if (CS%id_kd_conv > 0) & | |
| 288 | ! call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, unscale=US%Z2_T_to_m2_s) | |
| 289 | ! if (CS%id_kv_conv > 0) & | |
| 290 | ! call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, unscale=US%Z2_T_to_m2_s) | |
| 291 | 0 | call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s) |
| 292 | 0 | call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s) |
| 293 | endif | |
| 294 | ||
| 295 | ! send diagnostics to post_data | |
| 296 | 0 | if (CS%id_N2 > 0) call post_data(CS%id_N2, N2_3d, CS%diag) |
| 297 | 0 | if (CS%id_kd_conv > 0) call post_data(CS%id_kd_conv, Kd_conv, CS%diag) |
| 298 | 0 | if (CS%id_kv_conv > 0) call post_data(CS%id_kv_conv, Kv_conv, CS%diag) |
| 299 | ||
| 300 | 0 | end subroutine calculate_CVMix_conv |
| 301 | ||
| 302 | !> Reads the parameter "USE_CVMix_CONVECTION" and returns state. | |
| 303 | !! This function allows other modules to know whether this parameterization will | |
| 304 | !! be used without needing to duplicate the log entry. | |
| 305 | 1 | logical function CVMix_conv_is_used(param_file) |
| 306 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 307 | call get_param(param_file, mdl, "USE_CVMix_CONVECTION", CVMix_conv_is_used, & | |
| 308 | 1 | default=.false., do_not_log=.true.) |
| 309 | ||
| 310 | 1 | end function CVMix_conv_is_used |
| 311 | ||
| 312 | 0 | end module MOM_CVMix_conv |