MOM6
Neverworld_initialization.F90
1 !> Initialization for the "Neverworld" configuration
2 module neverworld_initialization
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_sponge, only : sponge_cs, set_up_sponge_field, initialize_sponge
7 use mom_dyn_horgrid, only : dyn_horgrid_type
8 use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
9 use mom_file_parser, only : get_param, log_version, param_file_type
10 use mom_get_input, only : directories
11 use mom_grid, only : ocean_grid_type
12 use mom_tracer_registry, only : tracer_registry_type
13 use mom_unit_scaling, only : unit_scale_type
14 use mom_variables, only : thermo_var_ptrs
15 use mom_verticalgrid, only : verticalgrid_type
16 use mom_eos, only : calculate_density, calculate_density_derivs, eos_type
17 
18 use random_numbers_mod, only: initializerandomnumberstream, getrandomnumbers, randomnumberstream
19 
20 implicit none ; private
21 
22 #include <MOM_memory.h>
23 
24 public neverworld_initialize_topography
25 public neverworld_initialize_thickness
26 
27 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
28 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
29 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
30 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
31 
32 contains
33 
34 !> This subroutine sets up the Neverworld test case topography.
35 subroutine neverworld_initialize_topography(D, G, param_file, max_depth)
36  type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
37  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
38  intent(out) :: d !< Ocean bottom depth in the units of depth_max
39  type(param_file_type), intent(in) :: param_file !< Parameter file structure
40  real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units
41 
42  ! Local variables
43  real :: pi ! 3.1415926... calculated as 4*atan(1)
44  real :: d0 ! A constant to make the maximum !
45  ! basin depth MAXIMUM_DEPTH. !
46  real :: x, y
47  ! This include declares and sets the variable "version".
48 # include "version_variable.h"
49  character(len=40) :: mdl = "Neverworld_initialize_topography" ! This subroutine's name.
50  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
51  real :: nl_roughness_amp, nl_top_amp
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.", 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.", 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 
85 end subroutine neverworld_initialize_topography
86 
87 !> Returns the value of a cosine-bell function evaluated at x/L
88 real function cosbell(x, L)
89  real , intent(in) :: x !< non-dimensional position
90  real , intent(in) :: l !< non-dimensional width
91  real :: pi !< 3.1415926... calculated as 4*atan(1)
92 
93  pi = 4.0*atan(1.0)
94  cosbell = 0.5 * (1 + cos(pi*min(abs(x/l),1.0)))
95 end function cosbell
96 
97 !> Returns the value of a sin-spike function evaluated at x/L
98 real function spike(x, L)
99 
100  real , intent(in) :: x !< non-dimensional position
101  real , intent(in) :: l !< non-dimensional width
102  real :: pi !< 3.1415926... calculated as 4*atan(1)
103 
104  pi = 4.0*atan(1.0)
105  spike = (1 - sin(pi*min(abs(x/l),0.5)))
106 end 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.
110 !! If clip is present the top of the cone is cut off at "clip", which
111 !! effectively defaults to 1.
112 real function cone(x, x0, L, clip)
113  real, intent(in) :: x !< non-dimensional coordinate [nondim]
114  real, intent(in) :: x0 !< position of peak [nondim]
115  real, intent(in) :: l !< half-width of base of cone [nondim]
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)
120 end function cone
121 
122 !> Returns an s-curve s(x) s.t. s(x0)<=0, s(x0+L)>=1 and cubic in between.
123 real function scurve(x, x0, L)
124  real, intent(in) :: x !< non-dimensional coordinate [nondim]
125  real, intent(in) :: x0 !< position of peak [nondim]
126  real, intent(in) :: l !< half-width of base of cone [nondim]
127  real :: s
128 
129  s = max( 0., min( 1.,( x - x0 ) / l ) )
130  scurve = ( 3. - 2.*s ) * ( s * s )
131 end function scurve
132 
133 !> Returns a "coastal" profile.
134 real function cstprof(x, x0, L, lf, bf, sf, sh)
135  real, intent(in) :: x !< non-dimensional coordinate [nondim]
136  real, intent(in) :: x0 !< position of peak [nondim]
137  real, intent(in) :: l !< width of profile [nondim]
138  real, intent(in) :: lf !< fraction of width that is "land" [nondim]
139  real, intent(in) :: bf !< fraction of width that is "beach" [nondim]
140  real, intent(in) :: sf !< fraction of width that is "continental slope" [nondim]
141  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
142  real :: s
143 
144  s = max( 0., min( 1.,( x - x0 ) / l ) )
145  cstprof = sh * scurve(s-lf,0.,bf) + (1.-sh) * scurve(s - (1.-sf),0.,sf)
146 end function cstprof
147 
148 !> Distance between points x,y and a line segment (x0,y0) and (x0,y1).
149 real function dist_line_fixed_x(x, y, x0, y0, y1)
150  real, intent(in) :: x !< non-dimensional x-coordinate [nondim]
151  real, intent(in) :: y !< non-dimensional y-coordinate [nondim]
152  real, intent(in) :: x0 !< x-position of line segment [nondim]
153  real, intent(in) :: y0 !< y-position of line segment end[nondim]
154  real, intent(in) :: y1 !< y-position of line segment end[nondim]
155  real :: dx, yr, dy
156 
157  dx = x - x0
158  yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1
159  dy = y - yr ! =0 within y0<y<y1, =y0-y for y<y0, =y-y1 for y>y1
160  dist_line_fixed_x = sqrt( dx*dx + dy*dy )
161 end function dist_line_fixed_x
162 
163 !> Distance between points x,y and a line segment (x0,y0) and (x1,y0).
164 real function dist_line_fixed_y(x, y, x0, x1, y0)
165  real, intent(in) :: x !< non-dimensional x-coordinate [nondim]
166  real, intent(in) :: y !< non-dimensional y-coordinate [nondim]
167  real, intent(in) :: x0 !< x-position of line segment end[nondim]
168  real, intent(in) :: x1 !< x-position of line segment end[nondim]
169  real, intent(in) :: y0 !< y-position of line segment [nondim]
170  real :: dx, yr, dy
171 
172  dist_line_fixed_y = dist_line_fixed_x(y, x, y0, x0, x1)
173 end function dist_line_fixed_y
174 
175 !> A "coast profile" applied in an N-S line from lon0,lat0 to lon0,lat1.
176 real function ns_coast(lon, lat, lon0, lat0, lat1, dlon, sh)
177  real, intent(in) :: lon !< Longitude [degrees_E]
178  real, intent(in) :: lat !< Latitude [degrees_N]
179  real, intent(in) :: lon0 !< Longitude of coast [degrees_E]
180  real, intent(in) :: lat0 !< Latitude of coast end [degrees_N]
181  real, intent(in) :: lat1 !< Latitude of coast end [degrees_N]
182  real, intent(in) :: dlon !< "Radius" of coast profile [degrees]
183  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
184  real :: r
185 
186  r = dist_line_fixed_x( lon, lat, lon0, lat0, lat1 )
187  ns_coast = cstprof(r, 0., dlon, 0.125, 0.125, 0.5, sh)
188 end function ns_coast
189 
190 !> A "coast profile" applied in an E-W line from lon0,lat0 to lon1,lat0.
191 real function ew_coast(lon, lat, lon0, lon1, lat0, dlat, sh)
192  real, intent(in) :: lon !< Longitude [degrees_E]
193  real, intent(in) :: lat !< Latitude [degrees_N]
194  real, intent(in) :: lon0 !< Longitude of coast end [degrees_E]
195  real, intent(in) :: lon1 !< Longitude of coast end [degrees_E]
196  real, intent(in) :: lat0 !< Latitude of coast [degrees_N]
197  real, intent(in) :: dlat !< "Radius" of coast profile [degrees]
198  real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
199  real :: r
200 
201  r = dist_line_fixed_y( lon, lat, lon0, lon1, lat0 )
202  ew_coast = cstprof(r, 0., dlat, 0.125, 0.125, 0.5, sh)
203 end function ew_coast
204 
205 !> A NS ridge
206 real function ns_ridge(lon, lat, lon0, lat0, lat1, dlon, rh)
207  real, intent(in) :: lon !< Longitude [degrees_E]
208  real, intent(in) :: lat !< Latitude [degrees_N]
209  real, intent(in) :: lon0 !< Longitude of ridge center [degrees_E]
210  real, intent(in) :: lat0 !< Latitude of ridge end [degrees_N]
211  real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N]
212  real, intent(in) :: dlon !< "Radius" of ridge profile [degrees]
213  real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim]
214  real :: r
215 
216  r = dist_line_fixed_x( lon, lat, lon0, lat0, lat1 )
217  ns_ridge = 1. - rh * cone(r, 0., dlon)
218 end function ns_ridge
219 
220 
221 !> A circular ridge
222 real function circ_ridge(lon, lat, lon0, lat0, ring_radius, ring_thickness, ridge_height)
223  real, intent(in) :: lon !< Longitude [degrees_E]
224  real, intent(in) :: lat !< Latitude [degrees_N]
225  real, intent(in) :: lon0 !< Longitude of center of ring [degrees_E]
226  real, intent(in) :: lat0 !< Latitude of center of ring [degrees_N]
227  real, intent(in) :: ring_radius !< Radius of ring [degrees]
228  real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees]
229  real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim]
230  real :: r
231 
232  r = sqrt( (lon - lon0)**2 + (lat - lat0)**2 ) ! Pseudo-distance from a point
233  r = abs( r - ring_radius) ! Pseudo-distance from a circle
234  r = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height
235  circ_ridge = 1. - r ! nondim depths (1-frac_ridge_height) .. 1
236 end function circ_ridge
237 
238 !> This subroutine initializes layer thicknesses for the Neverworld test case,
239 !! by finding the depths of interfaces in a specified latitude-dependent
240 !! temperature profile with an exponentially decaying thermocline on top of a
241 !! linear stratification.
242 subroutine neverworld_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, P_ref)
243  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
244  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
245  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
246  real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being
247  !! initialized [H ~> m or kg m-2].
248  type(param_file_type), intent(in) :: param_file !< A structure indicating the open
249  !! file to parse for model
250  !! parameter values.
251  type(eos_type), pointer :: eqn_of_state !< Equation of state structure
252  real, intent(in) :: p_ref !< The coordinate-density
253  !! reference pressure [R L2 T-2 ~> Pa].
254  ! Local variables
255  real :: e0(szk_(g)+1) ! The resting interface heights, in depth units [Z ~> m],
256  ! usually negative because it is positive upward.
257  real, dimension(SZK_(G)) :: h_profile ! Vector of initial thickness profile [Z ~> m].
258  real :: e_interface ! Current interface position [Z ~> m].
259  real :: x,y,r1,r2 ! x,y and radial coordinates for computation of initial pert.
260  real :: pert_amp ! Amplitude of perturbations measured in Angstrom_H
261  real :: h_noise ! Amplitude of noise to scale h by
262  real :: noise ! Noise
263  type(randomnumberstream) :: rns ! Random numbers for stochastic tidal parameterization
264  character(len=40) :: mdl = "Neverworld_initialize_thickness" ! This subroutine's name.
265  integer :: i, j, k, k1, is, ie, js, je, nz, itt
266 
267  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
268 
269  call mom_mesg(" Neverworld_initialization.F90, Neverworld_initialize_thickness: setting thickness", 5)
270  call get_param(param_file, mdl, "INIT_THICKNESS_PROFILE", h_profile, &
271  "Profile of initial layer thicknesses.", units="m", scale=us%m_to_Z, &
272  fail_if_missing=.true.)
273  call get_param(param_file, mdl, "NL_THICKNESS_PERT_AMP", pert_amp, &
274  "Amplitude of finite scale perturbations as fraction of depth.", &
275  units="nondim", default=0.)
276  call get_param(param_file, mdl, "NL_THICKNESS_NOISE_AMP", h_noise, &
277  "Amplitude of noise to scale layer by.", units="nondim", default=0.)
278 
279  ! e0 is the notional position of interfaces
280  e0(1) = 0. ! The surface
281  do k=1,nz
282  e0(k+1) = e0(k) - h_profile(k)
283  enddo
284 
285  do j=js,je ; do i=is,ie
286  e_interface = -g%bathyT(i,j)
287  do k=nz,2,-1
288  h(i,j,k) = gv%Z_to_H * (e0(k) - e_interface) ! Nominal thickness
289  x=(g%geoLonT(i,j)-g%west_lon)/g%len_lon
290  y=(g%geoLatT(i,j)-g%south_lat)/g%len_lat
291  r1=sqrt((x-0.7)**2+(y-0.2)**2)
292  r2=sqrt((x-0.3)**2+(y-0.25)**2)
293  h(i,j,k) = h(i,j,k) + pert_amp * (e0(k) - e0(nz+1)) * gv%Z_to_H * &
294  (spike(r1,0.15)-spike(r2,0.15)) ! Prescribed perturbation
295  if (h_noise /= 0.) then
296  rns = initializerandomnumberstream( int( 4096*(x + (y+1.)) ) )
297  call getrandomnumbers(rns, noise) ! x will be in (0,1)
298  noise = h_noise * 2. * ( noise - 0.5 ) ! range -h_noise to h_noise
299  h(i,j,k) = ( 1. + noise ) * h(i,j,k)
300  endif
301  h(i,j,k) = max( gv%Angstrom_H, h(i,j,k) ) ! Limit to non-negative
302  e_interface = e_interface + gv%H_to_Z * h(i,j,k) ! Actual position of upper interface
303  enddo
304  h(i,j,1) = gv%Z_to_H * (e0(1) - e_interface) ! Nominal thickness
305  h(i,j,1) = max( gv%Angstrom_H, h(i,j,1) ) ! Limit to non-negative
306  enddo ; enddo
307 
308 end subroutine neverworld_initialize_thickness
309 
310 end module neverworld_initialization
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:80
mom_get_input::directories
Container for paths and parameter file names.
Definition: MOM_get_input.F90:20
mom_dyn_horgrid
Contains a shareable dynamic type for describing horizontal grids and metric data and utilty routines...
Definition: MOM_dyn_horgrid.F90:3
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_tracer_registry
This module contains the tracer_registry_type and the subroutines that handle registration of tracers...
Definition: MOM_tracer_registry.F90:5
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_get_input
Reads the only Fortran name list needed to boot-strap the model.
Definition: MOM_get_input.F90:6
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:108
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_sponge
Implements sponge regions in isopycnal mode.
Definition: MOM_sponge.F90:2
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:81
mom_tracer_registry::tracer_registry_type
Type to carry basic tracer information.
Definition: MOM_tracer_registry.F90:138
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_sponge::sponge_cs
This control structure holds memory and parameters for the MOM_sponge module.
Definition: MOM_sponge.F90:41
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
neverworld_initialization
Initialization for the "Neverworld" configuration.
Definition: Neverworld_initialization.F90:2
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_dyn_horgrid::dyn_horgrid_type
Describes the horizontal ocean grid with only dynamic memory arrays.
Definition: MOM_dyn_horgrid.F90:23
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:68