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