src/initialization/MOM_shared_initialization.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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.
7		module MOM_shared_initialization
8
9		use MOM_coms, only : max_across_PEs, reproducing_sum
10		use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast
11		use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID
12		use MOM_dyn_horgrid, only : dyn_horgrid_type
13		use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe
14		use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
15		use MOM_file_parser, only : get_param, log_param, param_file_type, log_version
16		use MOM_io, only : create_MOM_file, file_exists, field_size, get_filename_appendix
17		use MOM_io, only : MOM_infra_file, MOM_field
18		use MOM_io, only : MOM_read_data, MOM_read_vector, read_variable, stdout
19		use MOM_io, only : open_file_to_read, close_file_to_read, SINGLE_FILE, MULTIPLE
20		use MOM_io, only : slasher, vardesc, MOM_write_field, var_desc
21		use MOM_string_functions, only : uppercase
22		use MOM_unit_scaling, only : unit_scale_type
23
24		implicit none ; private
25
26		public MOM_shared_init_init
27		public MOM_initialize_rotation, MOM_calculate_grad_Coriolis
28		public initialize_topography_from_file, apply_topography_edits_from_file
29		public initialize_topography_named, limit_topography, diagnoseMaximumDepth
30		public set_rotation_planetary, set_rotation_beta_plane, initialize_grid_rotation_angle
31		public reset_face_lengths_named, reset_face_lengths_file, reset_face_lengths_list
32		public read_face_length_list, set_velocity_depth_max, set_velocity_depth_min
33		public set_subgrid_topo_at_vel_from_file
34		public compute_global_grid_integrals, write_ocean_geometry_file
35		public set_meanSL_from_file
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
42		contains
43
44		! -----------------------------------------------------------------------------
45		!> MOM_shared_init_init just writes the code version.
46	0	subroutine 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	0	"Sharable code to initialize time-invariant fields, like bathymetry and Coriolis parameters.")
56
57	0	end subroutine MOM_shared_init_init
58		! -----------------------------------------------------------------------------
59
60		!> MOM_initialize_rotation makes the appropriate call to set up the Coriolis parameter.
61	1	subroutine 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	1	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	1	default="2omegasinlat")
82	2	select case (trim(config))
83	1	case ("2omegasinlat"); call set_rotation_planetary(f, G, PF, US)
84	0	case ("beta"); call set_rotation_beta_plane(f, G, PF, US)
85	0	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	2	"Unrecognized rotation setup "//trim(config))
89		end select
90	1	call callTree_leave(trim(mdl)//'()')
91	1	end subroutine MOM_initialize_rotation
92
93		!> Calculates the components of grad f (Coriolis parameter)
94	1	subroutine 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	1	call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
107	1	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	0	dF_dx(:,:) = 0.0 ; dF_dy(:,:) = 0.0
111	0	return
112		endif
113
114	7261	do j=G%jsc, G%jec ; do i=G%isc, G%iec
115	7200	f1 = 0.5*( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) )
116	7200	f2 = 0.5*( G%CoriolisBu(I-1,J) + G%CoriolisBu(I-1,J-1) )
117	7200	dF_dx(i,j) = G%IdxT(i,j) * ( f1 - f2 )
118	7200	f1 = 0.5*( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) )
119	7200	f2 = 0.5*( G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J-1) )
120	7260	dF_dy(i,j) = G%IdyT(i,j) * ( f1 - f2 )
121		enddo ; enddo
122	1	call pass_vector(dF_dx, dF_dy, G%Domain, stagger=AGRID)
123	1	call callTree_leave(trim(mdl)//'()')
124
125		end subroutine MOM_calculate_grad_Coriolis
126
127		!> Return the global maximum ocean bottom depth in the same units as the input depth.
128	0	function 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	0	diagnoseMaximumDepth = D(G%isc,G%jsc)
136	0	do j=G%jsc, G%jec ; do i=G%isc, G%iec
137	0	diagnoseMaximumDepth = max(diagnoseMaximumDepth,D(i,j))
138		enddo ; enddo
139	0	call max_across_PEs(diagnoseMaximumDepth)
140	0	end function diagnoseMaximumDepth
141
142		!> Read time mean ocean sea level from a file
143	0	subroutine 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	0	call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
157
158	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
159	0	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	0	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	0	default="meanSL")
166	0	filename = trim(inputdir)//trim(file)
167	0	call log_param(param_file, mdl, "INPUTDIR/TOPO_FILE", filename)
168
169	0	if (.not.file_exists(filename, G%Domain)) &
170	0	call MOM_error(FATAL, " "//mdl//": Unable to open "//trim(filename))
171
172	0	call MOM_read_data(filename, trim(varname), meanSL, G%Domain, scale=US%m_to_Z)
173	0	call pass_var(meanSL, G%Domain)
174
175	0	call callTree_leave(trim(mdl)//'()')
176	0	end subroutine set_meanSL_from_file
177
178		!> Read gridded depths from file
179	0	subroutine 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	0	call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
191
192	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
193	0	inputdir = slasher(inputdir)
194		call get_param(param_file, mdl, "TOPO_FILE", topo_file, &
195		"The file from which the bathymetry is read.", &
196	0	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	0	default="depth")
200
201	0	filename = trim(inputdir)//trim(topo_file)
202	0	call log_param(param_file, mdl, "INPUTDIR/TOPO_FILE", filename)
203
204	0	if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
205	0	" initialize_topography_from_file: Unable to open "//trim(filename))
206
207	0	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	0	call MOM_read_data(filename, trim(topo_varname), D, G%Domain, scale=US%m_to_Z)
214
215	0	call apply_topography_edits_from_file(D, G, param_file, US)
216
217	0	call callTree_leave(trim(mdl)//'()')
218	0	end subroutine initialize_topography_from_file
219
220		!> Applies a list of topography overrides read from a netcdf file
221	0	subroutine 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	0	real, dimension(:), allocatable :: new_depth ! The new values of the depths [Z ~> m]
231	0	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	0	call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
240
241	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
242	0	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	0	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	0	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	0	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	0	units="m", default=-9999.0, scale=US%m_to_Z)
260	0	if (mask_depth == -9999.*US%m_to_Z) mask_depth = min_depth
261
262	0	if (len_trim(topo_edits_file)==0) return
263
264	0	topo_edits_file = trim(inputdir)//trim(topo_edits_file)
265	0	if (is_root_PE()) then
266	0	if (.not.file_exists(topo_edits_file, G%Domain)) &
267	0	call MOM_error(FATAL, trim(mdl)//': Unable to find file '//trim(topo_edits_file))
268	0	call open_file_to_read(topo_edits_file, ncid)
269		else
270	0	ncid = -1
271		endif
272
273		! Read and check the values of ni and nj in the file for consistency with this configuration.
274	0	call read_variable(topo_edits_file, 'ni', i_file, ncid_in=ncid)
275	0	call read_variable(topo_edits_file, 'nj', j_file, ncid_in=ncid)
276	0	if (i_file /= G%ieg) call MOM_error(FATAL, trim(mdl)//': Incompatible i-dimension of grid in '//&
277	0	trim(topo_edits_file))
278	0	if (j_file /= G%jeg) call MOM_error(FATAL, trim(mdl)//': Incompatible j-dimension of grid in '//&
279	0	trim(topo_edits_file))
280
281		! Get nEdits
282	0	call field_size(topo_edits_file, 'zEdit', sizes, ndims=ndims, ncid_in=ncid)
283	0	if (ndims /= 1) call MOM_error(FATAL, "The variable zEdit has an "//&
284	0	"unexpected number of dimensions in "//trim(topo_edits_file) )
285	0	n_edits = sizes(1)
286	0	allocate(ig(n_edits))
287	0	allocate(jg(n_edits))
288	0	allocate(new_depth(n_edits))
289
290		! Read iEdit, jEdit and zEdit
291	0	call read_variable(topo_edits_file, 'iEdit', ig, ncid_in=ncid)
292	0	call read_variable(topo_edits_file, 'jEdit', jg, ncid_in=ncid)
293	0	call read_variable(topo_edits_file, 'zEdit', new_depth, ncid_in=ncid, scale=US%m_to_Z)
294	0	call close_file_to_read(ncid, topo_edits_file)
295
296	0	do n = 1, n_edits
297	0	i = ig(n) - G%idg_offset + 1 ! +1 for python indexing
298	0	j = jg(n) - G%jdg_offset + 1
299	0	if (i>=G%isc .and. i<=G%iec .and. j>=G%jsc .and. j<=G%jec) then
300	0	if (new_depth(n) /= mask_depth) then
301		write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') &
302	0	'Ocean topography edit: ', n, ig(n), jg(n), D(i,j)US%Z_to_m, '->', abs(US%Z_to_mnew_depth(n)), i, j
303	0	D(i,j) = abs(new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
304		else
305	0	if (topo_edits_change_mask) then
306		write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') &
307	0	'Ocean topography edit: ',n,ig(n),jg(n),D(i,j)US%Z_to_m,'->',abs(US%Z_to_mnew_depth(n)),i,j
308	0	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	0	"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	0	deallocate( ig, jg, new_depth )
318
319	0	call callTree_leave(trim(mdl)//'()')
320	0	end subroutine apply_topography_edits_from_file
321
322		!> initialize the bathymetry based on one of several named idealized configurations
323	0	subroutine 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	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
344	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
345
346	0	call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
347		call MOM_mesg(" MOM_shared_initialization.F90, initialize_topography_named: "//&
348	0	"TOPO_CONFIG = "//trim(topog_config), 5)
349
350		call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
351	0	"The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z)
352	0	if (max_depth<=0.) call MOM_error(FATAL,"initialize_topography_named: "// &
353	0	"MAXIMUM_DEPTH has a non-sensical value! Was it set?")
354
355	0	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	0	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	0	"the named topographies.", units="m", default=400000.0, scale=US%m_to_L)
362		endif
363
364
365	0	PI = 4.0*atan(1.0)
366
367	0	if (trim(topog_config) == "flat") then
368	0	do j=js,je ; do i=is,ie ; D(i,j) = max_depth ; enddo ; enddo
369	0	elseif (trim(topog_config) == "spoon") then
370		D0 = (max_depth - Dedge) / &
371		((1.0 - exp(-0.5G%len_latG%Rad_Earth_LPI/(180.0 expdecay))) * &
372	0	(1.0 - exp(-0.5G%len_latG%Rad_Earth_LPI/(180.0 expdecay))))
373	0	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_LPI / &
379	0	(180.0*expdecay)) ))
380		enddo ; enddo
381	0	elseif (trim(topog_config) == "bowl") then
382		D0 = (max_depth - Dedge) / &
383		((1.0 - exp(-0.5G%len_latG%Rad_Earth_LPI/(180.0 expdecay))) * &
384	0	(1.0 - exp(-0.5G%len_latG%Rad_Earth_LPI/(180.0 expdecay))))
385
386		! This sets a bowl shaped (sort of) bottom topography, with a
387		! maximum depth of max_depth.
388	0	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_LPI/ &
392		(180.0expdecay))) &
393		(1.0 - exp((G%geoLatT(i,j) - (G%south_lat+G%len_lat))* &
394	0	G%Rad_Earth_LPI/(180.0expdecay)))))
395		enddo ; enddo
396	0	elseif (trim(topog_config) == "halfpipe") then
397	0	D0 = max_depth - Dedge
398	0	do j=js,je ; do i=is,ie
399	0	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	0	"Unrecognized topography name "//trim(topog_config))
404		endif
405
406		! This is here just for safety. Hopefully it doesn't do anything.
407	0	do j=js,je ; do i=is,ie
408	0	if (D(i,j) > max_depth) D(i,j) = max_depth
409	0	if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth
410		enddo ; enddo
411
412	0	call callTree_leave(trim(mdl)//'()')
413	0	end subroutine initialize_topography_named
414		! -----------------------------------------------------------------------------
415
416		! -----------------------------------------------------------------------------
417		!> limit_topography ensures that min_depth < D(x,y) < max_depth
418	1	subroutine 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	1	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	1	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	1	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	1	if (mask_depth == -9999.*US%m_to_Z) then
450	1	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	0	"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	8773	do j=G%jsd,G%jed ; do i=G%isd,G%ied
460	8772	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	0	do j=G%jsd,G%jed ; do i=G%isd,G%ied
468	0	if (D(i,j) > min(min_depth,mask_depth)) then
469	0	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	0	D(i,j) = min(min_depth,mask_depth)
474		endif
475		enddo ; enddo
476		endif
477
478	1	call callTree_leave(trim(mdl)//'()')
479	1	end subroutine limit_topography
480		! -----------------------------------------------------------------------------
481
482		! -----------------------------------------------------------------------------
483		!> This subroutine sets up the Coriolis parameter for a sphere
484	1	subroutine 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	1	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	1	units="s-1", default=7.2921e-5, scale=US%T_to_s)
502	1	PI = 4.0*atan(1.0)
503
504	9031	do I=G%IsdB,G%IedB ; do J=G%JsdB,G%JedB
505	9030	f(I,J) = ( 2.0 * omega ) * sin( ( PI * G%geoLatBu(I,J) ) / 180.)
506		enddo ; enddo
507
508	1	call callTree_leave(trim(mdl)//'()')
509	1	end subroutine set_rotation_planetary
510		! -----------------------------------------------------------------------------
511
512		! -----------------------------------------------------------------------------
513		!> This subroutine sets up the Coriolis parameter for a beta-plane or f-plane
514	0	subroutine 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	0	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	0	"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	0	"the betaplane option.", units="m-1 s-1", default=0.0, scale=US%T_to_s*US%L_to_m)
540	0	call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees")
541
542	0	PI = 4.0*atan(1.0)
543	0	y_scl = G%grid_unit_to_L
544	0	if (G%grid_unit_to_L <= 0.0) y_scl = PI * G%Rad_Earth_L / 180.
545
546	0	select case (axis_units(1:1))
547		case ("d")
548	0	beta_lat_ref_units = "degrees"
549		case ("k")
550	0	beta_lat_ref_units = "kilometers"
551		case ("m")
552	0	beta_lat_ref_units = "meters"
553		case default ; call MOM_error(FATAL, &
554	0	" 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	0	units=trim(beta_lat_ref_units), default=0.0)
560
561	0	do I=G%IsdB,G%IedB ; do J=G%JsdB,G%JedB
562	0	f(I,J) = f_0 + beta * ( (G%geoLatBu(I,J) - beta_lat_ref) * y_scl )
563		enddo ; enddo
564
565	0	call callTree_leave(trim(mdl)//'()')
566	0	end 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.
570	1	subroutine 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	1	default=.false.)
590
591	1	if (use_bugs) then
592	0	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	0	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	0	G%geoLatBu(I-1,J-1) - G%geoLatBu(I,J-1) )
599	0	G%sin_rot(i,j) = sin(angle) ! angle is the clockwise angle from lat/lon to ocean
600	0	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	0	call pass_var(G%cos_rot, G%Domain)
605	0	call pass_var(G%sin_rot, G%Domain)
606		else
607	1	pi_720deg = atan(1.0) / 180.0
608	1	len_lon = 360.0 ; if (G%len_lon > 0.0) len_lon = G%len_lon
609	7261	do j=G%jsc,G%jec ; do i=G%isc,G%iec
610	50400	do n=1,2 ; do m=1,2
611	43200	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	7200	(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	7200	(G%geoLatBu(I,J) - G%geoLatBu(I-1,J-1)) )
618	7200	G%sin_rot(i,j) = sin(angle) ! angle is the clockwise angle from lat/lon to ocean
619	7260	G%cos_rot(i,j) = cos(angle) ! grid (e.g. angle of ocean "north" from true north)
620		enddo ; enddo
621
622	1	call pass_vector(G%cos_rot, G%sin_rot, G%Domain, stagger=AGRID)
623		endif
624
625	1	end subroutine initialize_grid_rotation_angle
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.
630	28800	function 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	28800	if (Lx > 0.0) then
637	28800	x_mod = modulo(x - (xc - 0.5Lx), Lx) + (xc - 0.5Lx)
638		else
639	0	x_mod = x
640		endif
641	28800	end 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.
646	0	subroutine 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	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
661	0	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
662	0	pi_180 = (4.0*atan(1.0))/180.0
663
664	0	dx_2 = -1.0 ; dy_2 = -1.0
665	0	option = -1
666
667	0	select case ( trim(name) )
668	0	case ("global_1deg") ; option = 1 ; dx_2 = 0.5*1.0
669		case default ; call MOM_error(FATAL, "reset_face_lengths_named: "//&
670	0	"Unrecognized channel configuration name "//trim(name))
671		end select
672
673	0	if (option==1) then ! 1-degree settings.
674	0	do j=jsd,jed ; do I=IsdB,IedB ! Change any u-face lengths within this loop.
675	0	dy_2 = dx_2 * G%dyCu(I,j)G%IdxCu(I,j) cos(pi_180 * G%geoLatCu(I,j))
676
677	0	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	0	G%dy_Cu(I,j) = G%mask2dCu(I,j)12000.0US%m_to_L ! Gibraltar
680
681	0	if ((abs(G%geoLatCu(I,j)-12.5) < dy_2) .and. (abs(G%geoLonCu(I,j)-43.0) < dx_2)) &
682	0	G%dy_Cu(I,j) = G%mask2dCu(I,j)10000.0US%m_to_L ! Red Sea
683
684	0	if ((abs(G%geoLatCu(I,j)-40.5) < dy_2) .and. (abs(G%geoLonCu(I,j)-26.0) < dx_2)) &
685	0	G%dy_Cu(I,j) = G%mask2dCu(I,j)5000.0US%m_to_L ! Dardanelles
686
687	0	if ((abs(G%geoLatCu(I,j)-41.5) < dy_2) .and. (abs(G%geoLonCu(I,j)+220.0) < dx_2)) &
688	0	G%dy_Cu(I,j) = G%mask2dCu(I,j)35000.0US%m_to_L ! Tsugaru strait at 140.0e
689
690	0	if ((abs(G%geoLatCu(I,j)-45.5) < dy_2) .and. (abs(G%geoLonCu(I,j)+217.5) < 0.9)) &
691	0	G%dy_Cu(I,j) = G%mask2dCu(I,j)15000.0US%m_to_L ! Betw Hokkaido and Sakhalin at 217&218 = 142e
692
693		! Greater care needs to be taken in the tripolar region.
694	0	if ((abs(G%geoLatCu(I,j)-80.84) < 0.2) .and. (abs(G%geoLonCu(I,j)+64.9) < 0.8)) &
695	0	G%dy_Cu(I,j) = G%mask2dCu(I,j)38000.0US%m_to_L ! Smith Sound in Canadian Arch - tripolar region
696
697		enddo ; enddo
698
699	0	do J=JsdB,JedB ; do i=isd,ied ! Change any v-face lengths within this loop.
700	0	dy_2 = dx_2 * G%dyCv(i,J)G%IdxCv(i,J) cos(pi_180 * G%geoLatCv(i,J))
701	0	if ((abs(G%geoLatCv(i,J)-41.0) < dy_2) .and. (abs(G%geoLonCv(i,J)-28.5) < dx_2)) &
702	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)2500.0US%m_to_L ! Bosporus - should be 1000.0 m wide.
703
704	0	if ((abs(G%geoLatCv(i,J)-13.0) < dy_2) .and. (abs(G%geoLonCv(i,J)-42.5) < dx_2)) &
705	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)10000.0US%m_to_L ! Red Sea
706
707	0	if ((abs(G%geoLatCv(i,J)+2.8) < 0.8) .and. (abs(G%geoLonCv(i,J)+241.5) < dx_2)) &
708	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)40000.0US%m_to_L ! Makassar Straits at 241.5 W = 118.5 E
709
710	0	if ((abs(G%geoLatCv(i,J)-0.56) < 0.5) .and. (abs(G%geoLonCv(i,J)+240.5) < dx_2)) &
711	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)80000.0US%m_to_L ! entry to Makassar Straits at 240.5 W = 119.5 E
712
713	0	if ((abs(G%geoLatCv(i,J)-0.19) < 0.5) .and. (abs(G%geoLonCv(i,J)+230.5) < dx_2)) &
714	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)25000.0US%m_to_L ! Channel betw N Guinea and Halmahara 230.5 W = 129.5 E
715
716	0	if ((abs(G%geoLatCv(i,J)-0.19) < 0.5) .and. (abs(G%geoLonCv(i,J)+229.5) < dx_2)) &
717	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)25000.0US%m_to_L ! Channel betw N Guinea and Halmahara 229.5 W = 130.5 E
718
719	0	if ((abs(G%geoLatCv(i,J)-0.0) < 0.25) .and. (abs(G%geoLonCv(i,J)+228.5) < dx_2)) &
720	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)25000.0US%m_to_L ! Channel betw N Guinea and Halmahara 228.5 W = 131.5 E
721
722	0	if ((abs(G%geoLatCv(i,J)+8.5) < 0.5) .and. (abs(G%geoLonCv(i,J)+244.5) < dx_2)) &
723	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)20000.0US%m_to_L ! Lombok Straits at 244.5 W = 115.5 E
724
725	0	if ((abs(G%geoLatCv(i,J)+8.5) < 0.5) .and. (abs(G%geoLonCv(i,J)+235.5) < dx_2)) &
726	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)20000.0US%m_to_L ! Timor Straits at 235.5 W = 124.5 E
727
728	0	if ((abs(G%geoLatCv(i,J)-52.5) < dy_2) .and. (abs(G%geoLonCv(i,J)+218.5) < dx_2)) &
729	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)2500.0US%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	0	if ((abs(G%geoLatCv(i,J)-76.8) < 0.06) .and. (abs(G%geoLonCv(i,J)+88.7) < dx_2)) &
733	0	G%dx_Cv(i,J) = G%mask2dCv(i,J)8400.0US%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	0	do j=jsd,jed ; do I=IsdB,IedB
741	0	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	0	US%L_to_mG%dy_Cu(I,j), US%L_to_mG%dyCu(I,j), US%L_to_m*(G%dy_Cu(I,j)-G%dyCu(I,j)), &
745	0	G%geoLonCu(I,j), G%geoLatCu(I,j)
746	0	call MOM_error(FATAL,"reset_face_lengths_named "//mesg)
747		endif
748	0	G%areaCu(I,j) = G%dxCu(I,j) * G%dy_Cu(I,j)
749	0	G%IareaCu(I,j) = 0.0
750	0	if (G%areaCu(I,j) > 0.0) G%IareaCu(I,j) = G%mask2dCu(I,j) / (G%areaCu(I,j))
751		enddo ; enddo
752
753	0	do J=JsdB,JedB ; do i=isd,ied
754	0	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	0	US%L_to_mG%dx_Cv(i,J), US%L_to_mG%dxCv(i,J), US%L_to_m*(G%dx_Cv(i,J)-G%dxCv(i,J)), &
758	0	G%geoLonCv(i,J), G%geoLatCv(i,J)
759
760	0	call MOM_error(FATAL,"reset_face_lengths_named "//mesg)
761		endif
762	0	G%areaCv(i,J) = G%dyCv(i,J) * G%dx_Cv(i,J)
763	0	G%IareaCv(i,J) = 0.0
764	0	if (G%areaCv(i,J) > 0.0) G%IareaCv(i,J) = G%mask2dCv(i,J) / (G%areaCv(i,J))
765		enddo ; enddo
766
767	0	end 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.
773	0	subroutine 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	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
786	0	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
787		! These checks apply regardless of the chosen option.
788
789	0	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	0	default="ocean_geometry.nc")
794	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
795	0	inputdir = slasher(inputdir)
796	0	filename = trim(inputdir)//trim(chan_file)
797	0	call log_param(param_file, mdl, "INPUTDIR/CHANNEL_WIDTH_FILE", filename)
798
799	0	if (is_root_pe()) then ; if (.not.file_exists(filename)) &
800		call MOM_error(FATAL," reset_face_lengths_file: Unable to open "//&
801	0	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	0	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	0	default="dxCvo")
810
811	0	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	0	call pass_vector(G%dy_Cu, G%dx_Cv, G%Domain, To_All+SCALAR_PAIR, CGRID_NE)
813
814	0	do j=jsd,jed ; do I=IsdB,IedB
815	0	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	0	US%L_to_mG%dy_Cu(I,j), US%L_to_mG%dyCu(I,j), US%L_to_m*(G%dy_Cu(I,j)-G%dyCu(I,j)), &
819	0	G%geoLonCu(I,j), G%geoLatCu(I,j)
820	0	call MOM_error(FATAL,"reset_face_lengths_file "//mesg)
821		endif
822	0	G%areaCu(I,j) = G%dxCu(I,j) * G%dy_Cu(I,j)
823	0	G%IareaCu(I,j) = 0.0
824	0	if (G%areaCu(I,j) > 0.0) G%IareaCu(I,j) = G%mask2dCu(I,j) / (G%areaCu(I,j))
825		enddo ; enddo
826
827	0	do J=JsdB,JedB ; do i=isd,ied
828	0	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	0	US%L_to_mG%dx_Cv(i,J), US%L_to_mG%dxCv(i,J), US%L_to_m*(G%dx_Cv(i,J)-G%dxCv(i,J)), &
832	0	G%geoLonCv(i,J), G%geoLatCv(i,J)
833
834	0	call MOM_error(FATAL,"reset_face_lengths_file "//mesg)
835		endif
836	0	G%areaCv(i,J) = G%dyCv(i,J) * G%dx_Cv(i,J)
837	0	G%IareaCv(i,J) = 0.0
838	0	if (G%areaCv(i,J) > 0.0) G%IareaCv(i,J) = G%mask2dCv(i,J) / (G%areaCv(i,J))
839		enddo ; enddo
840
841	0	call callTree_leave(trim(mdl)//'()')
842	0	end 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.
848	0	subroutine 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	0	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	0	u_width, v_width ! The open width of faces [L ~> m]
862		integer, allocatable, dimension(:) :: &
863	0	u_line_no, v_line_no, & ! The line numbers in lines of u- and v-face lines
864	0	u_line_used, v_line_used ! The number of times each u- and v-line is used.
865		real, allocatable, dimension(:) :: &
866	0	Dmin_u, Dmax_u, Davg_u ! Porous barrier monomial fit params [Z ~> m]
867		real, allocatable, dimension(:) :: &
868	0	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	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
886	0	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
887
888	0	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	0	default="MOM_channel_list")
893	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
894	0	inputdir = slasher(inputdir)
895	0	filename = trim(inputdir)//trim(chan_file)
896	0	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	0	"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	0	default=.false.)
904
905	0	if (is_root_pe()) then
906		! Open the input file.
907	0	if (.not.file_exists(filename)) call MOM_error(FATAL, &
908	0	" reset_face_lengths_list: Unable to open "//trim(filename))
909
910		! Find an unused unit number.
911	0	do iounit=10,512
912	0	INQUIRE(iounit,OPENED=unit_in_use) ; if (.not.unit_in_use) exit
913		enddo
914	0	if (iounit >= 512) call MOM_error(FATAL, &
915	0	"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	0	form='FORMATTED', action='READ', position='REWIND', iostat=ios)
920	0	if (ios /= 0) call MOM_error(FATAL, &
921	0	"reset_face_lengths_list: Error opening "//trim(filename))
922
923		! Count the number of u_width and v_width entries.
924	0	call read_face_length_list(iounit, filename, num_lines, lines)
925		else
926	0	num_lines = 0
927		endif
928
929	0	len_lon = 360.0 ; if (G%len_lon > 0.0) len_lon = G%len_lon
930	0	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	0	call broadcast(num_lines, root_PE())
933	0	u_pt = 0 ; v_pt = 0
934	0	if (num_lines > 0) then
935	0	allocate(lines(num_lines))
936
937	0	allocate(u_lat(2,num_lines), source=-1e34)
938	0	allocate(u_lon(2,num_lines), source=-1e34)
939	0	allocate(u_width(num_lines), source=-1e34)
940	0	allocate(u_line_used(num_lines), source=0)
941	0	allocate(u_line_no(num_lines), source=0)
942
943	0	allocate(v_lat(2,num_lines), source=-1e34)
944	0	allocate(v_lon(2,num_lines), source=-1e34)
945	0	allocate(v_width(num_lines), source=-1e34)
946	0	allocate(v_line_used(num_lines), source=0)
947	0	allocate(v_line_no(num_lines), source=0)
948
949	0	allocate(Dmin_u(num_lines), source=0.0)
950	0	allocate(Dmax_u(num_lines), source=0.0)
951	0	allocate(Davg_u(num_lines), source=0.0)
952
953	0	allocate(Dmin_v(num_lines), source=0.0)
954	0	allocate(Dmax_v(num_lines), source=0.0)
955	0	allocate(Davg_v(num_lines), source=0.0)
956
957		! Actually read the lines.
958	0	if (is_root_pe()) then
959	0	call read_face_length_list(iounit, filename, nl_read, lines)
960	0	if (nl_read /= num_lines) &
961		call MOM_error(FATAL, 'reset_face_lengths_list : Found different '// &
962	0	'number of valid lines on second reading of '//trim(filename))
963	0	close(iounit) ; iounit = -1
964		endif
965
966		! Broadcast the lines.
967	0	call broadcast(lines, 120, root_PE())
968
969		! Populate the u_width, etc., data.
970	0	do ln=1,num_lines
971	0	line = lines(ln)
972		! Detect keywords
973	0	found_u = .false. ; found_v = .false.
974	0	found_u_por = .false. ; found_v_por = .false.
975	0	isu = index(uppercase(line), "U_WIDTH") ; if (isu > 0) found_u = .true.
976	0	isv = index(uppercase(line), "V_WIDTH") ; if (isv > 0) found_v = .true.
977	0	isu_por = index(uppercase(line), "U_WIDTH_POR") ; if (isu_por > 0) found_u_por = .true.
978	0	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	0	if (found_u) then
982	0	u_pt = u_pt + 1
983	0	if (found_u_por .eqv. .false.) then
984	0	read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
985	0	elseif (found_u_por) then
986	0	read(line(isu_por+12:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt), &
987	0	Dmin_u(u_pt), Dmax_u(u_pt), Davg_u(u_pt)
988		endif
989	0	u_width(u_pt) = US%m_to_L*u_width(u_pt) ! Rescale units equivalently to scale=US%m_to_L during read.
990	0	Dmin_u(u_pt) = US%m_to_Z*Dmin_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
991	0	Dmax_u(u_pt) = US%m_to_Z*Dmax_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
992	0	Davg_u(u_pt) = US%m_to_Z*Davg_u(u_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
993	0	u_line_no(u_pt) = ln
994	0	if (is_root_PE()) then
995	0	if (check_360) then
996	0	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	0	trim(filename))
1000	0	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	0	trim(filename))
1004		endif
1005	0	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	0	trim(filename))
1009	0	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	0	trim(filename))
1013	0	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	0	trim(filename))
1017	0	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	0	trim(filename))
1021		endif
1022	0	elseif (found_v) then
1023	0	v_pt = v_pt + 1
1024	0	if (found_v_por .eqv. .false.) then
1025	0	read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
1026	0	elseif (found_v_por) then
1027	0	read(line(isv+12:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt), &
1028	0	Dmin_v(v_pt), Dmax_v(v_pt), Davg_v(v_pt)
1029		endif
1030	0	v_width(v_pt) = US%m_to_L*v_width(v_pt) ! Rescale units equivalently to scale=US%m_to_L during read.
1031	0	Dmin_v(v_pt) = US%m_to_Z*Dmin_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1032	0	Dmax_v(v_pt) = US%m_to_Z*Dmax_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1033	0	Davg_v(v_pt) = US%m_to_Z*Davg_v(v_pt) ! Rescale units equivalently to scale=US%m_to_Z during read.
1034	0	v_line_no(v_pt) = ln
1035	0	if (is_root_PE()) then
1036	0	if (check_360) then
1037	0	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	0	trim(filename))
1041	0	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	0	trim(filename))
1045		endif
1046	0	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	0	trim(filename))
1050	0	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	0	trim(filename))
1054	0	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	0	trim(filename))
1058	0	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	0	trim(filename))
1062		endif
1063		endif
1064		enddo
1065
1066		endif
1067
1068	0	do j=jsd,jed ; do I=IsdB,IedB
1069	0	lat = G%geoLatCu(I,j) ; lon = G%geoLonCu(I,j)
1070	0	if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1071	0	else ; lon_p = lon ; lon_m = lon ; endif
1072
1073	0	do npt=1,u_pt
1074	0	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	0	((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then
1078
1079	0	G%dy_Cu(I,j) = G%mask2dCu(I,j) * min(G%dyCu(I,j), max(u_width(npt), 0.0))
1080	0	G%porous_DminU(I,j) = Dmin_u(npt)
1081	0	G%porous_DmaxU(I,j) = Dmax_u(npt)
1082	0	G%porous_DavgU(I,j) = Davg_u(npt)
1083
1084	0	if (j>=G%jsc .and. j<=G%jec .and. I>=G%isc .and. I<=G%iec) then ! Limit messages/checking to compute domain
1085	0	if ( G%mask2dCu(I,j) == 0.0 ) then
1086	0	write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",&
1087	0	u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") so grid metric is unmodified."
1088		else
1089	0	u_line_used(npt) = u_line_used(npt) + 1
1090		write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1091	0	"read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon," (",&
1092	0	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	0	"read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",G%porous_DminU(I,j),&
1095	0	G%porous_DmaxU(I,j), G%porous_DavgU(I,j),")m"
1096		endif
1097		endif
1098		endif
1099		enddo
1100
1101	0	G%areaCu(I,j) = G%dxCu(I,j) * G%dy_Cu(I,j)
1102	0	G%IareaCu(I,j) = 0.0
1103	0	if (G%areaCu(I,j) > 0.0) G%IareaCu(I,j) = G%mask2dCu(I,j) / (G%areaCu(I,j))
1104		enddo ; enddo
1105
1106	0	do J=JsdB,JedB ; do i=isd,ied
1107	0	lat = G%geoLatCv(i,J) ; lon = G%geoLonCv(i,J)
1108	0	if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1109	0	else ; lon_p = lon ; lon_m = lon ; endif
1110
1111	0	do npt=1,v_pt
1112	0	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	0	((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then
1116	0	G%dx_Cv(i,J) = G%mask2dCv(i,J) * min(G%dxCv(i,J), max(v_width(npt), 0.0))
1117	0	G%porous_DminV(i,J) = Dmin_v(npt)
1118	0	G%porous_DmaxV(i,J) = Dmax_v(npt)
1119	0	G%porous_DavgV(i,J) = Davg_v(npt)
1120
1121	0	if (i>=G%isc .and. i<=G%iec .and. J>=G%jsc .and. J<=G%jec) then ! Limit messages/checking to compute domain
1122	0	if ( G%mask2dCv(i,J) == 0.0 ) then
1123	0	write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",&
1124	0	v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") so grid metric is unmodified."
1125		else
1126	0	v_line_used(npt) = v_line_used(npt) + 1
1127		write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1128	0	"read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon," (",&
1129	0	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	0	"read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",G%porous_DminV(i,J),&
1132	0	G%porous_DmaxV(i,J), G%porous_DavgV(i,J),")m"
1133		endif
1134		endif
1135		endif
1136		enddo
1137
1138	0	G%areaCv(i,J) = G%dyCv(i,J) * G%dx_Cv(i,J)
1139	0	G%IareaCv(i,J) = 0.0
1140	0	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	0	unused = 0
1145	0	if (u_pt > 0) call sum_across_PEs(u_line_used, u_pt)
1146	0	if (v_pt > 0) call sum_across_PEs(v_line_used, v_pt)
1147	0	if (is_root_PE()) then
1148	0	unused = 0
1149	0	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	0	trim(lines(u_line_no(npt))) )
1152	0	unused = unused + 1
1153		endif ; enddo
1154	0	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	0	trim(lines(v_line_no(npt))) )
1157	0	unused = unused + 1
1158		endif ; enddo
1159	0	if (fatal_unused_lengths .and. (unused > 0)) call MOM_error(FATAL, &
1160	0	"reset_face_lengths_list causing MOM6 abort due to unused face length lines.")
1161		endif
1162
1163	0	if (num_lines > 0) then
1164	0	deallocate(lines)
1165	0	deallocate(u_line_used, v_line_used, u_line_no, v_line_no)
1166	0	deallocate(u_lat) ; deallocate(u_lon) ; deallocate(u_width)
1167	0	deallocate(v_lat) ; deallocate(v_lon) ; deallocate(v_width)
1168	0	deallocate(Dmin_u) ; deallocate(Dmax_u) ; deallocate(Davg_u)
1169	0	deallocate(Dmin_v) ; deallocate(Dmax_v) ; deallocate(Davg_v)
1170		endif
1171
1172	0	call callTree_leave(trim(mdl)//'()')
1173	0	end 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.
1178	0	subroutine 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	0	num_lines = 0
1192
1193	0	if (iounit <= 0) return
1194	0	rewind(iounit)
1195	0	do while(.true.)
1196	0	read(iounit, '(a)', end=8, err=9) line
1197	0	last = len_trim(line)
1198		! Eliminate either F90 or C comments from the line.
1199	0	icom = index(line(:last), "!") ; if (icom > 0) last = icom-1
1200	0	icom = index(line(:last), "/*") ; if (icom > 0) last = icom-1
1201	0	if (last < 1) cycle
1202
1203		! Detect keywords
1204	0	line_up = uppercase(line)
1205	0	found_u = .false. ; found_v = .false.
1206	0	isu = index(line_up(:last), "U_WIDTH") ; if (isu > 0) found_u = .true.
1207	0	isv = index(line_up(:last), "V_WIDTH") ; if (isv > 0) found_v = .true.
1208
1209	0	if (found_u .and. found_v) call MOM_error(FATAL, &
1210		"read_face_length_list : both U_WIDTH and V_WIDTH found when "//&
1211	0	"reading the line "//trim(line(:last))//" in file "//trim(filename))
1212	0	if (found_u .or. found_v) then
1213	0	num_lines = num_lines + 1
1214	0	if (associated(lines)) then
1215	0	lines(num_lines) = line(1:last)
1216		endif
1217		endif
1218		enddo ! while (.true.)
1219
1220		8 continue
1221	0	return
1222
1223		9 call MOM_error(FATAL, "read_face_length_list : "//&
1224	0	"Error while reading file "//trim(filename))
1225
1226	0	end 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
1239	0	subroutine 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	0	call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
1252
1253	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1254	0	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	0	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	0	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	0	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	0	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	0	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	0	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	0	default="depthv_av")
1278
1279	0	filename = trim(inputdir)//trim(topo_file)
1280	0	call log_param(param_file, mdl, "INPUTDIR/TOPO_AT_VEL_FILE", filename)
1281
1282	0	if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
1283	0	" 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	0	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	0	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	0	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	0	G%porous_DmaxU = -G%porous_DmaxU ; G%porous_DminU = -G%porous_DminU ; G%porous_DavgU = -G%porous_DavgU
1295	0	G%porous_DmaxV = -G%porous_DmaxV ; G%porous_DminV = -G%porous_DminV ; G%porous_DavgV = -G%porous_DavgV
1296
1297	0	call pass_vector(G%porous_DmaxU, G%porous_DmaxV, G%Domain, To_All+SCALAR_PAIR, CGRID_NE)
1298	0	call pass_vector(G%porous_DminU, G%porous_DminV, G%Domain, To_All+SCALAR_PAIR, CGRID_NE)
1299	0	call pass_vector(G%porous_DavgU, G%porous_DavgV, G%Domain, To_All+SCALAR_PAIR, CGRID_NE)
1300
1301	0	call callTree_leave(trim(mdl)//'()')
1302	0	end subroutine set_subgrid_topo_at_vel_from_file
1303		! -----------------------------------------------------------------------------
1304
1305		! -----------------------------------------------------------------------------
1306		!> Set the bathymetry at velocity points to be the maximum of the depths at the
1307		!! neighoring tracer points.
1308	0	subroutine 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	0	do I=G%isd,G%ied-1 ; do j=G%jsd,G%jed
1315	0	G%Dblock_u(I,j) = G%mask2dCu(I,j) * max(G%bathyT(i,j), G%bathyT(i+1,j))
1316	0	G%Dopen_u(I,j) = G%Dblock_u(I,j)
1317		enddo ; enddo
1318	0	do i=G%isd,G%ied ; do J=G%jsd,G%jed-1
1319	0	G%Dblock_v(I,J) = G%mask2dCv(i,J) * max(G%bathyT(i,j), G%bathyT(i,j+1))
1320	0	G%Dopen_v(I,J) = G%Dblock_v(I,J)
1321		enddo ; enddo
1322	0	end 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.
1328	0	subroutine 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	0	do I=G%isd,G%ied-1 ; do j=G%jsd,G%jed
1335	0	G%Dblock_u(I,j) = G%mask2dCu(I,j) * min(G%bathyT(i,j), G%bathyT(i+1,j))
1336	0	G%Dopen_u(I,j) = G%Dblock_u(I,j)
1337		enddo ; enddo
1338	0	do i=G%isd,G%ied ; do J=G%jsd,G%jed-1
1339	0	G%Dblock_v(I,J) = G%mask2dCv(i,J) * min(G%bathyT(i,j), G%bathyT(i,j+1))
1340	0	G%Dopen_v(I,J) = G%Dblock_v(I,J)
1341		enddo ; enddo
1342	0	end 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
1348	1	subroutine 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	2	real, dimension(G%isc:G%iec, G%jsc:G%jec) :: masked_area ! Masked cell areas [L2 ~> m2]
1354		integer :: i, j
1355
1356	7261	masked_area(:,:) = 0.
1357	1	G%areaT_global = 0.0 ; G%IareaT_global = 0.0
1358	7261	do j=G%jsc,G%jec ; do i=G%isc,G%iec
1359	7260	masked_area(i,j) = G%areaT(i,j) * G%mask2dT(i,j)
1360		enddo ; enddo
1361	1	G%areaT_global = reproducing_sum(masked_area, unscale=US%L_to_m**2)
1362
1363	1	if (G%areaT_global == 0.0) &
1364	0	call MOM_error(FATAL, "compute_global_grid_integrals: zero ocean area (check topography?)")
1365
1366	1	G%IareaT_global = 1.0 / G%areaT_global
1367	1	end 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.
1373	0	subroutine 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	0	vars ! Types with metadata about the variables and their staggering
1387		type(MOM_field), dimension(:), allocatable :: &
1388	0	fields ! Opaque types used by MOM_io to store variable metadata information
1389	0	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	0	call callTree_enter('write_ocean_geometry_file()')
1396
1397	0	nFlds = 19 ; if (G%bathymetry_at_vel) nFlds = 23
1398
1399	0	allocate(vars(nFlds))
1400	0	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	0	vars(1) = var_desc("geolatb","degree","latitude at corner (Bu) points",'q','1','1')
1413	0	vars(2) = var_desc("geolonb","degree","longitude at corner (Bu) points",'q','1','1')
1414	0	vars(3) = var_desc("geolat","degree", "latitude at tracer (T) points", 'h','1','1')
1415	0	vars(4) = var_desc("geolon","degree","longitude at tracer (T) points",'h','1','1')
1416	0	vars(5) = var_desc("D","meter","Basin Depth",'h','1','1')
1417	0	vars(6) = var_desc("f","s-1","Coriolis Parameter",'q','1','1')
1418	0	vars(7) = var_desc("dxCv","m","Zonal grid spacing at v points",'v','1','1')
1419	0	vars(8) = var_desc("dyCu","m","Meridional grid spacing at u points",'u','1','1')
1420	0	vars(9) = var_desc("dxCu","m","Zonal grid spacing at u points",'u','1','1')
1421	0	vars(10)= var_desc("dyCv","m","Meridional grid spacing at v points",'v','1','1')
1422	0	vars(11)= var_desc("dxT","m","Zonal grid spacing at h points",'h','1','1')
1423	0	vars(12)= var_desc("dyT","m","Meridional grid spacing at h points",'h','1','1')
1424	0	vars(13)= var_desc("dxBu","m","Zonal grid spacing at q points",'q','1','1')
1425	0	vars(14)= var_desc("dyBu","m","Meridional grid spacing at q points",'q','1','1')
1426	0	vars(15)= var_desc("Ah","m2","Area of h cells",'h','1','1')
1427	0	vars(16)= var_desc("Aq","m2","Area of q cells",'q','1','1')
1428
1429	0	vars(17)= var_desc("dxCvo","m","Open zonal grid spacing at v points",'v','1','1')
1430	0	vars(18)= var_desc("dyCuo","m","Open meridional grid spacing at u points",'u','1','1')
1431	0	vars(19)= var_desc("wet", "nondim", "land or ocean?", 'h','1','1')
1432
1433	0	if (G%bathymetry_at_vel) then
1434	0	vars(20) = var_desc("Dblock_u","m","Blocked depth at u points",'u','1','1')
1435	0	vars(21) = var_desc("Dopen_u","m","Open depth at u points",'u','1','1')
1436	0	vars(22) = var_desc("Dblock_v","m","Blocked depth at v points",'v','1','1')
1437	0	vars(23) = var_desc("Dopen_v","m","Open depth at v points",'v','1','1')
1438		endif
1439
1440	0	if (present(geom_file)) then
1441	0	filepath = trim(directory) // trim(geom_file)
1442		else
1443	0	filepath = trim(directory) // "ocean_geometry"
1444		endif
1445
1446		! Append ensemble run number to filename if it is an ensemble run
1447	0	call get_filename_appendix(filename_appendix)
1448	0	if (len_trim(filename_appendix) > 0) then
1449	0	geom_file_len = len_trim(filepath)
1450	0	if (filepath(geom_file_len-2:geom_file_len) == ".nc") then
1451	0	filepath = filepath(1:geom_file_len-3) // '.' // trim(filename_appendix) // ".nc"
1452		else
1453	0	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	0	default=.false.)
1462	0	file_threading = SINGLE_FILE
1463	0	if (multiple_files) file_threading = MULTIPLE
1464
1465		call create_MOM_file(IO_handle, trim(filepath), vars, nFlds, fields, &
1466	0	file_threading, dG=G)
1467
1468	0	call MOM_write_field(IO_handle, fields(1), G%Domain, G%geoLatBu)
1469	0	call MOM_write_field(IO_handle, fields(2), G%Domain, G%geoLonBu)
1470	0	call MOM_write_field(IO_handle, fields(3), G%Domain, G%geoLatT)
1471	0	call MOM_write_field(IO_handle, fields(4), G%Domain, G%geoLonT)
1472
1473	0	call MOM_write_field(IO_handle, fields(5), G%Domain, G%bathyT, unscale=US%Z_to_m)
1474	0	call MOM_write_field(IO_handle, fields(6), G%Domain, G%CoriolisBu, unscale=US%s_to_T)
1475
1476	0	call MOM_write_field(IO_handle, fields(7), G%Domain, G%dxCv, unscale=US%L_to_m)
1477	0	call MOM_write_field(IO_handle, fields(8), G%Domain, G%dyCu, unscale=US%L_to_m)
1478	0	call MOM_write_field(IO_handle, fields(9), G%Domain, G%dxCu, unscale=US%L_to_m)
1479	0	call MOM_write_field(IO_handle, fields(10), G%Domain, G%dyCv, unscale=US%L_to_m)
1480	0	call MOM_write_field(IO_handle, fields(11), G%Domain, G%dxT, unscale=US%L_to_m)
1481	0	call MOM_write_field(IO_handle, fields(12), G%Domain, G%dyT, unscale=US%L_to_m)
1482	0	call MOM_write_field(IO_handle, fields(13), G%Domain, G%dxBu, unscale=US%L_to_m)
1483	0	call MOM_write_field(IO_handle, fields(14), G%Domain, G%dyBu, unscale=US%L_to_m)
1484
1485	0	call MOM_write_field(IO_handle, fields(15), G%Domain, G%areaT, unscale=US%L_to_m**2)
1486	0	call MOM_write_field(IO_handle, fields(16), G%Domain, G%areaBu, unscale=US%L_to_m**2)
1487
1488	0	call MOM_write_field(IO_handle, fields(17), G%Domain, G%dx_Cv, unscale=US%L_to_m)
1489	0	call MOM_write_field(IO_handle, fields(18), G%Domain, G%dy_Cu, unscale=US%L_to_m)
1490	0	call MOM_write_field(IO_handle, fields(19), G%Domain, G%mask2dT)
1491
1492	0	if (G%bathymetry_at_vel) then
1493	0	call MOM_write_field(IO_handle, fields(20), G%Domain, G%Dblock_u, unscale=US%Z_to_m)
1494	0	call MOM_write_field(IO_handle, fields(21), G%Domain, G%Dopen_u, unscale=US%Z_to_m)
1495	0	call MOM_write_field(IO_handle, fields(22), G%Domain, G%Dblock_v, unscale=US%Z_to_m)
1496	0	call MOM_write_field(IO_handle, fields(23), G%Domain, G%Dopen_v, unscale=US%Z_to_m)
1497		endif
1498
1499	0	call IO_handle%close()
1500
1501	0	deallocate(vars, fields)
1502
1503	0	call callTree_leave('write_ocean_geometry_file()')
1504	0	end subroutine write_ocean_geometry_file
1505
1506		end module MOM_shared_initialization