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 double diffusion scheme. | |
| 6 | module MOM_CVMix_ddiff | |
| 7 | ||
| 8 | use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field | |
| 9 | use MOM_diag_mediator, only : post_data | |
| 10 | use MOM_EOS, only : calculate_density_derivs | |
| 11 | use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE | |
| 12 | use MOM_file_parser, only : openParameterBlock, closeParameterBlock | |
| 13 | use MOM_file_parser, only : get_param, log_version, param_file_type | |
| 14 | use MOM_debugging, only : hchksum | |
| 15 | use MOM_grid, only : ocean_grid_type | |
| 16 | use MOM_unit_scaling, only : unit_scale_type | |
| 17 | use MOM_variables, only : thermo_var_ptrs | |
| 18 | use MOM_verticalGrid, only : verticalGrid_type | |
| 19 | use cvmix_ddiff, only : cvmix_init_ddiff, CVMix_coeffs_ddiff | |
| 20 | use cvmix_kpp, only : CVmix_kpp_compute_kOBL_depth | |
| 21 | implicit none ; private | |
| 22 | ||
| 23 | #include <MOM_memory.h> | |
| 24 | ||
| 25 | public CVMix_ddiff_init, CVMix_ddiff_end, CVMix_ddiff_is_used, compute_ddiff_coeffs | |
| 26 | ||
| 27 | !> Control structure including parameters for CVMix double diffusion. | |
| 28 | type, public :: CVMix_ddiff_cs ; private | |
| 29 | ||
| 30 | ! Parameters | |
| 31 | real :: strat_param_max !< maximum value for the stratification parameter [nondim] | |
| 32 | real :: kappa_ddiff_s !< leading coefficient in formula for salt-fingering regime | |
| 33 | !! for salinity diffusion [Z2 T-1 ~> m2 s-1] | |
| 34 | real :: ddiff_exp1 !< interior exponent in salt-fingering regime formula [nondim] | |
| 35 | real :: ddiff_exp2 !< exterior exponent in salt-fingering regime formula [nondim] | |
| 36 | real :: mol_diff !< molecular diffusivity [Z2 T-1 ~> m2 s-1] | |
| 37 | real :: kappa_ddiff_param1 !< exterior coefficient in diffusive convection regime [nondim] | |
| 38 | real :: kappa_ddiff_param2 !< middle coefficient in diffusive convection regime [nondim] | |
| 39 | real :: kappa_ddiff_param3 !< interior coefficient in diffusive convection regime [nondim] | |
| 40 | real :: min_thickness !< Minimum thickness allowed [H ~> m or kg m-2] | |
| 41 | character(len=4) :: diff_conv_type !< type of diffusive convection to use. Options are Marmorino & | |
| 42 | !! Caldwell 1976 ("MC76"; default) and Kelley 1988, 1990 ("K90") | |
| 43 | logical :: debug !< If true, turn on debugging | |
| 44 | ||
| 45 | end type CVMix_ddiff_cs | |
| 46 | ||
| 47 | character(len=40) :: mdl = "MOM_CVMix_ddiff" !< This module's name. | |
| 48 | ||
| 49 | contains | |
| 50 | ||
| 51 | !> Initialized the CVMix double diffusion module. | |
| 52 | 1 | logical function CVMix_ddiff_init(Time, G, GV, US, param_file, diag, CS) |
| 53 | ||
| 54 | type(time_type), intent(in) :: Time !< The current time. | |
| 55 | type(ocean_grid_type), intent(in) :: G !< Grid structure. | |
| 56 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. | |
| 57 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 58 | type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle | |
| 59 | type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. | |
| 60 | type(CVMix_ddiff_cs), pointer :: CS !< This module's control structure. | |
| 61 | ||
| 62 | ! This include declares and sets the variable "version". | |
| 63 | # include "version_variable.h" | |
| 64 | ||
| 65 | 1 | if (associated(CS)) then |
| 66 | call MOM_error(WARNING, "CVMix_ddiff_init called with an associated "// & | |
| 67 | 0 | "control structure.") |
| 68 | 0 | return |
| 69 | endif | |
| 70 | ||
| 71 | ! Read parameters | |
| 72 | 1 | call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_init, default=.false., do_not_log=.true.) |
| 73 | call log_version(param_file, mdl, version, & | |
| 74 | "Parameterization of mixing due to double diffusion processes via CVMix", & | |
| 75 | 1 | all_default=.not.CVMix_ddiff_init) |
| 76 | call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_init, & | |
| 77 | "If true, turns on double diffusive processes via CVMix. "//& | |
| 78 | "Note that double diffusive processes on viscosity are ignored "//& | |
| 79 | "in CVMix, see http://cvmix.github.io/ for justification.", & | |
| 80 | 1 | default=.false.) |
| 81 | ||
| 82 | 1 | if (.not. CVMix_ddiff_init) return |
| 83 | 0 | allocate(CS) |
| 84 | ||
| 85 | 0 | call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) |
| 86 | ||
| 87 | call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, & | |
| 88 | 0 | units="m", scale=GV%m_to_H, default=0.001, do_not_log=.True.) |
| 89 | ||
| 90 | 0 | call openParameterBlock(param_file,'CVMIX_DDIFF') |
| 91 | ||
| 92 | call get_param(param_file, mdl, "STRAT_PARAM_MAX", CS%strat_param_max, & | |
| 93 | "The maximum value for the double dissusion stratification parameter", & | |
| 94 | 0 | units="nondim", default=2.55) |
| 95 | ||
| 96 | call get_param(param_file, mdl, "KAPPA_DDIFF_S", CS%kappa_ddiff_s, & | |
| 97 | "Leading coefficient in formula for salt-fingering regime for salinity diffusion.", & | |
| 98 | 0 | units="m2 s-1", default=1.0e-4, scale=US%m2_s_to_Z2_T) |
| 99 | ||
| 100 | call get_param(param_file, mdl, "DDIFF_EXP1", CS%ddiff_exp1, & | |
| 101 | "Interior exponent in salt-fingering regime formula.", & | |
| 102 | 0 | units="nondim", default=1.0) |
| 103 | ||
| 104 | call get_param(param_file, mdl, "DDIFF_EXP2", CS%ddiff_exp2, & | |
| 105 | "Exterior exponent in salt-fingering regime formula.", & | |
| 106 | 0 | units="nondim", default=3.0) |
| 107 | ||
| 108 | call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM1", CS%kappa_ddiff_param1, & | |
| 109 | "Exterior coefficient in diffusive convection regime.", & | |
| 110 | 0 | units="nondim", default=0.909) |
| 111 | ||
| 112 | call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM2", CS%kappa_ddiff_param2, & | |
| 113 | "Middle coefficient in diffusive convection regime.", & | |
| 114 | 0 | units="nondim", default=4.6) |
| 115 | ||
| 116 | call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM3", CS%kappa_ddiff_param3, & | |
| 117 | "Interior coefficient in diffusive convection regime.", & | |
| 118 | 0 | units="nondim", default=-0.54) |
| 119 | ||
| 120 | call get_param(param_file, mdl, "MOL_DIFF", CS%mol_diff, & | |
| 121 | "Molecular diffusivity used in CVMix double diffusion.", & | |
| 122 | 0 | units="m2 s-1", default=1.5e-6, scale=US%m2_s_to_Z2_T) |
| 123 | ||
| 124 | call get_param(param_file, mdl, "DIFF_CONV_TYPE", CS%diff_conv_type, & | |
| 125 | "type of diffusive convection to use. Options are Marmorino \n" //& | |
| 126 | "and Caldwell 1976 (MC76) and Kelley 1988, 1990 (K90).", & | |
| 127 | 0 | default="MC76") |
| 128 | ||
| 129 | 0 | call closeParameterBlock(param_file) |
| 130 | ||
| 131 | call cvmix_init_ddiff(strat_param_max=CS%strat_param_max, & | |
| 132 | kappa_ddiff_s=US%Z2_T_to_m2_s*CS%kappa_ddiff_s, & | |
| 133 | ddiff_exp1=CS%ddiff_exp1, & | |
| 134 | ddiff_exp2=CS%ddiff_exp2, & | |
| 135 | mol_diff=US%Z2_T_to_m2_s*CS%mol_diff, & | |
| 136 | kappa_ddiff_param1=CS%kappa_ddiff_param1, & | |
| 137 | kappa_ddiff_param2=CS%kappa_ddiff_param2, & | |
| 138 | kappa_ddiff_param3=CS%kappa_ddiff_param3, & | |
| 139 | 0 | diff_conv_type=CS%diff_conv_type) |
| 140 | ||
| 141 | 0 | end function CVMix_ddiff_init |
| 142 | ||
| 143 | !> Subroutine for computing vertical diffusion coefficients for the | |
| 144 | !! double diffusion mixing parameterization. | |
| 145 | 0 | subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho) |
| 146 | ||
| 147 | type(ocean_grid_type), intent(in) :: G !< Grid structure. | |
| 148 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. | |
| 149 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. | |
| 150 | type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. | |
| 151 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 152 | integer, intent(in) :: j !< Meridional grid index to work on. | |
| 153 | ! Kd_T and Kd_S are intent inout because only one j-row is set here, but they are essentially outputs. | |
| 154 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kd_T !< Interface double diffusion diapycnal | |
| 155 | !! diffusivity for temperature | |
| 156 | !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 157 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kd_S !< Interface double diffusion diapycnal | |
| 158 | !! diffusivity for salinity | |
| 159 | !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 160 | type(CVMix_ddiff_cs), pointer :: CS !< The control structure returned | |
| 161 | !! by a previous call to CVMix_ddiff_init. | |
| 162 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & | |
| 163 | optional, intent(inout) :: R_rho !< The density ratios at interfaces [nondim]. | |
| 164 | ||
| 165 | ! Local variables | |
| 166 | real, dimension(SZK_(GV)) :: & | |
| 167 | 0 | cellHeight, & !< Height of cell centers relative to the sea surface [H ~> m or kg m-2] |
| 168 | 0 | dRho_dT, & !< partial derivatives of density with temperature [R C-1 ~> kg m-3 degC-1] |
| 169 | 0 | dRho_dS, & !< partial derivatives of density with salinity [R S-1 ~> kg m-3 ppt-1] |
| 170 | 0 | pres_int, & !< pressure at each interface [R L2 T-2 ~> Pa] |
| 171 | 0 | temp_int, & !< temp and at interfaces [C ~> degC] |
| 172 | 0 | salt_int, & !< salt at at interfaces [S ~> ppt] |
| 173 | 0 | alpha_dT, & !< alpha*dT across interfaces [kg m-3] |
| 174 | 0 | beta_dS, & !< beta*dS across interfaces [kg m-3] |
| 175 | 0 | dT, & !< temperature difference between adjacent layers [C ~> degC] |
| 176 | 0 | dS !< salinity difference between adjacent layers [S ~> ppt] |
| 177 | real, dimension(SZK_(GV)+1) :: & | |
| 178 | 0 | Kd1_T, & !< Diapycanal diffusivity of temperature [m2 s-1]. |
| 179 | 0 | Kd1_S !< Diapycanal diffusivity of salinity [m2 s-1]. |
| 180 | ||
| 181 | 0 | real, dimension(SZK_(GV)+1) :: iFaceHeight !< Height of interfaces relative to the sea surface [H ~> m or kg m-2] |
| 182 | real :: dh, hcorr ! Limited thicknesses and a cumulative correction [H ~> m or kg m-2] | |
| 183 | integer :: i, k | |
| 184 | ||
| 185 | ! initialize dummy variables | |
| 186 | 0 | pres_int(:) = 0.0 ; temp_int(:) = 0.0 ; salt_int(:) = 0.0 |
| 187 | 0 | alpha_dT(:) = 0.0 ; beta_dS(:) = 0.0 ; dRho_dT(:) = 0.0 |
| 188 | 0 | dRho_dS(:) = 0.0 ; dT(:) = 0.0 ; dS(:) = 0.0 |
| 189 | ||
| 190 | ||
| 191 | ! GMM, I am leaving some code commented below. We need to pass BLD to | |
| 192 | ! this subroutine to avoid adding diffusivity above that. This needs | |
| 193 | ! to be done once we re-structure the order of the calls. | |
| 194 | !if (.not. associated(hbl)) then | |
| 195 | ! allocate(hbl(SZI_(G), SZJ_(G))) | |
| 196 | ! hbl(:,:) = 0.0 | |
| 197 | !endif | |
| 198 | ||
| 199 | 0 | do i = G%isc, G%iec |
| 200 | ||
| 201 | ! skip calling at land points | |
| 202 | 0 | if (G%mask2dT(i,j) == 0.) cycle |
| 203 | ||
| 204 | 0 | pres_int(1) = 0. ; if (associated(tv%p_surf)) pres_int(1) = tv%p_surf(i,j) |
| 205 | ! we don't have SST and SSS, so let's use values at top-most layer | |
| 206 | 0 | temp_int(1) = tv%T(i,j,1) ; salt_int(1) = tv%S(i,j,1) |
| 207 | 0 | do K=2,GV%ke |
| 208 | ! pressure at interface | |
| 209 | 0 | pres_int(K) = pres_int(K-1) + (GV%g_Earth * GV%H_to_RZ) * h(i,j,k-1) |
| 210 | ! temp and salt at interface | |
| 211 | ! for temp: (t1*h1 + t2*h2)/(h1+h2) | |
| 212 | 0 | temp_int(K) = (tv%T(i,j,k-1)*h(i,j,k-1) + tv%T(i,j,k)*h(i,j,k)) / (h(i,j,k-1)+h(i,j,k)) |
| 213 | 0 | salt_int(K) = (tv%S(i,j,k-1)*h(i,j,k-1) + tv%S(i,j,k)*h(i,j,k)) / (h(i,j,k-1)+h(i,j,k)) |
| 214 | ! dT and dS | |
| 215 | 0 | dT(K) = (tv%T(i,j,k-1)-tv%T(i,j,k)) |
| 216 | 0 | dS(K) = (tv%S(i,j,k-1)-tv%S(i,j,k)) |
| 217 | enddo ! k-loop finishes | |
| 218 | ||
| 219 | 0 | call calculate_density_derivs(temp_int, salt_int, pres_int, drho_dT, drho_dS, tv%eqn_of_state) |
| 220 | ||
| 221 | ! The "-1.0" below is needed so that the following criteria is satisfied: | |
| 222 | ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger" | |
| 223 | ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then "diffusive convection" | |
| 224 | 0 | do k=1,GV%ke |
| 225 | 0 | alpha_dT(k) = -1.0*US%R_to_kg_m3*drho_dT(k) * dT(k) |
| 226 | 0 | beta_dS(k) = US%R_to_kg_m3*drho_dS(k) * dS(k) |
| 227 | enddo | |
| 228 | ||
| 229 | 0 | if (present(R_rho)) then |
| 230 | 0 | do k=1,GV%ke |
| 231 | ! Set R_rho using Adcroft's rule of reciprocals. | |
| 232 | 0 | R_rho(i,j,k) = 0.0 ; if (abs(beta_dS(k)) > 0.0) R_rho(i,j,k) = alpha_dT(k) / beta_dS(k) |
| 233 | ! avoid NaN's again for safety, perhaps unnecessarily. | |
| 234 | 0 | if (R_rho(i,j,k) /= R_rho(i,j,k)) R_rho(i,j,k) = 0.0 |
| 235 | enddo | |
| 236 | endif | |
| 237 | ||
| 238 | 0 | iFaceHeight(1) = 0.0 ! BBL is all relative to the surface |
| 239 | 0 | hcorr = 0.0 |
| 240 | ! compute heights at cell center and interfaces | |
| 241 | 0 | do k=1,GV%ke |
| 242 | 0 | dh = h(i,j,k) ! Nominal thickness to use for increment, in height units |
| 243 | 0 | dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) |
| 244 | 0 | hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 |
| 245 | 0 | dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness |
| 246 | 0 | cellHeight(k) = iFaceHeight(k) - 0.5 * dh |
| 247 | 0 | iFaceHeight(k+1) = iFaceHeight(k) - dh |
| 248 | enddo | |
| 249 | ||
| 250 | ! gets index of the level and interface above hbl in [H ~> m or kg m-2] | |
| 251 | !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight, hbl(i,j)) | |
| 252 | ||
| 253 | 0 | Kd1_T(:) = 0.0 ; Kd1_S(:) = 0.0 |
| 254 | call CVMix_coeffs_ddiff(Tdiff_out=Kd1_T(:), & | |
| 255 | Sdiff_out=Kd1_S(:), & | |
| 256 | strat_param_num=alpha_dT(:), & | |
| 257 | strat_param_denom=beta_dS(:), & | |
| 258 | nlev=GV%ke, & | |
| 259 | 0 | max_nlev=GV%ke) |
| 260 | 0 | do K=1,GV%ke+1 |
| 261 | 0 | Kd_T(i,j,K) = GV%m2_s_to_HZ_T * Kd1_T(K) |
| 262 | 0 | Kd_S(i,j,K) = GV%m2_s_to_HZ_T * Kd1_S(K) |
| 263 | enddo | |
| 264 | ||
| 265 | ! Do not apply mixing due to convection within the boundary layer | |
| 266 | !do k=1,kOBL | |
| 267 | ! Kd_T(i,j,k) = 0.0 | |
| 268 | ! Kd_S(i,j,k) = 0.0 | |
| 269 | !enddo | |
| 270 | ||
| 271 | enddo ! i-loop | |
| 272 | ||
| 273 | 0 | end subroutine compute_ddiff_coeffs |
| 274 | ||
| 275 | !> Reads the parameter "USE_CVMIX_DDIFF" and returns state. | |
| 276 | !! This function allows other modules to know whether this parameterization will | |
| 277 | !! be used without needing to duplicate the log entry. | |
| 278 | 1 | logical function CVMix_ddiff_is_used(param_file) |
| 279 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 280 | call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_is_used, & | |
| 281 | 1 | default=.false., do_not_log=.true.) |
| 282 | ||
| 283 | 1 | end function CVMix_ddiff_is_used |
| 284 | ||
| 285 | !> Clear pointers and deallocate memory | |
| 286 | ! NOTE: Placeholder destructor | |
| 287 | 0 | subroutine CVMix_ddiff_end(CS) |
| 288 | type(CVMix_ddiff_cs), pointer :: CS !< Control structure for this module that | |
| 289 | !! will be deallocated in this subroutine | |
| 290 | 0 | end subroutine CVMix_ddiff_end |
| 291 | ||
| 292 | 0 | end module MOM_CVMix_ddiff |