Neverworld_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 "Neverworld" configuration
7
10use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
11use mom_file_parser, only : get_param, log_version, param_file_type
12use mom_get_input, only : directories
13use mom_grid, only : ocean_grid_type
14use mom_tracer_registry, only : tracer_registry_type
18
19use random_numbers_mod, only: initializerandomnumberstream, getrandomnumbers, randomnumberstream
20
21implicit none ; private
22
23#include <MOM_memory.h>
24
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
33contains
34
35!> This subroutine sets up the Neverworld test case topography.
36subroutine neverworld_initialize_topography(D, G, param_file, max_depth)
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 in the units of depth_max [A]
40 type(param_file_type), intent(in) :: param_file !< Parameter file structure
41 real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units [A]
42
43 ! Local variables
44 real :: pi ! 3.1415926... calculated as 4*atan(1) [nondim]
45 real :: x, y ! Lateral positions normalized by the domain size [nondim]
46 ! This include declares and sets the variable "version".
47# include "version_variable.h"
48 character(len=40) :: mdl = "Neverworld_initialize_topography" ! This subroutine's name.
49 real :: nl_top_amp ! Amplitude of large-scale topographic features as a fraction of the maximum depth [nondim]
50 real :: nl_roughness_amp ! Amplitude of topographic roughness as a fraction of the maximum depth [nondim]
51 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
52 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
53 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
54
55 call mom_mesg(" Neverworld_initialization.F90, Neverworld_initialize_topography: setting topography", 5)
56
57 call log_version(param_file, mdl, version, "")
58 call get_param(param_file, mdl, "NL_ROUGHNESS_AMP", nl_roughness_amp, &
59 "Amplitude of wavy signal in bathymetry.", units="nondim", default=0.05)
60 call get_param(param_file, mdl, "NL_CONTINENT_AMP", nl_top_amp, &
61 "Scale factor for topography - 0.0 for no continents.", units="nondim", default=1.0)
62
63 pi = 4.0*atan(1.0)
64
65! Calculate the depth of the bottom.
66 do j=js,je ; do i=is,ie
67 x = (g%geoLonT(i,j)-g%west_lon) / g%len_lon
68 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) = 1.0 - 1.1 * spike(y-1,0.12) - 1.1 * spike(y,0.12) - & !< The great northern wall and Antarctica
71 nl_top_amp*( &
72 (1.2 * spike(x,0.2) + 1.2 * spike(x-1.0,0.2)) * spike(min(0.0,y-0.3),0.2) & !< South America
73 + 1.2 * spike(x-0.5,0.2) * spike(min(0.0,y-0.55),0.2) & !< Africa
74 + 1.2 * (spike(x,0.12) + spike(x-1,0.12)) * spike(max(0.0,y-0.06),0.12) & !< Antarctic Peninsula
75 + 0.1 * (cosbell(x,0.1) + cosbell(x-1,0.1)) & !< Drake Passage ridge
76 + 0.5 * cosbell(x-0.16,0.05) * (cosbell(y-0.18,0.13)**0.4) & !< Scotia Arc East
77 + 0.4 * (cosbell(x-0.09,0.08)**0.4) * cosbell(y-0.26,0.05) & !< Scotia Arc North
78 + 0.4 * (cosbell(x-0.08,0.08)**0.4) * cosbell(y-0.1,0.05)) & !< Scotia Arc South
79 - nl_roughness_amp * cos(14*pi*x) * sin(14*pi*y) & !< roughness
80 - nl_roughness_amp * cos(20*pi*x) * cos(20*pi*y) !< roughness
81 if (d(i,j) < 0.0) d(i,j) = 0.0
82 d(i,j) = d(i,j) * max_depth
83 enddo ; enddo
84
86
87!> Returns the value of a cosine-bell function evaluated at x/L [nondim]
88real function cosbell(x, L)
89 real , intent(in) :: x !< Position in arbitrary units [A]
90 real , intent(in) :: l !< Width in arbitrary units [A]
91 real :: pi !< 3.1415926... calculated as 4*atan(1) [nondim]
92
93 pi = 4.0*atan(1.0)
94 cosbell = 0.5 * (1 + cos(pi*min(abs(x/l),1.0)))
95end function cosbell
96
97!> Returns the value of a sin-spike function evaluated at x/L [nondim]
98real function spike(x, L)
99
100 real , intent(in) :: x !< Position in arbitrary units [A]
101 real , intent(in) :: l !< Width in arbitrary units [A]
102 real :: pi !< 3.1415926... calculated as 4*atan(1) [nondim]
103
104 pi = 4.0*atan(1.0)
105 spike = (1 - sin(pi*min(abs(x/l),0.5)))
106end function spike
107
108!> Returns the value of a triangular function centered at x=x0 with value 1
109!! and linearly decreasing to 0 at x=x0+/-L, and 0 otherwise [nondim].
110!! If clip is present the top of the cone is cut off at "clip", which
111!! effectively defaults to 1.
112real function cone(x, x0, L, clip)
113 real, intent(in) :: x !< Coordinate in arbitrary units [A]
114 real, intent(in) :: x0 !< position of peak in arbitrary units [A]
115 real, intent(in) :: l !< half-width of base of cone in arbitrary units [A]
116 real, optional, intent(in) :: clip !< clipping height of cone [nondim]
117
118 cone = max( 0., 1. - abs(x - x0) / l )
119 if (present(clip)) cone = min(clip, cone)
120end function cone
121
122!> Returns an s-curve s(x) s.t. s(x0)<=0, s(x0+L)>=1 and cubic in between [nondim].
123real function scurve(x, x0, L)
124 real, intent(in) :: x !< Coordinate in arbitrary units [A]
125 real, intent(in) :: x0 !< position of peak in arbitrary units [A]
126 real, intent(in) :: l !< half-width of base of cone in arbitrary units [A]
127 real :: s ! A rescaled position [nondim]
128
129 s = max( 0., min( 1.,( x - x0 ) / l ) )
130 scurve = ( 3. - 2.*s ) * ( s * s )
131end function scurve
132
133! None of the following 7 functions appear to be used.
134
135!> Returns a "coastal" profile [nondim].
136real function cstprof(x, x0, L, lf, bf, sf, sh)
137 real, intent(in) :: x !< Coordinate in arbitrary units [A]
138 real, intent(in) :: x0 !< position of peak in arbitrary units [A]
139 real, intent(in) :: l !< width of profile in arbitrary units [A]
140 real, intent(in) :: lf !< fraction of width that is "land" [nondim]
141 real, intent(in) :: bf !< fraction of width that is "beach" [nondim]
142 real, intent(in) :: sf !< fraction of width that is "continental slope" [nondim]
143 real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
144 real :: s ! A rescaled position [nondim]
145
146 s = max( 0., min( 1.,( x - x0 ) / l ) )
147 cstprof = sh * scurve(s-lf,0.,bf) + (1.-sh) * scurve(s - (1.-sf),0.,sf)
148end function cstprof
149
150!> Distance between points x,y and a line segment (x0,y0) and (x0,y1) in arbitrary units [A].
151real function dist_line_fixed_x(x, y, x0, y0, y1)
152 real, intent(in) :: x !< X-coordinate in arbitrary units [A]
153 real, intent(in) :: y !< Y-coordinate in arbitrary units [A]
154 real, intent(in) :: x0 !< x-position of line segment in arbitrary units [A]
155 real, intent(in) :: y0 !< y-position of line segment end in arbitrary units [A]
156 real, intent(in) :: y1 !< y-position of line segment end in arbitrary units [A]
157 real :: dx, yr, dy ! Relative positions in arbitrary units [A]
158
159 dx = x - x0
160 yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1
161 dy = y - yr ! =0 within y0<y<y1, =y0-y for y<y0, =y-y1 for y>y1
162 dist_line_fixed_x = sqrt( (dx*dx) + (dy*dy) )
163end function dist_line_fixed_x
164
165!> Distance between points x,y and a line segment (x0,y0) and (x1,y0) in arbitrary units [A].
166real function dist_line_fixed_y(x, y, x0, x1, y0)
167 real, intent(in) :: x !< X-coordinate in arbitrary units [A]
168 real, intent(in) :: y !< Y-coordinate in arbitrary units [A]
169 real, intent(in) :: x0 !< x-position of line segment end in arbitrary units [A]
170 real, intent(in) :: x1 !< x-position of line segment end in arbitrary units [A]
171 real, intent(in) :: y0 !< y-position of line segment in arbitrary units [A]
172
173 dist_line_fixed_y = dist_line_fixed_x(y, x, y0, x0, x1)
174end function dist_line_fixed_y
175
176!> A "coast profile" applied in an N-S line from lon0,lat0 to lon0,lat1 [nondim].
177real function ns_coast(lon, lat, lon0, lat0, lat1, dlon, sh)
178 real, intent(in) :: lon !< Longitude [degrees_E]
179 real, intent(in) :: lat !< Latitude [degrees_N]
180 real, intent(in) :: lon0 !< Longitude of coast [degrees_E]
181 real, intent(in) :: lat0 !< Latitude of coast end [degrees_N]
182 real, intent(in) :: lat1 !< Latitude of coast end [degrees_N]
183 real, intent(in) :: dlon !< "Radius" of coast profile [degrees]
184 real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
185 real :: r ! A relative position [nondim]
186
187 r = dist_line_fixed_x( lon, lat, lon0, lat0, lat1 )
188 ns_coast = cstprof(r, 0., dlon, 0.125, 0.125, 0.5, sh)
189end function ns_coast
190
191!> A "coast profile" applied in an E-W line from lon0,lat0 to lon1,lat0 [nondim].
192real function ew_coast(lon, lat, lon0, lon1, lat0, dlat, sh)
193 real, intent(in) :: lon !< Longitude [degrees_E]
194 real, intent(in) :: lat !< Latitude [degrees_N]
195 real, intent(in) :: lon0 !< Longitude of coast end [degrees_E]
196 real, intent(in) :: lon1 !< Longitude of coast end [degrees_E]
197 real, intent(in) :: lat0 !< Latitude of coast [degrees_N]
198 real, intent(in) :: dlat !< "Radius" of coast profile [degrees]
199 real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
200 real :: r ! A relative position [nondim]
201
202 r = dist_line_fixed_y( lon, lat, lon0, lon1, lat0 )
203 ew_coast = cstprof(r, 0., dlat, 0.125, 0.125, 0.5, sh)
204end function ew_coast
205
206!> A NS ridge [nondim]
207real function ns_ridge(lon, lat, lon0, lat0, lat1, dlon, rh)
208 real, intent(in) :: lon !< Longitude [degrees_E]
209 real, intent(in) :: lat !< Latitude [degrees_N]
210 real, intent(in) :: lon0 !< Longitude of ridge center [degrees_E]
211 real, intent(in) :: lat0 !< Latitude of ridge end [degrees_N]
212 real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N]
213 real, intent(in) :: dlon !< "Radius" of ridge profile [degrees]
214 real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim]
215 real :: r ! A distance from a point [degrees]
216
217 r = dist_line_fixed_x( lon, lat, lon0, lat0, lat1 )
218 ns_ridge = 1. - rh * cone(r, 0., dlon)
219end function ns_ridge
220
221
222!> A circular ridge [nondim]
223real function circ_ridge(lon, lat, lon0, lat0, ring_radius, ring_thickness, ridge_height)
224 real, intent(in) :: lon !< Longitude [degrees_E]
225 real, intent(in) :: lat !< Latitude [degrees_N]
226 real, intent(in) :: lon0 !< Longitude of center of ring [degrees_E]
227 real, intent(in) :: lat0 !< Latitude of center of ring [degrees_N]
228 real, intent(in) :: ring_radius !< Radius of ring [degrees]
229 real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees]
230 real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim]
231 real :: r ! A relative position [degrees]
232 real :: frac_ht ! The fractional height of the topography [nondim]
233
234 r = sqrt( ((lon - lon0)**2) + ((lat - lat0)**2) ) ! Pseudo-distance from a point
235 r = abs( r - ring_radius) ! Pseudo-distance from a circle
236 frac_ht = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height
237 circ_ridge = 1. - frac_ht ! Fractional depths (1-frac_ridge_height) .. 1
238end function circ_ridge
239
240!> This subroutine initializes layer thicknesses for the Neverworld test case,
241!! by finding the depths of interfaces in a specified latitude-dependent
242!! temperature profile with an exponentially decaying thermocline on top of a
243!! linear stratification.
244subroutine neverworld_initialize_thickness(h, depth_tot, G, GV, US, param_file, P_ref)
245 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
246 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
247 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
248 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< The thickness that is being
249 !! initialized [Z ~> m]
250 real, dimension(SZI_(G),SZJ_(G)), &
251 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
252 type(param_file_type), intent(in) :: param_file !< A structure indicating the open
253 !! file to parse for model
254 !! parameter values.
255 real, intent(in) :: p_ref !< The coordinate-density
256 !! reference pressure [R L2 T-2 ~> Pa].
257 ! Local variables
258 real :: e0(szk_(gv)+1) ! The resting interface heights, in depth units [Z ~> m],
259 ! usually negative because it is positive upward.
260 real, dimension(SZK_(GV)) :: h_profile ! Vector of initial thickness profile [Z ~> m].
261 real :: e_interface ! Current interface position [Z ~> m].
262 real :: x, y ! horizontal coordinates for computation of the initial perturbation normalized
263 ! by the domain sizes [nondim]
264 real :: r1, r2 ! radial coordinates for computation of initial perturbation, normalized
265 ! by the domain sizes [nondim]
266 real :: pert_amp ! Amplitude of perturbations as a fraction of layer thicknesses [nondim]
267 real :: h_noise ! Amplitude of noise to scale h by [nondim]
268 real :: noise ! Fractional noise in the layer thicknesses [nondim]
269 type(randomnumberstream) :: rns ! Random numbers for stochastic tidal parameterization
270 character(len=40) :: mdl = "Neverworld_initialize_thickness" ! This subroutine's name.
271 integer :: i, j, k, is, ie, js, je, nz
272
273 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
274
275 call mom_mesg(" Neverworld_initialization.F90, Neverworld_initialize_thickness: setting thickness", 5)
276 call get_param(param_file, mdl, "INIT_THICKNESS_PROFILE", h_profile, &
277 "Profile of initial layer thicknesses.", units="m", scale=us%m_to_Z, &
278 fail_if_missing=.true.)
279 call get_param(param_file, mdl, "NL_THICKNESS_PERT_AMP", pert_amp, &
280 "Amplitude of finite scale perturbations as fraction of depth.", &
281 units="nondim", default=0.)
282 call get_param(param_file, mdl, "NL_THICKNESS_NOISE_AMP", h_noise, &
283 "Amplitude of noise to scale layer by.", units="nondim", default=0.)
284
285 ! e0 is the notional position of interfaces
286 e0(1) = 0. ! The surface
287 do k=1,nz
288 e0(k+1) = e0(k) - h_profile(k)
289 enddo
290
291 do j=js,je ; do i=is,ie
292 e_interface = -depth_tot(i,j)
293 do k=nz,2,-1
294 h(i,j,k) = e0(k) - e_interface ! Nominal thickness
295 x = (g%geoLonT(i,j)-g%west_lon)/g%len_lon
296 y = (g%geoLatT(i,j)-g%south_lat)/g%len_lat
297 r1 = sqrt(((x-0.7)**2) + ((y-0.2)**2))
298 r2 = sqrt(((x-0.3)**2) + ((y-0.25)**2))
299 h(i,j,k) = h(i,j,k) + pert_amp * (e0(k) - e0(nz+1)) * &
300 (spike(r1,0.15)-spike(r2,0.15)) ! Prescribed perturbation
301 if (h_noise /= 0.) then
302 rns = initializerandomnumberstream( int( 4096*(x + (y+1.)) ) )
303 call getrandomnumbers(rns, noise) ! x will be in (0,1)
304 noise = h_noise * 2. * ( noise - 0.5 ) ! range -h_noise to h_noise
305 h(i,j,k) = ( 1. + noise ) * h(i,j,k)
306 endif
307 h(i,j,k) = max( gv%Angstrom_Z, h(i,j,k) ) ! Limit to non-negative
308 e_interface = e_interface + h(i,j,k) ! Actual position of upper interface
309 enddo
310 h(i,j,1) = e0(1) - e_interface ! Nominal thickness
311 h(i,j,1) = max( gv%Angstrom_Z, h(i,j,1) ) ! Limit to non-negative
312 enddo ; enddo
313
315