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 | !> Initialization for the "bench mark" configuration | |
| 6 | module benchmark_initialization | |
| 7 | ||
| 8 | use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge | |
| 9 | use MOM_dyn_horgrid, only : dyn_horgrid_type | |
| 10 | use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe | |
| 11 | use MOM_file_parser, only : get_param, log_version, param_file_type | |
| 12 | use MOM_get_input, only : directories | |
| 13 | use MOM_grid, only : ocean_grid_type | |
| 14 | use MOM_tracer_registry, only : tracer_registry_type | |
| 15 | use MOM_unit_scaling, only : unit_scale_type | |
| 16 | use MOM_variables, only : thermo_var_ptrs | |
| 17 | use MOM_verticalGrid, only : verticalGrid_type | |
| 18 | use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type | |
| 19 | ||
| 20 | implicit none ; private | |
| 21 | ||
| 22 | #include <MOM_memory.h> | |
| 23 | ||
| 24 | public benchmark_initialize_topography | |
| 25 | public benchmark_initialize_thickness | |
| 26 | public benchmark_init_temperature_salinity | |
| 27 | ||
| 28 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 29 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 30 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 31 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 32 | ||
| 33 | contains | |
| 34 | ||
| 35 | !> This subroutine sets up the benchmark test case topography. | |
| 36 | 1 | subroutine benchmark_initialize_topography(D, G, param_file, max_depth, US) |
| 37 | type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type | |
| 38 | real, dimension(G%isd:G%ied,G%jsd:G%jed), & | |
| 39 | intent(out) :: D !< Ocean bottom depth [Z ~> m] | |
| 40 | type(param_file_type), intent(in) :: param_file !< Parameter file structure | |
| 41 | real, intent(in) :: max_depth !< Maximum model depth [Z ~> m] | |
| 42 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 43 | ||
| 44 | ! Local variables | |
| 45 | real :: min_depth ! The minimum basin depth [Z ~> m] | |
| 46 | real :: PI ! 3.1415926... calculated as 4*atan(1) [nondim] | |
| 47 | real :: D0 ! A constant to make the maximum basin depth MAXIMUM_DEPTH [Z ~> m] | |
| 48 | real :: x ! Longitude relative to the domain edge, normalized by its extent [nondim] | |
| 49 | real :: y ! Latitude relative to the domain edge, normalized by its extent [nondim] | |
| 50 | ! This include declares and sets the variable "version". | |
| 51 | # include "version_variable.h" | |
| 52 | character(len=40) :: mdl = "benchmark_initialize_topography" ! This subroutine's name. | |
| 53 | integer :: i, j, is, ie, js, je | |
| 54 | 1 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 55 | ||
| 56 | 1 | call MOM_mesg(" benchmark_initialization.F90, benchmark_initialize_topography: setting topography", 5) |
| 57 | ||
| 58 | 1 | call log_version(param_file, mdl, version, "") |
| 59 | call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & | |
| 60 | 1 | "The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z) |
| 61 | ||
| 62 | 1 | PI = 4.0*atan(1.0) |
| 63 | 1 | D0 = max_depth / 0.5 |
| 64 | ||
| 65 | ! Calculate the depth of the bottom. | |
| 66 | 64981 | do j=js,je ; do i=is,ie |
| 67 | 64800 | x = (G%geoLonT(i,j)-G%west_lon) / G%len_lon |
| 68 | 64800 | y = (G%geoLatT(i,j)-G%south_lat) / G%len_lat |
| 69 | ! This sets topography that has a reentrant channel to the south. | |
| 70 | D(i,j) = -D0 * ( y*(1.0 + 0.6*cos(4.0*PI*x)) & | |
| 71 | + 0.75*exp(-6.0*y) & | |
| 72 | 64800 | + 0.05*cos(10.0*PI*x) - 0.7 ) |
| 73 | 64800 | if (D(i,j) > max_depth) D(i,j) = max_depth |
| 74 | 64980 | if (D(i,j) < min_depth) D(i,j) = 0. |
| 75 | enddo ; enddo | |
| 76 | ||
| 77 | 1 | end subroutine benchmark_initialize_topography |
| 78 | ||
| 79 | !> Initializes layer thicknesses for the benchmark test case, | |
| 80 | !! by finding the depths of interfaces in a specified latitude-dependent | |
| 81 | !! temperature profile with an exponentially decaying thermocline on top of a | |
| 82 | !! linear stratification. | |
| 83 | 1 | subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, eqn_of_state, & |
| 84 | P_Ref, just_read) | |
| 85 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 86 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 87 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 88 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & | |
| 89 | intent(out) :: h !< The thickness that is being initialized [Z ~> m] | |
| 90 | real, dimension(SZI_(G),SZJ_(G)), & | |
| 91 | intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] | |
| 92 | type(param_file_type), intent(in) :: param_file !< A structure indicating the open file | |
| 93 | !! to parse for model parameter values. | |
| 94 | type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure | |
| 95 | real, intent(in) :: P_Ref !< The coordinate-density | |
| 96 | !! reference pressure [R L2 T-2 ~> Pa]. | |
| 97 | logical, intent(in) :: just_read !< If true, this call will | |
| 98 | !! only read parameters without changing h. | |
| 99 | ! Local variables | |
| 100 | 2 | real :: e0(SZK_(GV)+1) ! The resting interface heights, in depth units [Z ~> m], |
| 101 | ! usually negative because it is positive upward. | |
| 102 | 2 | real :: e_pert(SZK_(GV)+1) ! Interface height perturbations, positive upward, |
| 103 | ! in depth units [Z ~> m]. | |
| 104 | 2 | real :: eta1D(SZK_(GV)+1) ! Interface height relative to the sea surface |
| 105 | ! positive upward, in depth units [Z ~> m]. | |
| 106 | real :: SST ! The initial sea surface temperature [C ~> degC]. | |
| 107 | real :: S_ref ! A default value for salinities [S ~> ppt] | |
| 108 | real :: T_light ! A first guess at the temperature of the lightest layer [C ~> degC] | |
| 109 | real :: T_int ! The initial temperature of an interface [C ~> degC]. | |
| 110 | real :: ML_depth ! The specified initial mixed layer depth, in depth units [Z ~> m]. | |
| 111 | real :: thermocline_scale ! The e-folding scale of the thermocline, in depth units [Z ~> m]. | |
| 112 | real, dimension(SZK_(GV)) :: & | |
| 113 | 2 | T0, S0, & ! Profiles of temperature [C ~> degC] and salinity [S ~> ppt] |
| 114 | 2 | rho_guess, & ! Potential density at T0 & S0 [R ~> kg m-3]. |
| 115 | 2 | drho_dT, & ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1]. |
| 116 | 1 | drho_dS ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]. |
| 117 | 1 | real :: pres(SZK_(GV)) ! Reference pressure [R L2 T-2 ~> Pa]. |
| 118 | real :: a_exp ! The fraction of the overall stratification that is exponential [nondim] | |
| 119 | real :: I_ts, I_md ! Inverse lengthscales [Z-1 ~> m-1]. | |
| 120 | real :: T_frac ! A ratio of the interface temperature to the range | |
| 121 | ! between SST and the bottom temperature [nondim]. | |
| 122 | real :: err ! The normalized error between the profile's temperature and the | |
| 123 | ! interface temperature for a given z [nondim] | |
| 124 | real :: derr_dz ! The derivative of the normalized error between the profile's | |
| 125 | ! temperature and the interface temperature with z [Z-1 ~> m-1] | |
| 126 | real :: pi ! 3.1415926... calculated as 4*atan(1) [nondim] | |
| 127 | real :: z ! A work variable for the interface position [Z ~> m] | |
| 128 | ! This include declares and sets the variable "version". | |
| 129 | # include "version_variable.h" | |
| 130 | character(len=40) :: mdl = "benchmark_initialize_thickness" ! This subroutine's name. | |
| 131 | integer :: i, j, k, k1, is, ie, js, je, nz, itt | |
| 132 | ||
| 133 | 1 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 134 | ||
| 135 | 1 | if (.not.just_read) call log_version(param_file, mdl, version, "") |
| 136 | call get_param(param_file, mdl, "BENCHMARK_ML_DEPTH_IC", ML_depth, & | |
| 137 | "Initial mixed layer depth in the benchmark test case.", & | |
| 138 | 1 | units='m', default=50.0, scale=US%m_to_Z, do_not_log=just_read) |
| 139 | call get_param(param_file, mdl, "BENCHMARK_THERMOCLINE_SCALE", thermocline_scale, & | |
| 140 | "Initial thermocline depth scale in the benchmark test case.", & | |
| 141 | 1 | default=500.0, units="m", scale=US%m_to_Z, do_not_log=just_read) |
| 142 | call get_param(param_file, mdl, "BENCHMARK_T_LIGHT", T_light, & | |
| 143 | "A first guess at the temperature of the lightest layer in the benchmark test case.", & | |
| 144 | 1 | units="degC", default=29.0, scale=US%degC_to_C, do_not_log=just_read) |
| 145 | call get_param(param_file, mdl, "S_REF", S_ref, & | |
| 146 | "The uniform salinities used to initialize the benchmark test case.", & | |
| 147 | 1 | units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=just_read) |
| 148 | ||
| 149 | 1 | if (just_read) return ! This subroutine has no run-time parameters. |
| 150 | ||
| 151 | 1 | call MOM_mesg(" benchmark_initialization.F90, benchmark_initialize_thickness: setting thickness", 5) |
| 152 | ||
| 153 | 1 | k1 = GV%nk_rho_varies + 1 |
| 154 | ||
| 155 | 1 | a_exp = 0.9 |
| 156 | ||
| 157 | ! This block calculates T0(k) for the purpose of diagnosing where the | |
| 158 | ! interfaces will be found. | |
| 159 | 23 | do k=1,nz |
| 160 | 23 | pres(k) = P_Ref ; S0(k) = S_ref |
| 161 | enddo | |
| 162 | 1 | T0(k1) = T_light |
| 163 | 1 | call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state) |
| 164 | 1 | call calculate_density_derivs(T0(k1), S0(k1), pres(k1), drho_dT(k1), drho_dS(k1), eqn_of_state) |
| 165 | ||
| 166 | ! A first guess of the layers' temperatures. | |
| 167 | 23 | do k=1,nz |
| 168 | 23 | T0(k) = T0(k1) + (GV%Rlay(k) - rho_guess(k1)) / drho_dT(k1) |
| 169 | enddo | |
| 170 | ||
| 171 | ! Refine the guesses for each layer. | |
| 172 | 7 | do itt=1,6 |
| 173 | 6 | call calculate_density(T0, S0, pres, rho_guess, eqn_of_state) |
| 174 | 6 | call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state) |
| 175 | 139 | do k=1,nz |
| 176 | 138 | T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) |
| 177 | enddo | |
| 178 | enddo | |
| 179 | ||
| 180 | 1 | pi = 4.0*atan(1.0) |
| 181 | 1 | I_ts = 1.0 / thermocline_scale |
| 182 | 1 | I_md = 1.0 / G%max_depth |
| 183 | 64981 | do j=js,je ; do i=is,ie |
| 184 | SST = 0.5*(T0(k1)+T0(nz)) - 0.9*0.5*(T0(k1)-T0(nz)) * & | |
| 185 | 64800 | cos(pi*(G%geoLatT(i,j)-G%south_lat)/(G%len_lat)) |
| 186 | ||
| 187 | 1490400 | do k=1,nz ; e_pert(K) = 0.0 ; enddo |
| 188 | ||
| 189 | ! This sets the initial thickness (in [Z ~> m]) of the layers. The thicknesses | |
| 190 | ! are set to insure that: | |
| 191 | ! 1. each layer is at least GV%Angstrom_Z thick, and | |
| 192 | ! 2. the interfaces are where they should be based on the resting depths and | |
| 193 | ! interface height perturbations, as long at this doesn't interfere with 1. | |
| 194 | 64800 | eta1D(nz+1) = -depth_tot(i,j) |
| 195 | ||
| 196 | 1425600 | do k=nz,2,-1 |
| 197 | 1360800 | T_int = 0.5*(T0(k) + T0(k-1)) |
| 198 | 1360800 | T_frac = (T_int - T0(nz)) / (SST - T0(nz)) |
| 199 | ! Find the z such that T_frac = a exp(z/thermocline_scale) + (1-a) (z+D)/D | |
| 200 | 1360800 | z = 0.0 |
| 201 | 9525600 | do itt=1,6 |
| 202 | 8164800 | err = a_exp * exp(z*I_ts) + (1.0 - a_exp) * (z*I_md + 1.0) - T_frac |
| 203 | 8164800 | derr_dz = a_exp * I_ts * exp(z*I_ts) + (1.0 - a_exp) * I_md |
| 204 | 9525600 | z = z - err / derr_dz |
| 205 | enddo | |
| 206 | 1360800 | e0(K) = z |
| 207 | ! e0(K) = -ML_depth + thermocline_scale * log((T_int - T0(nz)) / (SST - T0(nz))) | |
| 208 | ||
| 209 | 1360800 | eta1D(K) = e0(K) + e_pert(K) |
| 210 | ||
| 211 | 1360800 | if (eta1D(K) > -ML_depth) eta1D(K) = -ML_depth |
| 212 | ||
| 213 | 1360800 | if (eta1D(K) < eta1D(K+1) + GV%Angstrom_Z) & |
| 214 | 1015628 | eta1D(K) = eta1D(K+1) + GV%Angstrom_Z |
| 215 | ||
| 216 | 1425600 | h(i,j,k) = max(eta1D(K) - eta1D(K+1), GV%Angstrom_Z) |
| 217 | enddo | |
| 218 | 64980 | h(i,j,1) = max(0.0 - eta1D(2), GV%Angstrom_Z) |
| 219 | ||
| 220 | enddo ; enddo | |
| 221 | ||
| 222 | end subroutine benchmark_initialize_thickness | |
| 223 | ||
| 224 | !> Initializes layer temperatures and salinities for benchmark | |
| 225 | 1 | subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, & |
| 226 | eqn_of_state, P_Ref, just_read) | |
| 227 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure | |
| 228 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 229 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature | |
| 230 | !! that is being initialized [C ~> degC] | |
| 231 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is being | |
| 232 | !! initialized [S ~> ppt] | |
| 233 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 234 | type(param_file_type), intent(in) :: param_file !< A structure indicating the | |
| 235 | !! open file to parse for | |
| 236 | !! model parameter values. | |
| 237 | type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure | |
| 238 | real, intent(in) :: P_Ref !< The coordinate-density | |
| 239 | !! reference pressure [R L2 T-2 ~> Pa] | |
| 240 | logical, intent(in) :: just_read !< If true, this call will only read | |
| 241 | !! parameters without changing T & S. | |
| 242 | ! Local variables | |
| 243 | 2 | real :: T0(SZK_(GV)) ! A profile of temperatures [C ~> degC] |
| 244 | 2 | real :: S0(SZK_(GV)) ! A profile of salinities [S ~> ppt] |
| 245 | real :: S_ref ! A default value for salinities [S ~> ppt] | |
| 246 | real :: T_light ! A first guess at the temperature of the lightest layer [C ~> degC] | |
| 247 | 2 | real :: pres(SZK_(GV)) ! Reference pressure [R L2 T-2 ~> Pa] |
| 248 | 2 | real :: drho_dT(SZK_(GV)) ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1] |
| 249 | 1 | real :: drho_dS(SZK_(GV)) ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1] |
| 250 | 1 | real :: rho_guess(SZK_(GV)) ! Potential density at T0 & S0 [R ~> kg m-3] |
| 251 | real :: PI ! 3.1415926... calculated as 4*atan(1) [nondim] | |
| 252 | real :: SST ! The initial sea surface temperature [C ~> degC] | |
| 253 | character(len=40) :: mdl = "benchmark_init_temperature_salinity" ! This subroutine's name. | |
| 254 | integer :: i, j, k, k1, is, ie, js, je, nz, itt | |
| 255 | ||
| 256 | 1 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 257 | ||
| 258 | call get_param(param_file, mdl, "S_REF", S_ref, & | |
| 259 | 1 | units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.) |
| 260 | call get_param(param_file, mdl, "BENCHMARK_T_LIGHT", T_light, & | |
| 261 | 1 | units="degC", default=29.0, scale=US%degC_to_C, do_not_log=.true.) |
| 262 | ||
| 263 | 1 | if (just_read) return ! All run-time parameters have been read, so return. |
| 264 | ||
| 265 | 1 | k1 = GV%nk_rho_varies + 1 |
| 266 | ||
| 267 | 23 | do k=1,nz |
| 268 | 23 | pres(k) = P_Ref ; S0(k) = S_ref |
| 269 | enddo | |
| 270 | ||
| 271 | 1 | T0(k1) = T_light |
| 272 | 1 | call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state) |
| 273 | 3 | call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, (/k1,k1/) ) |
| 274 | ||
| 275 | ! A first guess of the layers' temperatures. ! | |
| 276 | 23 | do k=1,nz |
| 277 | 23 | T0(k) = T0(k1) + (GV%Rlay(k) - rho_guess(k1)) / drho_dT(k1) |
| 278 | enddo | |
| 279 | ||
| 280 | ! Refine the guesses for each layer. ! | |
| 281 | 7 | do itt = 1,6 |
| 282 | 6 | call calculate_density(T0, S0, pres, rho_guess, eqn_of_state) |
| 283 | 6 | call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state) |
| 284 | 139 | do k=1,nz |
| 285 | 138 | T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) |
| 286 | enddo | |
| 287 | enddo | |
| 288 | ||
| 289 | 1429583 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 290 | 1425600 | T(i,j,k) = T0(k) |
| 291 | 1429560 | S(i,j,k) = S0(k) |
| 292 | enddo ; enddo ; enddo | |
| 293 | 1 | PI = 4.0*atan(1.0) |
| 294 | 64981 | do j=js,je ; do i=is,ie |
| 295 | SST = 0.5*(T0(k1)+T0(nz)) - 0.9*0.5*(T0(k1)-T0(nz)) * & | |
| 296 | 64800 | cos(PI*(G%geoLatT(i,j)-G%south_lat)/(G%len_lat)) |
| 297 | 324180 | do k=1,k1-1 |
| 298 | 324000 | T(i,j,k) = SST |
| 299 | enddo | |
| 300 | enddo ; enddo | |
| 301 | ||
| 302 | end subroutine benchmark_init_temperature_salinity | |
| 303 | ||
| 304 | end module benchmark_initialization |