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 :: found_u, found_v
876 logical :: unit_in_use
877 logical :: fatal_unused_lengths
878 integer :: unused
879 integer :: ios, iounit, isu, isv
880 integer :: num_lines, nl_read, ln, npt, u_pt, v_pt
881 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
882 integer :: isu_por, isv_por
883 logical :: found_u_por, found_v_por
884
885 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
886 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
887
888 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
889
890 call get_param(param_file, mdl, "CHANNEL_LIST_FILE", chan_file, &
891 "The file from which the list of narrowed channels is read.", &
892 default="MOM_channel_list")
893 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
894 inputdir = slasher(inputdir)
895 filename = trim(inputdir)//trim(chan_file)
896 call log_param(param_file, mdl, "INPUTDIR/CHANNEL_LIST_FILE", filename)
897 call get_param(param_file, mdl, "CHANNEL_LIST_360_LON_CHECK", check_360, &
898 "If true, the channel configuration list works for any "//&
899 "longitudes in the range of -360 to 360.", default=.true.)
900 call get_param(param_file, mdl, "FATAL_UNUSED_CHANNEL_WIDTHS", fatal_unused_lengths, &
901 "If true, trigger a fatal error if there are any channel widths in "//&
902 "CHANNEL_LIST_FILE that do not cause any open face widths to change.", &
903 default=.false.)
904
905 if (is_root_pe()) then
906 ! Open the input file.
907 if (.not.file_exists(filename)) call mom_error(fatal, &
908 " reset_face_lengths_list: Unable to open "//trim(filename))
909
910 ! Find an unused unit number.
911 do iounit=10,512
912 INQUIRE(iounit,opened=unit_in_use) ; if (.not.unit_in_use) exit
913 enddo
914 if (iounit >= 512) call mom_error(fatal, &
915 "reset_face_lengths_list: No unused file unit could be found.")
916
917 ! Open the parameter file.
918 open(iounit, file=trim(filename), access='SEQUENTIAL', &
919 form='FORMATTED', action='READ', position='REWIND', iostat=ios)
920 if (ios /= 0) call mom_error(fatal, &
921 "reset_face_lengths_list: Error opening "//trim(filename))
922
923 ! Count the number of u_width and v_width entries.
924 call read_face_length_list(iounit, filename, num_lines, lines)
925 else
926 num_lines = 0
927 endif
928
929 len_lon = 360.0 ; if (g%len_lon > 0.0) len_lon = g%len_lon
930 len_lat = 180.0 ; if (g%len_lat > 0.0) len_lat = g%len_lat
931 ! Broadcast the number of lines and allocate the required space.
932 call broadcast(num_lines, root_pe())
933 u_pt = 0 ; v_pt = 0
934 if (num_lines > 0) then
935 allocate(lines(num_lines))
936
937 allocate(u_lat(2,num_lines), source=-1e34)
938 allocate(u_lon(2,num_lines), source=-1e34)
939 allocate(u_width(num_lines), source=-1e34)
940 allocate(u_line_used(num_lines), source=0)
941 allocate(u_line_no(num_lines), source=0)
942
943 allocate(v_lat(2,num_lines), source=-1e34)
944 allocate(v_lon(2,num_lines), source=-1e34)
945 allocate(v_width(num_lines), source=-1e34)
946 allocate(v_line_used(num_lines), source=0)
947 allocate(v_line_no(num_lines), source=0)
948
949 allocate(dmin_u(num_lines), source=0.0)
950 allocate(dmax_u(num_lines), source=0.0)
951 allocate(davg_u(num_lines), source=0.0)
952
953 allocate(dmin_v(num_lines), source=0.0)
954 allocate(dmax_v(num_lines), source=0.0)
955 allocate(davg_v(num_lines), source=0.0)
956
957 ! Actually read the lines.
958 if (is_root_pe()) then
959 call read_face_length_list(iounit, filename, nl_read, lines)
960 if (nl_read /= num_lines) &
961 call mom_error(fatal, 'reset_face_lengths_list : Found different '// &
962 'number of valid lines on second reading of '//trim(filename))
963 close(iounit) ; iounit = -1
964 endif
965
966 ! Broadcast the lines.
967 call broadcast(lines, 120, root_pe())
968
969 ! Populate the u_width, etc., data.
970 do ln=1,num_lines
971 line = lines(ln)
972 ! Detect keywords
973 found_u = .false. ; found_v = .false.
974 found_u_por = .false. ; found_v_por = .false.
975 isu = index(uppercase(line), "U_WIDTH") ; if (isu > 0) found_u = .true.
976 isv = index(uppercase(line), "V_WIDTH") ; if (isv > 0) found_v = .true.
977 isu_por = index(uppercase(line), "U_WIDTH_POR") ; if (isu_por > 0) found_u_por = .true.
978 isv_por = index(uppercase(line), "V_WIDTH_POR") ; if (isv_por > 0) found_v_por = .true.
979
980 ! Store and check the relevant values.
981 if (found_u) then
982 u_pt = u_pt + 1
983 if (found_u_por .eqv. .false.) then
984 read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
985 elseif (found_u_por) then
986 read(line(isu_por+12:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt), &
987 dmin_u(u_pt), dmax_u(u_pt), davg_u(u_pt)
988 endif
989 u_width(u_pt) = us%m_to_L*u_width(u_pt) ! Rescale units equivalently to scale=US%m_to_L during read.
990 dmin_u(u_pt) = us%m_to_Z*dmin_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
991 dmax_u(u_pt) = us%m_to_Z*dmax_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
992 davg_u(u_pt) = us%m_to_Z*davg_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
993 u_line_no(u_pt) = ln
994 if (is_root_pe()) then
995 if (check_360) then
996 if ((abs(u_lon(1,u_pt)) > len_lon) .or. (abs(u_lon(2,u_pt)) > len_lon)) &
997 call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
998 "u-longitude found when reading line "//trim(line)//" from file "//&
999 trim(filename))
1000 if ((abs(u_lat(1,u_pt)) > len_lat) .or. (abs(u_lat(2,u_pt)) > len_lat)) &
1001 call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
1002 "u-latitude found when reading line "//trim(line)//" from file "//&
1003 trim(filename))
1004 endif
1005 if (u_lat(1,u_pt) > u_lat(2,u_pt)) &
1006 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1007 "u-face latitudes found when reading line "//trim(line)//" from file "//&
1008 trim(filename))
1009 if (u_lon(1,u_pt) > u_lon(2,u_pt)) &
1010 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1011 "u-face longitudes found when reading line "//trim(line)//" from file "//&
1012 trim(filename))
1013 if (u_width(u_pt) < 0.0) &
1014 call mom_error(warning, "reset_face_lengths_list : Negative "//&
1015 "u-width found when reading line "//trim(line)//" from file "//&
1016 trim(filename))
1017 if (dmin_u(u_pt) > dmax_u(u_pt)) &
1018 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1019 "topographical min/max found when reading line "//trim(line)//" from file "//&
1020 trim(filename))
1021 endif
1022 elseif (found_v) then
1023 v_pt = v_pt + 1
1024 if (found_v_por .eqv. .false.) then
1025 read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
1026 elseif (found_v_por) then
1027 read(line(isv+12:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt), &
1028 dmin_v(v_pt), dmax_v(v_pt), davg_v(v_pt)
1029 endif
1030 v_width(v_pt) = us%m_to_L*v_width(v_pt) ! Rescale units equivalently to scale=US%m_to_L during read.
1031 dmin_v(v_pt) = us%m_to_Z*dmin_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1032 dmax_v(v_pt) = us%m_to_Z*dmax_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1033 davg_v(v_pt) = us%m_to_Z*davg_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1034 v_line_no(v_pt) = ln
1035 if (is_root_pe()) then
1036 if (check_360) then
1037 if ((abs(v_lon(1,v_pt)) > len_lon) .or. (abs(v_lon(2,v_pt)) > len_lon)) &
1038 call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
1039 "v-longitude found when reading line "//trim(line)//" from file "//&
1040 trim(filename))
1041 if ((abs(v_lat(1,v_pt)) > len_lat) .or. (abs(v_lat(2,v_pt)) > len_lat)) &
1042 call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
1043 "v-latitude found when reading line "//trim(line)//" from file "//&
1044 trim(filename))
1045 endif
1046 if (v_lat(1,v_pt) > v_lat(2,v_pt)) &
1047 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1048 "v-face latitudes found when reading line "//trim(line)//" from file "//&
1049 trim(filename))
1050 if (v_lon(1,v_pt) > v_lon(2,v_pt)) &
1051 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1052 "v-face longitudes found when reading line "//trim(line)//" from file "//&
1053 trim(filename))
1054 if (v_width(v_pt) < 0.0) &
1055 call mom_error(warning, "reset_face_lengths_list : Negative "//&
1056 "v-width found when reading line "//trim(line)//" from file "//&
1057 trim(filename))
1058 if (dmin_v(v_pt) > dmax_v(v_pt)) &
1059 call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
1060 "topographical min/max found when reading line "//trim(line)//" from file "//&
1061 trim(filename))
1062 endif
1063 endif
1064 enddo
1065
1066 endif
1067
1068 do j=jsd,jed ; do i=isdb,iedb
1069 lat = g%geoLatCu(i,j) ; lon = g%geoLonCu(i,j)
1070 if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1071 else ; lon_p = lon ; lon_m = lon ; endif
1072
1073 do npt=1,u_pt
1074 if (((lat >= u_lat(1,npt)) .and. (lat <= u_lat(2,npt))) .and. &
1075 (((lon >= u_lon(1,npt)) .and. (lon <= u_lon(2,npt))) .or. &
1076 ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. &
1077 ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then
1078
1079 g%dy_Cu(i,j) = g%mask2dCu(i,j) * min(g%dyCu(i,j), max(u_width(npt), 0.0))
1080 g%porous_DminU(i,j) = dmin_u(npt)
1081 g%porous_DmaxU(i,j) = dmax_u(npt)
1082 g%porous_DavgU(i,j) = davg_u(npt)
1083
1084 if (j>=g%jsc .and. j<=g%jec .and. i>=g%isc .and. i<=g%iec) then ! Limit messages/checking to compute domain
1085 if ( g%mask2dCu(i,j) == 0.0 ) then
1086 write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",&
1087 u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") so grid metric is unmodified."
1088 else
1089 u_line_used(npt) = u_line_used(npt) + 1
1090 write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1091 "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon," (",&
1092 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"
1093 write(stdout,'(A,3F8.2,A)') &
1094 "read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",g%porous_DminU(i,j),&
1095 g%porous_DmaxU(i,j), g%porous_DavgU(i,j),")m"
1096 endif
1097 endif
1098 endif
1099 enddo
1100
1101 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
1102 g%IareaCu(i,j) = 0.0
1103 if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
1104 enddo ; enddo
1105
1106 do j=jsdb,jedb ; do i=isd,ied
1107 lat = g%geoLatCv(i,j) ; lon = g%geoLonCv(i,j)
1108 if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1109 else ; lon_p = lon ; lon_m = lon ; endif
1110
1111 do npt=1,v_pt
1112 if (((lat >= v_lat(1,npt)) .and. (lat <= v_lat(2,npt))) .and. &
1113 (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. &
1114 ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. &
1115 ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then
1116 g%dx_Cv(i,j) = g%mask2dCv(i,j) * min(g%dxCv(i,j), max(v_width(npt), 0.0))
1117 g%porous_DminV(i,j) = dmin_v(npt)
1118 g%porous_DmaxV(i,j) = dmax_v(npt)
1119 g%porous_DavgV(i,j) = davg_v(npt)
1120
1121 if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec) then ! Limit messages/checking to compute domain
1122 if ( g%mask2dCv(i,j) == 0.0 ) then
1123 write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",&
1124 v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") so grid metric is unmodified."
1125 else
1126 v_line_used(npt) = v_line_used(npt) + 1
1127 write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1128 "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon," (",&
1129 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"
1130 write(stdout,'(A,3F8.2,A)') &
1131 "read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",g%porous_DminV(i,j),&
1132 g%porous_DmaxV(i,j), g%porous_DavgV(i,j),")m"
1133 endif
1134 endif
1135 endif
1136 enddo
1137
1138 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
1139 g%IareaCv(i,j) = 0.0
1140 if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
1141 enddo ; enddo
1142
1143 ! Verify that all channel widths have been used
1144 unused = 0
1145 if (u_pt > 0) call sum_across_pes(u_line_used, u_pt)
1146 if (v_pt > 0) call sum_across_pes(v_line_used, v_pt)
1147 if (is_root_pe()) then
1148 unused = 0
1149 do npt=1,u_pt ; if (u_line_used(npt) == 0) then
1150 call mom_error(warning, "reset_face_lengths_list unused u-face line: "//&
1151 trim(lines(u_line_no(npt))) )
1152 unused = unused + 1
1153 endif ; enddo
1154 do npt=1,v_pt ; if (v_line_used(npt) == 0) then
1155 call mom_error(warning, "reset_face_lengths_list unused v-face line: "//&
1156 trim(lines(v_line_no(npt))) )
1157 unused = unused + 1
1158 endif ; enddo
1159 if (fatal_unused_lengths .and. (unused > 0)) call mom_error(fatal, &
1160 "reset_face_lengths_list causing MOM6 abort due to unused face length lines.")
1161 endif
1162
1163 if (num_lines > 0) then
1164 deallocate(lines)
1165 deallocate(u_line_used, v_line_used, u_line_no, v_line_no)
1166 deallocate(u_lat) ; deallocate(u_lon) ; deallocate(u_width)
1167 deallocate(v_lat) ; deallocate(v_lon) ; deallocate(v_width)
1168 deallocate(dmin_u) ; deallocate(dmax_u) ; deallocate(davg_u)
1169 deallocate(dmin_v) ; deallocate(dmax_v) ; deallocate(davg_v)
1170 endif
1171
1172 call calltree_leave(trim(mdl)//'()')
1173end subroutine reset_face_lengths_list
1174! -----------------------------------------------------------------------------
1175
1176! -----------------------------------------------------------------------------
1177!> This subroutine reads and counts the non-blank lines in the face length list file, after removing comments.
1178subroutine read_face_length_list(iounit, filename, num_lines, lines)
1179 integer, intent(in) :: iounit !< An open I/O unit number for the file
1180 character(len=*), intent(in) :: filename !< The name of the face-length file to read
1181 integer, intent(out) :: num_lines !< The number of non-blank lines in the file
1182 character(len=120), dimension(:), pointer :: lines !< The non-blank lines, after removing comments
1183
1184 ! This subroutine reads and counts the non-blank lines in the face length
1185 ! list file, after removing comments.
1186 character(len=120) :: line, line_up
1187 logical :: found_u, found_v
1188 integer :: isu, isv, icom
1189 integer :: last
1190
1191 num_lines = 0
1192
1193 if (iounit <= 0) return
1194 rewind(iounit)
1195 do while(.true.)
1196 read(iounit, '(a)', end=8, err=9) line
1197 last = len_trim(line)
1198 ! Eliminate either F90 or C comments from the line.
1199 icom = index(line(:last), "!") ; if (icom > 0) last = icom-1
1200 icom = index(line(:last), "/*") ; if (icom > 0) last = icom-1
1201 if (last < 1) cycle
1202
1203 ! Detect keywords
1204 line_up = uppercase(line)
1205 found_u = .false. ; found_v = .false.
1206 isu = index(line_up(:last), "U_WIDTH") ; if (isu > 0) found_u = .true.
1207 isv = index(line_up(:last), "V_WIDTH") ; if (isv > 0) found_v = .true.
1208
1209 if (found_u .and. found_v) call mom_error(fatal, &
1210 "read_face_length_list : both U_WIDTH and V_WIDTH found when "//&
1211 "reading the line "//trim(line(:last))//" in file "//trim(filename))
1212 if (found_u .or. found_v) then
1213 num_lines = num_lines + 1
1214 if (associated(lines)) then
1215 lines(num_lines) = line(1:last)
1216 endif
1217 endif
1218 enddo ! while (.true.)
1219
12208 continue
1221 return
1222
12239 call mom_error(fatal, "read_face_length_list : "//&
1224 "Error while reading file "//trim(filename))
1225
1226end subroutine read_face_length_list
1227! -----------------------------------------------------------------------------
1228
1229! -----------------------------------------------------------------------------
1230!> Read from a file the maximum, minimum and average bathymetry at velocity points,
1231!! for the use of porous barrier.
1232!! Note that we assume the depth values in the sub-grid bathymetry file of the same
1233!! convention as in-cell bathymetry file, i.e. positive below the sea surface and
1234!! increasing downward; while in subroutine reset_face_lengths_list, it is implied
1235!! that read-in fields min_bathy, max_bathy and avg_bathy from the input file
1236!! CHANNEL_LIST_FILE all have negative values below the surface. Therefore, to ensure
1237!! backward compatibility, all signs of the variable are inverted here.
1238!! And porous_Dmax[UV] = shallowest point, porous_Dmin[UV] = deepest point
1239subroutine set_subgrid_topo_at_vel_from_file(G, param_file, US)
1240 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
1241 type(param_file_type), intent(in) :: param_file !< Parameter file structure
1242 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1243
1244 ! Local variables
1245 character(len=200) :: filename, topo_file, inputdir ! Strings for file/path
1246 character(len=200) :: varname_uhi, varname_ulo, varname_uav, &
1247 varname_vhi, varname_vlo, varname_vav ! Variable names in file
1248 character(len=40) :: mdl = "set_subgrid_topo_at_vel_from_file" ! This subroutine's name.
1249 integer :: i, j
1250
1251 call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
1252
1253 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1254 inputdir = slasher(inputdir)
1255 call get_param(param_file, mdl, "TOPO_AT_VEL_FILE", topo_file, &
1256 "The file from which the bathymetry parameters at the velocity points are read. "//&
1257 "While the names of the parameters reflect their physical locations, i.e. HIGH is above LOW, "//&
1258 "their signs follow the model's convention, which is positive below the sea surface", &
1259 default="topog_edge.nc")
1260 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_HIGH", varname_uhi, &
1261 "The variable name of the highest bathymetry at the u-cells in TOPO_AT_VEL_FILE.", &
1262 default="depthu_hi")
1263 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_LOW", varname_ulo, &
1264 "The variable name of the lowest bathymetry at the u-cells in TOPO_AT_VEL_FILE.", &
1265 default="depthu_lo")
1266 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_U_AVE", varname_uav, &
1267 "The variable name of the average bathymetry at the u-cells in TOPO_AT_VEL_FILE.", &
1268 default="depthu_av")
1269 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_HIGH", varname_vhi, &
1270 "The variable name of the highest bathymetry at the v-cells in TOPO_AT_VEL_FILE.", &
1271 default="depthv_hi")
1272 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_LOW", varname_vlo, &
1273 "The variable name of the lowest bathymetry at the v-cells in TOPO_AT_VEL_FILE.", &
1274 default="depthv_lo")
1275 call get_param(param_file, mdl, "TOPO_AT_VEL_VARNAME_V_AVE", varname_vav, &
1276 "The variable name of the average bathymetry at the v-cells in TOPO_AT_VEL_FILE.", &
1277 default="depthv_av")
1278
1279 filename = trim(inputdir)//trim(topo_file)
1280 call log_param(param_file, mdl, "INPUTDIR/TOPO_AT_VEL_FILE", filename)
1281
1282 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
1283 " set_subgrid_topo_at_vel_from_file: Unable to open "//trim(filename))
1284
1285 call mom_read_vector(filename, trim(varname_uhi), trim(varname_vhi), &
1286 g%porous_DmaxU, g%porous_DmaxV, g%Domain, stagger=cgrid_ne, scale=us%m_to_Z)
1287 call mom_read_vector(filename, trim(varname_ulo), trim(varname_vlo), &
1288 g%porous_DminU, g%porous_DminV, g%Domain, stagger=cgrid_ne, scale=us%m_to_Z)
1289 call mom_read_vector(filename, trim(varname_uav), trim(varname_vav), &
1290 g%porous_DavgU, g%porous_DavgV, g%Domain, stagger=cgrid_ne, scale=us%m_to_Z)
1291
1292 ! The signs of the depth parameters need to be inverted to be backward compatible with input files
1293 ! used by subroutine reset_face_lengths_list, which assumes depth is negative below the sea surface.
1294 g%porous_DmaxU = -g%porous_DmaxU ; g%porous_DminU = -g%porous_DminU ; g%porous_DavgU = -g%porous_DavgU
1295 g%porous_DmaxV = -g%porous_DmaxV ; g%porous_DminV = -g%porous_DminV ; g%porous_DavgV = -g%porous_DavgV
1296
1297 call pass_vector(g%porous_DmaxU, g%porous_DmaxV, g%Domain, to_all+scalar_pair, cgrid_ne)
1298 call pass_vector(g%porous_DminU, g%porous_DminV, g%Domain, to_all+scalar_pair, cgrid_ne)
1299 call pass_vector(g%porous_DavgU, g%porous_DavgV, g%Domain, to_all+scalar_pair, cgrid_ne)
1300
1301 call calltree_leave(trim(mdl)//'()')
1303! -----------------------------------------------------------------------------
1304
1305! -----------------------------------------------------------------------------
1306!> Set the bathymetry at velocity points to be the maximum of the depths at the
1307!! neighoring tracer points.
1308subroutine set_velocity_depth_max(G)
1309 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
1310 ! This subroutine sets the 4 bottom depths at velocity points to be the
1311 ! maximum of the adjacent depths.
1312 integer :: i, j
1313
1314 do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1315 g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1316 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1317 enddo ; enddo
1318 do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1319 g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1320 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1321 enddo ; enddo
1322end subroutine set_velocity_depth_max
1323! -----------------------------------------------------------------------------
1324
1325! -----------------------------------------------------------------------------
1326!> Set the bathymetry at velocity points to be the minimum of the depths at the
1327!! neighoring tracer points.
1328subroutine set_velocity_depth_min(G)
1329 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
1330 ! This subroutine sets the 4 bottom depths at velocity points to be the
1331 ! minimum of the adjacent depths.
1332 integer :: i, j
1333
1334 do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1335 g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1336 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1337 enddo ; enddo
1338 do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1339 g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1340 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1341 enddo ; enddo
1342end subroutine set_velocity_depth_min
1343! -----------------------------------------------------------------------------
1344
1345! -----------------------------------------------------------------------------
1346!> Pre-compute global integrals of grid quantities (like masked ocean area) for
1347!! later use in reporting diagnostics
1348subroutine compute_global_grid_integrals(G, US)
1349 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
1350 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1351
1352 ! Local variables
1353 real, dimension(G%isc:G%iec, G%jsc:G%jec) :: masked_area ! Masked cell areas [L2 ~> m2]
1354 integer :: i, j
1355
1356 masked_area(:,:) = 0.
1357 g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1358 do j=g%jsc,g%jec ; do i=g%isc,g%iec
1359 masked_area(i,j) = g%areaT(i,j) * g%mask2dT(i,j)
1360 enddo ; enddo
1361 g%areaT_global = reproducing_sum(masked_area, unscale=us%L_to_m**2)
1362
1363 if (g%areaT_global == 0.0) &
1364 call mom_error(fatal, "compute_global_grid_integrals: zero ocean area (check topography?)")
1365
1366 g%IareaT_global = 1.0 / g%areaT_global
1367end subroutine compute_global_grid_integrals
1368! -----------------------------------------------------------------------------
1369
1370! -----------------------------------------------------------------------------
1371!> Write out a file describing the topography, Coriolis parameter, grid locations
1372!! and various other fixed fields from the grid.
1373subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file)
1374 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid
1375 type(param_file_type), intent(in) :: param_file !< Parameter file structure
1376 character(len=*), intent(in) :: directory !< The directory into which to place the geometry file.
1377 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1378 character(len=*), optional, intent(in) :: geom_file !< If present, the name of the geometry file
1379 !! (otherwise the file is "ocean_geometry")
1380
1381 ! Local variables.
1382 character(len=240) :: filepath ! The full path to the file to write
1383 character(len=40) :: mdl = "write_ocean_geometry_file"
1384 character(len=32) :: filename_appendix = '' ! Appendix to geom filename for ensemble runs
1385 type(vardesc), dimension(:), allocatable :: &
1386 vars ! Types with metadata about the variables and their staggering
1387 type(mom_field), dimension(:), allocatable :: &
1388 fields ! Opaque types used by MOM_io to store variable metadata information
1389 type(mom_infra_file) :: io_handle ! The I/O handle of the fileset
1390 integer :: nflds ! The number of variables in this file
1391 integer :: file_threading
1392 integer :: geom_file_len ! geometry file name length
1393 logical :: multiple_files
1394
1395 call calltree_enter('write_ocean_geometry_file()')
1396
1397 nflds = 19 ; if (g%bathymetry_at_vel) nflds = 23
1398
1399 allocate(vars(nflds))
1400 allocate(fields(nflds))
1401
1402 ! var_desc populates a type defined in MOM_io.F90. The arguments, in order, are:
1403 ! (1) the variable name for the NetCDF file
1404 ! (2) the units of the variable when output
1405 ! (3) the variable's long name
1406 ! (4) a character indicating the horizontal grid, which may be '1' (column),
1407 ! 'h', 'q', 'u', or 'v', for the corresponding C-grid variable
1408 ! (5) a character indicating the vertical grid, which may be 'L' (layer),
1409 ! 'i' (interface), or '1' (no vertical location)
1410 ! (6) a character indicating the time levels of the field, which may be
1411 ! 's' (snap-shot), 'p' (periodic), or '1' (no time variation)
1412 vars(1) = var_desc("geolatb","degree","latitude at corner (Bu) points",'q','1','1')
1413 vars(2) = var_desc("geolonb","degree","longitude at corner (Bu) points",'q','1','1')
1414 vars(3) = var_desc("geolat","degree", "latitude at tracer (T) points", 'h','1','1')
1415 vars(4) = var_desc("geolon","degree","longitude at tracer (T) points",'h','1','1')
1416 vars(5) = var_desc("D","meter","Basin Depth",'h','1','1')
1417 vars(6) = var_desc("f","s-1","Coriolis Parameter",'q','1','1')
1418 vars(7) = var_desc("dxCv","m","Zonal grid spacing at v points",'v','1','1')
1419 vars(8) = var_desc("dyCu","m","Meridional grid spacing at u points",'u','1','1')
1420 vars(9) = var_desc("dxCu","m","Zonal grid spacing at u points",'u','1','1')
1421 vars(10)= var_desc("dyCv","m","Meridional grid spacing at v points",'v','1','1')
1422 vars(11)= var_desc("dxT","m","Zonal grid spacing at h points",'h','1','1')
1423 vars(12)= var_desc("dyT","m","Meridional grid spacing at h points",'h','1','1')
1424 vars(13)= var_desc("dxBu","m","Zonal grid spacing at q points",'q','1','1')
1425 vars(14)= var_desc("dyBu","m","Meridional grid spacing at q points",'q','1','1')
1426 vars(15)= var_desc("Ah","m2","Area of h cells",'h','1','1')
1427 vars(16)= var_desc("Aq","m2","Area of q cells",'q','1','1')
1428
1429 vars(17)= var_desc("dxCvo","m","Open zonal grid spacing at v points",'v','1','1')
1430 vars(18)= var_desc("dyCuo","m","Open meridional grid spacing at u points",'u','1','1')
1431 vars(19)= var_desc("wet", "nondim", "land or ocean?", 'h','1','1')
1432
1433 if (g%bathymetry_at_vel) then
1434 vars(20) = var_desc("Dblock_u","m","Blocked depth at u points",'u','1','1')
1435 vars(21) = var_desc("Dopen_u","m","Open depth at u points",'u','1','1')
1436 vars(22) = var_desc("Dblock_v","m","Blocked depth at v points",'v','1','1')
1437 vars(23) = var_desc("Dopen_v","m","Open depth at v points",'v','1','1')
1438 endif
1439
1440 if (present(geom_file)) then
1441 filepath = trim(directory) // trim(geom_file)
1442 else
1443 filepath = trim(directory) // "ocean_geometry"
1444 endif
1445
1446 ! Append ensemble run number to filename if it is an ensemble run
1447 call get_filename_appendix(filename_appendix)
1448 if (len_trim(filename_appendix) > 0) then
1449 geom_file_len = len_trim(filepath)
1450 if (filepath(geom_file_len-2:geom_file_len) == ".nc") then
1451 filepath = filepath(1:geom_file_len-3) // '.' // trim(filename_appendix) // ".nc"
1452 else
1453 filepath = filepath // '.' // trim(filename_appendix)
1454 endif
1455 endif
1456
1457 call get_param(param_file, mdl, "PARALLEL_RESTARTFILES", multiple_files, &
1458 "If true, the IO layout is used to group processors that write to the same "//&
1459 "restart file or each processor writes its own (numbered) restart file. "//&
1460 "If false, a single restart file is generated combining output from all PEs.", &
1461 default=.false.)
1462 file_threading = single_file
1463 if (multiple_files) file_threading = multiple
1464
1465 call create_mom_file(io_handle, trim(filepath), vars, nflds, fields, &
1466 file_threading, dg=g)
1467
1468 call mom_write_field(io_handle, fields(1), g%Domain, g%geoLatBu)
1469 call mom_write_field(io_handle, fields(2), g%Domain, g%geoLonBu)
1470 call mom_write_field(io_handle, fields(3), g%Domain, g%geoLatT)
1471 call mom_write_field(io_handle, fields(4), g%Domain, g%geoLonT)
1472
1473 call mom_write_field(io_handle, fields(5), g%Domain, g%bathyT, unscale=us%Z_to_m)
1474 call mom_write_field(io_handle, fields(6), g%Domain, g%CoriolisBu, unscale=us%s_to_T)
1475
1476 call mom_write_field(io_handle, fields(7), g%Domain, g%dxCv, unscale=us%L_to_m)
1477 call mom_write_field(io_handle, fields(8), g%Domain, g%dyCu, unscale=us%L_to_m)
1478 call mom_write_field(io_handle, fields(9), g%Domain, g%dxCu, unscale=us%L_to_m)
1479 call mom_write_field(io_handle, fields(10), g%Domain, g%dyCv, unscale=us%L_to_m)
1480 call mom_write_field(io_handle, fields(11), g%Domain, g%dxT, unscale=us%L_to_m)
1481 call mom_write_field(io_handle, fields(12), g%Domain, g%dyT, unscale=us%L_to_m)
1482 call mom_write_field(io_handle, fields(13), g%Domain, g%dxBu, unscale=us%L_to_m)
1483 call mom_write_field(io_handle, fields(14), g%Domain, g%dyBu, unscale=us%L_to_m)
1484
1485 call mom_write_field(io_handle, fields(15), g%Domain, g%areaT, unscale=us%L_to_m**2)
1486 call mom_write_field(io_handle, fields(16), g%Domain, g%areaBu, unscale=us%L_to_m**2)
1487
1488 call mom_write_field(io_handle, fields(17), g%Domain, g%dx_Cv, unscale=us%L_to_m)
1489 call mom_write_field(io_handle, fields(18), g%Domain, g%dy_Cu, unscale=us%L_to_m)
1490 call mom_write_field(io_handle, fields(19), g%Domain, g%mask2dT)
1491
1492 if (g%bathymetry_at_vel) then
1493 call mom_write_field(io_handle, fields(20), g%Domain, g%Dblock_u, unscale=us%Z_to_m)
1494 call mom_write_field(io_handle, fields(21), g%Domain, g%Dopen_u, unscale=us%Z_to_m)
1495 call mom_write_field(io_handle, fields(22), g%Domain, g%Dblock_v, unscale=us%Z_to_m)
1496 call mom_write_field(io_handle, fields(23), g%Domain, g%Dopen_v, unscale=us%Z_to_m)
1497 endif
1498
1499 call io_handle%close()
1500
1501 deallocate(vars, fields)
1502
1503 call calltree_leave('write_ocean_geometry_file()')
1504end subroutine write_ocean_geometry_file
1505