Phillips_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!> Initialization for the "Phillips" channel configuration
7
8use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
10use mom_file_parser, only : get_param, log_version, param_file_type
11use mom_get_input, only : directories
12use mom_grid, only : ocean_grid_type
14use mom_tracer_registry, only : tracer_registry_type
18
19implicit none ; private
20
21#include <MOM_memory.h>
22
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! This include declares and sets the variable "version".
34#include "version_variable.h"
35
36contains
37
38!> Initialize the thickness field for the Phillips model test case.
39subroutine phillips_initialize_thickness(h, depth_tot, G, GV, US, param_file, just_read)
40 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
41 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
42 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
43 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
44 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
45 real, dimension(SZI_(G),SZJ_(G)), &
46 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
47 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
48 !! to parse for model parameter values.
49 logical, intent(in) :: just_read !< If true, this call will only read
50 !! parameters without changing h.
51
52 real :: eta0(szk_(gv)+1) ! The 1-d nominal positions of the interfaces [Z ~> m]
53 real :: eta_im(szj_(g),szk_(gv)+1) ! A temporary array for zonal-mean eta [Z ~> m]
54 real :: eta1d(szk_(gv)+1) ! Interface height relative to the sea surface, positive upward [Z ~> m]
55 real :: jet_width ! The width of the zonal-mean jet in the same units as geolat, often [km]
56 real :: jet_height ! The interface height scale associated with the zonal-mean jet [Z ~> m]
57 real :: y_2 ! The y-position relative to the center of the domain in the same units as
58 ! geolat, often [km]
59 real :: half_strat ! The fractional depth where the stratification is centered [nondim]
60 real :: half_depth ! The depth where the stratification is centered [Z ~> m]
61 real :: km_to_grid_unit ! The conversion factor from km to the units of latitude, often 1 [nondim],
62 ! but this could be 1000 [m km-1]
63 logical :: reentrant_y ! If true, model is re-entrant in the y direction
64 character(len=40) :: mdl = "Phillips_initialize_thickness" ! This subroutine's name.
65 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
66 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
67
68 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
69 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
70
71 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "Phillips_initialization: "//&
72 "Phillips_initialize_thickness is only set to work with Cartesian axis units.")
73 if (abs(g%grid_unit_to_L*us%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km.
74 km_to_grid_unit = 1.0
75 elseif (abs(g%grid_unit_to_L*us%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m.
76 km_to_grid_unit = 1000.0
77 else
78 call mom_error(fatal, "Phillips_initialization: "//&
79 "Phillips_initialize_thickness is not recognizing the value of G%grid_unit_to_L.")
80 endif
81
82 eta_im(:,:) = 0.0
83
84 if (.not.just_read) call log_version(param_file, mdl, version)
85 call get_param(param_file, mdl, "HALF_STRAT_DEPTH", half_strat, &
86 "The fractional depth where the stratification is centered.", &
87 units="nondim", default=0.5, do_not_log=just_read)
88 call get_param(param_file, mdl, "JET_WIDTH", jet_width, &
89 "The width of the zonal-mean jet.", units="km", scale=km_to_grid_unit, &
90 fail_if_missing=.not.just_read, do_not_log=just_read)
91 call get_param(param_file, mdl, "JET_HEIGHT", jet_height, &
92 "The interface height scale associated with the zonal-mean jet.", &
93 units="m", scale=us%m_to_Z, fail_if_missing=.not.just_read, do_not_log=just_read)
94
95 ! If re-entrant in the Y direction, we use a sine function instead of a
96 ! tanh. The ratio len_lat/jet_width should be an integer in this case.
97 call get_param(param_file, mdl, "REENTRANT_Y", reentrant_y, &
98 default=.false., do_not_log=.true.)
99
100 if (just_read) return ! All run-time parameters have been read, so return.
101
102 half_depth = g%max_depth*half_strat
103 eta0(1) = 0.0 ; eta0(nz+1) = -g%max_depth
104 do k=2,1+nz/2 ; eta0(k) = -half_depth*(2.0*(k-1)/real(nz)) ; enddo
105 do k=2+nz/2,nz+1
106 eta0(k) = -g%max_depth - 2.0*(g%max_depth-half_depth) * ((k-(nz+1))/real(nz))
107 enddo
108 pi = 4.0*atan(1.0)
109
110 do j=js,je
111 eta_im(j,1) = 0.0 ; eta_im(j,nz+1) = -g%max_depth
112 enddo
113 do k=2,nz ; do j=js,je
114 y_2 = g%geoLatT(is,j) - g%south_lat - 0.5*g%len_lat
115 eta_im(j,k) = eta0(k) + jet_height * tanh(y_2 / jet_width)
116 ! or ... + jet_height * atan(y_2 / jet_width)
117 if (reentrant_y) then
118 y_2 = 2.*pi*y_2
119 eta_im(j,k) = eta0(k) + jet_height * sin(y_2 / jet_width)
120 endif
121 if (eta_im(j,k) > 0.0) eta_im(j,k) = 0.0
122 if (eta_im(j,k) < -g%max_depth) eta_im(j,k) = -g%max_depth
123 enddo ; enddo
124
125 do j=js,je ; do i=is,ie
126 ! This sets the initial thickness in [H ~> m or kg m-2] of the layers. The
127 ! thicknesses are set to insure that: 1. each layer is at least an Angstrom thick, and
128 ! 2. the interfaces are where they should be based on the resting depths and interface
129 ! height perturbations, as long at this doesn't interfere with 1.
130 eta1d(nz+1) = -depth_tot(i,j)
131 do k=nz,1,-1
132 eta1d(k) = eta_im(j,k)
133 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
134 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
135 h(i,j,k) = gv%Angstrom_Z
136 else
137 h(i,j,k) = eta1d(k) - eta1d(k+1)
138 endif
139 enddo
140 enddo ; enddo
141
143
144!> Initialize the velocity fields for the Phillips model test case
145subroutine phillips_initialize_velocity(u, v, G, GV, US, param_file, just_read)
146 type(ocean_grid_type), intent(in) :: g !< Grid structure
147 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
148 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
149 intent(out) :: u !< i-component of velocity [L T-1 ~> m s-1]
150 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
151 intent(out) :: v !< j-component of velocity [L T-1 ~> m s-1]
152 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
153 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
154 !! parse for modelparameter values.
155 logical, intent(in) :: just_read !< If true, this call will only read
156 !! parameters without changing u & v.
157
158 real :: jet_width_grid ! The width of the zonal-mean jet in the same units as geolat, often [km]
159 real :: jet_width_l ! The width of the zonal-mean jet [L ~> m]
160 real :: i_jet_width ! The inverse of the width of the zonal-mean jet [L-1 ~> m-1]
161 real :: jet_height ! The interface height scale associated with the zonal-mean jet [Z ~> m]
162 real :: x_2 ! The x-position relative to the center of the domain normalized by the
163 ! domain width [nondim]
164 real :: y_2_grid ! The y-position relative to the center of the domain in the same units
165 ! as geolat, often [km]
166 real :: y_2_l ! The y-position relative to the center of the domain [L ~> m]
167 real :: y_2_norm ! The y-position relative to the center of the domain normalized by the
168 ! domain width [nondim]
169 real :: velocity_amplitude ! The amplitude of velocity perturbations [L T-1 ~> m s-1]
170 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
171 real :: km_to_grid_unit ! The conversion factor from km to the units of latitude, often 1 [nondim],
172 ! but this could be 1000 [m km-1]
173 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
174 integer :: answer_date ! The vintage of the expressions in the Phillips_initialization code.
175 ! Values below 20250101 recover the answers from the end of 2018, while
176 ! higher values use mathematically equivalent expressions that are fully
177 ! rescalable.
178 integer :: i, j, k, is, ie, js, je, nz, m
179 logical :: reentrant_y ! If true, model is re-entrant in the y direction
180 character(len=40) :: mdl = "Phillips_initialize_velocity" ! This subroutine's name.
181 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
182
183 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "Phillips_initialization: "//&
184 "Phillips_initialize_velocity is only set to work with Cartesian axis units.")
185 if (abs(g%grid_unit_to_L*us%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km.
186 km_to_grid_unit = 1.0
187 elseif (abs(g%grid_unit_to_L*us%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m.
188 km_to_grid_unit = 1000.0
189 else
190 call mom_error(fatal, "Phillips_initialization: "//&
191 "Phillips_initialize_velocity is not recognizing the value of G%grid_unit_to_L.")
192 endif
193
194 if (.not.just_read) call log_version(param_file, mdl, version)
195 call get_param(param_file, mdl, "VELOCITY_IC_PERTURB_AMP", velocity_amplitude, &
196 "The magnitude of the initial velocity perturbation.", &
197 units="m s-1", default=0.001, scale=us%m_s_to_L_T, do_not_log=just_read)
198 call get_param(param_file, mdl, "JET_WIDTH", jet_width_l, &
199 "The width of the zonal-mean jet.", units="km", scale=1000.0*us%m_to_L, &
200 fail_if_missing=.not.just_read, do_not_log=just_read)
201 call get_param(param_file, mdl, "JET_WIDTH", jet_width_grid, &
202 "The width of the zonal-mean jet.", units="km", scale=km_to_grid_unit, &
203 fail_if_missing=.not.just_read, do_not_log=just_read)
204 call get_param(param_file, mdl, "JET_HEIGHT", jet_height, &
205 "The interface height scale associated with the zonal-mean jet.", &
206 units="m", scale=us%m_to_Z, fail_if_missing=.not.just_read, do_not_log=just_read)
207 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
208 "This sets the default value for the various _ANSWER_DATE parameters.", &
209 default=99991231)
210 call get_param(param_file, mdl, "PHILLIPS_ANSWER_DATE", answer_date, &
211 "The vintage of the expressions in the Phillips_initialization code. Values "//&
212 "below 20250101 recover the answers from the end of 2018, while higher "//&
213 "values use mathematically equivalent expressions that are fully rescalable.", &
214 default=min(20241201,default_answer_date)) !### Change this to default=default_answer_date)
215 ! If re-entrant in the Y direction, we use a sine function instead of a
216 ! tanh. The ratio len_lat/jet_width_grid should be an integer in this case.
217 call get_param(param_file, mdl, "REENTRANT_Y", reentrant_y, &
218 default=.false., do_not_log=.true.)
219
220 if (just_read) return ! All run-time parameters have been read, so return.
221
222 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, 'Phillips_initialization.F90: '// &
223 "Phillips_initialize_velocity() is only set to work with Cartesian axis units.")
224
225 u(:,:,:) = 0.0
226 v(:,:,:) = 0.0
227
228 pi = 4.0*atan(1.0)
229
230 ! Use thermal wind shear to give a geostrophically balanced flow.
231 if (answer_date < 20250101) then
232 do k=nz-1,1 ; do j=js,je ; do i=is-1,ie
233 y_2_grid = g%geoLatCu(i,j) - g%south_lat - 0.5*g%len_lat
234 if (reentrant_y) then
235 y_2_grid = 2.*pi*y_2_grid
236 u(i,j,k) = u(i,j,k+1) + (1.e-3 * (jet_height / (us%m_to_L*jet_width_grid)) * &
237 cos(y_2_grid/jet_width_grid) )
238 else
239 ! This uses d/d y_2 atan(y_2 / jet_width)
240 ! u(I,j,k) = u(I,j,k+1) + ( jet_height / &
241 ! (1.0e3*US%m_to_L*jet_width_grid * (1.0 + (y_2_grid / jet_width_grid)**2))) * &
242 ! (2.0 * GV%g_prime(K+1) / (G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1)))
243 ! This uses d/d y_2 tanh(y_2 / jet_width)
244 u(i,j,k) = u(i,j,k+1) + (1e-3 * (jet_height / (us%m_to_L*jet_width_grid)) * &
245 (sech(y_2_grid / jet_width_grid))**2 ) * &
246 (2.0 * gv%g_prime(k+1) / (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)))
247 endif
248 enddo ; enddo ; enddo
249 else
250 i_jet_width = 1.0 / jet_width_l
251 do k=nz-1,1 ; do j=js,je ; do i=is-1,ie
252 y_2_l = (g%geoLatCu(i,j) - (g%south_lat + 0.5*g%len_lat)) * g%grid_unit_to_L
253 if (reentrant_y) then
254 u(i,j,k) = u(i,j,k+1) + ((jet_height * i_jet_width) * cos(2.*pi*(y_2_l*i_jet_width)) )
255 else
256 ! This uses d/d y_2 atan(y_2 / jet_width)
257 ! u(I,j,k) = u(I,j,k+1) + ( (jet_height*I_jet_width) / (1.0 + (y_2_L*I_jet_width)**2)) * &
258 ! (2.0 * GV%g_prime(K+1) / (G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1)))
259 ! This uses d/d y_2_L tanh(y_2_L*I_jet_width)
260 u(i,j,k) = u(i,j,k+1) + ((jet_height * i_jet_width) * (sech(y_2_l*i_jet_width))**2 ) * &
261 (2.0 * gv%g_prime(k+1) / (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)))
262 endif
263 enddo ; enddo ; enddo
264 endif
265
266 do k=1,nz ; do j=js,je ; do i=is-1,ie
267 y_2_norm = (g%geoLatCu(i,j) - g%south_lat - 0.5*g%len_lat) / g%len_lat
268 x_2 = (g%geoLonCu(i,j) - g%west_lon - 0.5*g%len_lon) / g%len_lon
269 if (g%geoLonCu(i,j) == g%west_lon) then
270 ! This modification is required so that the perturbations are identical for
271 ! symmetric and non-symmetric memory. It is exactly equivalent to
272 ! taking the longitude at the eastern edge of the domain, so that x_2 ~= 0.5.
273 x_2 = ((g%west_lon + g%len_lon*real(g%ieg-(g%isg-1))/real(g%Domain%niglobal)) - &
274 g%west_lon - 0.5*g%len_lon) / g%len_lon
275 endif
276 u(i,j,k) = u(i,j,k) + velocity_amplitude * ((real(k)-0.5)/real(nz)) * &
277 (0.5 - abs(2.0*x_2) + 0.1*abs(cos(10.0*pi*x_2)) - abs(sin(5.0*pi*y_2_norm)))
278 do m=1,10
279 u(i,j,k) = u(i,j,k) + 0.2*velocity_amplitude * ((real(k)-0.5)/real(nz)) * &
280 cos(2.0*m*pi*x_2 + 2*m) * cos(6.0*pi*y_2_norm)
281 enddo
282 enddo ; enddo ; enddo
283
284end subroutine phillips_initialize_velocity
285
286!> Sets up the inverse restoration time (Idamp), and the values towards which the interface
287!! heights and an arbitrary number of tracers should be restored within each sponge for the Phillips
288!! model test case
289subroutine phillips_initialize_sponges(G, GV, US, tv, param_file, CSp, h)
290 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
291 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
292 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
293 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers
294 !! to any available thermodynamic
295 !! fields, potential temperature and
296 !! salinity or mixed layer density.
297 !! Absent fields have NULL ptrs.
298 type(param_file_type), intent(in) :: param_file !< A structure indicating the
299 !! open file to parse for model
300 !! parameter values.
301 type(sponge_cs), pointer :: csp !< A pointer that is set to point to
302 !! the control structure for the
303 !! sponge module.
304 real, intent(in), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< Thickness field [H ~> m or kg m-2].
305
306 ! Local variables
307 real :: eta0(szk_(gv)+1) ! The 1-d nominal positions of the interfaces [Z ~> m]
308 real :: eta(szi_(g),szj_(g),szk_(gv)+1) ! A temporary array for interface heights [Z ~> m].
309 real :: temp(szi_(g),szj_(g),szk_(gv)) ! A temporary array for other variables [various]
310 real :: idamp(szi_(g),szj_(g)) ! The sponge damping rate [T-1 ~> s-1]
311 real :: eta_im(szj_(g),szk_(gv)+1) ! A temporary array for zonal-mean eta [Z ~> m].
312 real :: idamp_im(szj_(g)) ! The inverse zonal-mean damping rate [T-1 ~> s-1].
313 real :: damp_rate ! The inverse zonal-mean damping rate [T-1 ~> s-1].
314 real :: jet_width ! The width of the zonal mean jet in the same units as geolat, often [km]
315 real :: jet_height ! The interface height scale associated with the zonal-mean jet [Z ~> m].
316 real :: y_2 ! The y-position relative to the channel center in the same units as
317 ! geolat, often [km]
318 real :: half_strat ! The fractional depth where the straficiation is centered [nondim].
319 real :: half_depth ! The depth where the stratification is centered [Z ~> m].
320 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
321 real :: km_to_grid_unit ! The conversion factor from km to the units of latitude, often 1 [nondim],
322 ! but this could be 1000 [m km-1]
323 logical :: reentrant_y ! If true, model is re-entrant in the y direction
324 character(len=40) :: mdl = "Phillips_initialize_sponges" ! This subroutine's name.
325
326 integer :: j, k, is, ie, js, je, isd, ied, jsd, jed, nz
327 logical, save :: first_call = .true.
328
329 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
330 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
331
332 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "Phillips_initialization: "//&
333 "Phillips_initialize_sponges is only set to work with Cartesian axis units.")
334 if (abs(g%grid_unit_to_L*us%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km.
335 km_to_grid_unit = 1.0
336 elseif (abs(g%grid_unit_to_L*us%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m.
337 km_to_grid_unit = 1000.0
338 else
339 call mom_error(fatal, "Phillips_initialization: "//&
340 "Phillips_initialize_sponges is not recognizing the value of G%grid_unit_to_L.")
341 endif
342
343 eta(:,:,:) = 0.0 ; temp(:,:,:) = 0.0 ; idamp(:,:) = 0.0
344 eta_im(:,:) = 0.0 ; idamp_im(:) = 0.0
345
346 if (first_call) call log_version(param_file, mdl, version)
347 first_call = .false.
348 call get_param(param_file, mdl, "HALF_STRAT_DEPTH", half_strat, &
349 "The fractional depth where the stratificaiton is centered.", &
350 units="nondim", default=0.5)
351 call get_param(param_file, mdl, "SPONGE_RATE", damp_rate, &
352 "The rate at which the zonal-mean sponges damp.", &
353 units="s-1", default=1.0/(10.0*86400.0), scale=us%T_to_s)
354
355 call get_param(param_file, mdl, "JET_WIDTH", jet_width, &
356 "The width of the zonal-mean jet.", units="km", scale=km_to_grid_unit, &
357 fail_if_missing=.true.)
358 call get_param(param_file, mdl, "JET_HEIGHT", jet_height, &
359 "The interface height scale associated with the zonal-mean jet.", &
360 units="m", scale=us%m_to_Z, fail_if_missing=.true.)
361 ! If re-entrant in the Y direction, we use a sine function instead of a
362 ! tanh. The ratio len_lat/jet_width should be an integer in this case.
363 call get_param(param_file, mdl, "REENTRANT_Y", reentrant_y, &
364 default=.false., do_not_log=.true.)
365
366 half_depth = g%max_depth*half_strat
367 eta0(1) = 0.0 ; eta0(nz+1) = -g%max_depth
368 do k=2,1+nz/2 ; eta0(k) = -half_depth*(2.0*(k-1)/real(nz)) ; enddo
369 do k=2+nz/2,nz+1
370 eta0(k) = -g%max_depth - 2.0*(g%max_depth-half_depth) * ((k-(nz+1))/real(nz))
371 enddo
372 pi = 4.0*atan(1.0)
373
374 do j=js,je
375 idamp_im(j) = damp_rate
376 eta_im(j,1) = 0.0 ; eta_im(j,nz+1) = -g%max_depth
377 enddo
378 do k=2,nz ; do j=js,je
379 y_2 = g%geoLatT(is,j) - g%south_lat - 0.5*g%len_lat
380 eta_im(j,k) = eta0(k) + jet_height * tanh(y_2 / jet_width)
381 if (reentrant_y) then
382 y_2 = 2.*pi*y_2
383 eta_im(j,k) = eta0(k) + jet_height * sin(y_2 / jet_width)
384 endif
385 if (eta_im(j,k) > 0.0) eta_im(j,k) = 0.0
386 if (eta_im(j,k) < -g%max_depth) eta_im(j,k) = -g%max_depth
387 enddo ; enddo
388
389 call initialize_sponge(idamp, eta, g, param_file, csp, gv, idamp_im, eta_im)
390
391end subroutine phillips_initialize_sponges
392
393!> sech calculates the hyperbolic secant.
394function sech(x)
395 real, intent(in) :: x !< Input value [nondim].
396 real :: sech !< Result [nondim].
397
398 ! This is here to prevent overflows or underflows.
399 if (abs(x) > 228.) then
400 sech = 0.0
401 else
402 sech = 2.0 / (exp(x) + exp(-x))
403 endif
404end function sech
405
406!> Initialize topography.
407subroutine phillips_initialize_topography(D, G, param_file, max_depth, US)
408 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
409 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
410 intent(out) :: d !< Ocean bottom depth [Z ~> m]
411 type(param_file_type), intent(in) :: param_file !< Parameter file structure
412 real, intent(in) :: max_depth !< Maximum model depth [Z ~> m]
413 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
414
415 ! Local variables
416 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
417 real :: htop ! The maximum height of the topography above max_depth [Z ~> m]
418 real :: wtop ! meridional width of topographic features [km]
419 real :: ltop ! zonal width of topographic features [km]
420 real :: offset ! meridional offset from the center of topographic features [km]
421 real :: dist ! zonal width of topographic features [km]
422 real :: x1, x2, x3, x4, y1, y2 ! Various positions in the domain [km]
423 integer :: i, j, is, ie, js, je
424 character(len=40) :: mdl = "Phillips_initialize_topography" ! This subroutine's name.
425
426 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
427
428 pi = 4.0*atan(1.0)
429
430 call get_param(param_file, mdl, "PHILLIPS_HTOP", htop, &
431 "The maximum height of the topography.", units="m", scale=us%m_to_Z, &
432 fail_if_missing=.true.)
433! Htop=0.375*max_depth ! max height of topog. above max_depth
434 wtop = 0.5*g%len_lat ! meridional width of drake and mount
435 ltop = 0.25*g%len_lon ! zonal width of topographic features
436 offset = 0.1*g%len_lat ! meridional offset from center
437 dist = 0.333*g%len_lon ! distance between drake and mount, this should be longer than Ltop/2
438
439 y1 = g%south_lat+0.5*g%len_lat+offset-0.5*wtop ; y2 = y1+wtop
440 x1 = g%west_lon+0.1*g%len_lon ; x2 = x1+ltop ; x3 = x1+dist ; x4 = x3+3.0/2.0*ltop
441
442 do j=js,je ; do i=is,ie
443 d(i,j)=0.0
444 if (g%geoLonT(i,j)>x1 .and. g%geoLonT(i,j)<x2) then
445 d(i,j) = htop*sin(pi*(g%geoLonT(i,j)-x1)/(x2-x1))**2
446 if (g%geoLatT(i,j)>y1 .and. g%geoLatT(i,j)<y2) then
447 d(i,j) = d(i,j)*(1-sin(pi*(g%geoLatT(i,j)-y1)/(y2-y1))**2)
448 endif
449 elseif (g%geoLonT(i,j)>x3 .and. g%geoLonT(i,j)<x4 .and. &
450 g%geoLatT(i,j)>y1 .and. g%geoLatT(i,j)<y2) then
451 d(i,j) = 2.0/3.0*htop*sin(pi*(g%geoLonT(i,j)-x3)/(x4-x3))**2 &
452 *sin(pi*(g%geoLatT(i,j)-y1)/(y2-y1))**2
453 endif
454 d(i,j) = max_depth - d(i,j)
455 enddo ; enddo
456
458
459!> \namespace phillips_initialization
460!!
461!! By Robert Hallberg, April 1994 - June 2002
462!!
463!! This subroutine initializes the fields for the simulations.
464!! The one argument passed to initialize, Time, is set to the
465!! current time of the simulation. The fields which are initialized
466!! here are:
467!! u - Zonal velocity [L T-1 ~> m s-1].
468!! v - Meridional velocity [L T-1 ~> m s-1].
469!! h - Layer thickness [H ~> m or kg m-2] (must be positive)
470!! D - Basin depth [Z ~> m] (positive downward)
471!! f - The Coriolis parameter [T-1 ~> s-1].
472!! If ENABLE_THERMODYNAMICS is defined:
473!! T - Temperature [C ~> degC].
474!! S - Salinity [S ~> ppt].
475!! If SPONGE is defined:
476!! A series of subroutine calls are made to set up the damping
477!! rates and reference profiles for all variables that are damped
478!! in the sponge.
479!! Any user provided tracer code is also first linked through this
480!! subroutine.
481!!
482!! Forcing-related fields (taux, tauy, buoy, ustar, etc.) are set
483!! in MOM_surface_forcing.F90.
484!!
485!! These variables are all set in the set of subroutines (in this
486!! file) Phillips_initialize_thickness, Phillips_initialize_velocity,
487!! Phillips_initialize_topography and Phillips_initialize_sponges
488!! that seet up fields that are specific to the Phillips instability
489!! test case.
490