seamount_initialization.F90

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.
7
8use mom_domains, only : sum_across_pes
10use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
11use mom_file_parser, only : get_param, param_file_type
12use mom_get_input, only : directories
13use mom_grid, only : ocean_grid_type
15use mom_tracer_registry, only : tracer_registry_type
19use regrid_consts, only : coordinatemode, default_coordinate_mode
20use regrid_consts, only : regridding_layer, regridding_zstar
21use regrid_consts, only : regridding_rho, regridding_sigma
22
23implicit none ; private
24
25#include <MOM_memory.h>
26
27character(len=40) :: mdl = "seamount_initialization" !< This module's name.
28
29! The following routines are visible to the outside world
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
39contains
40
41!> Initialization of topography.
42subroutine 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 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 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 units=g%y_ax_unit_short, default=0.)
67
68 lx = lx / g%len_lon
69 ly = ly / g%len_lat
70 rlx = 0. ; if (lx>0.) rlx = 1. / lx
71 rly = 0. ; if (ly>0.) rly = 1. / ly
72
73 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 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon - 0.5
76 y = ( g%geoLatT(i,j) - g%south_lat ) / g%len_lat - 0.5
77 d(i,j) = g%max_depth * ( 1.0 - delta * exp(-((rlx*x)**2) - ((rly*y)**2)) )
78 enddo ; enddo
79
81
82!> Initialization of thicknesses.
83!! This subroutine initializes the layer thicknesses to be uniform.
84subroutine 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 real :: e0(szk_(gv)+1) ! The resting interface heights [Z ~> m], usually
98 ! negative because it is positive upward.
99 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
108
109 if (.not.just_read) &
110 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 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 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 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 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 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 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 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 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 default=2048.0, units="m-1", scale=us%Z_to_m, do_not_log=just_read)
146 if (just_read) return ! All run-time parameters have been read, so return.
147
148 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 ( (real(k)-1.5) / real(nz-1) ) ) / s_range
157 ! Force round numbers ... the above expression has irrational factors ...
158 if (eta_ic_quanta > 0.0) &
159 e0(k) = nint(eta_ic_quanta*e0(k)) / eta_ic_quanta
160 e0(k) = min(real(1-k)*gv%Angstrom_Z, e0(k)) ! Bound by surface
161 e0(k) = max(-g%max_depth, e0(k)) ! Bound by bottom
162 enddo
163 do j=js,je ; do i=is,ie
164 eta1d(nz+1) = -depth_tot(i,j)
165 do k=nz,1,-1
166 eta1d(k) = e0(k)
167 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
168 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
169 h(i,j,k) = gv%Angstrom_Z
170 else
171 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 if (just_read) return ! All run-time parameters have been read, so return.
178 do j=js,je ; do i=is,ie
179 eta1d(nz+1) = -depth_tot(i,j)
180 do k=nz,1,-1
181 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
182 if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
183 eta1d(k) = eta1d(k+1) + min_thickness
184 h(i,j,k) = min_thickness
185 else
186 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 if (just_read) return ! All run-time parameters have been read, so return.
193 do j=js,je ; do i=is,ie
194 h(i,j,:) = depth_tot(i,j) / real(nz)
195 enddo ; enddo
196
197end select
198
200
201!> Initial values for temperature and salinity
202subroutine 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 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 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 '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 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 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 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 units="degC", default=0., scale=us%degC_to_C, do_not_log=just_read)
248
249 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 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 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 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 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 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 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 units="nondim", default=1.0, do_not_log=.true.)
266 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 k_light = gv%nk_rho_varies + 1
270 do j=js,je ; do i=is,ie
271 t(i,j,k_light) = t_light ; s(i,j,k_light) = s_light
272 enddo ; enddo
273 a1 = 2.0 * res_rat / (1.0 + res_rat)
274 do k=k_light+1,nz
275 k_frac = real(k-k_light)/real(nz-k_light)
276 frac_dense = a1 * k_frac + (1.0 - a1) * k_frac**2
277 do j=js,je ; do i=is,ie
278 t(i,j,k) = frac_dense * (t_dense - t_light) + t_light
279 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 if (just_read) return ! All run-time parameters have been read, so return.
284 do j=js,je ; do i=is,ie
285 xi0 = 0.0
286 do k = 1,nz
287 xi1 = xi0 + h(i,j,k) / g%max_depth
288 select case ( trim(density_profile) )
289 case ('linear')
290 !S(i,j,k) = S_surf + S_range * 0.5 * (xi0 + xi1)
291 s(i,j,k) = s_surf + ( 0.5 * s_range ) * (xi0 + xi1) ! Coded this way to reproduce old hard-coded answers
292 t(i,j,k) = t_surf + t_range * 0.5 * (xi0 + xi1)
293 case ('parabolic')
294 s(i,j,k) = s_surf + s_range * (2.0 / 3.0) * (xi1**3 - xi0**3) / (xi1 - xi0)
295 t(i,j,k) = t_surf + t_range * (2.0 / 3.0) * (xi1**3 - xi0**3) / (xi1 - xi0)
296 case ('exponential')
297 r = 0.8 ! small values give sharp profiles
298 s(i,j,k) = s_surf + s_range * (exp(xi1/r)-exp(xi0/r)) / (xi1 - xi0)
299 t(i,j,k) = t_surf + t_range * (exp(xi1/r)-exp(xi0/r)) / (xi1 - xi0)
300 case default
301 call mom_error(fatal, 'Unknown value for "INITIAL_DENSITY_PROFILE"')
302 end select
303 xi0 = xi1
304 enddo
305 enddo ; enddo
306 end select
307
309