MOM_grid_initialize.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!> Initializes horizontal grid
7
8use mom_checksums, only : hchksum, bchksum, uvchksum, hchksum_pair, bchksum_pair
9use mom_domains, only : pass_var, pass_vector, pe_here, root_pe, broadcast
10use mom_domains, only : agrid, bgrid_ne, cgrid_ne, to_all, scalar_pair
11use mom_domains, only : to_north, to_south, to_east, to_west
12use mom_domains, only : mom_domain_type, clone_mom_domain, deallocate_mom_domain
14use mom_error_handler, only : mom_error, mom_mesg, fatal, is_root_pe
16use mom_file_parser, only : get_param, log_param, log_version, param_file_type
17use mom_io, only : mom_read_data, slasher, file_exists, stdout
18use mom_io, only : corner, north_face, east_face
20
21implicit none ; private
22
24
25! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
26! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
27! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
28! vary with the Boussinesq approximation, the Boussinesq variant is given first.
29
30!> Global positioning system (aka container for information to describe the grid)
31type, private :: gps ; private
32 real :: len_lon !< The longitudinal or x-direction length of the domain [degrees_E] or [km] or [m].
33 real :: len_lat !< The latitudinal or y-direction length of the domain [degrees_N] or [km] or [m].
34 real :: west_lon !< The western longitude of the domain or the equivalent
35 !! starting value for the x-axis [degrees_E] or [km] or [m].
36 real :: south_lat !< The southern latitude of the domain or the equivalent
37 !! starting value for the y-axis [degrees_N] or [km] or [m].
38 real :: rad_earth_l !< The radius of the Earth in rescaled units [L ~> m]
39 real :: lat_enhance_factor !< The amount by which the meridional resolution
40 !! is enhanced within LAT_EQ_ENHANCE of the equator [nondim]
41 real :: lat_eq_enhance !< The latitude range to the north and south of the equator
42 !! over which the resolution is enhanced [degrees_N]
43 logical :: isotropic !< If true, an isotropic grid on a sphere (also known as a Mercator grid)
44 !! is used. With an isotropic grid, the meridional extent of the domain
45 !! (LENLAT), the zonal extent (LENLON), and the number of grid points in each
46 !! direction are _not_ independent. In MOM the meridional extent is determined
47 !! to fit the zonal extent and the number of grid points, while grid is
48 !! perfectly isotropic.
49 logical :: equator_reference !< If true, the grid is defined to have the equator at the
50 !! nearest q- or h- grid point to (-LOWLAT*NJGLOBAL/LENLAT).
51 integer :: niglobal !< The number of i-points in the global grid computational domain
52 integer :: njglobal !< The number of j-points in the global grid computational domain
53end type gps
54
55contains
56
57!> set_grid_metrics is used to set the primary values in the model's horizontal
58!! grid. The bathymetry, land-sea mask and any restricted channel widths are
59!! not known yet, so these are set later.
60subroutine set_grid_metrics(G, param_file, US)
61 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
62 type(param_file_type), intent(in) :: param_file !< Parameter file structure
63 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
64
65 ! Local variables
66 ! This include declares and sets the variable "version".
67# include "version_variable.h"
68 logical :: debug
69 character(len=256) :: config
70
71 call calltree_enter("set_grid_metrics(), MOM_grid_initialize.F90")
72 call log_version(param_file, "MOM_grid_init", version, "")
73 call get_param(param_file, "MOM_grid_init", "GRID_CONFIG", config, &
74 "A character string that determines the method for "//&
75 "defining the horizontal grid. Current options are: \n"//&
76 " \t mosaic - read the grid from a mosaic (supergrid) \n"//&
77 " \t file set by GRID_FILE.\n"//&
78 " \t cartesian - use a (flat) Cartesian grid.\n"//&
79 " \t spherical - use a simple spherical grid.\n"//&
80 " \t mercator - use a Mercator spherical grid.", &
81 fail_if_missing=.true.)
82 call get_param(param_file, "MOM_grid_init", "DEBUG", debug, &
83 "If true, write out verbose debugging data.", &
84 default=.false., debuggingparam=.true.)
85
86 ! These are defaults that may be changed in the next select block.
87 g%x_axis_units = "degrees_east" ; g%y_axis_units = "degrees_north"
88 g%x_ax_unit_short = "degrees_E" ; g%y_ax_unit_short = "degrees_N"
89 g%grid_unit_to_L = 0.0
90 g%Rad_Earth_L = -1.0*us%m_to_L ; g%len_lat = 0.0 ; g%len_lon = 0.0
91 select case (trim(config))
92 case ("mosaic"); call set_grid_metrics_from_mosaic(g, param_file, us)
93 case ("cartesian"); call set_grid_metrics_cartesian(g, param_file, us)
94 case ("spherical"); call set_grid_metrics_spherical(g, param_file, us)
95 case ("mercator"); call set_grid_metrics_mercator(g, param_file, us)
96 case ("file"); call mom_error(fatal, "MOM_grid_init: set_grid_metrics "//&
97 'GRID_CONFIG "file" is no longer a supported option. Use a '//&
98 'mosaic file ("mosaic") or one of the analytic forms instead.')
99 case default ; call mom_error(fatal, "MOM_grid_init: set_grid_metrics "//&
100 "Unrecognized grid configuration "//trim(config))
101 end select
102 if (g%Rad_Earth_L <= 0.0) then
103 ! The grid metrics were set with an option that does not explicitly initialize Rad_Earth.
104 call get_param(param_file, "MOM_grid_init", "RAD_EARTH", g%Rad_Earth_L, &
105 "The radius of the Earth.", units="m", default=6.378e6, scale=us%m_to_L)
106 endif
107
108 ! Calculate derived metrics (i.e. reciprocals and products)
109 call calltree_enter("set_derived_metrics(), MOM_grid_initialize.F90")
110 call set_derived_dyn_horgrid(g, us)
111 call calltree_leave("set_derived_metrics()")
112
113 if (debug) call grid_metrics_chksum('MOM_grid_init/set_grid_metrics', g, us)
114
115 call calltree_leave("set_grid_metrics()")
116end subroutine set_grid_metrics
117
118! ------------------------------------------------------------------------------
119
120!> grid_metrics_chksum performs a set of checksums on metrics on the grid for
121!! debugging.
122subroutine grid_metrics_chksum(parent, G, US)
123 character(len=*), intent(in) :: parent !< A string identifying the caller
124 type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
125 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
126
127 integer :: halo
128
129 halo = min(g%ied-g%iec, g%jed-g%jec, 1)
130
131 call hchksum_pair(trim(parent)//': d[xy]T', g%dxT, g%dyT, g%HI, &
132 haloshift=halo, unscale=us%L_to_m, scalar_pair=.true.)
133
134 call uvchksum(trim(parent)//': dxC[uv]', g%dxCu, g%dyCv, g%HI, haloshift=halo, unscale=us%L_to_m)
135
136 call uvchksum(trim(parent)//': dxC[uv]', g%dyCu, g%dxCv, g%HI, haloshift=halo, unscale=us%L_to_m)
137
138 call bchksum_pair(trim(parent)//': dxB[uv]', g%dxBu, g%dyBu, g%HI, haloshift=halo, unscale=us%L_to_m)
139
140 call hchksum_pair(trim(parent)//': Id[xy]T', g%IdxT, g%IdyT, g%HI, &
141 haloshift=halo, unscale=us%m_to_L, scalar_pair=.true.)
142
143 call uvchksum(trim(parent)//': Id[xy]C[uv]', g%IdxCu, g%IdyCv, g%HI, haloshift=halo, unscale=us%m_to_L)
144
145 call uvchksum(trim(parent)//': Id[xy]C[uv]', g%IdyCu, g%IdxCv, g%HI, haloshift=halo, unscale=us%m_to_L)
146
147 call bchksum_pair(trim(parent)//': Id[xy]B[uv]', g%IdxBu, g%IdyBu, g%HI, haloshift=halo, unscale=us%m_to_L)
148
149 call hchksum(g%areaT, trim(parent)//': areaT',g%HI, haloshift=halo, unscale=us%L_to_m**2)
150 call bchksum(g%areaBu, trim(parent)//': areaBu',g%HI, haloshift=halo, unscale=us%L_to_m**2)
151
152 call hchksum(g%IareaT, trim(parent)//': IareaT',g%HI, haloshift=halo, unscale=us%m_to_L**2)
153 call bchksum(g%IareaBu, trim(parent)//': IareaBu',g%HI, haloshift=halo, unscale=us%m_to_L**2)
154
155 call hchksum(g%geoLonT,trim(parent)//': geoLonT',g%HI, haloshift=halo)
156 call hchksum(g%geoLatT,trim(parent)//': geoLatT',g%HI, haloshift=halo)
157
158 call bchksum(g%geoLonBu, trim(parent)//': geoLonBu',g%HI, haloshift=halo)
159 call bchksum(g%geoLatBu, trim(parent)//': geoLatBu',g%HI, haloshift=halo)
160
161 call uvchksum(trim(parent)//': geoLonC[uv]', g%geoLonCu, g%geoLonCv, g%HI, haloshift=halo)
162
163 call uvchksum(trim(parent)//': geoLatC[uv]', g%geoLatCu, g%geoLatCv, g%HI, haloshift=halo)
164
165end subroutine grid_metrics_chksum
166
167! ------------------------------------------------------------------------------
168
169!> Sets the grid metrics from a mosaic file.
170subroutine set_grid_metrics_from_mosaic(G, param_file, US)
171 type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
172 type(param_file_type), intent(in) :: param_file !< Parameter file structure
173 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
174
175 ! Local variables
176 ! These are symmetric arrays, corresponding to the data in the mosaic file
177 real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpT ! Areas [L2 ~> m2]
178 real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpU ! East face supergrid spacing [L ~> m]
179 real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpV ! North face supergrid spacing [L ~> m]
180 real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpZ ! Corner latitudes [degrees_N] or
181 ! longitudes [degrees_E]
182 real, dimension(:,:), allocatable :: tmpGlbl ! A global array of axis labels [degrees_N] or [km] or [m]
183 character(len=200) :: filename, grid_file, inputdir
184 character(len=64) :: mdl = "MOM_grid_init set_grid_metrics_from_mosaic"
185 type(mom_domain_type), pointer :: SGdom => null() ! Supergrid domain
186 logical :: lon_bug ! If true use an older buggy answer in the tripolar longitude.
187 integer :: i, j, i2, j2, ni, nj
188 integer :: start(4), nread(4)
189
190 call calltree_enter("set_grid_metrics_from_mosaic(), MOM_grid_initialize.F90")
191
192 call get_param(param_file, mdl, "GRID_FILE", grid_file, &
193 "Name of the file from which to read horizontal grid data.", &
194 fail_if_missing=.true.)
195 call get_param(param_file, mdl, "USE_TRIPOLAR_GEOLONB_BUG", lon_bug, &
196 "If true, use older code that incorrectly sets the longitude "//&
197 "in some points along the tripolar fold to be off by 360 degrees.", &
198 default=.false.)
199 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
200 inputdir = slasher(inputdir)
201 filename = trim(adjustl(inputdir)) // trim(adjustl(grid_file))
202 call log_param(param_file, mdl, "INPUTDIR/GRID_FILE", filename)
203 if (.not.file_exists(filename)) &
204 call mom_error(fatal," set_grid_metrics_from_mosaic: Unable to open "//&
205 trim(filename))
206
207 !<MISSING CODE TO READ REFINEMENT LEVEL>
208
209 call clone_mom_domain(g%domain, sgdom, symmetric=.true., domain_name="MOM_MOSAIC", &
210 refine=2, extra_halo=1)
211
212 ! Read X from the supergrid
213 tmpz(:,:) = 999.
214 call mom_read_data(filename, 'x', tmpz, sgdom, position=corner)
215
216 if (lon_bug) then
217 call pass_var(tmpz, sgdom, position=corner)
218 else
219 call pass_var(tmpz, sgdom, position=corner, inner_halo=0)
220 endif
221 call extrapolate_metric(tmpz, 2*(g%jsc-g%jsd)+2, missing=999.)
222 do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
223 g%geoLonT(i,j) = tmpz(i2-1,j2-1)
224 enddo ; enddo
225 do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
226 g%geoLonBu(i,j) = tmpz(i2,j2)
227 enddo ; enddo
228 do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
229 g%geoLonCu(i,j) = tmpz(i2,j2-1)
230 enddo ; enddo
231 do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
232 g%geoLonCv(i,j) = tmpz(i2-1,j2)
233 enddo ; enddo
234 ! For some reason, this messes up the solution...
235 ! call pass_var(G%geoLonBu, G%domain, position=CORNER)
236
237 ! Read Y from the supergrid
238 tmpz(:,:) = 999.
239 call mom_read_data(filename, 'y', tmpz, sgdom, position=corner)
240
241 call pass_var(tmpz, sgdom, position=corner)
242 call extrapolate_metric(tmpz, 2*(g%jsc-g%jsd)+2, missing=999.)
243 do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
244 g%geoLatT(i,j) = tmpz(i2-1,j2-1)
245 enddo ; enddo
246 do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
247 g%geoLatBu(i,j) = tmpz(i2,j2)
248 enddo ; enddo
249 do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
250 g%geoLatCu(i,j) = tmpz(i2,j2-1)
251 enddo ; enddo
252 do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
253 g%geoLatCv(i,j) = tmpz(i2-1,j2)
254 enddo ; enddo
255
256 ! This routine could be modified to support the use of a mosaic using Cartesian grid coordinates,
257 ! in which case the values of G%x_axis_units, G%y_axis_units and G%grid_unit_to_L would need to be
258 ! reset appropriately here, but this option has not yet been implemented, and the grid coordinates
259 ! are assumed to be degrees of longitude and latitude.
260
261 ! Read DX,DY from the supergrid
262 tmpu(:,:) = 0. ; tmpv(:,:) = 0.
263 call mom_read_data(filename, 'dx', tmpv, sgdom, position=north_face, scale=us%m_to_L)
264 call mom_read_data(filename, 'dy', tmpu, sgdom, position=east_face, scale=us%m_to_L)
265 call pass_vector(tmpu, tmpv, sgdom, to_all+scalar_pair, cgrid_ne)
266 call extrapolate_metric(tmpv, 2*(g%jsc-g%jsd)+2, missing=0.)
267 call extrapolate_metric(tmpu, 2*(g%jsc-g%jsd)+2, missing=0.)
268
269 do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
270 g%dxT(i,j) = tmpv(i2-1,j2-1) + tmpv(i2,j2-1)
271 g%dyT(i,j) = tmpu(i2-1,j2-1) + tmpu(i2-1,j2)
272 enddo ; enddo
273
274 do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
275 g%dxCu(i,j) = tmpv(i2,j2-1) + tmpv(i2+1,j2-1)
276 g%dyCu(i,j) = tmpu(i2,j2-1) + tmpu(i2,j2)
277 enddo ; enddo
278
279 do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
280 g%dxCv(i,j) = tmpv(i2-1,j2) + tmpv(i2,j2)
281 g%dyCv(i,j) = tmpu(i2-1,j2) + tmpu(i2-1,j2+1)
282 enddo ; enddo
283
284 do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
285 g%dxBu(i,j) = tmpv(i2,j2) + tmpv(i2+1,j2)
286 g%dyBu(i,j) = tmpu(i2,j2) + tmpu(i2,j2+1)
287 enddo ; enddo
288
289 ! Read AREA from the supergrid
290 tmpt(:,:) = 0.
291 call mom_read_data(filename, 'area', tmpt, sgdom, scale=us%m_to_L**2)
292 call pass_var(tmpt, sgdom)
293 call extrapolate_metric(tmpt, 2*(g%jsc-g%jsd)+2, missing=0.)
294
295 do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
296 g%areaT(i,j) = (tmpt(i2-1,j2-1) + tmpt(i2,j2)) + &
297 (tmpt(i2-1,j2) + tmpt(i2,j2-1))
298 enddo ; enddo
299 do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
300 g%areaBu(i,j) = (tmpt(i2,j2) + tmpt(i2+1,j2+1)) + &
301 (tmpt(i2,j2+1) + tmpt(i2+1,j2))
302 enddo ; enddo
303
304 ni = sgdom%niglobal
305 nj = sgdom%njglobal
306 call deallocate_mom_domain(sgdom)
307
308 call pass_vector(g%dyCu, g%dxCv, g%Domain, to_all+scalar_pair, cgrid_ne)
309 call pass_vector(g%dxCu, g%dyCv, g%Domain, to_all+scalar_pair, cgrid_ne)
310 call pass_vector(g%dxBu, g%dyBu, g%Domain, to_all+scalar_pair, bgrid_ne)
311 call pass_var(g%areaT, g%Domain)
312 call pass_var(g%areaBu, g%Domain, position=corner)
313
314 ! Construct axes for diagnostic output (only necessary because "ferret" uses
315 ! broken convention for interpretting netCDF files).
316 start(:) = 1 ; nread(:) = 1
317 start(2) = 2 ; nread(1) = ni+1 ; nread(2) = 2
318 allocate( tmpglbl(ni+1,2) )
319 if (is_root_pe()) &
320 call mom_read_data(filename, "x", tmpglbl, start, nread, &
321 no_domain=.true., turns=g%HI%turns)
322 call broadcast(tmpglbl, 2*(ni+1), root_pe())
323
324 ! I don't know why the second axis is 1 or 2 here. -RWH
325 do i=g%isg,g%ieg
326 g%gridLonT(i) = tmpglbl(2*(i-g%isg)+2,2)
327 enddo
328 ! Note that the dynamic grid always uses symmetric memory for the global
329 ! arrays G%gridLatB and G%gridLonB.
330 do i=g%isg-1,g%ieg
331 g%gridLonB(i) = tmpglbl(2*(i-g%isg)+3,1)
332 enddo
333 deallocate( tmpglbl )
334
335 allocate( tmpglbl(1, nj+1) )
336 start(:) = 1 ; nread(:) = 1
337 start(1) = int(ni/4)+1 ; nread(2) = nj+1
338 if (is_root_pe()) &
339 call mom_read_data(filename, "y", tmpglbl, start, nread, &
340 no_domain=.true., turns=g%HI%turns)
341 call broadcast(tmpglbl, nj+1, root_pe())
342
343 do j=g%jsg,g%jeg
344 g%gridLatT(j) = tmpglbl(1,2*(j-g%jsg)+2)
345 enddo
346 do j=g%jsg-1,g%jeg
347 g%gridLatB(j) = tmpglbl(1,2*(j-g%jsg)+3)
348 enddo
349 deallocate( tmpglbl )
350
351 call calltree_leave("set_grid_metrics_from_mosaic()")
352end subroutine set_grid_metrics_from_mosaic
353
354
355! ------------------------------------------------------------------------------
356
357!> Calculate the values of the metric terms for a Cartesian grid that
358!! might be used and save them in arrays.
359!!
360!! Within this subroutine, the x- and y- grid spacings and their
361!! inverses and the cell areas centered on h, q, u, and v points are
362!! calculated, as are the geographic locations of each of these 4
363!! sets of points.
364subroutine set_grid_metrics_cartesian(G, param_file, US)
365 type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
366 type(param_file_type), intent(in) :: param_file !< Parameter file structure
367 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
368 ! Local variables
369 integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, I1off, J1off
370 integer :: niglobal, njglobal
371 real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) ! Axis labels [degrees_N] or [km] or [m]
372 real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) ! Axis labels [degrees_E] or [km] or [m]
373 real :: dx_everywhere, dy_everywhere ! Grid spacings [L ~> m].
374 real :: I_dx, I_dy ! Inverse grid spacings [L-1 ~> m-1].
375 real :: PI ! The ratio of the circumference of a circle to its diameter [nondim]
376 character(len=80) :: units_temp
377 character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_cartesian"
378
379 niglobal = g%Domain%niglobal ; njglobal = g%Domain%njglobal
380 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
381 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
382 i1off = g%idg_offset ; j1off = g%jdg_offset
383
384 call calltree_enter("set_grid_metrics_cartesian(), MOM_grid_initialize.F90")
385
386 pi = 4.0*atan(1.0)
387
388 call get_param(param_file, mdl, "AXIS_UNITS", units_temp, &
389 "The units for the Cartesian axes. Valid entries are: \n"//&
390 " \t degrees - degrees of latitude and longitude \n"//&
391 " \t m or meter(s) - meters \n"//&
392 " \t k or km or kilometer(s) - kilometers", default="degrees")
393 if (trim(units_temp) == "k") units_temp = "km"
394
395 call get_param(param_file, mdl, "SOUTHLAT", g%south_lat, &
396 "The southern latitude of the domain or the equivalent "//&
397 "starting value for the y-axis.", units=units_temp, &
398 fail_if_missing=.true.)
399 call get_param(param_file, mdl, "LENLAT", g%len_lat, &
400 "The latitudinal or y-direction length of the domain.", &
401 units=units_temp, fail_if_missing=.true.)
402 call get_param(param_file, mdl, "WESTLON", g%west_lon, &
403 "The western longitude of the domain or the equivalent "//&
404 "starting value for the x-axis.", units=units_temp, &
405 default=0.0)
406 call get_param(param_file, mdl, "LENLON", g%len_lon, &
407 "The longitudinal or x-direction length of the domain.", &
408 units=units_temp, fail_if_missing=.true.)
409 call get_param(param_file, mdl, "RAD_EARTH", g%Rad_Earth_L, &
410 "The radius of the Earth.", units="m", default=6.378e6, scale=us%m_to_L)
411
412 if (units_temp(1:1) == 'k') then
413 g%x_axis_units = "kilometers" ; g%y_axis_units = "kilometers"
414 g%x_ax_unit_short = "km" ; g%y_ax_unit_short = "km"
415 elseif (units_temp(1:1) == 'm') then
416 g%x_axis_units = "meters" ; g%y_axis_units = "meters"
417 g%x_ax_unit_short = "m" ; g%y_ax_unit_short = "m"
418 endif
419 call log_param(param_file, mdl, "explicit AXIS_UNITS", g%x_axis_units)
420
421 ! Note that the dynamic grid always uses symmetric memory for the global
422 ! arrays G%gridLatB and G%gridLonB.
423 do j=g%jsg-1,g%jeg
424 g%gridLatB(j) = g%south_lat + g%len_lat*real(j-(g%jsg-1))/real(njglobal)
425 enddo
426 do j=g%jsg,g%jeg
427 g%gridLatT(j) = g%south_lat + g%len_lat*(real(j-g%jsg)+0.5)/real(njglobal)
428 enddo
429 do i=g%isg-1,g%ieg
430 g%gridLonB(i) = g%west_lon + g%len_lon*real(i-(g%isg-1))/real(niglobal)
431 enddo
432 do i=g%isg,g%ieg
433 g%gridLonT(i) = g%west_lon + g%len_lon*(real(i-g%isg)+0.5)/real(niglobal)
434 enddo
435
436 do j=jsdb,jedb
437 grid_latb(j) = g%south_lat + g%len_lat*real(j+j1off-(g%jsg-1))/real(njglobal)
438 enddo
439 do j=jsd,jed
440 grid_latt(j) = g%south_lat + g%len_lat*(real(j+j1off-g%jsg)+0.5)/real(njglobal)
441 enddo
442 do i=isdb,iedb
443 grid_lonb(i) = g%west_lon + g%len_lon*real(i+i1off-(g%isg-1))/real(niglobal)
444 enddo
445 do i=isd,ied
446 grid_lont(i) = g%west_lon + g%len_lon*(real(i+i1off-g%isg)+0.5)/real(niglobal)
447 enddo
448
449 if (units_temp(1:1) == 'k') then ! Axes are measured in km.
450 g%grid_unit_to_L = 1000.0*us%m_to_L
451 dx_everywhere = 1000.0*us%m_to_L * g%len_lon / (real(niglobal))
452 dy_everywhere = 1000.0*us%m_to_L * g%len_lat / (real(njglobal))
453 elseif (units_temp(1:1) == 'm') then ! Axes are measured in m.
454 g%grid_unit_to_L = us%m_to_L
455 dx_everywhere = us%m_to_L*g%len_lon / (real(niglobal))
456 dy_everywhere = us%m_to_L*g%len_lat / (real(njglobal))
457 else ! Axes are measured in degrees of latitude and longitude.
458 dx_everywhere = g%Rad_Earth_L * g%len_lon * pi / (180.0 * niglobal)
459 dy_everywhere = g%Rad_Earth_L * g%len_lat * pi / (180.0 * njglobal)
460 endif
461
462 i_dx = 1.0 / dx_everywhere ; i_dy = 1.0 / dy_everywhere
463
464 do j=jsdb,jedb ; do i=isdb,iedb
465 g%geoLonBu(i,j) = grid_lonb(i) ; g%geoLatBu(i,j) = grid_latb(j)
466
467 g%dxBu(i,j) = dx_everywhere ; g%IdxBu(i,j) = i_dx
468 g%dyBu(i,j) = dy_everywhere ; g%IdyBu(i,j) = i_dy
469 g%areaBu(i,j) = dx_everywhere * dy_everywhere ; g%IareaBu(i,j) = i_dx * i_dy
470 enddo ; enddo
471
472 do j=jsd,jed ; do i=isd,ied
473 g%geoLonT(i,j) = grid_lont(i) ; g%geoLatT(i,j) = grid_latt(j)
474 g%dxT(i,j) = dx_everywhere ; g%IdxT(i,j) = i_dx
475 g%dyT(i,j) = dy_everywhere ; g%IdyT(i,j) = i_dy
476 g%areaT(i,j) = dx_everywhere * dy_everywhere ; g%IareaT(i,j) = i_dx * i_dy
477 enddo ; enddo
478
479 do j=jsd,jed ; do i=isdb,iedb
480 g%geoLonCu(i,j) = grid_lonb(i) ; g%geoLatCu(i,j) = grid_latt(j)
481
482 g%dxCu(i,j) = dx_everywhere ; g%IdxCu(i,j) = i_dx
483 g%dyCu(i,j) = dy_everywhere ; g%IdyCu(i,j) = i_dy
484 enddo ; enddo
485
486 do j=jsdb,jedb ; do i=isd,ied
487 g%geoLonCv(i,j) = grid_lont(i) ; g%geoLatCv(i,j) = grid_latb(j)
488
489 g%dxCv(i,j) = dx_everywhere ; g%IdxCv(i,j) = i_dx
490 g%dyCv(i,j) = dy_everywhere ; g%IdyCv(i,j) = i_dy
491 enddo ; enddo
492
493 call calltree_leave("set_grid_metrics_cartesian()")
494end subroutine set_grid_metrics_cartesian
495
496! ------------------------------------------------------------------------------
497
498!> Calculate the values of the metric terms that might be used
499!! and save them in arrays.
500!!
501!! Within this subroutine, the x- and y- grid spacings and their
502!! inverses and the cell areas centered on h, q, u, and v points are
503!! calculated, as are the geographic locations of each of these 4
504!! sets of points.
505subroutine set_grid_metrics_spherical(G, param_file, US)
506 type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
507 type(param_file_type), intent(in) :: param_file !< Parameter file structure
508 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
509 ! Local variables
510 real :: PI ! PI = 3.1415926... as 4*atan(1) [nondim]
511 real :: PI_180 ! The conversion factor from degrees to radians [radians degree-1]
512 integer :: i, j, isd, ied, jsd, jed
513 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
514 integer :: i_offset, j_offset
515 real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) ! Axis labels [degrees_N]
516 real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) ! Axis labels [degrees_E]
517 real :: dLon ! The change in longitude between successive grid points [degrees_E]
518 real :: dLat ! The change in latitude between successive grid points [degrees_N]
519 real :: dL_di ! dLon rescaled from degrees to radians [radians]
520 real :: latitude ! The latitude of a grid point [degrees_N]
521 character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_spherical"
522
523 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
524 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
525 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
526 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
527 i_offset = g%idg_offset ; j_offset = g%jdg_offset
528
529 call calltree_enter("set_grid_metrics_spherical(), MOM_grid_initialize.F90")
530
531! Calculate the values of the metric terms that might be used
532! and save them in arrays.
533 pi = 4.0*atan(1.0) ; pi_180 = atan(1.0)/45.
534
535 call get_param(param_file, mdl, "SOUTHLAT", g%south_lat, &
536 "The southern latitude of the domain.", units="degrees_N", &
537 fail_if_missing=.true.)
538 call get_param(param_file, mdl, "LENLAT", g%len_lat, &
539 "The latitudinal length of the domain.", units="degrees_N", &
540 fail_if_missing=.true.)
541 call get_param(param_file, mdl, "WESTLON", g%west_lon, &
542 "The western longitude of the domain.", units="degrees_E", &
543 default=0.0)
544 call get_param(param_file, mdl, "LENLON", g%len_lon, &
545 "The longitudinal length of the domain.", units="degrees_E", &
546 fail_if_missing=.true.)
547 call get_param(param_file, mdl, "RAD_EARTH", g%Rad_Earth_L, &
548 "The radius of the Earth.", units="m", default=6.378e6, scale=us%m_to_L)
549
550 dlon = g%len_lon/g%Domain%niglobal
551 dlat = g%len_lat/g%Domain%njglobal
552
553 ! Note that the dynamic grid always uses symmetric memory for the global
554 ! arrays G%gridLatB and G%gridLonB.
555 do j=g%jsg-1,g%jeg
556 latitude = g%south_lat + dlat*(real(j-(g%jsg-1)))
557 g%gridLatB(j) = min(max(latitude,-90.),90.)
558 enddo
559 do j=g%jsg,g%jeg
560 latitude = g%south_lat + dlat*(real(j-g%jsg)+0.5)
561 g%gridLatT(j) = min(max(latitude,-90.),90.)
562 enddo
563 do i=g%isg-1,g%ieg
564 g%gridLonB(i) = g%west_lon + dlon*(real(i-(g%isg-1)))
565 enddo
566 do i=g%isg,g%ieg
567 g%gridLonT(i) = g%west_lon + dlon*(real(i-g%isg)+0.5)
568 enddo
569
570 do j=jsdb,jedb
571 latitude = g%south_lat + dlat* real(j+j_offset-(g%jsg-1))
572 grid_latb(j) = min(max(latitude,-90.),90.)
573 enddo
574 do j=jsd,jed
575 latitude = g%south_lat + dlat*(real(j+j_offset-g%jsg)+0.5)
576 grid_latt(j) = min(max(latitude,-90.),90.)
577 enddo
578 do i=isdb,iedb
579 grid_lonb(i) = g%west_lon + dlon*real(i+i_offset-(g%isg-1))
580 enddo
581 do i=isd,ied
582 grid_lont(i) = g%west_lon + dlon*(real(i+i_offset-g%isg)+0.5)
583 enddo
584
585 dl_di = (g%len_lon * 4.0*atan(1.0)) / (180.0 * g%Domain%niglobal)
586 do j=jsdb,jedb ; do i=isdb,iedb
587 g%geoLonBu(i,j) = grid_lonb(i)
588 g%geoLatBu(i,j) = grid_latb(j)
589
590 ! The following line is needed to reproduce the solution from
591 ! set_grid_metrics_mercator when used to generate a simple spherical grid.
592 g%dxBu(i,j) = g%Rad_Earth_L * cos( g%geoLatBu(i,j)*pi_180 ) * dl_di
593! G%dxBu(I,J) = G%Rad_Earth_L * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 )
594 g%dyBu(i,j) = g%Rad_Earth_L * dlat*pi_180
595 g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
596 enddo ; enddo
597
598 do j=jsdb,jedb ; do i=isd,ied
599 g%geoLonCv(i,j) = grid_lont(i)
600 g%geoLatCv(i,j) = grid_latb(j)
601
602 ! The following line is needed to reproduce the solution from
603 ! set_grid_metrics_mercator when used to generate a simple spherical grid.
604 g%dxCv(i,j) = g%Rad_Earth_L * cos( g%geoLatCv(i,j)*pi_180 ) * dl_di
605! G%dxCv(i,J) = G%Rad_Earth_L * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 )
606 g%dyCv(i,j) = g%Rad_Earth_L * dlat*pi_180
607 enddo ; enddo
608
609 do j=jsd,jed ; do i=isdb,iedb
610 g%geoLonCu(i,j) = grid_lonb(i)
611 g%geoLatCu(i,j) = grid_latt(j)
612
613 ! The following line is needed to reproduce the solution from
614 ! set_grid_metrics_mercator when used to generate a simple spherical grid.
615 g%dxCu(i,j) = g%Rad_Earth_L * cos( g%geoLatCu(i,j)*pi_180 ) * dl_di
616! G%dxCu(I,j) = G%Rad_Earth_L * dLon*PI_180 * COS( latitude )
617 g%dyCu(i,j) = g%Rad_Earth_L * dlat*pi_180
618 enddo ; enddo
619
620 do j=jsd,jed ; do i=isd,ied
621 g%geoLonT(i,j) = grid_lont(i)
622 g%geoLatT(i,j) = grid_latt(j)
623
624 ! The following line is needed to reproduce the solution from
625 ! set_grid_metrics_mercator when used to generate a simple spherical grid.
626 g%dxT(i,j) = g%Rad_Earth_L * cos( g%geoLatT(i,j)*pi_180 ) * dl_di
627! G%dxT(i,j) = G%Rad_Earth_L * dLon*PI_180 * COS( latitude )
628 g%dyT(i,j) = g%Rad_Earth_L * dlat*pi_180
629
630! latitude = G%geoLatCv(i,J)*PI_180 ! In radians
631! dL_di = G%geoLatCv(i,max(jsd,J-1))*PI_180 ! In radians
632! G%areaT(i,j) = Rad_Earth_L**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di))
633 g%areaT(i,j) = g%dxT(i,j) * g%dyT(i,j)
634 enddo ; enddo
635
636 call calltree_leave("set_grid_metrics_spherical()")
637end subroutine set_grid_metrics_spherical
638
639!> Calculate the values of the metric terms that might be used
640!! and save them in arrays.
641!!
642!! Within this subroutine, the x- and y- grid spacings and their
643!! inverses and the cell areas centered on h, q, u, and v points are
644!! calculated, as are the geographic locations of each of these 4
645!! sets of points.
646subroutine set_grid_metrics_mercator(G, param_file, US)
647 type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
648 type(param_file_type), intent(in) :: param_file !< Parameter file structure
649 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
650 ! Local variables
651 integer :: i, j, isd, ied, jsd, jed
652 integer :: I_off, J_off
653 type(gps) :: GP
654 character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_mercator"
655 real :: PI, PI_2 ! PI = 3.1415926... as 4*atan(1), PI_2 = (PI) /2.0 [nondim]
656 real :: y_q, y_h ! Latitudes of a point [radians]
657 real :: id ! The i-grid space positions whose longitude is being sought [gridpoints]
658 real :: jd ! The j-grid space positions whose latitude is being sought [gridpoints]
659 real :: x_q, x_h ! Longitudes of a point [radians]
660 real, dimension(G%isd:G%ied,G%jsd:G%jed) :: &
661 xh, yh ! Latitude and longitude of h points in radians [radians]
662 real, dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: &
663 xu, yu ! Latitude and longitude of u points in radians [radians]
664 real, dimension(G%isd:G%ied,G%JsdB:G%JedB) :: &
665 xv, yv ! Latitude and longitude of v points in radians [radians]
666 real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: &
667 xq, yq ! Latitude and longitude of q points in radians [radians]
668 real :: fnRef ! fnRef is the value of Int_dj_dy or
669 ! Int_dj_dy at a latitude or longitude that is
670 ! being set to be at grid index jRef or iRef [gridpoints]
671 real :: jRef, iRef ! The grid index at which fnRef is evaluated [gridpoints]
672 integer :: itt1, itt2
673 logical, parameter :: simple_area = .true.
674 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
675
676 ! All of the metric terms should be defined over the domain from
677 ! isd to ied. Outside of the physical domain, both the metrics
678 ! and their inverses may be set to zero.
679 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
680 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
681 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
682 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
683 i_off = g%idg_offset ; j_off = g%jdg_offset
684
685 gp%niglobal = g%Domain%niglobal
686 gp%njglobal = g%Domain%njglobal
687
688 call calltree_enter("set_grid_metrics_mercator(), MOM_grid_initialize.F90")
689
690 ! Calculate the values of the metric terms that might be used
691 ! and save them in arrays.
692 pi = 4.0*atan(1.0) ; pi_2 = 0.5*pi
693
694 call get_param(param_file, mdl, "SOUTHLAT", gp%south_lat, &
695 "The southern latitude of the domain.", units="degrees_N", &
696 fail_if_missing=.true.)
697 call get_param(param_file, mdl, "LENLAT", gp%len_lat, &
698 "The latitudinal length of the domain.", units="degrees_N", &
699 fail_if_missing=.true.)
700 call get_param(param_file, mdl, "WESTLON", gp%west_lon, &
701 "The western longitude of the domain.", units="degrees_E", &
702 default=0.0)
703 call get_param(param_file, mdl, "LENLON", gp%len_lon, &
704 "The longitudinal length of the domain.", units="degrees_E", &
705 fail_if_missing=.true.)
706 call get_param(param_file, mdl, "RAD_EARTH", gp%Rad_Earth_L, &
707 "The radius of the Earth.", units="m", default=6.378e6, scale=us%m_to_L)
708 g%south_lat = gp%south_lat ; g%len_lat = gp%len_lat
709 g%west_lon = gp%west_lon ; g%len_lon = gp%len_lon
710 g%Rad_Earth_L = gp%Rad_Earth_L
711
712 call get_param(param_file, mdl, "ISOTROPIC", gp%isotropic, &
713 "If true, an isotropic grid on a sphere (also known as "//&
714 "a Mercator grid) is used. With an isotropic grid, the "//&
715 "meridional extent of the domain (LENLAT), the zonal "//&
716 "extent (LENLON), and the number of grid points in each "//&
717 "direction are _not_ independent. In MOM the meridional "//&
718 "extent is determined to fit the zonal extent and the "//&
719 "number of grid points, while grid is perfectly isotropic.", &
720 default=.false.)
721 call get_param(param_file, mdl, "EQUATOR_REFERENCE", gp%equator_reference, &
722 "If true, the grid is defined to have the equator at the "//&
723 "nearest q- or h- grid point to (-LOWLAT*NJGLOBAL/LENLAT).", &
724 default=.true.)
725 call get_param(param_file, mdl, "LAT_ENHANCE_FACTOR", gp%Lat_enhance_factor, &
726 "The amount by which the meridional resolution is "//&
727 "enhanced within LAT_EQ_ENHANCE of the equator.", &
728 units="nondim", default=1.0)
729 call get_param(param_file, mdl, "LAT_EQ_ENHANCE", gp%Lat_eq_enhance, &
730 "The latitude range to the north and south of the equator "//&
731 "over which the resolution is enhanced.", units="degrees_N", &
732 default=0.0)
733
734 ! With an isotropic grid, the north-south extent of the domain,
735 ! the east-west extent, and the number of grid points in each
736 ! direction are _not_ independent. Here the north-south extent
737 ! will be determined to fit the east-west extent and the number of
738 ! grid points. The grid is perfectly isotropic.
739 if (gp%equator_reference) then
740 ! With the following expression, the equator will always be placed
741 ! on either h or q points, in a position consistent with the ratio
742 ! GP%south_lat to GP%len_lat.
743 jref = (g%jsg-1) + 0.5*floor(gp%njglobal*((-1.0*gp%south_lat*2.0)/gp%len_lat)+0.5)
744 fnref = int_dj_dy(0.0, gp)
745 else
746 ! The following line sets the reference latitude GP%south_lat at j=js-1 (or -2?)
747 jref = (g%jsg-1)
748 fnref = int_dj_dy((gp%south_lat*pi/180.0), gp)
749 endif
750
751 ! These calculations no longer depend on the order in which they
752 ! are performed because they all use the same (poor) starting guess and
753 ! iterate to convergence.
754 ! Note that the dynamic grid always uses symmetric memory for the global
755 ! arrays G%gridLatB and G%gridLonB.
756 do j=g%jsg-1,g%jeg
757 jd = fnref + (j - jref)
758 y_q = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt2)
759 g%gridLatB(j) = y_q*180.0/pi
760 ! if (is_root_pe()) &
761 ! write(stdout, '("J, y_q = ",I0,", ",ES14.4," itts = ",I0)') j, y_q, itt2
762 enddo
763 do j=g%jsg,g%jeg
764 jd = fnref + (j - jref) - 0.5
765 y_h = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt1)
766 g%gridLatT(j) = y_h*180.0/pi
767 ! if (is_root_pe()) &
768 ! write(stdout, '("j, y_h = ",I0,", ",ES14.4," itts = ",I0)') j, y_h, itt1
769 enddo
770 do j=jsdb+j_off,jedb+j_off
771 jd = fnref + (j - jref)
772 y_q = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt2)
773 do i=isdb,iedb ; yq(i,j-j_off) = y_q ; enddo
774 do i=isd,ied ; yv(i,j-j_off) = y_q ; enddo
775 enddo
776 do j=jsd+j_off,jed+j_off
777 jd = fnref + (j - jref) - 0.5
778 y_h = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt1)
779 if ((j >= jsd+j_off) .and. (j <= jed+j_off)) then
780 do i=isd,ied ; yh(i,j-j_off) = y_h ; enddo
781 do i=isdb,iedb ; yu(i,j-j_off) = y_h ; enddo
782 endif
783 enddo
784
785 ! Determine the longitudes of the various points.
786
787 ! These two lines place the western edge of the domain at GP%west_lon.
788 iref = (g%isg-1) + gp%niglobal
789 fnref = int_di_dx(((gp%west_lon+gp%len_lon)*pi/180.0), gp)
790
791 ! These calculations no longer depend on the order in which they
792 ! are performed because they all use the same (poor) starting guess and
793 ! iterate to convergence.
794 do i=g%isg-1,g%ieg
795 id = fnref + (i - iref)
796 x_q = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt2)
797 g%gridLonB(i) = x_q*180.0/pi
798 enddo
799 do i=g%isg,g%ieg
800 id = fnref + (i - iref) - 0.5
801 x_h = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt1)
802 g%gridLonT(i) = x_h*180.0/pi
803 enddo
804 do i=isdb+i_off,iedb+i_off
805 id = fnref + (i - iref)
806 x_q = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt2)
807 do j=jsdb,jedb ; xq(i-i_off,j) = x_q ; enddo
808 do j=jsd,jed ; xu(i-i_off,j) = x_q ; enddo
809 enddo
810 do i=isd+i_off,ied+i_off
811 id = fnref + (i - iref) - 0.5
812 x_h = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt1)
813 do j=jsd,jed ; xh(i-i_off,j) = x_h ; enddo
814 do j=jsdb,jedb ; xv(i-i_off,j) = x_h ; enddo
815 enddo
816
817 do j=jsdb,jedb ; do i=isdb,iedb
818 g%geoLonBu(i,j) = xq(i,j)*180.0/pi
819 g%geoLatBu(i,j) = yq(i,j)*180.0/pi
820 g%dxBu(i,j) = ds_di(xq(i,j), yq(i,j), gp)
821 g%dyBu(i,j) = ds_dj(xq(i,j), yq(i,j), gp)
822
823 g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
824 g%IareaBu(i,j) = 1.0 / (g%areaBu(i,j))
825 enddo ; enddo
826
827 do j=jsd,jed ; do i=isd,ied
828 g%geoLonT(i,j) = xh(i,j)*180.0/pi
829 g%geoLatT(i,j) = yh(i,j)*180.0/pi
830 g%dxT(i,j) = ds_di(xh(i,j), yh(i,j), gp)
831 g%dyT(i,j) = ds_dj(xh(i,j), yh(i,j), gp)
832
833 g%areaT(i,j) = g%dxT(i,j)*g%dyT(i,j)
834 g%IareaT(i,j) = 1.0 / (g%areaT(i,j))
835 enddo ; enddo
836
837 do j=jsd,jed ; do i=isdb,iedb
838 g%geoLonCu(i,j) = xu(i,j)*180.0/pi
839 g%geoLatCu(i,j) = yu(i,j)*180.0/pi
840 g%dxCu(i,j) = ds_di(xu(i,j), yu(i,j), gp)
841 g%dyCu(i,j) = ds_dj(xu(i,j), yu(i,j), gp)
842 enddo ; enddo
843
844 do j=jsdb,jedb ; do i=isd,ied
845 g%geoLonCv(i,j) = xv(i,j)*180.0/pi
846 g%geoLatCv(i,j) = yv(i,j)*180.0/pi
847 g%dxCv(i,j) = ds_di(xv(i,j), yv(i,j), gp)
848 g%dyCv(i,j) = ds_dj(xv(i,j), yv(i,j), gp)
849 enddo ; enddo
850
851 if (.not.simple_area) then
852 do j=jsdb+1,jed ; do i=isdb+1,ied
853 g%areaT(i,j) = gp%Rad_Earth_L**2 * &
854 (dl(xq(i-1,j-1),xq(i-1,j),yq(i-1,j-1),yq(i-1,j)) + &
855 (dl(xq(i-1,j),xq(i,j),yq(i-1,j),yq(i,j)) + &
856 (dl(xq(i,j),xq(i,j-1),yq(i,j),yq(i,j-1)) + &
857 dl(xq(i,j-1),xq(i-1,j-1),yq(i,j-1),yq(i-1,j-1)))))
858 enddo ; enddo
859 if ((isdb == isd) .or. (jsdb == jsq)) then
860 ! Fill in row and column 1 to calculate the area in the southernmost
861 ! and westernmost land cells when we are not using symmetric memory.
862 ! The pass_var call updates these values if they are not land cells.
863 g%areaT(isd+1,jsd) = g%areaT(isd+1,jsd+1)
864 do j=jsd,jed ; g%areaT(isd,j) = g%areaT(isd+1,j) ; enddo
865 do i=isd,ied ; g%areaT(i,jsd) = g%areaT(i,jsd+1) ; enddo
866 ! Now replace the data in the halos, if value values exist.
867 call pass_var(g%areaT,g%Domain)
868 endif
869 do j=jsd,jed ; do i=isd,ied
870 g%IareaT(i,j) = 1.0 / (g%areaT(i,j))
871 enddo ; enddo
872 endif
873
874 call calltree_leave("set_grid_metrics_mercator()")
875end subroutine set_grid_metrics_mercator
876
877
878!> This function returns the grid spacing in the logical x direction in [L ~> m].
879function ds_di(x, y, GP)
880 real, intent(in) :: x !< The longitude in question [radians]
881 real, intent(in) :: y !< The latitude in question [radians]
882 type(gps), intent(in) :: gp !< A structure of grid parameters
883
884 real :: ds_di ! The returned grid spacing [L ~> m]
885
886 ds_di = gp%Rad_Earth_L * cos(y) * dx_di(x,gp)
887 ! In general, this might be...
888 ! ds_di = GP%Rad_Earth_L * sqrt( cos(y)*cos(y) * dx_di(x,y,GP)*dx_di(x,y,GP) + &
889 ! dy_di(x,y,GP)*dy_di(x,y,GP))
890end function ds_di
891
892!> This function returns the grid spacing in the logical y direction in [L ~> m].
893function ds_dj(x, y, GP)
894 real, intent(in) :: x !< The longitude in question [radians]
895 real, intent(in) :: y !< The latitude in question [radians]
896 type(gps), intent(in) :: gp !< A structure of grid parameters
897
898 real :: ds_dj ! The returned grid spacing [L ~> m]
899
900 ds_dj = gp%Rad_Earth_L * dy_dj(y,gp)
901 ! In general, this might be...
902 ! ds_dj = GP%Rad_Earth_L * sqrt( cos(y)*cos(y) * dx_dj(x,y,GP)*dx_dj(x,y,GP) + &
903 ! dy_dj(x,y,GP)*dy_dj(x,y,GP))
904end function ds_dj
905
906!> This function returns the contribution from the line integral along one of the four sides of a
907!! cell face to the area of a cell, in [radians2], assuming that the sides follow a linear path in
908!! latitude and longitude (i.e., on a Mercator grid).
909function dl(x1, x2, y1, y2)
910 real, intent(in) :: x1 !< Segment starting longitude [radians]
911 real, intent(in) :: x2 !< Segment ending longitude [radians]
912 real, intent(in) :: y1 !< Segment starting latitude [radians]
913 real, intent(in) :: y2 !< Segment ending latitude [radians]
914 ! Local variables
915 real :: dl ! A contribution to the spanned area the surface of the sphere [radian2]
916 real :: r ! A contribution from the range of latitudes, including trigonometric factors [radians]
917 real :: dy ! The spanned range of latitudes [radians]
918
919 dy = y2 - y1
920
921 if (abs(dy) > 2.5e-8) then
922 r = ((1.0 - cos(dy))*cos(y1) + sin(dy)*sin(y1)) / dy
923 else
924 r = (0.5*dy*cos(y1) + sin(y1))
925 endif
926 dl = r * (x2 - x1)
927
928end function dl
929
930!> This subroutine finds and returns the value of y at which the monotonically increasing
931!! function fn takes the value fnval, also returning in ittmax the number of iterations of
932!! Newton's method that were used to polish the root.
933function find_root( fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax)
934 real :: find_root !< The value of y where fn(y) = fnval that will be returned [radians]
935 real, external :: fn !< The external function whose root is being sought [gridpoints]
936 real, external :: dy_df !< The inverse of the derivative of that function [radian gridpoint-1]
937 type(gps), intent(in) :: gp !< A structure of grid parameters
938 real, intent(in) :: fnval !< The value of fn being sought [gridpoints]
939 real, intent(in) :: y1 !< A first guess for y [radians]
940 real, intent(in) :: ymin !< The minimum permitted value of y [radians]
941 real, intent(in) :: ymax !< The maximum permitted value of y [radians]
942 integer, intent(out) :: ittmax !< The number of iterations used to polish the root
943 ! Local variables
944 real :: y, y_next ! Successive guesses at the root position [radians]
945 real :: ybot, ytop ! Brackets bounding the root [radians]
946 real :: fnbot, fntop ! Values of fn at the bounding values of y [gridpoints]
947 real :: dy_dfn ! The inverse of the local derivative of fn with y [radian gridpoint-1]
948 real :: dy ! The jump to the next guess of y [radians]
949 real :: fny ! The difference between fn(y) and the target value [gridpoints]
950 integer :: itt
951 character(len=256) :: warnmesg
952
953! Bracket the root. Do not use the bounding values because the value at the
954! function at the bounds could be infinite, as is the case for the Mercator
955! grid recursion relation. (I.e., this is a search on an open interval.)
956 ybot = y1
957 fnbot = fn(ybot,gp) - fnval ; itt = 0
958 do while (fnbot > 0.0)
959 if ((ybot - 2.0*dy_df(ybot,gp)) < (0.5*(ybot+ymin))) then
960 ! Go twice as far as the secant method would normally go.
961 ybot = ybot - 2.0*dy_df(ybot,gp)
962 else ! But stay within the open interval!
963 ybot = 0.5*(ybot+ymin) ; itt = itt + 1
964 endif
965 fnbot = fn(ybot,gp) - fnval
966
967 if ((itt > 50) .and. (fnbot > 0.0)) then
968 write(warnmesg, '("PE ",I0," unable to find bottom bound for grid function. &
969 &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4,&
970 &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') &
971 pe_here(),ybot,ymin,fn(ybot,gp),dy_df(ybot,gp),fnval, fnbot
972 call mom_error(fatal,warnmesg)
973 endif
974 enddo
975
976 ytop = y1
977 fntop = fn(ytop,gp) - fnval ; itt = 0
978 do while (fntop < 0.0)
979 if ((ytop + 2.0*dy_df(ytop,gp)) < (0.5*(ytop+ymax))) then
980 ! Go twice as far as the secant method would normally go.
981 ytop = ytop + 2.0*dy_df(ytop,gp)
982 else ! But stay within the open interval!
983 ytop = 0.5*(ytop+ymax) ; itt = itt + 1
984 endif
985 fntop = fn(ytop,gp) - fnval
986
987 if ((itt > 50) .and. (fntop < 0.0)) then
988 write(warnmesg, '("PE ",I0," unable to find top bound for grid function. &
989 &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4, &
990 &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') &
991 pe_here(),ytop,ymax,fn(ytop,gp),dy_df(ytop,gp),fnval,fntop
992 call mom_error(fatal,warnmesg)
993 endif
994 enddo
995
996 ! Find the root using a bracketed variant of Newton's method, starting
997 ! with a false-positon method first guess.
998 if ((fntop < 0.0) .or. (fnbot > 0.0) .or. (ytop < ybot)) then
999 write(warnmesg, '("PE ",I0," find_root failed to bracket function. y = ",&
1000 &2ES10.4,", fn = ",2ES10.4,".")') pe_here(),ybot,ytop,fnbot,fntop
1001 call mom_error(fatal, warnmesg)
1002 endif
1003
1004 if (fntop == 0.0) then ; y = ytop ; fny = fntop
1005 elseif (fnbot == 0.0) then ; y = ybot ; fny = fnbot
1006 else
1007 y = (ybot*fntop - ytop*fnbot) / (fntop - fnbot)
1008 fny = fn(y,gp) - fnval
1009 if (fny < 0.0) then ; fnbot = fny ; ybot = y
1010 else ; fntop = fny ; ytop = y ; endif
1011 endif
1012
1013 do itt=1,50
1014 dy_dfn = dy_df(y,gp)
1015
1016 dy = -1.0* fny * dy_dfn
1017 y_next = y + dy
1018 if ((y_next >= ytop) .or. (y_next <= ybot)) then
1019 ! The Newton's method estimate has escaped bracketing, so use the
1020 ! false-position method instead. The complicated test is to properly
1021 ! handle the case where the iteration is down to roundoff level differences.
1022 y_next = y
1023 if (abs(fntop - fnbot) > epsilon(y) * (abs(fntop) + abs(fnbot))) &
1024 y_next = (ybot*fntop - ytop*fnbot) / (fntop - fnbot)
1025 endif
1026
1027 dy = y_next - y
1028 if (abs(dy) < (2.0*epsilon(y)*(abs(y) + abs(y_next)) + 1.0e-20)) then
1029 y = y_next ; exit
1030 endif
1031 y = y_next
1032
1033 fny = fn(y,gp) - fnval
1034 if (fny > 0.0) then ; ytop = y ; fntop = fny
1035 elseif (fny < 0.0) then ; ybot = y ; fnbot = fny
1036 else ; exit ; endif
1037
1038 enddo
1039 if (abs(y) < 1e-12) y = 0.0
1040
1041 ittmax = itt
1042 find_root = y
1043end function find_root
1044
1045!> This function calculates and returns the value of dx/di in [radian gridpoint-1],
1046!! where x is the longitude in Radians, and i is the integral east-west grid index.
1047function dx_di(x, GP)
1048 real, intent(in) :: x !< The longitude in question [radians]
1049 type(gps), intent(in) :: gp !< A structure of grid parameters
1050 real :: dx_di ! The derivative of zonal position with the grid index [radian gridpoint-1]
1051
1052 dx_di = (gp%len_lon * 4.0*atan(1.0)) / (180.0 * gp%niglobal)
1053
1054end function dx_di
1055
1056!> This function calculates and returns the integral of the inverse
1057!! of dx/di to the point x, in radians [gridpoints]
1058function int_di_dx(x, GP)
1059 real, intent(in) :: x !< The longitude in question [radians]
1060 type(gps), intent(in) :: gp !< A structure of grid parameters
1061 real :: int_di_dx ! A position in the global i-index space [gridpoints]
1062
1063 int_di_dx = x * ((180.0 * gp%niglobal) / (gp%len_lon * 4.0*atan(1.0)))
1064
1065end function int_di_dx
1066
1067!> This subroutine calculates and returns the value of dy/dj in [radian gridpoint-1],
1068!! where y is the latitude in Radians, and j is the integral north-south grid index.
1069function dy_dj(y, GP)
1070 real, intent(in) :: y !< The latitude in question [radians]
1071 type(gps), intent(in) :: gp !< A structure of grid parameters
1072 real :: dy_dj ! The derivative of meridional position with the grid index [radian gridpoint-1]
1073 ! Local variables
1074 real :: pi ! 3.1415926... calculated as 4*atan(1) [nondim]
1075 real :: c0 ! The constant that converts the nominal y-spacing in
1076 ! gridpoints to the nominal spacing in Radians [radian gridpoint-1]
1077 real :: y_eq_enhance ! The latitude in radians within which the resolution
1078 ! is enhanced [radians]
1079 pi = 4.0*atan(1.0)
1080 if (gp%isotropic) then
1081 c0 = (gp%len_lon * pi) / (180.0 * gp%niglobal)
1082 y_eq_enhance = pi*abs(gp%lat_eq_enhance)/180.0
1083 if (abs(y) < y_eq_enhance) then
1084 dy_dj = c0 * (cos(y) / (1.0 + 0.5*cos(y) * (gp%lat_enhance_factor - 1.0) * &
1085 (1.0+cos(pi*y/y_eq_enhance)) ))
1086 else
1087 dy_dj = c0 * cos(y)
1088 endif
1089 else
1090 c0 = (gp%len_lat * pi) / (180.0 * gp%njglobal)
1091 dy_dj = c0
1092 endif
1093
1094end function dy_dj
1095
1096!> This subroutine calculates and returns the integral of the inverse
1097!! of dy/dj to the point y in radians [gridpoints]
1098function int_dj_dy(y, GP)
1099 real, intent(in) :: y !< The latitude in question [radians]
1100 type(gps), intent(in) :: gp !< A structure of grid parameters
1101 real :: int_dj_dy ! The grid position of latitude y [gridpoints]
1102 ! Local variables
1103 real :: i_c0 ! The inverse of the constant that converts the
1104 ! nominal spacing in gridpoints to the nominal
1105 ! spacing in Radians [gridpoint radian-1]
1106 real :: pi ! 3.1415926... calculated as 4*atan(1) [nondim]
1107 real :: y_eq_enhance ! The latitude in radians from from the equator within which the meridional
1108 ! grid spacing is enhanced by a factor of GP%lat_enhance_factor [radians]
1109 real :: r ! The y grid position in the global index space [gridpoints]
1110
1111 pi = 4.0*atan(1.0)
1112 if (gp%isotropic) then
1113 i_c0 = (180.0 * gp%niglobal) / (gp%len_lon * pi)
1114 y_eq_enhance = pi*abs(gp%lat_eq_enhance)/180.0
1115
1116 if (y >= 0.0) then
1117 r = i_c0 * log((1.0 + sin(y))/cos(y))
1118 else
1119 r = -1.0 * i_c0 * log((1.0 - sin(y))/cos(y))
1120 endif
1121
1122 if (y >= y_eq_enhance) then
1123 r = r + i_c0*0.5*(gp%lat_enhance_factor - 1.0)*y_eq_enhance
1124 elseif (y <= -y_eq_enhance) then
1125 r = r - i_c0*0.5*(gp%lat_enhance_factor - 1.0)*y_eq_enhance
1126 else
1127 r = r + i_c0*0.5*(gp%lat_enhance_factor - 1.0) * &
1128 (y + (y_eq_enhance/pi)*sin(pi*y/y_eq_enhance))
1129 endif
1130 else
1131 i_c0 = (180.0 * gp%njglobal) / (gp%len_lat * pi)
1132 r = i_c0 * y
1133 endif
1134
1135 int_dj_dy = r
1136end function int_dj_dy
1137
1138!> Extrapolates missing metric data into all the halo regions.
1139subroutine extrapolate_metric(var, jh, missing)
1140 real, dimension(:,:), intent(inout) :: var !< The array in which to fill in halos in arbitrary units [A]
1141 integer, intent(in) :: jh !< The size of the halos to be filled
1142 real, optional, intent(in) :: missing !< The missing data fill value, 0 by default [A]
1143 ! Local variables
1144 real :: badval ! A bad data value [A]
1145 integer :: i, j
1146
1147 badval = 0.0 ; if (present(missing)) badval = missing
1148
1149 ! Fill in southern halo by extrapolating from the computational domain
1150 do j=lbound(var,2)+jh,lbound(var,2),-1 ; do i=lbound(var,1),ubound(var,1)
1151 if (var(i,j)==badval) var(i,j) = 2.0*var(i,j+1)-var(i,j+2)
1152 enddo ; enddo
1153
1154 ! Fill in northern halo by extrapolating from the computational domain
1155 do j=ubound(var,2)-jh,ubound(var,2) ; do i=lbound(var,1),ubound(var,1)
1156 if (var(i,j)==badval) var(i,j) = 2.0*var(i,j-1)-var(i,j-2)
1157 enddo ; enddo
1158
1159 ! Fill in western halo by extrapolating from the computational domain
1160 do j=lbound(var,2),ubound(var,2) ; do i=lbound(var,1)+jh,lbound(var,1),-1
1161 if (var(i,j)==badval) var(i,j) = 2.0*var(i+1,j)-var(i+2,j)
1162 enddo ; enddo
1163
1164 ! Fill in eastern halo by extrapolating from the computational domain
1165 do j=lbound(var,2),ubound(var,2) ; do i=ubound(var,1)-jh,ubound(var,1)
1166 if (var(i,j)==badval) var(i,j) = 2.0*var(i-1,j)-var(i-2,j)
1167 enddo ; enddo
1168
1169end subroutine extrapolate_metric
1170
1171!> This function implements Adcroft's rule for reciprocals, namely that
1172!! Adcroft_Inv(x) = 1/x for |x|>0 or 0 for x=0.
1173function adcroft_reciprocal(val) result(I_val)
1174 real, intent(in) :: val !< The value being inverted in arbitrary units [A]
1175 real :: i_val !< The Adcroft reciprocal of val [A-1]
1176
1177 i_val = 0.0
1178 if (val /= 0.0) i_val = 1.0/val
1179end function adcroft_reciprocal
1180
1181!> Initializes the grid masks and any metrics that come with masks already applied.
1182!!
1183!! Initialize_masks sets mask2dT, mask2dCu, mask2dCv, and mask2dBu to mask out
1184!! flow over any points which are shallower than Dmask and permit an
1185!! appropriate treatment of the boundary conditions. mask2dCu and mask2dCv
1186!! are 0.0 at any points adjacent to a land point. mask2dBu is 0.0 at
1187!! any land or boundary point. For points in the ocean interior or at open boundary
1188!! condition points, mask2dCu, mask2dCv, and mask2dBu are all 1.0.
1189subroutine initialize_masks(G, PF, US, OBC_dir_u, OBC_dir_v, open_corner_OBCs, maskT)
1190 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
1191 type(param_file_type), intent(in) :: pf !< Parameter file structure
1192 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1193 integer, dimension(G%IsdB:G%IedB,G%jsd:G%jed), &
1194 optional, intent(in) :: obc_dir_u !< Trinary values that indicate whether there
1195 !! is an open boundary condition at zonal velocity
1196 !! faces and their orientation, with 0 for no OBC,
1197 !! a positive value for an Eastern OBC and
1198 !! a negative value for a Western OBC.
1199 integer, dimension(G%isd:G%ied,G%JsdB:G%JedB), &
1200 optional, intent(in) :: obc_dir_v !< Trinary values that indicate whether there
1201 !! is an open boundary condition at zonal velocity
1202 !! faces and their orientation, with 0 for no OBC,
1203 !! a positive value for a Northern OBC and
1204 !! a negative value for a Southern OBC.
1205 logical, optional, intent(in) :: open_corner_obcs !< If present and true, the bay-like corner
1206 !! between two orthogonal open boundary segments is open,
1207 !! otherwise it is closed.
1208 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
1209 optional, intent(in) :: maskt !< If present, this array is used to set the
1210 !! the mask at tracer points instead of using the
1211 !! bathymetry to determine the masks [nondim]
1212
1213 ! Local variables
1214 real :: dmask ! The depth for masking in the same units as G%bathyT [Z ~> m].
1215 real :: min_depth ! The minimum ocean depth in the same units as G%bathyT [Z ~> m].
1216 real :: mask_depth ! The depth shallower than which to mask a point as land [Z ~> m].
1217 logical :: open_corners ! If true, the bay-like corner between two orthogonal open boundary segments is open
1218 character(len=40) :: mdl = "MOM_grid_init initialize_masks"
1219 integer :: i, j
1220
1221 call calltree_enter("initialize_masks(), MOM_grid_initialize.F90")
1222
1223 call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
1224 "If MASKING_DEPTH is unspecified, then anything shallower than "//&
1225 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
1226 "If MASKING_DEPTH is specified, then all depths shallower than "//&
1227 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
1228 units="m", default=0.0, scale=us%m_to_Z)
1229 call get_param(pf, mdl, "MASKING_DEPTH", mask_depth, &
1230 "The depth below which to mask points as land points, for which all "//&
1231 "fluxes are zeroed out. MASKING_DEPTH is ignored if it has the special "//&
1232 "default value.", &
1233 units="m", default=-9999.0, scale=us%m_to_Z, do_not_log=present(maskt))
1234
1235 dmask = mask_depth
1236 if (mask_depth == -9999.0*us%m_to_Z) dmask = min_depth
1237
1238 open_corners = .false. ; if (present(open_corner_obcs)) open_corners = open_corner_obcs
1239
1240 g%mask2dT(:,:) = 0.0 ; g%mask2dCu(:,:) = 0.0 ; g%mask2dCv(:,:) = 0.0 ; g%mask2dBu(:,:) = 0.0
1241
1242 ! Construct the h-point or T-point mask
1243 if (present(maskt)) then
1244 do j=g%jsd,g%jed ; do i=g%isd,g%ied
1245 g%mask2dT(i,j) = max(min(maskt(i,j), 1.0), 0.0)
1246 enddo ; enddo
1247 else
1248 do j=g%jsd,g%jed ; do i=g%isd,g%ied
1249 if (g%bathyT(i,j) <= dmask) then
1250 g%mask2dT(i,j) = 0.0
1251 else
1252 g%mask2dT(i,j) = 1.0
1253 endif
1254 enddo ; enddo
1255 endif
1256
1257 call pass_var(g%mask2dT, g%Domain)
1258
1259 do j=g%jsd,g%jed ; do i=g%isd,g%ied-1
1260 g%mask2dCu(i,j) = g%mask2dT(i,j) * g%mask2dT(i+1,j)
1261 enddo ; enddo
1262
1263 if (present(obc_dir_u)) then
1264 do j=g%jsd,g%jed ; do i=g%isd,g%ied-1
1265 if (obc_dir_u(i,j) > 0) g%mask2dCu(i,j) = g%mask2dT(i,j)
1266 if (obc_dir_u(i,j) < 0) g%mask2dCu(i,j) = g%mask2dT(i+1,j)
1267 enddo ; enddo
1268 endif
1269
1270 do j=g%jsd,g%jed ; do i=g%isd,g%ied-1
1271 ! This mask may be revised later after the open boundary positions are specified.
1272 g%OBCmaskCu(i,j) = g%mask2dCu(i,j)
1273 enddo ; enddo
1274
1275 do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied
1276 g%mask2dCv(i,j) = g%mask2dT(i,j) * g%mask2dT(i,j+1)
1277 enddo ; enddo
1278
1279 if (present(obc_dir_v)) then
1280 do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied
1281 if (obc_dir_v(i,j) > 0) g%mask2dCv(i,j) = g%mask2dT(i,j)
1282 if (obc_dir_v(i,j) < 0) g%mask2dCv(i,j) = g%mask2dT(i,j+1)
1283 enddo ; enddo
1284 endif
1285
1286 do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied
1287 ! This mask may be revised later after the open boundary positions are specified.
1288 g%OBCmaskCv(i,j) = g%mask2dCv(i,j)
1289 enddo ; enddo
1290
1291 ! The mask at the vertex can be determined from the masks at the faces.
1292 ! This works at interior ocean points or at convex OBC points.
1293 do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied-1
1294 g%mask2dBu(i,j) = (g%mask2dCu(i,j) * g%mask2dCu(i,j+1)) * (g%mask2dCv(i,j) * g%mask2dCv(i+1,j))
1295 enddo ; enddo
1296
1297 ! This block resets masks at the vertices when there are OBCs. The right logic is that if there
1298 ! are 2 or more unmasked OBCs, this point should be open, but to recreate the previous answers,
1299 if (present(obc_dir_u)) then
1300 do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied-1
1301 ! These are conditions to set open vertex points on a straight north-south coastline
1302 if ((g%mask2dCu(i,j) * obc_dir_u(i,j)) * (g%mask2dCu(i,j+1) * obc_dir_u(i,j+1)) > 0.) &
1303 g%mask2dBu(i,j) = 1.0
1304 enddo ; enddo
1305 endif
1306 if (present(obc_dir_v)) then
1307 do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied-1
1308 ! These are conditions to set open vertex points on a straight east-west coastline
1309 if ((g%mask2dCv(i,j) * obc_dir_v(i,j)) * (g%mask2dCv(i+1,j) * obc_dir_v(i+1,j)) > 0.) &
1310 g%mask2dBu(i,j) = 1.0
1311 enddo ; enddo
1312 endif
1313 if (open_corners .and. present(obc_dir_u) .and. present(obc_dir_v)) then
1314 do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied-1
1315 ! These are the 4 conditions to set an open point in a concave (bay-like) corner
1316 if ((g%mask2dCu(i,j+1) * obc_dir_u(i,j+1) < 0.) .and. (g%mask2dCv(i+1,j) * obc_dir_v(i+1,j) < 0.)) &
1317 g%mask2dBu(i,j) = 1.0 ! Southwestern corner
1318 if ((g%mask2dCu(i,j+1) * obc_dir_u(i,j+1) > 0.) .and. (g%mask2dCv(i,j) * obc_dir_v(i,j) < 0.)) &
1319 g%mask2dBu(i,j) = 1.0 ! Southeastern corner
1320 if ((g%mask2dCu(i,j) * obc_dir_u(i,j) < 0.) .and. (g%mask2dCv(i+1,j) * obc_dir_v(i+1,j) > 0.)) &
1321 g%mask2dBu(i,j) = 1.0 ! Northwestern corner
1322 if ((g%mask2dCu(i,j) * obc_dir_u(i,j) > 0.) .and. (g%mask2dCv(i,j) * obc_dir_v(i,j) > 0.)) &
1323 g%mask2dBu(i,j) = 1.0 ! Northeastern corner
1324 enddo ; enddo
1325 endif
1326
1327 call pass_var(g%mask2dBu, g%Domain, position=corner)
1328 call pass_vector(g%mask2dCu, g%mask2dCv, g%Domain, to_all+scalar_pair, cgrid_ne)
1329
1330 do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
1331 ! This open face length may be revised later.
1332 g%dy_Cu(i,j) = g%mask2dCu(i,j) * g%dyCu(i,j)
1333 g%IdxCu_OBCmask(i,j) = g%OBCmaskCu(i,j) * g%IdxCu(i,j)
1334 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
1335 g%IareaCu(i,j) = g%mask2dCu(i,j) * adcroft_reciprocal(g%areaCu(i,j))
1336 enddo ; enddo
1337
1338 do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
1339 ! This open face length may be revised later.
1340 g%dx_Cv(i,j) = g%mask2dCv(i,j) * g%dxCv(i,j)
1341 g%IdyCv_OBCmask(i,j) = g%OBCmaskCv(i,j) * g%IdyCv(i,j)
1342 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
1343 g%IareaCv(i,j) = g%mask2dCv(i,j) * adcroft_reciprocal(g%areaCv(i,j))
1344 enddo ; enddo
1345
1346 call calltree_leave("initialize_masks()")
1347end subroutine initialize_masks
1348
1349!> \namespace mom_grid_initialize
1350!!
1351!! The metric terms have the form Dzp, IDzp, or DXDYp, where z can
1352!! be X or Y, and p can be q, u, v, or h. z describes the direction
1353!! of the metric, while p describes the location. IDzp is the
1354!! inverse of Dzp, while DXDYp is the product of DXp and DYp except
1355!! that areaT is calculated analytically from the latitudes and
1356!! longitudes of the surrounding q points.
1357!!
1358!! On a sphere, a variety of grids can be implemented by defining
1359!! analytic expressions for dx_di, dy_dj (where x and y are latitude
1360!! and longitude, and i and j are grid indices) and the expressions
1361!! for the integrals of their inverses in the four subroutines
1362!! dy_dj, Int_dj_dy, dx_di, and Int_di_dx.
1363!!
1364!! initialize_masks sets up land masks based on the depth field.
1365!! The one argument is the minimum ocean depth. Depths that are
1366!! less than this are interpreted as land points.
1367
1368end module mom_grid_initialize