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 | !> Configures the model for the idealized seamount test case. | |
| 6 | module seamount_initialization | |
| 7 | ||
| 8 | use MOM_domains, only : sum_across_PEs | |
| 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, param_file_type | |
| 12 | use MOM_get_input, only : directories | |
| 13 | use MOM_grid, only : ocean_grid_type | |
| 14 | use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS | |
| 15 | use MOM_tracer_registry, only : tracer_registry_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 regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE | |
| 20 | use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR | |
| 21 | use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA | |
| 22 | ||
| 23 | implicit none ; private | |
| 24 | ||
| 25 | #include <MOM_memory.h> | |
| 26 | ||
| 27 | character(len=40) :: mdl = "seamount_initialization" !< This module's name. | |
| 28 | ||
| 29 | ! The following routines are visible to the outside world | |
| 30 | public seamount_initialize_topography | |
| 31 | public seamount_initialize_thickness | |
| 32 | public seamount_initialize_temperature_salinity | |
| 33 | ||
| 34 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 35 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 36 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 37 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 38 | ||
| 39 | contains | |
| 40 | ||
| 41 | !> Initialization of topography. | |
| 42 | 0 | subroutine seamount_initialize_topography( D, G, param_file, max_depth ) |
| 43 | type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type | |
| 44 | real, dimension(G%isd:G%ied,G%jsd:G%jed), & | |
| 45 | intent(out) :: D !< Ocean bottom depth [Z ~> m] | |
| 46 | type(param_file_type), intent(in) :: param_file !< Parameter file structure | |
| 47 | real, intent(in) :: max_depth !< Maximum ocean depth [Z ~> m] | |
| 48 | ||
| 49 | ! Local variables | |
| 50 | real :: delta ! Height of the seamount as a fraction of the maximum ocean depth [nondim] | |
| 51 | real :: x, y ! Normalized positions relative to the domain center [nondim] | |
| 52 | real :: Lx, Ly ! Seamount length scales normalized by the relevant domain sizes [nondim] | |
| 53 | real :: rLx, rLy ! The Adcroft reciprocals of Lx and Ly [nondim] | |
| 54 | integer :: i, j | |
| 55 | ||
| 56 | call get_param(param_file, mdl,"SEAMOUNT_DELTA", delta, & | |
| 57 | "Non-dimensional height of seamount.", & | |
| 58 | 0 | units="nondim", default=0.5) |
| 59 | call get_param(param_file, mdl,"SEAMOUNT_X_LENGTH_SCALE", Lx, & | |
| 60 | "Length scale of seamount in x-direction. "//& | |
| 61 | "Set to zero make topography uniform in the x-direction.", & | |
| 62 | 0 | units=G%x_ax_unit_short, default=20.) |
| 63 | call get_param(param_file, mdl,"SEAMOUNT_Y_LENGTH_SCALE", Ly, & | |
| 64 | "Length scale of seamount in y-direction. "//& | |
| 65 | "Set to zero make topography uniform in the y-direction.", & | |
| 66 | 0 | units=G%y_ax_unit_short, default=0.) |
| 67 | ||
| 68 | 0 | Lx = Lx / G%len_lon |
| 69 | 0 | Ly = Ly / G%len_lat |
| 70 | 0 | rLx = 0. ; if (Lx>0.) rLx = 1. / Lx |
| 71 | 0 | rLy = 0. ; if (Ly>0.) rLy = 1. / Ly |
| 72 | ||
| 73 | 0 | do j=G%jsc,G%jec ; do i=G%isc,G%iec |
| 74 | ! Compute normalized zonal coordinates (x,y=0 at center of domain) | |
| 75 | 0 | x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5 |
| 76 | 0 | y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5 |
| 77 | 0 | D(i,j) = G%max_depth * ( 1.0 - delta * exp(-((rLx*x)**2) - ((rLy*y)**2)) ) |
| 78 | enddo ; enddo | |
| 79 | ||
| 80 | 0 | end subroutine seamount_initialize_topography |
| 81 | ||
| 82 | !> Initialization of thicknesses. | |
| 83 | !! This subroutine initializes the layer thicknesses to be uniform. | |
| 84 | 0 | subroutine seamount_initialize_thickness (h, depth_tot, G, GV, US, param_file, 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 | logical, intent(in) :: just_read !< If true, this call will only read | |
| 95 | !! parameters without changing h. | |
| 96 | ||
| 97 | 0 | real :: e0(SZK_(GV)+1) ! The resting interface heights [Z ~> m], usually |
| 98 | ! negative because it is positive upward. | |
| 99 | 0 | real :: eta1D(SZK_(GV)+1) ! Interface height relative to the sea surface, positive upward [Z ~> m] |
| 100 | real :: min_thickness ! The minimum layer thicknesses [Z ~> m]. | |
| 101 | real :: S_ref ! A default value for salinities [S ~> ppt]. | |
| 102 | real :: S_surf, S_range, S_light, S_dense ! Various salinities [S ~> ppt]. | |
| 103 | real :: eta_IC_quanta ! The granularity of quantization of intial interface heights [Z-1 ~> m-1]. | |
| 104 | character(len=20) :: verticalCoordinate | |
| 105 | integer :: i, j, k, is, ie, js, je, nz | |
| 106 | ||
| 107 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 108 | ||
| 109 | 0 | if (.not.just_read) & |
| 110 | 0 | call MOM_mesg("seamount_initialization.F90, seamount_initialize_thickness: setting thickness") |
| 111 | ||
| 112 | call get_param(param_file, mdl,"MIN_THICKNESS",min_thickness, & | |
| 113 | 'Minimum thickness for layer', & | |
| 114 | 0 | units='m', default=1.0e-3, do_not_log=just_read, scale=US%m_to_Z) |
| 115 | call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE",verticalCoordinate, & | |
| 116 | 0 | default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) |
| 117 | ||
| 118 | ! WARNING: this routine specifies the interface heights so that the last layer | |
| 119 | ! is vanished, even at maximum depth. In order to have a uniform | |
| 120 | ! layer distribution, use this line of code within the loop: | |
| 121 | ! e0(k) = -G%max_depth * real(k-1) / real(nz) | |
| 122 | ! To obtain a thickness distribution where the last layer is | |
| 123 | ! vanished and the other thicknesses uniformly distributed, use: | |
| 124 | ! e0(k) = -G%max_depth * real(k-1) / real(nz-1) | |
| 125 | !do k=1,nz+1 | |
| 126 | ! e0(k) = -G%max_depth * real(k-1) / real(nz) | |
| 127 | !enddo | |
| 128 | ||
| 129 | 0 | select case ( coordinateMode(verticalCoordinate) ) |
| 130 | ||
| 131 | case ( REGRIDDING_LAYER, REGRIDDING_RHO ) ! Initial thicknesses for isopycnal coordinates | |
| 132 | call get_param(param_file, mdl,"INITIAL_SSS", S_surf, & | |
| 133 | 0 | units="ppt", default=34., scale=US%ppt_to_S, do_not_log=.true.) |
| 134 | call get_param(param_file, mdl,"INITIAL_S_RANGE", S_range, & | |
| 135 | 0 | units="ppt", default=2., scale=US%ppt_to_S, do_not_log=.true.) |
| 136 | call get_param(param_file, mdl, "S_REF", S_ref, & | |
| 137 | 0 | units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.) |
| 138 | call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", S_light, & | |
| 139 | 0 | units="ppt", default=US%S_to_ppt*S_Ref, scale=US%ppt_to_S, do_not_log=.true.) |
| 140 | call get_param(param_file, mdl, "TS_RANGE_S_DENSE", S_dense, & | |
| 141 | 0 | units="ppt", default=US%S_to_ppt*S_Ref, scale=US%ppt_to_S, do_not_log=.true.) |
| 142 | call get_param(param_file, mdl, "INTERFACE_IC_QUANTA", eta_IC_quanta, & | |
| 143 | "The granularity of initial interface height values "//& | |
| 144 | "per meter, to avoid sensivity to order-of-arithmetic changes.", & | |
| 145 | 0 | default=2048.0, units="m-1", scale=US%Z_to_m, do_not_log=just_read) |
| 146 | 0 | if (just_read) return ! All run-time parameters have been read, so return. |
| 147 | ||
| 148 | 0 | do K=1,nz+1 |
| 149 | ! Salinity of layer k is S_light + (k-1)/(nz-1) * (S_dense - S_light) | |
| 150 | ! Salinity of interface K is S_light + (K-3/2)/(nz-1) * (S_dense - S_light) | |
| 151 | ! Salinity at depth z should be S(z) = S_surf - S_range * z/max_depth | |
| 152 | ! Equating: S_surf - S_range * z/max_depth = S_light + (K-3/2)/(nz-1) * (S_dense - S_light) | |
| 153 | ! Equating: - S_range * z/max_depth = S_light - S_surf + (K-3/2)/(nz-1) * (S_dense - S_light) | |
| 154 | ! Equating: z/max_depth = - ( S_light - S_surf + (K-3/2)/(nz-1) * (S_dense - S_light) ) / S_range | |
| 155 | e0(K) = - G%max_depth * ( ( S_light - S_surf ) + ( S_dense - S_light ) * & | |
| 156 | 0 | ( (real(K)-1.5) / real(nz-1) ) ) / S_range |
| 157 | ! Force round numbers ... the above expression has irrational factors ... | |
| 158 | 0 | if (eta_IC_quanta > 0.0) & |
| 159 | 0 | e0(K) = nint(eta_IC_quanta*e0(K)) / eta_IC_quanta |
| 160 | 0 | e0(K) = min(real(1-K)*GV%Angstrom_Z, e0(K)) ! Bound by surface |
| 161 | 0 | e0(K) = max(-G%max_depth, e0(K)) ! Bound by bottom |
| 162 | enddo | |
| 163 | 0 | do j=js,je ; do i=is,ie |
| 164 | 0 | eta1D(nz+1) = -depth_tot(i,j) |
| 165 | 0 | do k=nz,1,-1 |
| 166 | 0 | eta1D(k) = e0(k) |
| 167 | 0 | if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then |
| 168 | 0 | eta1D(k) = eta1D(k+1) + GV%Angstrom_Z |
| 169 | 0 | h(i,j,k) = GV%Angstrom_Z |
| 170 | else | |
| 171 | 0 | h(i,j,k) = eta1D(k) - eta1D(k+1) |
| 172 | endif | |
| 173 | enddo | |
| 174 | enddo ; enddo | |
| 175 | ||
| 176 | case ( REGRIDDING_ZSTAR ) ! Initial thicknesses for z coordinates | |
| 177 | 0 | if (just_read) return ! All run-time parameters have been read, so return. |
| 178 | 0 | do j=js,je ; do i=is,ie |
| 179 | 0 | eta1D(nz+1) = -depth_tot(i,j) |
| 180 | 0 | do k=nz,1,-1 |
| 181 | 0 | eta1D(k) = -G%max_depth * real(k-1) / real(nz) |
| 182 | 0 | if (eta1D(k) < (eta1D(k+1) + min_thickness)) then |
| 183 | 0 | eta1D(k) = eta1D(k+1) + min_thickness |
| 184 | 0 | h(i,j,k) = min_thickness |
| 185 | else | |
| 186 | 0 | h(i,j,k) = eta1D(k) - eta1D(k+1) |
| 187 | endif | |
| 188 | enddo | |
| 189 | enddo ; enddo | |
| 190 | ||
| 191 | case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates | |
| 192 | 0 | if (just_read) return ! All run-time parameters have been read, so return. |
| 193 | 0 | do j=js,je ; do i=is,ie |
| 194 | 0 | h(i,j,:) = depth_tot(i,j) / real(nz) |
| 195 | enddo ; enddo | |
| 196 | ||
| 197 | end select | |
| 198 | ||
| 199 | end subroutine seamount_initialize_thickness | |
| 200 | ||
| 201 | !> Initial values for temperature and salinity | |
| 202 | 1 | subroutine seamount_initialize_temperature_salinity(T, S, h, G, GV, US, param_file, just_read) |
| 203 | type(ocean_grid_type), intent(in) :: G !< Ocean grid structure | |
| 204 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure | |
| 205 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC] | |
| 206 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt] | |
| 207 | real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m] | |
| 208 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 209 | type(param_file_type), intent(in) :: param_file !< Parameter file structure | |
| 210 | logical, intent(in) :: just_read !< If true, this call will | |
| 211 | !! only read parameters without changing T & S. | |
| 212 | ||
| 213 | ! Local variables | |
| 214 | real :: xi0, xi1 ! Fractional positions within the depth range [nondim] | |
| 215 | real :: r ! A nondimensional sharpness parameter with an exponetial profile [nondim] | |
| 216 | real :: S_Ref ! Default salinity range parameters [S ~> ppt]. | |
| 217 | real :: T_Ref ! Default temperature range parameters [C ~> degC]. | |
| 218 | real :: S_Light, S_Dense, S_surf, S_range ! Salinity range parameters [S ~> ppt]. | |
| 219 | real :: T_Light, T_Dense, T_surf, T_range ! Temperature range parameters [C ~> degC]. | |
| 220 | real :: res_rat ! The ratio of density space resolution in the denser part | |
| 221 | ! of the range to that in the lighter part of the range. | |
| 222 | ! Setting this greater than 1 increases the resolution for | |
| 223 | ! the denser water [nondim]. | |
| 224 | real :: a1, frac_dense, k_frac ! Nondimensional temporary variables [nondim] | |
| 225 | integer :: i, j, k, is, ie, js, je, nz, k_light | |
| 226 | ||
| 227 | character(len=20) :: verticalCoordinate, density_profile | |
| 228 | ||
| 229 | 1 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 230 | ||
| 231 | call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalCoordinate, & | |
| 232 | 1 | default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) |
| 233 | call get_param(param_file, mdl,"INITIAL_DENSITY_PROFILE", density_profile, & | |
| 234 | 'Initial profile shape. Valid values are "linear", "parabolic" '//& | |
| 235 | 1 | 'and "exponential".', default='linear', do_not_log=just_read) |
| 236 | call get_param(param_file, mdl,"INITIAL_SSS", S_surf, & | |
| 237 | 'Initial surface salinity', & | |
| 238 | 1 | units="ppt", default=34., scale=US%ppt_to_S, do_not_log=just_read) |
| 239 | call get_param(param_file, mdl,"INITIAL_SST", T_surf, & | |
| 240 | 'Initial surface temperature', & | |
| 241 | 1 | units="degC", default=0., scale=US%degC_to_C, do_not_log=just_read) |
| 242 | call get_param(param_file, mdl,"INITIAL_S_RANGE", S_range, & | |
| 243 | 'Initial salinity range (bottom - surface)', & | |
| 244 | 1 | units="ppt", default=2., scale=US%ppt_to_S, do_not_log=just_read) |
| 245 | call get_param(param_file, mdl,"INITIAL_T_RANGE", T_range, & | |
| 246 | 'Initial temperature range (bottom - surface)', & | |
| 247 | 1 | units="degC", default=0., scale=US%degC_to_C, do_not_log=just_read) |
| 248 | ||
| 249 | 0 | select case ( coordinateMode(verticalCoordinate) ) |
| 250 | case ( REGRIDDING_LAYER ) ! Initial thicknesses for layer isopycnal coordinates | |
| 251 | ! These parameters are used in MOM_fixed_initialization.F90 when CONFIG_COORD="ts_range" | |
| 252 | call get_param(param_file, mdl, "T_REF", T_ref, & | |
| 253 | 0 | units="degC", default=10.0, scale=US%degC_to_C, do_not_log=.true.) |
| 254 | call get_param(param_file, mdl, "TS_RANGE_T_LIGHT", T_light, & | |
| 255 | 0 | units="degC", default=US%C_to_degC*T_Ref, scale=US%degC_to_C, do_not_log=.true.) |
| 256 | call get_param(param_file, mdl, "TS_RANGE_T_DENSE", T_dense, & | |
| 257 | 0 | units="degC", default=US%C_to_degC*T_Ref, scale=US%degC_to_C, do_not_log=.true.) |
| 258 | call get_param(param_file, mdl, "S_REF", S_ref, & | |
| 259 | 0 | units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.) |
| 260 | call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", S_light, & | |
| 261 | 0 | units="ppt", default=US%S_to_ppt*S_Ref, scale=US%ppt_to_S, do_not_log=.true.) |
| 262 | call get_param(param_file, mdl, "TS_RANGE_S_DENSE", S_dense, & | |
| 263 | 0 | units="ppt", default=US%S_to_ppt*S_Ref, scale=US%ppt_to_S, do_not_log=.true.) |
| 264 | call get_param(param_file, mdl, "TS_RANGE_RESOLN_RATIO", res_rat, & | |
| 265 | 0 | units="nondim", default=1.0, do_not_log=.true.) |
| 266 | 0 | if (just_read) return ! All run-time parameters have been read, so return. |
| 267 | ||
| 268 | ! Emulate the T,S used in the "ts_range" coordinate configuration code | |
| 269 | 0 | k_light = GV%nk_rho_varies + 1 |
| 270 | 0 | do j=js,je ; do i=is,ie |
| 271 | 0 | T(i,j,k_light) = T_light ; S(i,j,k_light) = S_light |
| 272 | enddo ; enddo | |
| 273 | 0 | a1 = 2.0 * res_rat / (1.0 + res_rat) |
| 274 | 0 | do k=k_light+1,nz |
| 275 | 0 | k_frac = real(k-k_light)/real(nz-k_light) |
| 276 | 0 | frac_dense = a1 * k_frac + (1.0 - a1) * k_frac**2 |
| 277 | 0 | do j=js,je ; do i=is,ie |
| 278 | 0 | T(i,j,k) = frac_dense * (T_Dense - T_Light) + T_Light |
| 279 | 0 | S(i,j,k) = frac_dense * (S_Dense - S_Light) + S_Light |
| 280 | enddo ; enddo | |
| 281 | enddo | |
| 282 | case ( REGRIDDING_SIGMA, REGRIDDING_ZSTAR, REGRIDDING_RHO ) ! All other coordinate use FV initialization | |
| 283 | 1 | if (just_read) return ! All run-time parameters have been read, so return. |
| 284 | 7262 | do j=js,je ; do i=is,ie |
| 285 | 7200 | xi0 = 0.0 |
| 286 | 547260 | do k = 1,nz |
| 287 | 540000 | xi1 = xi0 + h(i,j,k) / G%max_depth |
| 288 | 1080000 | select case ( trim(density_profile) ) |
| 289 | case ('linear') | |
| 290 | !S(i,j,k) = S_surf + S_range * 0.5 * (xi0 + xi1) | |
| 291 | 0 | S(i,j,k) = S_surf + ( 0.5 * S_range ) * (xi0 + xi1) ! Coded this way to reproduce old hard-coded answers |
| 292 | 0 | T(i,j,k) = T_surf + T_range * 0.5 * (xi0 + xi1) |
| 293 | case ('parabolic') | |
| 294 | 540000 | S(i,j,k) = S_surf + S_range * (2.0 / 3.0) * (xi1**3 - xi0**3) / (xi1 - xi0) |
| 295 | 540000 | T(i,j,k) = T_surf + T_range * (2.0 / 3.0) * (xi1**3 - xi0**3) / (xi1 - xi0) |
| 296 | case ('exponential') | |
| 297 | 0 | r = 0.8 ! small values give sharp profiles |
| 298 | 0 | S(i,j,k) = S_surf + S_range * (exp(xi1/r)-exp(xi0/r)) / (xi1 - xi0) |
| 299 | 0 | T(i,j,k) = T_surf + T_range * (exp(xi1/r)-exp(xi0/r)) / (xi1 - xi0) |
| 300 | case default | |
| 301 | 1080000 | call MOM_error(FATAL, 'Unknown value for "INITIAL_DENSITY_PROFILE"') |
| 302 | end select | |
| 303 | 547200 | xi0 = xi1 |
| 304 | enddo | |
| 305 | enddo ; enddo | |
| 306 | end select | |
| 307 | ||
| 308 | end subroutine seamount_initialize_temperature_salinity | |
| 309 | ||
| 310 | end module seamount_initialization |