MOM_shared_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!> Code that initializes fixed aspects of the model grid, such as horizontal
6!! grid metrics, topography and Coriolis, and can be shared between components.
8
9use mom_coms, only : max_across_pes, reproducing_sum
10use mom_domains, only : pass_var, pass_vector, sum_across_pes, broadcast
11use mom_domains, only : root_pe, to_all, scalar_pair, cgrid_ne, agrid
13use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
15use mom_file_parser, only : get_param, log_param, param_file_type, log_version
17use mom_io, only : mom_infra_file, mom_field
18use mom_io, only : mom_read_data, mom_read_vector, read_variable, stdout
19use mom_io, only : open_file_to_read, close_file_to_read, single_file, multiple
20use mom_io, only : slasher, vardesc, mom_write_field, var_desc
23
24implicit none ; private
25
36
37! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
38! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
39! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
40! vary with the Boussinesq approximation, the Boussinesq variant is given first.
41
42contains
43
44! -----------------------------------------------------------------------------
45!> MOM_shared_init_init just writes the code version.
46subroutine mom_shared_init_init(PF)
47 type(param_file_type), intent(in) :: pf !< A structure indicating the open file
48 !! to parse for model parameter values.
49
50 character(len=40) :: mdl = "MOM_shared_initialization" ! This module's name.
51
52! This include declares and sets the variable "version".
53#include "version_variable.h"
54 call log_version(pf, mdl, version, &
55 "Sharable code to initialize time-invariant fields, like bathymetry and Coriolis parameters.")
56
57end subroutine mom_shared_init_init
58! -----------------------------------------------------------------------------
59
60!> MOM_initialize_rotation makes the appropriate call to set up the Coriolis parameter.
61subroutine mom_initialize_rotation(f, G, PF, US)
62 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
63 real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f !< The Coriolis parameter [T-1 ~> s-1]
64 type(param_file_type), intent(in) :: pf !< Parameter file structure
65 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
66
67! This subroutine makes the appropriate call to set up the Coriolis parameter.
68! This is a separate subroutine so that it can be made public and shared with
69! the ice-sheet code or other components.
70! Set up the Coriolis parameter, f, either analytically or from file.
71 character(len=40) :: mdl = "MOM_initialize_rotation" ! This subroutine's name.
72 character(len=200) :: config
73
74 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
75 call get_param(pf, mdl, "ROTATION", config, &
76 "This specifies how the Coriolis parameter is specified: \n"//&
77 " \t 2omegasinlat - Use twice the planetary rotation rate \n"//&
78 " \t\t times the sine of latitude.\n"//&
79 " \t betaplane - Use a beta-plane or f-plane.\n"//&
80 " \t USER - call a user modified routine.", &
81 default="2omegasinlat")
82 select case (trim(config))
83 case ("2omegasinlat"); call set_rotation_planetary(f, g, pf, us)
84 case ("beta"); call set_rotation_beta_plane(f, g, pf, us)
85 case ("betaplane"); call set_rotation_beta_plane(f, g, pf, us)
86 !case ("nonrotating") ! Note from AJA: Missing case?
87 case default ; call mom_error(fatal,"MOM_initialize: "// &
88 "Unrecognized rotation setup "//trim(config))
89 end select
90 call calltree_leave(trim(mdl)//'()')
91end subroutine mom_initialize_rotation
92
93!> Calculates the components of grad f (Coriolis parameter)
94subroutine mom_calculate_grad_coriolis(dF_dx, dF_dy, G, US)
95 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
96 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
97 intent(out) :: df_dx !< x-component of grad f [T-1 L-1 ~> s-1 m-1]
98 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
99 intent(out) :: df_dy !< y-component of grad f [T-1 L-1 ~> s-1 m-1]
100 type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type
101 ! Local variables
102 character(len=40) :: mdl = "MOM_calculate_grad_Coriolis" ! This subroutine's name.
103 integer :: i,j
104 real :: f1, f2 ! Average of adjacent Coriolis parameters [T-1 ~> s-1]
105
106 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
107 if ((lbound(g%CoriolisBu,1) > g%isc-1) .or. &
108 (lbound(g%CoriolisBu,2) > g%jsc-1)) then
109 ! The gradient of the Coriolis parameter can not be calculated with this grid.
110 df_dx(:,:) = 0.0 ; df_dy(:,:) = 0.0
111 return
112 endif
113
114 do j=g%jsc, g%jec ; do i=g%isc, g%iec
115 f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
116 f2 = 0.5*( g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1) )
117 df_dx(i,j) = g%IdxT(i,j) * ( f1 - f2 )
118 f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
119 f2 = 0.5*( g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1) )
120 df_dy(i,j) = g%IdyT(i,j) * ( f1 - f2 )
121 enddo ; enddo
122 call pass_vector(df_dx, df_dy, g%Domain, stagger=agrid)
123 call calltree_leave(trim(mdl)//'()')
124
125end subroutine mom_calculate_grad_coriolis
126
127!> Return the global maximum ocean bottom depth in the same units as the input depth.
128function diagnosemaximumdepth(D, G)
129 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
130 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
131 intent(in) :: d !< Ocean bottom depth in [m] or [Z ~> m]
132 real :: diagnosemaximumdepth !< The global maximum ocean bottom depth in [m] or [Z ~> m]
133 ! Local variables
134 integer :: i,j
135 diagnosemaximumdepth = d(g%isc,g%jsc)
136 do j=g%jsc, g%jec ; do i=g%isc, g%iec
138 enddo ; enddo
139 call max_across_pes(diagnosemaximumdepth)
140end function diagnosemaximumdepth
141
142!> Read time mean ocean sea level from a file
143subroutine set_meansl_from_file(meanSL, G, param_file, US)
144 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
145 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
146 intent(out) :: meansl !< Mean sea level referenced to a zero
147 !! reference height at tracer points [Z ~> m].
148 type(param_file_type), intent(in) :: param_file !< Parameter file structure
149 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
150 ! Local variables
151 character(len=200) :: filename, file, inputdir ! Strings for file/path
152 character(len=200) :: varname ! Variable name in file
153 character(len=40) :: mdl = "set_meanSL_from_file" ! This subroutine's name.
154 integer :: i, j
155
156 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
157
158 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
159 inputdir = slasher(inputdir)
160 call get_param(param_file, mdl, "MEAN_SEA_LEVEL_FILE", file, &
161 "The file from which the mean sea level is read.", &
162 default="mean_sea_level.nc")
163 call get_param(param_file, mdl, "MEAN_SEA_LEVEL_VARNAME", varname, &
164 "The name of the mean sea level variable in MEAN_SEA_LEVEL_FILE.", &
165 default="meanSL")
166 filename = trim(inputdir)//trim(file)
167 call log_param(param_file, mdl, "INPUTDIR/TOPO_FILE", filename)
168
169 if (.not.file_exists(filename, g%Domain)) &
170 call mom_error(fatal, " "//mdl//": Unable to open "//trim(filename))
171
172 call mom_read_data(filename, trim(varname), meansl, g%Domain, scale=us%m_to_Z)
173 call pass_var(meansl, g%Domain)
174
175 call calltree_leave(trim(mdl)//'()')
176end subroutine set_meansl_from_file
177
178!> Read gridded depths from file
179subroutine initialize_topography_from_file(D, G, param_file, US)
180 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
181 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
182 intent(out) :: d !< Ocean bottom depth [Z ~> m]
183 type(param_file_type), intent(in) :: param_file !< Parameter file structure
184 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
185 ! Local variables
186 character(len=200) :: filename, topo_file, inputdir ! Strings for file/path
187 character(len=200) :: topo_varname ! Variable name in file
188 character(len=40) :: mdl = "initialize_topography_from_file" ! This subroutine's name.
189
190 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
191
192 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
193 inputdir = slasher(inputdir)
194 call get_param(param_file, mdl, "TOPO_FILE", topo_file, &
195 "The file from which the bathymetry is read.", &
196 default="topog.nc")
197 call get_param(param_file, mdl, "TOPO_VARNAME", topo_varname, &
198 "The name of the bathymetry variable in TOPO_FILE.", &
199 default="depth")
200
201 filename = trim(inputdir)//trim(topo_file)
202 call log_param(param_file, mdl, "INPUTDIR/TOPO_FILE", filename)
203
204 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
205 " initialize_topography_from_file: Unable to open "//trim(filename))
206
207 d(:,:) = -9.0e30*us%m_to_Z ! Initializing to a very large negative depth (tall mountains) everywhere
208 ! before reading from a file should do nothing. However, in the instance of
209 ! masked-out PEs, halo regions are not updated when a processor does not
210 ! exist. We need to ensure the depth in masked-out PEs appears to be that
211 ! of land so this line does that in the halo regions. For non-masked PEs
212 ! the halo region is filled properly with a later pass_var().
213 call mom_read_data(filename, trim(topo_varname), d, g%Domain, scale=us%m_to_Z)
214
215 call apply_topography_edits_from_file(d, g, param_file, us)
216
217 call calltree_leave(trim(mdl)//'()')
219
220!> Applies a list of topography overrides read from a netcdf file
221subroutine apply_topography_edits_from_file(D, G, param_file, US)
222 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
223 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
224 intent(inout) :: d !< Ocean bottom depth [m] or [Z ~> m] if
225 !! US is present
226 type(param_file_type), intent(in) :: param_file !< Parameter file structure
227 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
228
229 ! Local variables
230 real, dimension(:), allocatable :: new_depth ! The new values of the depths [Z ~> m]
231 integer, dimension(:), allocatable :: ig, jg ! The global indicies of the points to modify
232 character(len=200) :: topo_edits_file, inputdir ! Strings for file/path
233 character(len=40) :: mdl = "apply_topography_edits_from_file" ! This subroutine's name.
234 integer :: i, j, n, ncid, n_edits, i_file, j_file, ndims, sizes(8)
235 logical :: topo_edits_change_mask
236 real :: min_depth ! The shallowest value of wet points [Z ~> m]
237 real :: mask_depth ! The depth defining the land-sea boundary [Z ~> m]
238
239 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
240
241 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
242 inputdir = slasher(inputdir)
243 call get_param(param_file, mdl, "TOPO_EDITS_FILE", topo_edits_file, &
244 "The file from which to read a list of i,j,z topography overrides.", &
245 default="")
246 call get_param(param_file, mdl, "ALLOW_LANDMASK_CHANGES", topo_edits_change_mask, &
247 "If true, allow topography overrides to change land mask.", &
248 default=.false.)
249 call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
250 "If MASKING_DEPTH is unspecified, then anything shallower than "//&
251 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
252 "If MASKING_DEPTH is specified, then all depths shallower than "//&
253 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
254 units="m", default=0.0, scale=us%m_to_Z)
255 call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, &
256 "The depth below which to mask points as land points, for which all "//&
257 "fluxes are zeroed out. MASKING_DEPTH is ignored if it has the special "//&
258 "default value.", &
259 units="m", default=-9999.0, scale=us%m_to_Z)
260 if (mask_depth == -9999.*us%m_to_Z) mask_depth = min_depth
261
262 if (len_trim(topo_edits_file)==0) return
263
264 topo_edits_file = trim(inputdir)//trim(topo_edits_file)
265 if (is_root_pe()) then
266 if (.not.file_exists(topo_edits_file, g%Domain)) &
267 call mom_error(fatal, trim(mdl)//': Unable to find file '//trim(topo_edits_file))
268 call open_file_to_read(topo_edits_file, ncid)
269 else
270 ncid = -1
271 endif
272
273 ! Read and check the values of ni and nj in the file for consistency with this configuration.
274 call read_variable(topo_edits_file, 'ni', i_file, ncid_in=ncid)
275 call read_variable(topo_edits_file, 'nj', j_file, ncid_in=ncid)
276 if (i_file /= g%ieg) call mom_error(fatal, trim(mdl)//': Incompatible i-dimension of grid in '//&
277 trim(topo_edits_file))
278 if (j_file /= g%jeg) call mom_error(fatal, trim(mdl)//': Incompatible j-dimension of grid in '//&
279 trim(topo_edits_file))
280
281 ! Get nEdits
282 call field_size(topo_edits_file, 'zEdit', sizes, ndims=ndims, ncid_in=ncid)
283 if (ndims /= 1) call mom_error(fatal, "The variable zEdit has an "//&
284 "unexpected number of dimensions in "//trim(topo_edits_file) )
285 n_edits = sizes(1)
286 allocate(ig(n_edits))
287 allocate(jg(n_edits))
288 allocate(new_depth(n_edits))
289
290 ! Read iEdit, jEdit and zEdit
291 call read_variable(topo_edits_file, 'iEdit', ig, ncid_in=ncid)
292 call read_variable(topo_edits_file, 'jEdit', jg, ncid_in=ncid)
293 call read_variable(topo_edits_file, 'zEdit', new_depth, ncid_in=ncid, scale=us%m_to_Z)
294 call close_file_to_read(ncid, topo_edits_file)
295
296 do n = 1, n_edits
297 i = ig(n) - g%idg_offset + 1 ! +1 for python indexing
298 j = jg(n) - g%jdg_offset + 1
299 if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec) then
300 if (new_depth(n) /= mask_depth) then
301 write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') &
302 'Ocean topography edit: ', n, ig(n), jg(n), d(i,j)*us%Z_to_m, '->', abs(us%Z_to_m*new_depth(n)), i, j
303 d(i,j) = abs(new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
304 else
305 if (topo_edits_change_mask) then
306 write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') &
307 'Ocean topography edit: ',n,ig(n),jg(n),d(i,j)*us%Z_to_m,'->',abs(us%Z_to_m*new_depth(n)),i,j
308 d(i,j) = abs(new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
309 else
310 call mom_error(fatal, ' apply_topography_edits_from_file: '//&
311 "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file))
312 endif
313 endif
314 endif
315 enddo
316
317 deallocate( ig, jg, new_depth )
318
319 call calltree_leave(trim(mdl)//'()')
321
322!> initialize the bathymetry based on one of several named idealized configurations
323subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth, US)
324 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
325 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
326 intent(out) :: d !< Ocean bottom depth [Z ~> m]
327 type(param_file_type), intent(in) :: param_file !< Parameter file structure
328 character(len=*), intent(in) :: topog_config !< The name of an idealized
329 !! topographic configuration
330 real, intent(in) :: max_depth !< Maximum depth [Z ~> m]
331 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
332
333 ! This subroutine places the bottom depth in m into D(:,:), shaped according to the named config.
334
335 ! Local variables
336 real :: min_depth ! The minimum depth [Z ~> m].
337 real :: pi ! 3.1415926... calculated as 4*atan(1) [nondim]
338 real :: d0 ! A constant to make the maximum basin depth MAXIMUM_DEPTH [Z ~> m]
339 real :: expdecay ! A decay scale of associated with the sloping boundaries [L ~> m]
340 real :: dedge ! The depth at the basin edge [Z ~> m]
341 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
342 character(len=40) :: mdl = "initialize_topography_named" ! This subroutine's name.
343 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
344 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
345
346 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
347 call mom_mesg(" MOM_shared_initialization.F90, initialize_topography_named: "//&
348 "TOPO_CONFIG = "//trim(topog_config), 5)
349
350 call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
351 "The minimum depth of the ocean.", units="m", default=0.0, scale=us%m_to_Z)
352 if (max_depth<=0.) call mom_error(fatal,"initialize_topography_named: "// &
353 "MAXIMUM_DEPTH has a non-sensical value! Was it set?")
354
355 if (trim(topog_config) /= "flat") then
356 call get_param(param_file, mdl, "EDGE_DEPTH", dedge, &
357 "The depth at the edge of one of the named topographies.", &
358 units="m", default=100.0, scale=us%m_to_Z)
359 call get_param(param_file, mdl, "TOPOG_SLOPE_SCALE", expdecay, &
360 "The exponential decay scale used in defining some of "//&
361 "the named topographies.", units="m", default=400000.0, scale=us%m_to_L)
362 endif
363
364
365 pi = 4.0*atan(1.0)
366
367 if (trim(topog_config) == "flat") then
368 do j=js,je ; do i=is,ie ; d(i,j) = max_depth ; enddo ; enddo
369 elseif (trim(topog_config) == "spoon") then
370 d0 = (max_depth - dedge) / &
371 ((1.0 - exp(-0.5*g%len_lat*g%Rad_Earth_L*pi/(180.0 *expdecay))) * &
372 (1.0 - exp(-0.5*g%len_lat*g%Rad_Earth_L*pi/(180.0 *expdecay))))
373 do j=js,je ; do i=is,ie
374 ! This sets a bowl shaped (sort of) bottom topography, with a !
375 ! maximum depth of max_depth. !
376 d(i,j) = dedge + d0 * &
377 (sin(pi * (g%geoLonT(i,j) - (g%west_lon)) / g%len_lon) * &
378 (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))*g%Rad_Earth_L*pi / &
379 (180.0*expdecay)) ))
380 enddo ; enddo
381 elseif (trim(topog_config) == "bowl") then
382 d0 = (max_depth - dedge) / &
383 ((1.0 - exp(-0.5*g%len_lat*g%Rad_Earth_L*pi/(180.0 *expdecay))) * &
384 (1.0 - exp(-0.5*g%len_lat*g%Rad_Earth_L*pi/(180.0 *expdecay))))
385
386 ! This sets a bowl shaped (sort of) bottom topography, with a
387 ! maximum depth of max_depth.
388 do j=js,je ; do i=is,ie
389 d(i,j) = dedge + d0 * &
390 (sin(pi * (g%geoLonT(i,j) - g%west_lon) / g%len_lon) * &
391 ((1.0 - exp(-(g%geoLatT(i,j) - g%south_lat)*g%Rad_Earth_L*pi/ &
392 (180.0*expdecay))) * &
393 (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))* &
394 g%Rad_Earth_L*pi/(180.0*expdecay)))))
395 enddo ; enddo
396 elseif (trim(topog_config) == "halfpipe") then
397 d0 = max_depth - dedge
398 do j=js,je ; do i=is,ie
399 d(i,j) = dedge + d0 * abs(sin(pi*(g%geoLatT(i,j) - g%south_lat)/g%len_lat))
400 enddo ; enddo
401 else
402 call mom_error(fatal,"initialize_topography_named: "// &
403 "Unrecognized topography name "//trim(topog_config))
404 endif
405
406 ! This is here just for safety. Hopefully it doesn't do anything.
407 do j=js,je ; do i=is,ie
408 if (d(i,j) > max_depth) d(i,j) = max_depth
409 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
410 enddo ; enddo
411
412 call calltree_leave(trim(mdl)//'()')
413end subroutine initialize_topography_named
414! -----------------------------------------------------------------------------
415
416! -----------------------------------------------------------------------------
417!> limit_topography ensures that min_depth < D(x,y) < max_depth
418subroutine limit_topography(D, G, param_file, max_depth, US)
419 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
420 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
421 intent(inout) :: d !< Ocean bottom depth [Z ~> m]
422 type(param_file_type), intent(in) :: param_file !< Parameter file structure
423 real, intent(in) :: max_depth !< Maximum depth of model [Z ~> m]
424 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
425
426 ! Local variables
427 integer :: i, j
428 character(len=40) :: mdl = "limit_topography" ! This subroutine's name.
429 real :: min_depth ! The shallowest value of wet points [Z ~> m]
430 real :: mask_depth ! The depth defining the land-sea boundary [Z ~> m]
431
432 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
433
434 call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
435 "If MASKING_DEPTH is unspecified, then anything shallower than "//&
436 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
437 "If MASKING_DEPTH is specified, then all depths shallower than "//&
438 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
439 units="m", default=0.0, scale=us%m_to_Z)
440 call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, &
441 "The depth below which to mask points as land points, for which all "//&
442 "fluxes are zeroed out. MASKING_DEPTH is ignored if it has the special "//&
443 "default value.", &
444 units="m", default=-9999.0, scale=us%m_to_Z, do_not_log=.true.)
445
446 ! Make sure that min_depth < D(x,y) < max_depth for ocean points
447 ! TBD: The following f.p. equivalence uses a special value. Originally, any negative value
448 ! indicated the branch. We should create a logical flag to indicate this branch.
449 if (mask_depth == -9999.*us%m_to_Z) then
450 if (min_depth<0.) then
451 call mom_error(fatal, trim(mdl)//": MINIMUM_DEPTH<0 does not work as expected "//&
452 "unless MASKING_DEPTH has been set appropriately. Set a meaningful "//&
453 "MASKING_DEPTH to enabled negative depths (land elevations) and to "//&
454 "enable flooding.")
455 endif
456 ! This is the old path way. The 0.5*min_depth is obscure and is retained to be
457 ! backward reproducible. If you are looking at the following line you should probably
458 ! set MASKING_DEPTH. This path way does not work for negative depths, i.e. flooding.
459 do j=g%jsd,g%jed ; do i=g%isd,g%ied
460 d(i,j) = min( max( d(i,j), 0.5*min_depth ), max_depth )
461 enddo ; enddo
462 else
463 ! This is the preferred path way.
464 ! mask_depth has a meaningful value; anything shallower than mask_depth is land.
465 ! If min_depth<mask_depth (which happens when using positive depths and not changing
466 ! MINIMUM_DEPTH) then the shallower is used to modify and determine values on land points.
467 do j=g%jsd,g%jed ; do i=g%isd,g%ied
468 if (d(i,j) > min(min_depth,mask_depth)) then
469 d(i,j) = min( max( d(i,j), min_depth ), max_depth )
470 else
471 ! This statement is required for cases with masked-out PEs over the land,
472 ! to remove the large initialized values (-9e30) from the halos.
473 d(i,j) = min(min_depth,mask_depth)
474 endif
475 enddo ; enddo
476 endif
477
478 call calltree_leave(trim(mdl)//'()')
479end subroutine limit_topography
480! -----------------------------------------------------------------------------
481
482! -----------------------------------------------------------------------------
483!> This subroutine sets up the Coriolis parameter for a sphere
484subroutine set_rotation_planetary(f, G, param_file, US)
485 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid
486 real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
487 intent(out) :: f !< Coriolis parameter (vertical component) [T-1 ~> s-1]
488 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
489 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
490
491! This subroutine sets up the Coriolis parameter for a sphere
492 character(len=30) :: mdl = "set_rotation_planetary" ! This subroutine's name.
493 integer :: i, j
494 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
495 real :: omega ! The planetary rotation rate [T-1 ~> s-1]
496
497 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
498
499 call get_param(param_file, "set_rotation_planetary", "OMEGA", omega, &
500 "The rotation rate of the earth.", &
501 units="s-1", default=7.2921e-5, scale=us%T_to_s)
502 pi = 4.0*atan(1.0)
503
504 do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
505 f(i,j) = ( 2.0 * omega ) * sin( ( pi * g%geoLatBu(i,j) ) / 180.)
506 enddo ; enddo
507
508 call calltree_leave(trim(mdl)//'()')
509end subroutine set_rotation_planetary
510! -----------------------------------------------------------------------------
511
512! -----------------------------------------------------------------------------
513!> This subroutine sets up the Coriolis parameter for a beta-plane or f-plane
514subroutine set_rotation_beta_plane(f, G, param_file, US)
515 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid
516 real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
517 intent(out) :: f !< Coriolis parameter (vertical component) [T-1 ~> s-1]
518 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
519 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
520
521! This subroutine sets up the Coriolis parameter for a beta-plane
522 integer :: i, j
523 real :: f_0 ! The reference value of the Coriolis parameter [T-1 ~> s-1]
524 real :: beta ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1]
525 real :: beta_lat_ref ! The reference latitude for the beta plane [degrees_N] or [km] or [m]
526 real :: y_scl ! A scaling factor from the units of latitude [L lat-1 ~> m lat-1]
527 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
528 character(len=40) :: mdl = "set_rotation_beta_plane" ! This subroutine's name.
529 character(len=200) :: axis_units
530 character(len=40) :: beta_lat_ref_units
531
532 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
533
534 call get_param(param_file, mdl, "F_0", f_0, &
535 "The reference value of the Coriolis parameter with the "//&
536 "betaplane option.", units="s-1", default=0.0, scale=us%T_to_s)
537 call get_param(param_file, mdl, "BETA", beta, &
538 "The northward gradient of the Coriolis parameter with "//&
539 "the betaplane option.", units="m-1 s-1", default=0.0, scale=us%T_to_s*us%L_to_m)
540 call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees")
541
542 pi = 4.0*atan(1.0)
543 y_scl = g%grid_unit_to_L
544 if (g%grid_unit_to_L <= 0.0) y_scl = pi * g%Rad_Earth_L / 180.
545
546 select case (axis_units(1:1))
547 case ("d")
548 beta_lat_ref_units = "degrees"
549 case ("k")
550 beta_lat_ref_units = "kilometers"
551 case ("m")
552 beta_lat_ref_units = "meters"
553 case default ; call mom_error(fatal, &
554 " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units))
555 end select
556
557 call get_param(param_file, mdl, "BETA_LAT_REF", beta_lat_ref, &
558 "The reference latitude (origin) of the beta-plane", &
559 units=trim(beta_lat_ref_units), default=0.0)
560
561 do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
562 f(i,j) = f_0 + beta * ( (g%geoLatBu(i,j) - beta_lat_ref) * y_scl )
563 enddo ; enddo
564
565 call calltree_leave(trim(mdl)//'()')
566end subroutine set_rotation_beta_plane
567
568!> initialize_grid_rotation_angle initializes the arrays with the sine and
569!! cosine of the angle between logical north on the grid and true north.
570subroutine initialize_grid_rotation_angle(G, PF)
571 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
572 type(param_file_type), intent(in) :: pf !< A structure indicating the open file
573 !! to parse for model parameter values.
574
575 real :: angle ! The clockwise angle of the grid relative to true north [degrees]
576 real :: lon_scale ! The trigonometric scaling factor converting changes in longitude
577 ! to equivalent distances in latitudes [nondim]
578 real :: len_lon ! The periodic range of longitudes, usually 360 degrees [degrees_E].
579 real :: pi_720deg ! One quarter the conversion factor from degrees to radians [radian degree-1]
580 real :: lonb(2,2) ! The longitude of a point, shifted to have about the same value [degrees_E].
581 character(len=40) :: mdl = "initialize_grid_rotation_angle" ! This subroutine's name.
582 logical :: use_bugs
583 integer :: i, j, m, n
584
585 call get_param(pf, mdl, "GRID_ROTATION_ANGLE_BUGS", use_bugs, &
586 "If true, use an older algorithm to calculate the sine and "//&
587 "cosines needed rotate between grid-oriented directions and "//&
588 "true north and east. Differences arise at the tripolar fold.", &
589 default=.false.)
590
591 if (use_bugs) then
592 do j=g%jsc,g%jec ; do i=g%isc,g%iec
593 lon_scale = cos((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j-1 ) + &
594 g%geoLatBu(i-1,j) + g%geoLatBu(i,j)) * atan(1.0)/180)
595 angle = atan2((g%geoLonBu(i-1,j) + g%geoLonBu(i,j) - &
596 g%geoLonBu(i-1,j-1) - g%geoLonBu(i,j-1))*lon_scale, &
597 g%geoLatBu(i-1,j) + g%geoLatBu(i,j) - &
598 g%geoLatBu(i-1,j-1) - g%geoLatBu(i,j-1) )
599 g%sin_rot(i,j) = sin(angle) ! angle is the clockwise angle from lat/lon to ocean
600 g%cos_rot(i,j) = cos(angle) ! grid (e.g. angle of ocean "north" from true north)
601 enddo ; enddo
602
603 ! This is not right at a tripolar or cubed-sphere fold.
604 call pass_var(g%cos_rot, g%Domain)
605 call pass_var(g%sin_rot, g%Domain)
606 else
607 pi_720deg = atan(1.0) / 180.0
608 len_lon = 360.0 ; if (g%len_lon > 0.0) len_lon = g%len_lon
609 do j=g%jsc,g%jec ; do i=g%isc,g%iec
610 do n=1,2 ; do m=1,2
611 lonb(m,n) = modulo_around_point(g%geoLonBu(i+m-2,j+n-2), g%geoLonT(i,j), len_lon)
612 enddo ; enddo
613 lon_scale = cos(pi_720deg*((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j)) + &
614 (g%geoLatBu(i,j-1) + g%geoLatBu(i-1,j)) ) )
615 angle = atan2(lon_scale*((lonb(1,2) - lonb(2,1)) + (lonb(2,2) - lonb(1,1))), &
616 (g%geoLatBu(i-1,j) - g%geoLatBu(i,j-1)) + &
617 (g%geoLatBu(i,j) - g%geoLatBu(i-1,j-1)) )
618 g%sin_rot(i,j) = sin(angle) ! angle is the clockwise angle from lat/lon to ocean
619 g%cos_rot(i,j) = cos(angle) ! grid (e.g. angle of ocean "north" from true north)
620 enddo ; enddo
621
622 call pass_vector(g%cos_rot, g%sin_rot, g%Domain, stagger=agrid)
623 endif
624
626
627! -----------------------------------------------------------------------------
628!> Return the modulo value of x in an interval [xc-(Lx/2) xc+(Lx/2)]
629!! If Lx<=0, then it returns x without applying modulo arithmetic.
630function modulo_around_point(x, xc, Lx) result(x_mod)
631 real, intent(in) :: x !< Value to which to apply modulo arithmetic [A]
632 real, intent(in) :: xc !< Center of modulo range [A]
633 real, intent(in) :: lx !< Modulo range width [A]
634 real :: x_mod !< x shifted by an integer multiple of Lx to be close to xc [A].
635
636 if (lx > 0.0) then
637 x_mod = modulo(x - (xc - 0.5*lx), lx) + (xc - 0.5*lx)
638 else
639 x_mod = x
640 endif
641end function modulo_around_point
642
643! -----------------------------------------------------------------------------
644!> This subroutine sets the open face lengths at selected points to restrict
645!! passages to their observed widths based on a named set of sizes.
646subroutine reset_face_lengths_named(G, param_file, name, US)
647 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
648 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
649 character(len=*), intent(in) :: name !< The name for the set of face lengths. Only "global_1deg"
650 !! is currently implemented.
651 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
652
653 ! Local variables
654 character(len=256) :: mesg ! Message for error messages.
655 real :: dx_2 ! Half the local zonal grid spacing [degrees_E]
656 real :: dy_2 ! Half the local meridional grid spacing [degrees_N]
657 real :: pi_180 ! Conversion factor from degrees to radians [nondim]
658 integer :: option
659 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
660 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
661 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
662 pi_180 = (4.0*atan(1.0))/180.0
663
664 dx_2 = -1.0 ; dy_2 = -1.0
665 option = -1
666
667 select case ( trim(name) )
668 case ("global_1deg") ; option = 1 ; dx_2 = 0.5*1.0
669 case default ; call mom_error(fatal, "reset_face_lengths_named: "//&
670 "Unrecognized channel configuration name "//trim(name))
671 end select
672
673 if (option==1) then ! 1-degree settings.
674 do j=jsd,jed ; do i=isdb,iedb ! Change any u-face lengths within this loop.
675 dy_2 = dx_2 * g%dyCu(i,j)*g%IdxCu(i,j) * cos(pi_180 * g%geoLatCu(i,j))
676
677 if ((abs(g%geoLatCu(i,j)-35.5) < dy_2) .and. (g%geoLonCu(i,j) < -4.5) .and. &
678 (g%geoLonCu(i,j) > -6.5)) &
679 g%dy_Cu(i,j) = g%mask2dCu(i,j)*12000.0*us%m_to_L ! Gibraltar
680
681 if ((abs(g%geoLatCu(i,j)-12.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-43.0) < dx_2)) &
682 g%dy_Cu(i,j) = g%mask2dCu(i,j)*10000.0*us%m_to_L ! Red Sea
683
684 if ((abs(g%geoLatCu(i,j)-40.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-26.0) < dx_2)) &
685 g%dy_Cu(i,j) = g%mask2dCu(i,j)*5000.0*us%m_to_L ! Dardanelles
686
687 if ((abs(g%geoLatCu(i,j)-41.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+220.0) < dx_2)) &
688 g%dy_Cu(i,j) = g%mask2dCu(i,j)*35000.0*us%m_to_L ! Tsugaru strait at 140.0e
689
690 if ((abs(g%geoLatCu(i,j)-45.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+217.5) < 0.9)) &
691 g%dy_Cu(i,j) = g%mask2dCu(i,j)*15000.0*us%m_to_L ! Betw Hokkaido and Sakhalin at 217&218 = 142e
692
693 ! Greater care needs to be taken in the tripolar region.
694 if ((abs(g%geoLatCu(i,j)-80.84) < 0.2) .and. (abs(g%geoLonCu(i,j)+64.9) < 0.8)) &
695 g%dy_Cu(i,j) = g%mask2dCu(i,j)*38000.0*us%m_to_L ! Smith Sound in Canadian Arch - tripolar region
696
697 enddo ; enddo
698
699 do j=jsdb,jedb ; do i=isd,ied ! Change any v-face lengths within this loop.
700 dy_2 = dx_2 * g%dyCv(i,j)*g%IdxCv(i,j) * cos(pi_180 * g%geoLatCv(i,j))
701 if ((abs(g%geoLatCv(i,j)-41.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-28.5) < dx_2)) &
702 g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0*us%m_to_L ! Bosporus - should be 1000.0 m wide.
703
704 if ((abs(g%geoLatCv(i,j)-13.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-42.5) < dx_2)) &
705 g%dx_Cv(i,j) = g%mask2dCv(i,j)*10000.0*us%m_to_L ! Red Sea
706
707 if ((abs(g%geoLatCv(i,j)+2.8) < 0.8) .and. (abs(g%geoLonCv(i,j)+241.5) < dx_2)) &
708 g%dx_Cv(i,j) = g%mask2dCv(i,j)*40000.0*us%m_to_L ! Makassar Straits at 241.5 W = 118.5 E
709
710 if ((abs(g%geoLatCv(i,j)-0.56) < 0.5) .and. (abs(g%geoLonCv(i,j)+240.5) < dx_2)) &
711 g%dx_Cv(i,j) = g%mask2dCv(i,j)*80000.0*us%m_to_L ! entry to Makassar Straits at 240.5 W = 119.5 E
712
713 if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+230.5) < dx_2)) &
714 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*us%m_to_L ! Channel betw N Guinea and Halmahara 230.5 W = 129.5 E
715
716 if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+229.5) < dx_2)) &
717 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*us%m_to_L ! Channel betw N Guinea and Halmahara 229.5 W = 130.5 E
718
719 if ((abs(g%geoLatCv(i,j)-0.0) < 0.25) .and. (abs(g%geoLonCv(i,j)+228.5) < dx_2)) &
720 g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*us%m_to_L ! Channel betw N Guinea and Halmahara 228.5 W = 131.5 E
721
722 if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+244.5) < dx_2)) &
723 g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0*us%m_to_L ! Lombok Straits at 244.5 W = 115.5 E
724
725 if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+235.5) < dx_2)) &
726 g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0*us%m_to_L ! Timor Straits at 235.5 W = 124.5 E
727
728 if ((abs(g%geoLatCv(i,j)-52.5) < dy_2) .and. (abs(g%geoLonCv(i,j)+218.5) < dx_2)) &
729 g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0*us%m_to_L ! Russia and Sakhalin Straits at 218.5 W = 141.5 E
730
731 ! Greater care needs to be taken in the tripolar region.
732 if ((abs(g%geoLatCv(i,j)-76.8) < 0.06) .and. (abs(g%geoLonCv(i,j)+88.7) < dx_2)) &
733 g%dx_Cv(i,j) = g%mask2dCv(i,j)*8400.0*us%m_to_L ! Jones Sound in Canadian Arch - tripolar region
734
735 enddo ; enddo
736 endif
737
738 ! These checks apply regardless of the chosen option.
739
740 do j=jsd,jed ; do i=isdb,iedb
741 if (g%dy_Cu(i,j) > g%dyCu(i,j)) then
742 write(mesg,'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
743 &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
744 us%L_to_m*g%dy_Cu(i,j), us%L_to_m*g%dyCu(i,j), us%L_to_m*(g%dy_Cu(i,j)-g%dyCu(i,j)), &
745 g%geoLonCu(i,j), g%geoLatCu(i,j)
746 call mom_error(fatal,"reset_face_lengths_named "//mesg)
747 endif
748 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
749 g%IareaCu(i,j) = 0.0
750 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
751 enddo ; enddo
752
753 do j=jsdb,jedb ; do i=isd,ied
754 if (g%dx_Cv(i,j) > g%dxCv(i,j)) then
755 write(mesg,'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
756 &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
757 us%L_to_m*g%dx_Cv(i,j), us%L_to_m*g%dxCv(i,j), us%L_to_m*(g%dx_Cv(i,j)-g%dxCv(i,j)), &
758 g%geoLonCv(i,j), g%geoLatCv(i,j)
759
760 call mom_error(fatal,"reset_face_lengths_named "//mesg)
761 endif
762 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
763 g%IareaCv(i,j) = 0.0
764 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
765 enddo ; enddo
766
767end subroutine reset_face_lengths_named
768! -----------------------------------------------------------------------------
769
770! -----------------------------------------------------------------------------
771!> This subroutine sets the open face lengths at selected points to restrict
772!! passages to their observed widths from a arrays read from a file.
773subroutine reset_face_lengths_file(G, param_file, US)
774 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
775 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
776 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
777
778 ! Local variables
779 character(len=40) :: mdl = "reset_face_lengths_file" ! This subroutine's name.
780 character(len=256) :: mesg ! Message for error messages.
781 character(len=200) :: filename, chan_file, inputdir ! Strings for file/path
782 character(len=64) :: dxcv_open_var, dycu_open_var ! Open face length names in files
783 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
784
785 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
786 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
787 ! These checks apply regardless of the chosen option.
788
789 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
790
791 call get_param(param_file, mdl, "CHANNEL_WIDTH_FILE", chan_file, &
792 "The file from which the list of narrowed channels is read.", &
793 default="ocean_geometry.nc")
794 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
795 inputdir = slasher(inputdir)
796 filename = trim(inputdir)//trim(chan_file)
797 call log_param(param_file, mdl, "INPUTDIR/CHANNEL_WIDTH_FILE", filename)
798
799 if (is_root_pe()) then ; if (.not.file_exists(filename)) &
800 call mom_error(fatal," reset_face_lengths_file: Unable to open "//&
801 trim(filename))
802 endif
803
804 call get_param(param_file, mdl, "OPEN_DY_CU_VAR", dycu_open_var, &
805 "The u-face open face length variable in CHANNEL_WIDTH_FILE.", &
806 default="dyCuo")
807 call get_param(param_file, mdl, "OPEN_DX_CV_VAR", dxcv_open_var, &
808 "The v-face open face length variable in CHANNEL_WIDTH_FILE.", &
809 default="dxCvo")
810
811 call mom_read_vector(filename, dycu_open_var, dxcv_open_var, g%dy_Cu, g%dx_Cv, g%Domain, scale=us%m_to_L)
812 call pass_vector(g%dy_Cu, g%dx_Cv, g%Domain, to_all+scalar_pair, cgrid_ne)
813
814 do j=jsd,jed ; do i=isdb,iedb
815 if (g%dy_Cu(i,j) > g%dyCu(i,j)) then
816 write(mesg,'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
817 &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
818 us%L_to_m*g%dy_Cu(i,j), us%L_to_m*g%dyCu(i,j), us%L_to_m*(g%dy_Cu(i,j)-g%dyCu(i,j)), &
819 g%geoLonCu(i,j), g%geoLatCu(i,j)
820 call mom_error(fatal,"reset_face_lengths_file "//mesg)
821 endif
822 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
823 g%IareaCu(i,j) = 0.0
824 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
825 enddo ; enddo
826
827 do j=jsdb,jedb ; do i=isd,ied
828 if (g%dx_Cv(i,j) > g%dxCv(i,j)) then
829 write(mesg,'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
830 &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
831 us%L_to_m*g%dx_Cv(i,j), us%L_to_m*g%dxCv(i,j), us%L_to_m*(g%dx_Cv(i,j)-g%dxCv(i,j)), &
832 g%geoLonCv(i,j), g%geoLatCv(i,j)
833
834 call mom_error(fatal,"reset_face_lengths_file "//mesg)
835 endif
836 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
837 g%IareaCv(i,j) = 0.0
838 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
839 enddo ; enddo
840
841 call calltree_leave(trim(mdl)//'()')
842end subroutine reset_face_lengths_file
843! -----------------------------------------------------------------------------
844
845! -----------------------------------------------------------------------------
846!> This subroutine sets the open face lengths at selected points to restrict
847!! passages to their observed widths from a list read from a file.
848subroutine reset_face_lengths_list(G, param_file, US)
849 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
850 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
851 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
852
853 ! Local variables
854 character(len=120), pointer, dimension(:) :: lines => null()
855 character(len=120) :: line
856 character(len=200) :: filename, chan_file, inputdir ! Strings for file/path
857 character(len=40) :: mdl = "reset_face_lengths_list" ! This subroutine's name.
858 real, allocatable, dimension(:,:) :: &
859 u_lat, u_lon, v_lat, v_lon ! The latitude and longitude ranges of faces [degrees_N] or [degrees_E]
860 real, allocatable, dimension(:) :: &
861 u_width, v_width ! The open width of faces [L ~> m]
862 integer, allocatable, dimension(:) :: &
863 u_line_no, v_line_no, & ! The line numbers in lines of u- and v-face lines
864 u_line_used, v_line_used ! The number of times each u- and v-line is used.
865 real, allocatable, dimension(:) :: &
866 dmin_u, dmax_u, davg_u ! Porous barrier monomial fit params [Z ~> m]
867 real, allocatable, dimension(:) :: &
868 dmin_v, dmax_v, davg_v ! Porous barrier monomial fit params [Z ~> m]
869 real :: lat, lon ! The latitude and longitude of a point [degrees_N] and [degrees_E].
870 real :: len_lon ! The periodic range of longitudes, usually 360 degrees [degrees_E].
871 real :: len_lat ! The range of latitudes, usually 180 degrees [degrees_N].
872 real :: lon_p, lon_m ! The longitude of a point shifted by 360 degrees [degrees_E].
873 logical :: check_360 ! If true, check for longitudes that are shifted by
874 ! +/- 360 degrees from the specified range of values.
875 logical :: unit_in_use
876 logical :: fatal_unused_lengths
877 integer :: unused
878 integer :: ios, iounit, isu, isv, isu_por, isv_por
879 integer :: num_lines, nl_read, ln, npt, u_pt, v_pt
880 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
881
882 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
883 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
884
885 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
886
887 call get_param(param_file, mdl, "CHANNEL_LIST_FILE", chan_file, &
888 "The file from which the list of narrowed channels is read.", &
889 default="MOM_channel_list")
890 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
891 inputdir = slasher(inputdir)
892 filename = trim(inputdir)//trim(chan_file)
893 call log_param(param_file, mdl, "INPUTDIR/CHANNEL_LIST_FILE", filename)
894 call get_param(param_file, mdl, "CHANNEL_LIST_360_LON_CHECK", check_360, &
895 "If true, the channel configuration list works for any "//&
896 "longitudes in the range of -360 to 360.", default=.true.)
897 call get_param(param_file, mdl, "FATAL_UNUSED_CHANNEL_WIDTHS", fatal_unused_lengths, &
898 "If true, trigger a fatal error if there are any channel widths in "//&
899 "CHANNEL_LIST_FILE that do not cause any open face widths to change.", &
900 default=.false.)
901
902 if (is_root_pe()) then
903 ! Open the input file.
904 if (.not.file_exists(filename)) call mom_error(fatal, &
905 " reset_face_lengths_list: Unable to open "//trim(filename))
906
907 ! Find an unused unit number.
908 do iounit=10,512
909 INQUIRE(iounit,opened=unit_in_use) ; if (.not.unit_in_use) exit
910 enddo
911 if (iounit >= 512) call mom_error(fatal, &
912 "reset_face_lengths_list: No unused file unit could be found.")
913
914 ! Open the parameter file.
915 open(iounit, file=trim(filename), access='SEQUENTIAL', &
916 form='FORMATTED', action='READ', position='REWIND', iostat=ios)
917 if (ios /= 0) call mom_error(fatal, &
918 "reset_face_lengths_list: Error opening "//trim(filename))
919
920 ! Count the number of u_width and v_width entries.
921 call read_face_length_list(iounit, filename, num_lines, lines)
922 else
923 num_lines = 0
924 endif
925
926 len_lon = 360.0 ; if (g%len_lon > 0.0) len_lon = g%len_lon
927 len_lat = 180.0 ; if (g%len_lat > 0.0) len_lat = g%len_lat
928 ! Broadcast the number of lines and allocate the required space.
929 call broadcast(num_lines, root_pe())
930 u_pt = 0 ; v_pt = 0
931 if (num_lines > 0) then
932 allocate(lines(num_lines))
933
934 allocate(u_lat(2,num_lines), source=-1e34)
935 allocate(u_lon(2,num_lines), source=-1e34)
936 allocate(u_width(num_lines), source=-1e34)
937 allocate(u_line_used(num_lines), source=0)
938 allocate(u_line_no(num_lines), source=0)
939
940 allocate(v_lat(2,num_lines), source=-1e34)
941 allocate(v_lon(2,num_lines), source=-1e34)
942 allocate(v_width(num_lines), source=-1e34)
943 allocate(v_line_used(num_lines), source=0)
944 allocate(v_line_no(num_lines), source=0)
945
946 allocate(dmin_u(num_lines), source=0.0)
947 allocate(dmax_u(num_lines), source=0.0)
948 allocate(davg_u(num_lines), source=0.0)
949
950 allocate(dmin_v(num_lines), source=0.0)
951 allocate(dmax_v(num_lines), source=0.0)
952 allocate(davg_v(num_lines), source=0.0)
953
954 ! Actually read the lines.
955 if (is_root_pe()) then
956 call read_face_length_list(iounit, filename, nl_read, lines)
957 if (nl_read /= num_lines) &
958 call mom_error(fatal, 'reset_face_lengths_list : Found different '// &
959 'number of valid lines on second reading of '//trim(filename))
960 close(iounit) ; iounit = -1
961 endif
962
963 ! Broadcast the lines.
964 call broadcast(lines, 120, root_pe())
965
966 ! Populate the u_width, etc., data.
967 do ln=1,num_lines
968 line = lines(ln)
969 ! Detect keywords
970 isu = index(uppercase(line), "U_WIDTH")
971 isv = index(uppercase(line), "V_WIDTH")
972 isu_por = index(uppercase(line), "U_WIDTH_POR")
973 isv_por = index(uppercase(line), "V_WIDTH_POR")
974
975 ! Store and check the relevant values.
976 if (isu > 0) then ! This line includes "U_WIDTH".
977 u_pt = u_pt + 1
978 if (isu_por <= 0) then ! This line sets "U_WIDTH"
979 read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
980 else ! This line sets "U_WIDTH_POR"
981 read(line(isu_por+12:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt), &
982 dmin_u(u_pt), dmax_u(u_pt), davg_u(u_pt)
983 endif
984 u_width(u_pt) = us%m_to_L*u_width(u_pt) ! Rescale units equivalently to scale=US%m_to_L during read.
985 dmin_u(u_pt) = us%m_to_Z*dmin_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
986 dmax_u(u_pt) = us%m_to_Z*dmax_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
987 davg_u(u_pt) = us%m_to_Z*davg_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
988 u_line_no(u_pt) = ln
989 if (is_root_pe()) then
990 if (check_360) then
991 if ((abs(u_lon(1,u_pt)) > len_lon) .or. (abs(u_lon(2,u_pt)) > len_lon)) &
992 call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
993 "u-longitude found when reading line "//trim(line)//" from file "//&
994 trim(filename))
995 if ((abs(u_lat(1,u_pt)) > len_lat) .or. (abs(u_lat(2,u_pt)) > len_lat)) &
996 call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
997 "u-latitude found when reading line "//trim(line)//" from file "//&
998 trim(filename))
999 endif
1000 if (u_lat(1,u_pt) > u_lat(2,u_pt)) &
1001 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1002 "u-face latitudes found when reading line "//trim(line)//" from file "//&
1003 trim(filename))
1004 if (u_lon(1,u_pt) > u_lon(2,u_pt)) &
1005 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1006 "u-face longitudes found when reading line "//trim(line)//" from file "//&
1007 trim(filename))
1008 if (u_width(u_pt) < 0.0) &
1009 call mom_error(warning, "reset_face_lengths_list : Negative "//&
1010 "u-width found when reading line "//trim(line)//" from file "//&
1011 trim(filename))
1012 if (dmin_u(u_pt) > dmax_u(u_pt)) &
1013 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1014 "topographical min/max found when reading line "//trim(line)//" from file "//&
1015 trim(filename))
1016 endif
1017 elseif (isv > 0) then ! This line includes "V_WIDTH".
1018 v_pt = v_pt + 1
1019 if (isv_por <= 0) then ! This line sets "V_WIDTH"
1020 read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
1021 else ! This line sets "V_WIDTH_POR"
1022 read(line(isv_por+12:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt), &
1023 dmin_v(v_pt), dmax_v(v_pt), davg_v(v_pt)
1024 endif
1025 v_width(v_pt) = us%m_to_L*v_width(v_pt) ! Rescale units equivalently to scale=US%m_to_L during read.
1026 dmin_v(v_pt) = us%m_to_Z*dmin_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1027 dmax_v(v_pt) = us%m_to_Z*dmax_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1028 davg_v(v_pt) = us%m_to_Z*davg_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1029 v_line_no(v_pt) = ln
1030 if (is_root_pe()) then
1031 if (check_360) then
1032 if ((abs(v_lon(1,v_pt)) > len_lon) .or. (abs(v_lon(2,v_pt)) > len_lon)) &
1033 call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
1034 "v-longitude found when reading line "//trim(line)//" from file "//&
1035 trim(filename))
1036 if ((abs(v_lat(1,v_pt)) > len_lat) .or. (abs(v_lat(2,v_pt)) > len_lat)) &
1037 call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
1038 "v-latitude found when reading line "//trim(line)//" from file "//&
1039 trim(filename))
1040 endif
1041 if (v_lat(1,v_pt) > v_lat(2,v_pt)) &
1042 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1043 "v-face latitudes found when reading line "//trim(line)//" from file "//&
1044 trim(filename))
1045 if (v_lon(1,v_pt) > v_lon(2,v_pt)) &
1046 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1047 "v-face longitudes found when reading line "//trim(line)//" from file "//&
1048 trim(filename))
1049 if (v_width(v_pt) < 0.0) &
1050 call mom_error(warning, "reset_face_lengths_list : Negative "//&
1051 "v-width found when reading line "//trim(line)//" from file "//&
1052 trim(filename))
1053 if (dmin_v(v_pt) > dmax_v(v_pt)) &
1054 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1055 "topographical min/max found when reading line "//trim(line)//" from file "//&
1056 trim(filename))
1057 endif
1058 endif
1059 enddo
1060
1061 endif
1062
1063 do j=jsd,jed ; do i=isdb,iedb
1064 lat = g%geoLatCu(i,j) ; lon = g%geoLonCu(i,j)
1065 if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1066 else ; lon_p = lon ; lon_m = lon ; endif
1067
1068 do npt=1,u_pt
1069 if (((lat >= u_lat(1,npt)) .and. (lat <= u_lat(2,npt))) .and. &
1070 (((lon >= u_lon(1,npt)) .and. (lon <= u_lon(2,npt))) .or. &
1071 ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. &
1072 ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then
1073
1074 g%dy_Cu(i,j) = g%mask2dCu(i,j) * min(g%dyCu(i,j), max(u_width(npt), 0.0))
1075 g%porous_DminU(i,j) = dmin_u(npt)
1076 g%porous_DmaxU(i,j) = dmax_u(npt)
1077 g%porous_DavgU(i,j) = davg_u(npt)
1078
1079 if (j>=g%jsc .and. j<=g%jec .and. i>=g%isc .and. i<=g%iec) then ! Limit messages/checking to compute domain
1080 if ( g%mask2dCu(i,j) == 0.0 ) then
1081 write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",&
1082 u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") so grid metric is unmodified."
1083 else
1084 u_line_used(npt) = u_line_used(npt) + 1
1085 write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1086 "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon," (",&
1087 u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") to ",us%L_to_m*g%dy_Cu(i,j),"m"
1088 write(stdout,'(A,3F8.2,A)') &
1089 "read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",g%porous_DminU(i,j),&
1090 g%porous_DmaxU(i,j), g%porous_DavgU(i,j),")m"
1091 endif
1092 endif
1093 endif
1094 enddo
1095
1096 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
1097 g%IareaCu(i,j) = 0.0
1098 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
1099 enddo ; enddo
1100
1101 do j=jsdb,jedb ; do i=isd,ied
1102 lat = g%geoLatCv(i,j) ; lon = g%geoLonCv(i,j)
1103 if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1104 else ; lon_p = lon ; lon_m = lon ; endif
1105
1106 do npt=1,v_pt
1107 if (((lat >= v_lat(1,npt)) .and. (lat <= v_lat(2,npt))) .and. &
1108 (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. &
1109 ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. &
1110 ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then
1111 g%dx_Cv(i,j) = g%mask2dCv(i,j) * min(g%dxCv(i,j), max(v_width(npt), 0.0))
1112 g%porous_DminV(i,j) = dmin_v(npt)
1113 g%porous_DmaxV(i,j) = dmax_v(npt)
1114 g%porous_DavgV(i,j) = davg_v(npt)
1115
1116 if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec) then ! Limit messages/checking to compute domain
1117 if ( g%mask2dCv(i,j) == 0.0 ) then
1118 write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",&
1119 v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") so grid metric is unmodified."
1120 else
1121 v_line_used(npt) = v_line_used(npt) + 1
1122 write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1123 "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon," (",&
1124 v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") to ",us%L_to_m*g%dx_Cv(i,j),"m"
1125 write(stdout,'(A,3F8.2,A)') &
1126 "read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",g%porous_DminV(i,j),&
1127 g%porous_DmaxV(i,j), g%porous_DavgV(i,j),")m"
1128 endif
1129 endif
1130 endif
1131 enddo
1132
1133 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
1134 g%IareaCv(i,j) = 0.0
1135 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
1136 enddo ; enddo
1137
1138 ! Verify that all channel widths have been used
1139 unused = 0
1140 if (u_pt > 0) call sum_across_pes(u_line_used, u_pt)
1141 if (v_pt > 0) call sum_across_pes(v_line_used, v_pt)
1142 if (is_root_pe()) then
1143 unused = 0
1144 do npt=1,u_pt ; if (u_line_used(npt) == 0) then
1145 call mom_error(warning, "reset_face_lengths_list unused u-face line: "//&
1146 trim(lines(u_line_no(npt))) )
1147 unused = unused + 1
1148 endif ; enddo
1149 do npt=1,v_pt ; if (v_line_used(npt) == 0) then
1150 call mom_error(warning, "reset_face_lengths_list unused v-face line: "//&
1151 trim(lines(v_line_no(npt))) )
1152 unused = unused + 1
1153 endif ; enddo
1154 if (fatal_unused_lengths .and. (unused > 0)) call mom_error(fatal, &
1155 "reset_face_lengths_list causing MOM6 abort due to unused face length lines.")
1156 endif
1157
1158 if (num_lines > 0) then
1159 deallocate(lines)
1160 deallocate(u_line_used, v_line_used, u_line_no, v_line_no)
1161 deallocate(u_lat) ; deallocate(u_lon) ; deallocate(u_width)
1162 deallocate(v_lat) ; deallocate(v_lon) ; deallocate(v_width)
1163 deallocate(dmin_u) ; deallocate(dmax_u) ; deallocate(davg_u)
1164 deallocate(dmin_v) ; deallocate(dmax_v) ; deallocate(davg_v)
1165 endif
1166
1167 call calltree_leave(trim(mdl)//'()')
1168end subroutine reset_face_lengths_list
1169! -----------------------------------------------------------------------------
1170
1171! -----------------------------------------------------------------------------
1172!> This subroutine reads and counts the non-blank lines in the face length list file, after removing comments.
1173subroutine read_face_length_list(iounit, filename, num_lines, lines)
1174 integer, intent(in) :: iounit !< An open I/O unit number for the file
1175 character(len=*), intent(in) :: filename !< The name of the face-length file to read
1176 integer, intent(out) :: num_lines !< The number of non-blank lines in the file
1177 character(len=120), dimension(:), pointer :: lines !< The non-blank lines, after removing comments
1178
1179 ! This subroutine reads and counts the non-blank lines in the face length
1180 ! list file, after removing comments.
1181 character(len=120) :: line, line_up
1182 logical :: found_u, found_v
1183 integer :: icom
1184 integer :: last
1185
1186 num_lines = 0
1187
1188 if (iounit <= 0) return
1189 rewind(iounit)
1190 do while(.true.)
1191 read(iounit, '(a)', end=8, err=9) line
1192 last = len_trim(line)
1193 ! Eliminate either F90 or C comments from the line.
1194 icom = index(line(:last), "!") ; if (icom > 0) last = icom-1
1195 icom = index(line(:last), "/*") ; if (icom > 0) last = icom-1
1196 if (last < 1) cycle
1197
1198 ! Detect keywords
1199 line_up = uppercase(line)
1200 found_u = (index(line_up(:last), "U_WIDTH") > 0)
1201 found_v = (index(line_up(:last), "V_WIDTH") > 0)
1202
1203 if (found_u .and. found_v) call mom_error(fatal, &
1204 "read_face_length_list : both U_WIDTH and V_WIDTH found when "//&
1205 "reading the line "//trim(line(:last))//" in file "//trim(filename))
1206 if (found_u .or. found_v) then
1207 num_lines = num_lines + 1
1208 if (associated(lines)) then
1209 lines(num_lines) = line(1:last)
1210 endif
1211 endif
1212 enddo ! while (.true.)
1213
12148 continue
1215 return
1216
12179 call mom_error(fatal, "read_face_length_list : "//&
1218 "Error while reading file "//trim(filename))
1219
1220end subroutine read_face_length_list
1221! -----------------------------------------------------------------------------
1222
1223! -----------------------------------------------------------------------------
1224!> Read from a file the maximum, minimum and average bathymetry at velocity points,
1225!! for the use of porous barrier.
1226!! Note that we assume the depth values in the sub-grid bathymetry file of the same
1227!! convention as in-cell bathymetry file, i.e. positive below the sea surface and
1228!! increasing downward; while in subroutine reset_face_lengths_list, it is implied
1229!! that read-in fields min_bathy, max_bathy and avg_bathy from the input file
1230!! CHANNEL_LIST_FILE all have negative values below the surface. Therefore, to ensure
1231!! backward compatibility, all signs of the variable are inverted here.
1232!! And porous_Dmax[UV] = shallowest point, porous_Dmin[UV] = deepest point
1233subroutine set_subgrid_topo_at_vel_from_file(G, param_file, US)
1234 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
1235 type(param_file_type), intent(in) :: param_file !< Parameter file structure
1236 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1237
1238 ! Local variables
1239 character(len=200) :: filename, topo_file, inputdir ! Strings for file/path
1240 character(len=200) :: varname_uhi, varname_ulo, varname_uav, &
1241 varname_vhi, varname_vlo, varname_vav ! Variable names in file
1242 character(len=40) :: mdl = "set_subgrid_topo_at_vel_from_file" ! This subroutine's name.
1243 integer :: i, j
1244
1245 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
1246
1247 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1248 inputdir = slasher(inputdir)
1249 call get_param(param_file, mdl, "TOPO_AT_VEL_FILE", topo_file, &
1250 "The file from which the bathymetry parameters at the velocity points are read. "//&
1251 "While the names of the parameters reflect their physical locations, i.e. HIGH is above LOW, "//&
1252 "their signs follow the model's convention, which is positive below the sea surface", &
1253 default="topog_edge.nc")
1254 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_HIGH", varname_uhi, &
1255 "The variable name of the highest bathymetry at the u-cells in TOPO_AT_VEL_FILE.", &
1256 default="depthu_hi")
1257 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_LOW", varname_ulo, &
1258 "The variable name of the lowest bathymetry at the u-cells in TOPO_AT_VEL_FILE.", &
1259 default="depthu_lo")
1260 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_AVE", varname_uav, &
1261 "The variable name of the average bathymetry at the u-cells in TOPO_AT_VEL_FILE.", &
1262 default="depthu_av")
1263 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_HIGH", varname_vhi, &
1264 "The variable name of the highest bathymetry at the v-cells in TOPO_AT_VEL_FILE.", &
1265 default="depthv_hi")
1266 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_LOW", varname_vlo, &
1267 "The variable name of the lowest bathymetry at the v-cells in TOPO_AT_VEL_FILE.", &
1268 default="depthv_lo")
1269 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_AVE", varname_vav, &
1270 "The variable name of the average bathymetry at the v-cells in TOPO_AT_VEL_FILE.", &
1271 default="depthv_av")
1272
1273 filename = trim(inputdir)//trim(topo_file)
1274 call log_param(param_file, mdl, "INPUTDIR/TOPO_AT_VEL_FILE", filename)
1275
1276 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
1277 " set_subgrid_topo_at_vel_from_file: Unable to open "//trim(filename))
1278
1279 call mom_read_vector(filename, trim(varname_uhi), trim(varname_vhi), &
1280 g%porous_DmaxU, g%porous_DmaxV, g%Domain, stagger=cgrid_ne, scale=us%m_to_Z)
1281 call mom_read_vector(filename, trim(varname_ulo), trim(varname_vlo), &
1282 g%porous_DminU, g%porous_DminV, g%Domain, stagger=cgrid_ne, scale=us%m_to_Z)
1283 call mom_read_vector(filename, trim(varname_uav), trim(varname_vav), &
1284 g%porous_DavgU, g%porous_DavgV, g%Domain, stagger=cgrid_ne, scale=us%m_to_Z)
1285
1286 ! The signs of the depth parameters need to be inverted to be backward compatible with input files
1287 ! used by subroutine reset_face_lengths_list, which assumes depth is negative below the sea surface.
1288 g%porous_DmaxU = -g%porous_DmaxU ; g%porous_DminU = -g%porous_DminU ; g%porous_DavgU = -g%porous_DavgU
1289 g%porous_DmaxV = -g%porous_DmaxV ; g%porous_DminV = -g%porous_DminV ; g%porous_DavgV = -g%porous_DavgV
1290
1291 call pass_vector(g%porous_DmaxU, g%porous_DmaxV, g%Domain, to_all+scalar_pair, cgrid_ne)
1292 call pass_vector(g%porous_DminU, g%porous_DminV, g%Domain, to_all+scalar_pair, cgrid_ne)
1293 call pass_vector(g%porous_DavgU, g%porous_DavgV, g%Domain, to_all+scalar_pair, cgrid_ne)
1294
1295 call calltree_leave(trim(mdl)//'()')
1297! -----------------------------------------------------------------------------
1298
1299! -----------------------------------------------------------------------------
1300!> Set the bathymetry at velocity points to be the maximum of the depths at the
1301!! neighoring tracer points.
1302subroutine set_velocity_depth_max(G)
1303 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
1304 ! This subroutine sets the 4 bottom depths at velocity points to be the
1305 ! maximum of the adjacent depths.
1306 integer :: i, j
1307
1308 do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1309 g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1310 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1311 enddo ; enddo
1312 do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1313 g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1314 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1315 enddo ; enddo
1316end subroutine set_velocity_depth_max
1317! -----------------------------------------------------------------------------
1318
1319! -----------------------------------------------------------------------------
1320!> Set the bathymetry at velocity points to be the minimum of the depths at the
1321!! neighoring tracer points.
1322subroutine set_velocity_depth_min(G)
1323 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
1324 ! This subroutine sets the 4 bottom depths at velocity points to be the
1325 ! minimum of the adjacent depths.
1326 integer :: i, j
1327
1328 do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1329 g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1330 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1331 enddo ; enddo
1332 do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1333 g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1334 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1335 enddo ; enddo
1336end subroutine set_velocity_depth_min
1337! -----------------------------------------------------------------------------
1338
1339! -----------------------------------------------------------------------------
1340!> Pre-compute global integrals of grid quantities (like masked ocean area) for
1341!! later use in reporting diagnostics
1342subroutine compute_global_grid_integrals(G, US)
1343 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
1344 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1345
1346 ! Local variables
1347 real, dimension(G%isc:G%iec, G%jsc:G%jec) :: masked_area ! Masked cell areas [L2 ~> m2]
1348 integer :: i, j
1349
1350 masked_area(:,:) = 0.
1351 g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1352 do j=g%jsc,g%jec ; do i=g%isc,g%iec
1353 masked_area(i,j) = g%areaT(i,j) * g%mask2dT(i,j)
1354 enddo ; enddo
1355 g%areaT_global = reproducing_sum(masked_area, unscale=us%L_to_m**2)
1356
1357 if (g%areaT_global == 0.0) &
1358 call mom_error(fatal, "compute_global_grid_integrals: zero ocean area (check topography?)")
1359
1360 g%IareaT_global = 1.0 / g%areaT_global
1361end subroutine compute_global_grid_integrals
1362! -----------------------------------------------------------------------------
1363
1364! -----------------------------------------------------------------------------
1365!> Write out a file describing the topography, Coriolis parameter, grid locations
1366!! and various other fixed fields from the grid.
1367subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file)
1368 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
1369 type(param_file_type), intent(in) :: param_file !< Parameter file structure
1370 character(len=*), intent(in) :: directory !< The directory into which to place the geometry file.
1371 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1372 character(len=*), optional, intent(in) :: geom_file !< If present, the name of the geometry file
1373 !! (otherwise the file is "ocean_geometry")
1374
1375 ! Local variables.
1376 character(len=240) :: filepath ! The full path to the file to write
1377 character(len=40) :: mdl = "write_ocean_geometry_file"
1378 character(len=32) :: filename_appendix = '' ! Appendix to geom filename for ensemble runs
1379 type(vardesc), dimension(:), allocatable :: &
1380 vars ! Types with metadata about the variables and their staggering
1381 type(mom_field), dimension(:), allocatable :: &
1382 fields ! Opaque types used by MOM_io to store variable metadata information
1383 type(mom_infra_file) :: io_handle ! The I/O handle of the fileset
1384 integer :: nflds ! The number of variables in this file
1385 integer :: file_threading
1386 integer :: geom_file_len ! geometry file name length
1387 logical :: multiple_files
1388
1389 call calltree_enter('write_ocean_geometry_file()')
1390
1391 nflds = 19 ; if (g%bathymetry_at_vel) nflds = 23
1392
1393 allocate(vars(nflds))
1394 allocate(fields(nflds))
1395
1396 ! var_desc populates a type defined in MOM_io.F90. The arguments, in order, are:
1397 ! (1) the variable name for the NetCDF file
1398 ! (2) the units of the variable when output
1399 ! (3) the variable's long name
1400 ! (4) a character indicating the horizontal grid, which may be '1' (column),
1401 ! 'h', 'q', 'u', or 'v', for the corresponding C-grid variable
1402 ! (5) a character indicating the vertical grid, which may be 'L' (layer),
1403 ! 'i' (interface), or '1' (no vertical location)
1404 ! (6) a character indicating the time levels of the field, which may be
1405 ! 's' (snap-shot), 'p' (periodic), or '1' (no time variation)
1406 vars(1) = var_desc("geolatb","degree","latitude at corner (Bu) points",'q','1','1')
1407 vars(2) = var_desc("geolonb","degree","longitude at corner (Bu) points",'q','1','1')
1408 vars(3) = var_desc("geolat","degree", "latitude at tracer (T) points", 'h','1','1')
1409 vars(4) = var_desc("geolon","degree","longitude at tracer (T) points",'h','1','1')
1410 vars(5) = var_desc("D","meter","Basin Depth",'h','1','1')
1411 vars(6) = var_desc("f","s-1","Coriolis Parameter",'q','1','1')
1412 vars(7) = var_desc("dxCv","m","Zonal grid spacing at v points",'v','1','1')
1413 vars(8) = var_desc("dyCu","m","Meridional grid spacing at u points",'u','1','1')
1414 vars(9) = var_desc("dxCu","m","Zonal grid spacing at u points",'u','1','1')
1415 vars(10)= var_desc("dyCv","m","Meridional grid spacing at v points",'v','1','1')
1416 vars(11)= var_desc("dxT","m","Zonal grid spacing at h points",'h','1','1')
1417 vars(12)= var_desc("dyT","m","Meridional grid spacing at h points",'h','1','1')
1418 vars(13)= var_desc("dxBu","m","Zonal grid spacing at q points",'q','1','1')
1419 vars(14)= var_desc("dyBu","m","Meridional grid spacing at q points",'q','1','1')
1420 vars(15)= var_desc("Ah","m2","Area of h cells",'h','1','1')
1421 vars(16)= var_desc("Aq","m2","Area of q cells",'q','1','1')
1422
1423 vars(17)= var_desc("dxCvo","m","Open zonal grid spacing at v points",'v','1','1')
1424 vars(18)= var_desc("dyCuo","m","Open meridional grid spacing at u points",'u','1','1')
1425 vars(19)= var_desc("wet", "nondim", "land or ocean?", 'h','1','1')
1426
1427 if (g%bathymetry_at_vel) then
1428 vars(20) = var_desc("Dblock_u","m","Blocked depth at u points",'u','1','1')
1429 vars(21) = var_desc("Dopen_u","m","Open depth at u points",'u','1','1')
1430 vars(22) = var_desc("Dblock_v","m","Blocked depth at v points",'v','1','1')
1431 vars(23) = var_desc("Dopen_v","m","Open depth at v points",'v','1','1')
1432 endif
1433
1434 if (present(geom_file)) then
1435 filepath = trim(directory) // trim(geom_file)
1436 else
1437 filepath = trim(directory) // "ocean_geometry"
1438 endif
1439
1440 ! Append ensemble run number to filename if it is an ensemble run
1441 call get_filename_appendix(filename_appendix)
1442 if (len_trim(filename_appendix) > 0) then
1443 geom_file_len = len_trim(filepath)
1444 if (filepath(geom_file_len-2:geom_file_len) == ".nc") then
1445 filepath = filepath(1:geom_file_len-3) // '.' // trim(filename_appendix) // ".nc"
1446 else
1447 filepath = filepath // '.' // trim(filename_appendix)
1448 endif
1449 endif
1450
1451 call get_param(param_file, mdl, "PARALLEL_RESTARTFILES", multiple_files, &
1452 "If true, the IO layout is used to group processors that write to the same "//&
1453 "restart file or each processor writes its own (numbered) restart file. "//&
1454 "If false, a single restart file is generated combining output from all PEs.", &
1455 default=.false.)
1456 file_threading = single_file
1457 if (multiple_files) file_threading = multiple
1458
1459 call create_mom_file(io_handle, trim(filepath), vars, nflds, fields, &
1460 file_threading, dg=g)
1461
1462 call mom_write_field(io_handle, fields(1), g%Domain, g%geoLatBu)
1463 call mom_write_field(io_handle, fields(2), g%Domain, g%geoLonBu)
1464 call mom_write_field(io_handle, fields(3), g%Domain, g%geoLatT)
1465 call mom_write_field(io_handle, fields(4), g%Domain, g%geoLonT)
1466
1467 call mom_write_field(io_handle, fields(5), g%Domain, g%bathyT, unscale=us%Z_to_m)
1468 call mom_write_field(io_handle, fields(6), g%Domain, g%CoriolisBu, unscale=us%s_to_T)
1469
1470 call mom_write_field(io_handle, fields(7), g%Domain, g%dxCv, unscale=us%L_to_m)
1471 call mom_write_field(io_handle, fields(8), g%Domain, g%dyCu, unscale=us%L_to_m)
1472 call mom_write_field(io_handle, fields(9), g%Domain, g%dxCu, unscale=us%L_to_m)
1473 call mom_write_field(io_handle, fields(10), g%Domain, g%dyCv, unscale=us%L_to_m)
1474 call mom_write_field(io_handle, fields(11), g%Domain, g%dxT, unscale=us%L_to_m)
1475 call mom_write_field(io_handle, fields(12), g%Domain, g%dyT, unscale=us%L_to_m)
1476 call mom_write_field(io_handle, fields(13), g%Domain, g%dxBu, unscale=us%L_to_m)
1477 call mom_write_field(io_handle, fields(14), g%Domain, g%dyBu, unscale=us%L_to_m)
1478
1479 call mom_write_field(io_handle, fields(15), g%Domain, g%areaT, unscale=us%L_to_m**2)
1480 call mom_write_field(io_handle, fields(16), g%Domain, g%areaBu, unscale=us%L_to_m**2)
1481
1482 call mom_write_field(io_handle, fields(17), g%Domain, g%dx_Cv, unscale=us%L_to_m)
1483 call mom_write_field(io_handle, fields(18), g%Domain, g%dy_Cu, unscale=us%L_to_m)
1484 call mom_write_field(io_handle, fields(19), g%Domain, g%mask2dT)
1485
1486 if (g%bathymetry_at_vel) then
1487 call mom_write_field(io_handle, fields(20), g%Domain, g%Dblock_u, unscale=us%Z_to_m)
1488 call mom_write_field(io_handle, fields(21), g%Domain, g%Dopen_u, unscale=us%Z_to_m)
1489 call mom_write_field(io_handle, fields(22), g%Domain, g%Dblock_v, unscale=us%Z_to_m)
1490 call mom_write_field(io_handle, fields(23), g%Domain, g%Dopen_v, unscale=us%Z_to_m)
1491 endif
1492
1493 call io_handle%close()
1494
1495 deallocate(vars, fields)
1496
1497 call calltree_leave('write_ocean_geometry_file()')
1498end subroutine write_ocean_geometry_file
1499