basin_builder.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!> An idealized topography building system
6module basin_builder
7
9use 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
15
16implicit none ; private
17
18#include <MOM_memory.h>
19
21
22! This include declares and sets the variable "version".
23# include "version_variable.h"
24character(len=40) :: mdl = "basin_builder" !< This module's name.
25
26contains
27
28!> Constructs idealized topography from simple functions
29subroutine basin_builder_topography(D, G, param_file, max_depth)
30 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
31 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
32 intent(out) :: d !< Ocean bottom depth in the units of depth_max [A]
33 type(param_file_type), intent(in) :: param_file !< Parameter file structure
34 real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units [A]
35 ! Local variables
36 character(len=17) :: pname1, pname2 ! For construction of parameter names
37 character(len=20) :: funcs ! Basin build function
38 real, dimension(20) :: pars ! Parameters for each function [various]
39 real :: lon ! Longitude [degrees_E]
40 real :: lat ! Latitude [degrees_N]
41 integer :: i, j, n, n_funcs
42
43 call mom_mesg(" basin_builder.F90, basin_builder_topography: setting topography", 5)
44 call log_version(param_file, mdl, version, "")
45
46 do j=g%jsc,g%jec ; do i=g%isc,g%iec
47 d(i,j) = 1.0
48 enddo ; enddo
49
50 call get_param(param_file, mdl, "BBUILDER_N", n_funcs, &
51 "Number of pieces of topography to use.", fail_if_missing=.true.)
52
53 do n=1,n_funcs
54 write( pname1, "('BBUILDER_',i3.3,'_FUNC')" ) n
55 write( pname2, "('BBUILDER_',i3.3,'_PARS')" ) n
56 call get_param(param_file, mdl, pname1, funcs, &
57 "The basin builder function to apply with parameters "//&
58 trim(pname2)//". Choices are: NS_COAST, EW_COAST, "//&
59 "CIRC_CONIC_RIDGE, NS_CONIC_RIDGE, CIRC_SCURVE_RIDGE, "//&
60 "NS_SCURVE_RIDGE.", &
61 fail_if_missing=.true.)
62 pars(:) = 0.
63 if (trim(lowercase(funcs)) == 'ns_coast') then
64 call get_param(param_file, mdl, pname2, pars(1:5), &
65 "NS_COAST parameters: longitude, starting latitude, "//&
66 "ending latitude, footprint radius, shelf depth.", &
67 units="degrees_E,degrees_N,degrees_N,degrees,m", &
68 fail_if_missing=.true.)
69 pars(5) = pars(5) / max_depth
70 do j=g%jsc,g%jec ; do i=g%isc,g%iec
71 lon = g%geoLonT(i,j)
72 lat = g%geoLatT(i,j)
73 d(i,j) = min( d(i,j), ns_coast(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
74 enddo ; enddo
75 elseif (trim(lowercase(funcs)) == 'ns_conic_ridge') then
76 call get_param(param_file, mdl, pname2, pars(1:5), &
77 "NS_CONIC_RIDGE parameters: longitude, starting latitude, "//&
78 "ending latitude, footprint radius, ridge height.", &
79 units="degrees_E,degrees_N,degrees_N,degrees,m", &
80 fail_if_missing=.true.)
81 pars(5) = pars(5) / max_depth
82 do j=g%jsc,g%jec ; do i=g%isc,g%iec
83 lon = g%geoLonT(i,j)
84 lat = g%geoLatT(i,j)
85 d(i,j) = min( d(i,j), ns_conic_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
86 enddo ; enddo
87 elseif (trim(lowercase(funcs)) == 'ns_scurve_ridge') then
88 call get_param(param_file, mdl, pname2, pars(1:5), &
89 "NS_SCURVE_RIDGE parameters: longitude, starting latitude, "//&
90 "ending latitude, footprint radius, ridge height.", &
91 units="degrees_E,degrees_N,degrees_N,degrees,m", &
92 fail_if_missing=.true.)
93 pars(5) = pars(5) / max_depth
94 do j=g%jsc,g%jec ; do i=g%isc,g%iec
95 lon = g%geoLonT(i,j)
96 lat = g%geoLatT(i,j)
97 d(i,j) = min( d(i,j), ns_scurve_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
98 enddo ; enddo
99 elseif (trim(lowercase(funcs)) == 'angled_coast') then
100 call get_param(param_file, mdl, pname2, pars(1:4), &
101 "ANGLED_COAST parameters: longitude intersection with Equator, "//&
102 "latitude intersection with Prime Meridian, footprint radius, shelf depth.", &
103 units="degrees_E,degrees_N,degrees,m", &
104 fail_if_missing=.true.)
105 pars(4) = pars(4) / max_depth
106 do j=g%jsc,g%jec ; do i=g%isc,g%iec
107 lon = g%geoLonT(i,j)
108 lat = g%geoLatT(i,j)
109 d(i,j) = min( d(i,j), angled_coast(lon, lat, pars(1), pars(2), pars(3), pars(4)) )
110 enddo ; enddo
111 elseif (trim(lowercase(funcs)) == 'ew_coast') then
112 call get_param(param_file, mdl, pname2, pars(1:5), &
113 "EW_COAST parameters: latitude, starting longitude, "//&
114 "ending longitude, footprint radius, shelf depth.", &
115 units="degrees_N,degrees_E,degrees_E,degrees,m", &
116 fail_if_missing=.true.)
117 pars(5) = pars(5) / max_depth
118 do j=g%jsc,g%jec ; do i=g%isc,g%iec
119 lon = g%geoLonT(i,j)
120 lat = g%geoLatT(i,j)
121 d(i,j) = min( d(i,j), ew_coast(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
122 enddo ; enddo
123 elseif (trim(lowercase(funcs)) == 'circ_conic_ridge') then
124 call get_param(param_file, mdl, pname2, pars(1:5), &
125 "CIRC_CONIC_RIDGE parameters: center longitude, center latitude, "//&
126 "ring radius, footprint radius, ridge height.", &
127 units="degrees_E,degrees_N,degrees,degrees,m", &
128 fail_if_missing=.true.)
129 pars(5) = pars(5) / max_depth
130 do j=g%jsc,g%jec ; do i=g%isc,g%iec
131 lon = g%geoLonT(i,j)
132 lat = g%geoLatT(i,j)
133 d(i,j) = min( d(i,j), circ_conic_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
134 enddo ; enddo
135 elseif (trim(lowercase(funcs)) == 'circ_scurve_ridge') then
136 call get_param(param_file, mdl, pname2, pars(1:5), &
137 "CIRC_SCURVe_RIDGE parameters: center longitude, center latitude, "//&
138 "ring radius, footprint radius, ridge height.", &
139 units="degrees_E,degrees_N,degrees,degrees,m", &
140 fail_if_missing=.true.)
141 pars(5) = pars(5) / max_depth
142 do j=g%jsc,g%jec ; do i=g%isc,g%iec
143 lon = g%geoLonT(i,j)
144 lat = g%geoLatT(i,j)
145 d(i,j) = min( d(i,j), circ_scurve_ridge(lon, lat, pars(1), pars(2), pars(3), pars(4), pars(5)) )
146 enddo ; enddo
147 else
148 call mom_error(fatal, "basin_builder.F90, basin_builer_topography:\n"//&
149 "Unrecognized function "//trim(funcs))
150 endif
151
152 enddo ! n
153
154 do j=g%jsc,g%jec ; do i=g%isc,g%iec
155 ! Dimensionalize by scaling 1 to max_depth
156 d(i,j) = d(i,j) * max_depth
157 enddo ; enddo
158
159end subroutine basin_builder_topography
160
161!> Returns the value of a triangular function centered at x=x0 with value 1
162!! and linearly decreasing to 0 at x=x0+/-L, and 0 otherwise [nondim].
163!! If clip is present the top of the cone is cut off at "clip", which
164!! effectively defaults to 1.
165real function cone(x, x0, L, clip)
166 real, intent(in) :: x !< Coordinate in arbitrary units [A]
167 real, intent(in) :: x0 !< position of peak in arbitrary units [A]
168 real, intent(in) :: l !< half-width of base of cone in arbitrary units [A]
169 real, optional, intent(in) :: clip !< clipping height of cone [nondim]
170
171 cone = max( 0., 1. - abs(x - x0) / l )
172 if (present(clip)) cone = min(clip, cone)
173end function cone
174
175!> Returns an s-curve s(x) s.t. s(x0)<=0, s(x0+L)>=1 and cubic in between [nondim].
176real function scurve(x, x0, L)
177 real, intent(in) :: x !< Coordinate in arbitrary units [A]
178 real, intent(in) :: x0 !< position of peak in arbitrary units [A]
179 real, intent(in) :: l !< half-width of base of cone in arbitrary units [A]
180 real :: s ! A rescaled position [nondim]
181
182 s = max( 0., min( 1.,( x - x0 ) / l ) )
183 scurve = ( 3. - 2.*s ) * ( s * s )
184end function scurve
185
186!> Returns a "coastal" profile [nondim].
187real function cstprof(x, x0, L, lf, bf, sf, sh)
188 real, intent(in) :: x !< Coordinate in arbitrary units [A]
189 real, intent(in) :: x0 !< position of peak in arbitrary units [A]
190 real, intent(in) :: l !< width of profile in arbitrary units [A]
191 real, intent(in) :: lf !< fraction of width that is "land" [nondim]
192 real, intent(in) :: bf !< fraction of width that is "beach" [nondim]
193 real, intent(in) :: sf !< fraction of width that is "continental slope" [nondim]
194 real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
195 real :: s ! A rescaled position [nondim]
196
197 s = max( 0., min( 1.,( x - x0 ) / l ) )
198 cstprof = sh * scurve(s-lf,0.,bf) + (1.-sh) * scurve(s - (1.-sf),0.,sf)
199end function cstprof
200
201!> Distance between points x,y and a line segment (x0,y0) and (x0,y1) in arbitrary units [A].
202real function dist_line_fixed_x(x, y, x0, y0, y1)
203 real, intent(in) :: x !< X-coordinate in arbitrary units [A]
204 real, intent(in) :: y !< Y-coordinate in arbitrary units [A]
205 real, intent(in) :: x0 !< x-position of line segment in arbitrary units [A]
206 real, intent(in) :: y0 !< y-position of line segment end in arbitrary units [A]
207 real, intent(in) :: y1 !< y-position of line segment end in arbitrary units [A]
208 real :: dx, yr, dy ! Relative positions in arbitrary units [A]
209
210 dx = x - x0
211 yr = min( max(y0,y1), max( min(y0,y1), y ) ) ! bound y by y0,y1
212 dy = y - yr ! =0 within y0<y<y1, =y0-y for y<y0, =y-y1 for y>y1
213 dist_line_fixed_x = sqrt( (dx*dx) + (dy*dy) )
214end function dist_line_fixed_x
215
216!> Distance between points x,y and a line segment (x0,y0) and (x1,y0) in arbitrary units [A].
217real function dist_line_fixed_y(x, y, x0, x1, y0)
218 real, intent(in) :: x !< X-coordinate in arbitrary units [A]
219 real, intent(in) :: y !< Y-coordinate in arbitrary units [A]
220 real, intent(in) :: x0 !< x-position of line segment end in arbitrary units [A]
221 real, intent(in) :: x1 !< x-position of line segment end in arbitrary units [A]
222 real, intent(in) :: y0 !< y-position of line segment in arbitrary units [A]
223
224 dist_line_fixed_y = dist_line_fixed_x(y, x, y0, x0, x1)
225end function dist_line_fixed_y
226
227!> An "angled coast profile" [nondim].
228real function angled_coast(lon, lat, lon_eq, lat_mer, dr, sh)
229 real, intent(in) :: lon !< Longitude [degrees_E]
230 real, intent(in) :: lat !< Latitude [degrees_N]
231 real, intent(in) :: lon_eq !< Longitude intersection with Equator [degrees_E]
232 real, intent(in) :: lat_mer !< Latitude intersection with Prime Meridian [degrees_N]
233 real, intent(in) :: dr !< "Radius" of coast profile [degrees]
234 real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
235 real :: r ! A relative position [degrees]
236 real :: i_dr ! The inverse of a distance [degrees-1]
237
238 i_dr = 1/sqrt( lat_mer*lat_mer + lon_eq*lon_eq )
239 r = i_dr * ( lat_mer*lon + lon_eq*lat - lon_eq*lat_mer)
240 angled_coast = cstprof(r, 0., dr, 0.125, 0.125, 0.5, sh)
241end function angled_coast
242
243!> A "coast profile" applied in an N-S line from lonC,lat0 to lonC,lat1 [nondim].
244real function ns_coast(lon, lat, lonC, lat0, lat1, dlon, sh)
245 real, intent(in) :: lon !< Longitude [degrees_E]
246 real, intent(in) :: lat !< Latitude [degrees_N]
247 real, intent(in) :: lonc !< Longitude of coast [degrees_E]
248 real, intent(in) :: lat0 !< Latitude of coast end [degrees_N]
249 real, intent(in) :: lat1 !< Latitude of coast end [degrees_N]
250 real, intent(in) :: dlon !< "Radius" of coast profile [degrees]
251 real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
252 real :: r ! A relative position [degrees]
253
254 r = dist_line_fixed_x( lon, lat, lonc, lat0, lat1 )
255 ns_coast = cstprof(r, 0., dlon, 0.125, 0.125, 0.5, sh)
256end function ns_coast
257
258!> A "coast profile" applied in an E-W line from lon0,latC to lon1,latC [nondim].
259real function ew_coast(lon, lat, latC, lon0, lon1, dlat, sh)
260 real, intent(in) :: lon !< Longitude [degrees_E]
261 real, intent(in) :: lat !< Latitude [degrees_N]
262 real, intent(in) :: latc !< Latitude of coast [degrees_N]
263 real, intent(in) :: lon0 !< Longitude of coast end [degrees_E]
264 real, intent(in) :: lon1 !< Longitude of coast end [degrees_E]
265 real, intent(in) :: dlat !< "Radius" of coast profile [degrees]
266 real, intent(in) :: sh !< depth of shelf as fraction of full depth [nondim]
267 real :: r ! A relative position [degrees]
268
269 r = dist_line_fixed_y( lon, lat, lon0, lon1, latc )
270 ew_coast = cstprof(r, 0., dlat, 0.125, 0.125, 0.5, sh)
271end function ew_coast
272
273!> A NS ridge with a cone profile [nondim]
274real function ns_conic_ridge(lon, lat, lonC, lat0, lat1, dlon, rh)
275 real, intent(in) :: lon !< Longitude [degrees_E]
276 real, intent(in) :: lat !< Latitude [degrees_N]
277 real, intent(in) :: lonc !< Longitude of ridge center [degrees_E]
278 real, intent(in) :: lat0 !< Latitude of ridge end [degrees_N]
279 real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N]
280 real, intent(in) :: dlon !< "Radius" of ridge profile [degrees]
281 real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim]
282 real :: r ! A relative position [degrees]
283
284 r = dist_line_fixed_x( lon, lat, lonc, lat0, lat1 )
285 ns_conic_ridge = 1. - rh * cone(r, 0., dlon)
286end function ns_conic_ridge
287
288!> A NS ridge with an scurve profile [nondim]
289real function ns_scurve_ridge(lon, lat, lonC, lat0, lat1, dlon, rh)
290 real, intent(in) :: lon !< Longitude [degrees_E]
291 real, intent(in) :: lat !< Latitude [degrees_N]
292 real, intent(in) :: lonc !< Longitude of ridge center [degrees_E]
293 real, intent(in) :: lat0 !< Latitude of ridge end [degrees_N]
294 real, intent(in) :: lat1 !< Latitude of ridge end [degrees_N]
295 real, intent(in) :: dlon !< "Radius" of ridge profile [degrees]
296 real, intent(in) :: rh !< depth of ridge as fraction of full depth [nondim]
297 real :: r ! A relative position [degrees]
298
299 r = dist_line_fixed_x( lon, lat, lonc, lat0, lat1 )
300 ns_scurve_ridge = 1. - rh * (1. - scurve(r, 0., dlon) )
301end function ns_scurve_ridge
302
303!> A circular ridge with cutoff conic profile [nondim]
304real function circ_conic_ridge(lon, lat, lon0, lat0, ring_radius, ring_thickness, ridge_height)
305 real, intent(in) :: lon !< Longitude [degrees_E]
306 real, intent(in) :: lat !< Latitude [degrees_N]
307 real, intent(in) :: lon0 !< Longitude of center of ring [degrees_E]
308 real, intent(in) :: lat0 !< Latitude of center of ring [degrees_N]
309 real, intent(in) :: ring_radius !< Radius of ring [degrees]
310 real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees]
311 real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim]
312 real :: r ! A relative position [degrees]
313 real :: frac_ht ! The fractional height of the topography [nondim]
314
315 r = sqrt( ((lon - lon0)**2) + ((lat - lat0)**2) ) ! Pseudo-distance from a point
316 r = abs( r - ring_radius) ! Pseudo-distance from a circle
317 frac_ht = cone(r, 0., ring_thickness, ridge_height) ! 0 .. frac_ridge_height
318 circ_conic_ridge = 1. - frac_ht ! nondim depths (1-frac_ridge_height) .. 1
319end function circ_conic_ridge
320
321!> A circular ridge with cutoff scurve profile [nondim]
322real function circ_scurve_ridge(lon, lat, lon0, lat0, ring_radius, ring_thickness, ridge_height)
323 real, intent(in) :: lon !< Longitude [degrees_E]
324 real, intent(in) :: lat !< Latitude [degrees_N]
325 real, intent(in) :: lon0 !< Longitude of center of ring [degrees_E]
326 real, intent(in) :: lat0 !< Latitude of center of ring [degrees_N]
327 real, intent(in) :: ring_radius !< Radius of ring [degrees]
328 real, intent(in) :: ring_thickness !< Radial thickness of ring [degrees]
329 real, intent(in) :: ridge_height !< Ridge height as fraction of full depth [nondim]
330 real :: r ! A relative position [degrees]
331 real :: s ! A function of the normalized position [nondim]
332 real :: frac_ht ! The fractional height of the topography [nondim]
333
334 r = sqrt( ((lon - lon0)**2) + ((lat - lat0)**2) ) ! Pseudo-distance from a point
335 r = abs( r - ring_radius) ! Pseudo-distance from a circle
336 s = 1. - scurve(r, 0., ring_thickness) ! 0 .. 1
337 frac_ht = s * ridge_height ! 0 .. frac_ridge_height
338 circ_scurve_ridge = 1. - frac_ht ! nondim depths (1-frac_ridge_height) .. 1
339end function circ_scurve_ridge
340
341end module basin_builder