MOM6
MOM_grid.F90
Go to the documentation of this file.
1 !> Provides the ocean grid type
2 module mom_grid
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
7 use mom_domains, only : mom_domain_type, get_domain_extent, compute_block_extent
9 use mom_error_handler, only : mom_error, mom_mesg, fatal
12 
13 implicit none ; private
14 
15 #include <MOM_memory.h>
16 
19 
20 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
21 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
22 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
23 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
24 
25 !> Ocean grid type. See mom_grid for details.
26 type, public :: ocean_grid_type
27  type(mom_domain_type), pointer :: domain => null() !< Ocean model domain
28  type(mom_domain_type), pointer :: domain_aux => null() !< A non-symmetric auxiliary domain type.
29  type(hor_index_type) :: hi !< Horizontal index ranges
30  type(hor_index_type) :: hid2 !< Horizontal index ranges for level-2-downsampling
31 
32  integer :: isc !< The start i-index of cell centers within the computational domain
33  integer :: iec !< The end i-index of cell centers within the computational domain
34  integer :: jsc !< The start j-index of cell centers within the computational domain
35  integer :: jec !< The end j-index of cell centers within the computational domain
36 
37  integer :: isd !< The start i-index of cell centers within the data domain
38  integer :: ied !< The end i-index of cell centers within the data domain
39  integer :: jsd !< The start j-index of cell centers within the data domain
40  integer :: jed !< The end j-index of cell centers within the data domain
41 
42  integer :: isg !< The start i-index of cell centers within the global domain
43  integer :: ieg !< The end i-index of cell centers within the global domain
44  integer :: jsg !< The start j-index of cell centers within the global domain
45  integer :: jeg !< The end j-index of cell centers within the global domain
46 
47  integer :: iscb !< The start i-index of cell vertices within the computational domain
48  integer :: iecb !< The end i-index of cell vertices within the computational domain
49  integer :: jscb !< The start j-index of cell vertices within the computational domain
50  integer :: jecb !< The end j-index of cell vertices within the computational domain
51 
52  integer :: isdb !< The start i-index of cell vertices within the data domain
53  integer :: iedb !< The end i-index of cell vertices within the data domain
54  integer :: jsdb !< The start j-index of cell vertices within the data domain
55  integer :: jedb !< The end j-index of cell vertices within the data domain
56 
57  integer :: isgb !< The start i-index of cell vertices within the global domain
58  integer :: iegb !< The end i-index of cell vertices within the global domain
59  integer :: jsgb !< The start j-index of cell vertices within the global domain
60  integer :: jegb !< The end j-index of cell vertices within the global domain
61 
62  integer :: isd_global !< The value of isd in the global index space (decompoistion invariant).
63  integer :: jsd_global !< The value of isd in the global index space (decompoistion invariant).
64  integer :: idg_offset !< The offset between the corresponding global and local i-indices.
65  integer :: jdg_offset !< The offset between the corresponding global and local j-indices.
66  integer :: ke !< The number of layers in the vertical.
67  logical :: symmetric !< True if symmetric memory is used.
68  logical :: nonblocking_updates !< If true, non-blocking halo updates are
69  !! allowed. The default is .false. (for now).
70  integer :: first_direction !< An integer that indicates which direction is
71  !! to be updated first in directionally split
72  !! parts of the calculation. This can be altered
73  !! during the course of the run via calls to
74  !! set_first_direction.
75 
76  real allocable_, dimension(NIMEM_,NJMEM_) :: &
77  mask2dt, & !< 0 for land points and 1 for ocean points on the h-grid [nondim].
78  geolatt, & !< The geographic latitude at q points in degrees of latitude or m.
79  geolont, & !< The geographic longitude at q points in degrees of longitude or m.
80  dxt, & !< dxT is delta x at h points [L ~> m].
81  idxt, & !< 1/dxT [L-1 ~> m-1].
82  dyt, & !< dyT is delta y at h points [L ~> m].
83  idyt, & !< IdyT is 1/dyT [L-1 ~> m-1].
84  areat, & !< The area of an h-cell [L2 ~> m2].
85  iareat, & !< 1/areaT [L-2 ~> m-2].
86  sin_rot, & !< The sine of the angular rotation between the local model grid's northward
87  !! and the true northward directions.
88  cos_rot !< The cosine of the angular rotation between the local model grid's northward
89  !! and the true northward directions.
90 
91  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
92  mask2dcu, & !< 0 for boundary points and 1 for ocean points on the u grid [nondim].
93  geolatcu, & !< The geographic latitude at u points in degrees of latitude or m.
94  geoloncu, & !< The geographic longitude at u points in degrees of longitude or m.
95  dxcu, & !< dxCu is delta x at u points [L ~> m].
96  idxcu, & !< 1/dxCu [L-1 ~> m-1].
97  dycu, & !< dyCu is delta y at u points [L ~> m].
98  idycu, & !< 1/dyCu [L-1 ~> m-1].
99  dy_cu, & !< The unblocked lengths of the u-faces of the h-cell [L ~> m].
100  iareacu, & !< The masked inverse areas of u-grid cells [L-2 ~> m-2].
101  areacu !< The areas of the u-grid cells [L2 ~> m2].
102 
103  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
104  mask2dcv, & !< 0 for boundary points and 1 for ocean points on the v grid [nondim].
105  geolatcv, & !< The geographic latitude at v points in degrees of latitude or m.
106  geoloncv, & !< The geographic longitude at v points in degrees of longitude or m.
107  dxcv, & !< dxCv is delta x at v points [L ~> m].
108  idxcv, & !< 1/dxCv [L-1 ~> m-1].
109  dycv, & !< dyCv is delta y at v points [L ~> m].
110  idycv, & !< 1/dyCv [L-1 ~> m-1].
111  dx_cv, & !< The unblocked lengths of the v-faces of the h-cell [L ~> m].
112  iareacv, & !< The masked inverse areas of v-grid cells [L-2 ~> m-2].
113  areacv !< The areas of the v-grid cells [L2 ~> m2].
114 
115  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
116  mask2dbu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim].
117  geolatbu, & !< The geographic latitude at q points in degrees of latitude or m.
118  geolonbu, & !< The geographic longitude at q points in degrees of longitude or m.
119  dxbu, & !< dxBu is delta x at q points [L ~> m].
120  idxbu, & !< 1/dxBu [L-1 ~> m-1].
121  dybu, & !< dyBu is delta y at q points [L ~> m].
122  idybu, & !< 1/dyBu [L-1 ~> m-1].
123  areabu, & !< areaBu is the area of a q-cell [L2 ~> m2]
124  iareabu !< IareaBu = 1/areaBu [L-2 ~> m-2].
125 
126  real, pointer, dimension(:) :: &
127  gridlatt => null(), & !< The latitude of T points for the purpose of labeling the output axes.
128  !! On many grids this is the same as geoLatT.
129  gridlatb => null() !< The latitude of B points for the purpose of labeling the output axes.
130  !! On many grids this is the same as geoLatBu.
131  real, pointer, dimension(:) :: &
132  gridlont => null(), & !< The longitude of T points for the purpose of labeling the output axes.
133  !! On many grids this is the same as geoLonT.
134  gridlonb => null() !< The longitude of B points for the purpose of labeling the output axes.
135  !! On many grids this is the same as geoLonBu.
136  character(len=40) :: &
137  x_axis_units, & !< The units that are used in labeling the x coordinate axes.
138  y_axis_units !< The units that are used in labeling the y coordinate axes.
139 
140  real allocable_, dimension(NIMEM_,NJMEM_) :: &
141  bathyt !< Ocean bottom depth at tracer points, in depth units [Z ~> m].
142 
143  logical :: bathymetry_at_vel !< If true, there are separate values for the
144  !! basin depths at velocity points. Otherwise the effects of
145  !! of topography are entirely determined from thickness points.
146  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
147  dblock_u, & !< Topographic depths at u-points at which the flow is blocked [Z ~> m].
148  dopen_u !< Topographic depths at u-points at which the flow is open at width dy_Cu [Z ~> m].
149  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
150  dblock_v, & !< Topographic depths at v-points at which the flow is blocked [Z ~> m].
151  dopen_v !< Topographic depths at v-points at which the flow is open at width dx_Cv [Z ~> m].
152  real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
153  coriolisbu !< The Coriolis parameter at corner points [T-1 ~> s-1].
154  real allocable_, dimension(NIMEM_,NJMEM_) :: &
155  df_dx, & !< Derivative d/dx f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
156  df_dy !< Derivative d/dy f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
157  real :: g_earth !< The gravitational acceleration [m2 Z-1 s-2 ~> m s-2].
158 
159  ! These variables are global sums that are useful for 1-d diagnostics and should not be rescaled.
160  real :: areat_global !< Global sum of h-cell area [m2]
161  real :: iareat_global !< Global sum of inverse h-cell area (1/areaT_global) [m-2].
162 
163  type(unit_scale_type), pointer :: us => null() !< A dimensional unit scaling type
164 
165 
166  ! These variables are for block structures.
167  integer :: nblocks !< The number of sub-PE blocks on this PE
168  type(hor_index_type), pointer :: block(:) => null() !< Index ranges for each block
169 
170  ! These parameters are run-time parameters that are used during some
171  ! initialization routines (but not all)
172  real :: south_lat !< The latitude (or y-coordinate) of the first v-line
173  real :: west_lon !< The longitude (or x-coordinate) of the first u-line
174  real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain
175  real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain
176  real :: rad_earth = 6.378e6 !< The radius of the planet [m].
177  real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m].
178 end type ocean_grid_type
179 
180 contains
181 
182 !> MOM_grid_init initializes the ocean grid array sizes and grid memory.
183 subroutine mom_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_vel)
184  type(ocean_grid_type), intent(inout) :: g !< The horizontal grid type
185  type(param_file_type), intent(in) :: param_file !< Parameter file handle
186  type(unit_scale_type), optional, pointer :: us !< A dimensional unit scaling type
187  type(hor_index_type), &
188  optional, intent(in) :: hi !< A hor_index_type for array extents
189  logical, optional, intent(in) :: global_indexing !< If true use global index
190  !! values instead of having the data domain on each
191  !! processor start at 1.
192  logical, optional, intent(in) :: bathymetry_at_vel !< If true, there are
193  !! separate values for the ocean bottom depths at
194  !! velocity points. Otherwise the effects of topography
195  !! are entirely determined from thickness points.
196 
197 ! This include declares and sets the variable "version".
198 #include "version_variable.h"
199  integer :: isd, ied, jsd, jed, nk
200  integer :: isdb, iedb, jsdb, jedb
201  integer :: ied_max, jed_max
202  integer :: niblock, njblock, nihalo, njhalo, nblocks, n, i, j
203  logical :: local_indexing ! If false use global index values instead of having
204  ! the data domain on each processor start at 1.
205 
206  integer, allocatable, dimension(:) :: ibegin, iend, jbegin, jend
207  character(len=40) :: mod_nm = "MOM_grid" ! This module's name.
208 
209 
210  ! Read all relevant parameters and write them to the model log.
211  call log_version(param_file, mod_nm, version, &
212  "Parameters providing information about the lateral grid.")
213 
214 
215  call get_param(param_file, mod_nm, "NIBLOCK", niblock, "The number of blocks "// &
216  "in the x-direction on each processor (for openmp).", default=1, &
217  layoutparam=.true.)
218  call get_param(param_file, mod_nm, "NJBLOCK", njblock, "The number of blocks "// &
219  "in the y-direction on each processor (for openmp).", default=1, &
220  layoutparam=.true.)
221 
222  if (present(us)) then ; if (associated(us)) g%US => us ; endif
223 
224  if (present(hi)) then
225  g%HI = hi
226 
227  g%isc = hi%isc ; g%iec = hi%iec ; g%jsc = hi%jsc ; g%jec = hi%jec
228  g%isd = hi%isd ; g%ied = hi%ied ; g%jsd = hi%jsd ; g%jed = hi%jed
229  g%isg = hi%isg ; g%ieg = hi%ieg ; g%jsg = hi%jsg ; g%jeg = hi%jeg
230 
231  g%IscB = hi%IscB ; g%IecB = hi%IecB ; g%JscB = hi%JscB ; g%JecB = hi%JecB
232  g%IsdB = hi%IsdB ; g%IedB = hi%IedB ; g%JsdB = hi%JsdB ; g%JedB = hi%JedB
233  g%IsgB = hi%IsgB ; g%IegB = hi%IegB ; g%JsgB = hi%JsgB ; g%JegB = hi%JegB
234 
235  g%idg_offset = hi%idg_offset ; g%jdg_offset = hi%jdg_offset
236  g%isd_global = g%isd + hi%idg_offset ; g%jsd_global = g%jsd + hi%jdg_offset
237  g%symmetric = hi%symmetric
238  else
239  local_indexing = .true.
240  if (present(global_indexing)) local_indexing = .not.global_indexing
241  call hor_index_init(g%Domain, g%HI, param_file, &
242  local_indexing=local_indexing)
243 
244  ! get_domain_extent ensures that domains start at 1 for compatibility between
245  ! static and dynamically allocated arrays, unless global_indexing is true.
246  call get_domain_extent(g%Domain, g%isc, g%iec, g%jsc, g%jec, &
247  g%isd, g%ied, g%jsd, g%jed, &
248  g%isg, g%ieg, g%jsg, g%jeg, &
249  g%idg_offset, g%jdg_offset, g%symmetric, &
250  local_indexing=local_indexing)
251  g%isd_global = g%isd+g%idg_offset ; g%jsd_global = g%jsd+g%jdg_offset
252  endif
253 
254  g%nonblocking_updates = g%Domain%nonblocking_updates
255 
256  ! Set array sizes for fields that are discretized at tracer cell boundaries.
257  g%IscB = g%isc ; g%JscB = g%jsc
258  g%IsdB = g%isd ; g%JsdB = g%jsd
259  g%IsgB = g%isg ; g%JsgB = g%jsg
260  if (g%symmetric) then
261  g%IscB = g%isc-1 ; g%JscB = g%jsc-1
262  g%IsdB = g%isd-1 ; g%JsdB = g%jsd-1
263  g%IsgB = g%isg-1 ; g%JsgB = g%jsg-1
264  endif
265  g%IecB = g%iec ; g%JecB = g%jec
266  g%IedB = g%ied ; g%JedB = g%jed
267  g%IegB = g%ieg ; g%JegB = g%jeg
268 
269  call mom_mesg(" MOM_grid.F90, MOM_grid_init: allocating metrics", 5)
270 
271  call allocate_metrics(g)
272 
273  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
274  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
275 
276  g%bathymetry_at_vel = .false.
277  if (present(bathymetry_at_vel)) g%bathymetry_at_vel = bathymetry_at_vel
278  if (g%bathymetry_at_vel) then
279  alloc_(g%Dblock_u(isdb:iedb, jsd:jed)) ; g%Dblock_u(:,:) = 0.0
280  alloc_(g%Dopen_u(isdb:iedb, jsd:jed)) ; g%Dopen_u(:,:) = 0.0
281  alloc_(g%Dblock_v(isd:ied, jsdb:jedb)) ; g%Dblock_v(:,:) = 0.0
282  alloc_(g%Dopen_v(isd:ied, jsdb:jedb)) ; g%Dopen_v(:,:) = 0.0
283  endif
284 
285 ! setup block indices.
286  nihalo = g%Domain%nihalo
287  njhalo = g%Domain%njhalo
288  nblocks = niblock * njblock
289  if (nblocks < 1) call mom_error(fatal, "MOM_grid_init: " // &
290  "nblocks(=NI_BLOCK*NJ_BLOCK) must be no less than 1")
291 
292  allocate(ibegin(niblock), iend(niblock), jbegin(njblock), jend(njblock))
293  call compute_block_extent(g%HI%isc,g%HI%iec,niblock,ibegin,iend)
294  call compute_block_extent(g%HI%jsc,g%HI%jec,njblock,jbegin,jend)
295  !-- make sure the last block is the largest.
296  do i = 1, niblock-1
297  if (iend(i)-ibegin(i) > iend(niblock)-ibegin(niblock) ) call mom_error(fatal, &
298  "MOM_grid_init: the last block size in x-direction is not the largest")
299  enddo
300  do j = 1, njblock-1
301  if (jend(j)-jbegin(j) > jend(njblock)-jbegin(njblock) ) call mom_error(fatal, &
302  "MOM_grid_init: the last block size in y-direction is not the largest")
303  enddo
304 
305  g%nblocks = nblocks
306  allocate(g%Block(nblocks))
307  ied_max = 1 ; jed_max = 1
308  do n = 1,nblocks
309  ! Copy all information from the array index type describing the local grid.
310  g%Block(n) = g%HI
311 
312  i = mod((n-1), niblock) + 1
313  j = (n-1)/niblock + 1
314  !--- isd and jsd are always 1 for each block to permit array reuse.
315  g%Block(n)%isd = 1 ; g%Block(n)%jsd = 1
316  g%Block(n)%isc = g%Block(n)%isd+nihalo
317  g%Block(n)%jsc = g%Block(n)%jsd+njhalo
318  g%Block(n)%iec = g%Block(n)%isc + iend(i) - ibegin(i)
319  g%Block(n)%jec = g%Block(n)%jsc + jend(j) - jbegin(j)
320  g%Block(n)%ied = g%Block(n)%iec + nihalo
321  g%Block(n)%jed = g%Block(n)%jec + njhalo
322  g%Block(n)%IscB = g%Block(n)%isc; g%Block(n)%IecB = g%Block(n)%iec
323  g%Block(n)%JscB = g%Block(n)%jsc; g%Block(n)%JecB = g%Block(n)%jec
324  ! For symmetric memory domains, the first block will have the extra point
325  ! at the lower boundary of its computational domain.
326  if (g%symmetric) then
327  if (i==1) g%Block(n)%IscB = g%Block(n)%IscB-1
328  if (j==1) g%Block(n)%JscB = g%Block(n)%JscB-1
329  endif
330  g%Block(n)%IsdB = g%Block(n)%isd; g%Block(n)%IedB = g%Block(n)%ied
331  g%Block(n)%JsdB = g%Block(n)%jsd; g%Block(n)%JedB = g%Block(n)%jed
332  !--- For symmetric memory domain, every block will have an extra point
333  !--- at the lower boundary of its data domain.
334  if (g%symmetric) then
335  g%Block(n)%IsdB = g%Block(n)%IsdB-1
336  g%Block(n)%JsdB = g%Block(n)%JsdB-1
337  endif
338  g%Block(n)%idg_offset = (ibegin(i) - g%Block(n)%isc) + g%HI%idg_offset
339  g%Block(n)%jdg_offset = (jbegin(j) - g%Block(n)%jsc) + g%HI%jdg_offset
340  ! Find the largest values of ied and jed so that all blocks will have the
341  ! same size in memory.
342  ied_max = max(ied_max, g%Block(n)%ied)
343  jed_max = max(jed_max, g%Block(n)%jed)
344  enddo
345 
346  ! Reset all of the data domain sizes to match the largest for array reuse,
347  ! recalling that all block have isd=jed=1 for array reuse.
348  do n = 1,nblocks
349  g%Block(n)%ied = ied_max ; g%Block(n)%IedB = ied_max
350  g%Block(n)%jed = jed_max ; g%Block(n)%JedB = jed_max
351  enddo
352 
353  !-- do some bounds error checking
354  if ( g%block(nblocks)%ied+g%block(nblocks)%idg_offset > g%HI%ied + g%HI%idg_offset ) &
355  call mom_error(fatal, "MOM_grid_init: G%ied_bk > G%ied")
356  if ( g%block(nblocks)%jed+g%block(nblocks)%jdg_offset > g%HI%jed + g%HI%jdg_offset ) &
357  call mom_error(fatal, "MOM_grid_init: G%jed_bk > G%jed")
358 
359  call get_domain_extent_dsamp2(g%Domain, g%HId2%isc, g%HId2%iec, g%HId2%jsc, g%HId2%jec,&
360  g%HId2%isd, g%HId2%ied, g%HId2%jsd, g%HId2%jed,&
361  g%HId2%isg, g%HId2%ieg, g%HId2%jsg, g%HId2%jeg)
362 
363  ! Set array sizes for fields that are discretized at tracer cell boundaries.
364  g%HId2%IscB = g%HId2%isc ; g%HId2%JscB = g%HId2%jsc
365  g%HId2%IsdB = g%HId2%isd ; g%HId2%JsdB = g%HId2%jsd
366  g%HId2%IsgB = g%HId2%isg ; g%HId2%JsgB = g%HId2%jsg
367  if (g%symmetric) then
368  g%HId2%IscB = g%HId2%isc-1 ; g%HId2%JscB = g%HId2%jsc-1
369  g%HId2%IsdB = g%HId2%isd-1 ; g%HId2%JsdB = g%HId2%jsd-1
370  g%HId2%IsgB = g%HId2%isg-1 ; g%HId2%JsgB = g%HId2%jsg-1
371  endif
372  g%HId2%IecB = g%HId2%iec ; g%HId2%JecB = g%HId2%jec
373  g%HId2%IedB = g%HId2%ied ; g%HId2%JedB = g%HId2%jed
374  g%HId2%IegB = g%HId2%ieg ; g%HId2%JegB = g%HId2%jeg
375 
376 end subroutine mom_grid_init
377 
378 !> rescale_grid_bathymetry permits a change in the internal units for the bathymetry on the grid,
379 !! both rescaling the depths and recording the new internal units.
380 subroutine rescale_grid_bathymetry(G, m_in_new_units)
381  type(ocean_grid_type), intent(inout) :: g !< The horizontal grid structure
382  real, intent(in) :: m_in_new_units !< The new internal representation of 1 m depth.
383 
384  ! Local variables
385  real :: rescale
386  integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
387 
388  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
389  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
390 
391  if (m_in_new_units == 1.0) return
392  if (m_in_new_units < 0.0) &
393  call mom_error(fatal, "rescale_grid_bathymetry: Negative depth units are not permitted.")
394  if (m_in_new_units == 0.0) &
395  call mom_error(fatal, "rescale_grid_bathymetry: Zero depth units are not permitted.")
396 
397  rescale = 1.0 / m_in_new_units
398  do j=jsd,jed ; do i=isd,ied
399  g%bathyT(i,j) = rescale*g%bathyT(i,j)
400  enddo ; enddo
401  if (g%bathymetry_at_vel) then ; do j=jsd,jed ; do i=isdb,iedb
402  g%Dblock_u(i,j) = rescale*g%Dblock_u(i,j) ; g%Dopen_u(i,j) = rescale*g%Dopen_u(i,j)
403  enddo ; enddo ; endif
404  if (g%bathymetry_at_vel) then ; do j=jsdb,jedb ; do i=isd,ied
405  g%Dblock_v(i,j) = rescale*g%Dblock_v(i,j) ; g%Dopen_v(i,j) = rescale*g%Dopen_v(i,j)
406  enddo ; enddo ; endif
407  g%max_depth = rescale*g%max_depth
408 
409 end subroutine rescale_grid_bathymetry
410 
411 !> set_derived_metrics calculates metric terms that are derived from other metrics.
412 subroutine set_derived_metrics(G, US)
413  type(ocean_grid_type), intent(inout) :: g !< The horizontal grid structure
414  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
415 ! Various inverse grid spacings and derived areas are calculated within this
416 ! subroutine.
417  integer :: i, j, isd, ied, jsd, jed
418  integer :: isdb, iedb, jsdb, jedb
419 
420  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
421  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
422 
423  do j=jsd,jed ; do i=isd,ied
424  if (g%dxT(i,j) < 0.0) g%dxT(i,j) = 0.0
425  if (g%dyT(i,j) < 0.0) g%dyT(i,j) = 0.0
426  g%IdxT(i,j) = adcroft_reciprocal(g%dxT(i,j))
427  g%IdyT(i,j) = adcroft_reciprocal(g%dyT(i,j))
428  g%IareaT(i,j) = adcroft_reciprocal(g%areaT(i,j))
429  enddo ; enddo
430 
431  do j=jsd,jed ; do i=isdb,iedb
432  if (g%dxCu(i,j) < 0.0) g%dxCu(i,j) = 0.0
433  if (g%dyCu(i,j) < 0.0) g%dyCu(i,j) = 0.0
434  g%IdxCu(i,j) = adcroft_reciprocal(g%dxCu(i,j))
435  g%IdyCu(i,j) = adcroft_reciprocal(g%dyCu(i,j))
436  enddo ; enddo
437 
438  do j=jsdb,jedb ; do i=isd,ied
439  if (g%dxCv(i,j) < 0.0) g%dxCv(i,j) = 0.0
440  if (g%dyCv(i,j) < 0.0) g%dyCv(i,j) = 0.0
441  g%IdxCv(i,j) = adcroft_reciprocal(g%dxCv(i,j))
442  g%IdyCv(i,j) = adcroft_reciprocal(g%dyCv(i,j))
443  enddo ; enddo
444 
445  do j=jsdb,jedb ; do i=isdb,iedb
446  if (g%dxBu(i,j) < 0.0) g%dxBu(i,j) = 0.0
447  if (g%dyBu(i,j) < 0.0) g%dyBu(i,j) = 0.0
448 
449  g%IdxBu(i,j) = adcroft_reciprocal(g%dxBu(i,j))
450  g%IdyBu(i,j) = adcroft_reciprocal(g%dyBu(i,j))
451  ! areaBu has usually been set to a positive area elsewhere.
452  if (g%areaBu(i,j) <= 0.0) g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
453  g%IareaBu(i,j) = adcroft_reciprocal(g%areaBu(i,j))
454  enddo ; enddo
455 end subroutine set_derived_metrics
456 
457 !> Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0.
458 function adcroft_reciprocal(val) result(I_val)
459  real, intent(in) :: val !< The value being inverted.
460  real :: i_val !< The Adcroft reciprocal of val.
461 
462  i_val = 0.0 ; if (val /= 0.0) i_val = 1.0/val
463 end function adcroft_reciprocal
464 
465 !> Returns true if the coordinates (x,y) are within the h-cell (i,j)
466 logical function ispointincell(G, i, j, x, y)
467  type(ocean_grid_type), intent(in) :: g !< Grid type
468  integer, intent(in) :: i !< i index of cell to test
469  integer, intent(in) :: j !< j index of cell to test
470  real, intent(in) :: x !< x coordinate of point
471  real, intent(in) :: y !< y coordinate of point
472  ! Local variables
473  real :: xne, xnw, xse, xsw, yne, ynw, yse, ysw
474  real :: p0, p1, p2, p3, l0, l1, l2, l3
475  ispointincell = .false.
476  xne = g%geoLonBu(i ,j ) ; yne = g%geoLatBu(i ,j )
477  xnw = g%geoLonBu(i-1,j ) ; ynw = g%geoLatBu(i-1,j )
478  xse = g%geoLonBu(i ,j-1) ; yse = g%geoLatBu(i ,j-1)
479  xsw = g%geoLonBu(i-1,j-1) ; ysw = g%geoLatBu(i-1,j-1)
480  ! This is a crude calculation that assume a geographic coordinate system
481  if (x<min(xne,xnw,xse,xsw) .or. x>max(xne,xnw,xse,xsw) .or. &
482  y<min(yne,ynw,yse,ysw) .or. y>max(yne,ynw,yse,ysw) ) then
483  return ! Avoid the more complicated calculation
484  endif
485  l0 = (x-xsw)*(yse-ysw) - (y-ysw)*(xse-xsw)
486  l1 = (x-xse)*(yne-yse) - (y-yse)*(xne-xse)
487  l2 = (x-xne)*(ynw-yne) - (y-yne)*(xnw-xne)
488  l3 = (x-xnw)*(ysw-ynw) - (y-ynw)*(xsw-xnw)
489 
490  p0 = sign(1., l0) ; if (l0 == 0.) p0=0.
491  p1 = sign(1., l1) ; if (l1 == 0.) p1=0.
492  p2 = sign(1., l2) ; if (l2 == 0.) p2=0.
493  p3 = sign(1., l3) ; if (l3 == 0.) p3=0.
494 
495  if ( (abs(p0)+abs(p2)) + (abs(p1)+abs(p3)) == abs((p0+p2) + (p1+p3)) ) then
496  ispointincell=.true.
497  endif
498 end function ispointincell
499 
500 !> Store an integer indicating which direction to work on first.
501 subroutine set_first_direction(G, y_first)
502  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
503  integer, intent(in) :: y_first !< The first direction to store
504 
505  g%first_direction = y_first
506 end subroutine set_first_direction
507 
508 !> Return global shape of horizontal grid
509 subroutine get_global_grid_size(G, niglobal, njglobal)
510  type(ocean_grid_type), intent(inout) :: g !< The horizontal grid type
511  integer, intent(out) :: niglobal !< i-index global size of grid
512  integer, intent(out) :: njglobal !< j-index global size of grid
513 
514  call get_global_shape(g%domain, niglobal, njglobal)
515 
516 end subroutine get_global_grid_size
517 
518 !> Allocate memory used by the ocean_grid_type and related structures.
519 subroutine allocate_metrics(G)
520  type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type
521  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isg, ieg, jsg, jeg
522 
523  ! This subroutine allocates the lateral elements of the ocean_grid_type that
524  ! are always used and zeros them out.
525 
526  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
527  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
528  isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
529 
530  alloc_(g%dxT(isd:ied,jsd:jed)) ; g%dxT(:,:) = 0.0
531  alloc_(g%dxCu(isdb:iedb,jsd:jed)) ; g%dxCu(:,:) = 0.0
532  alloc_(g%dxCv(isd:ied,jsdb:jedb)) ; g%dxCv(:,:) = 0.0
533  alloc_(g%dxBu(isdb:iedb,jsdb:jedb)) ; g%dxBu(:,:) = 0.0
534  alloc_(g%IdxT(isd:ied,jsd:jed)) ; g%IdxT(:,:) = 0.0
535  alloc_(g%IdxCu(isdb:iedb,jsd:jed)) ; g%IdxCu(:,:) = 0.0
536  alloc_(g%IdxCv(isd:ied,jsdb:jedb)) ; g%IdxCv(:,:) = 0.0
537  alloc_(g%IdxBu(isdb:iedb,jsdb:jedb)) ; g%IdxBu(:,:) = 0.0
538 
539  alloc_(g%dyT(isd:ied,jsd:jed)) ; g%dyT(:,:) = 0.0
540  alloc_(g%dyCu(isdb:iedb,jsd:jed)) ; g%dyCu(:,:) = 0.0
541  alloc_(g%dyCv(isd:ied,jsdb:jedb)) ; g%dyCv(:,:) = 0.0
542  alloc_(g%dyBu(isdb:iedb,jsdb:jedb)) ; g%dyBu(:,:) = 0.0
543  alloc_(g%IdyT(isd:ied,jsd:jed)) ; g%IdyT(:,:) = 0.0
544  alloc_(g%IdyCu(isdb:iedb,jsd:jed)) ; g%IdyCu(:,:) = 0.0
545  alloc_(g%IdyCv(isd:ied,jsdb:jedb)) ; g%IdyCv(:,:) = 0.0
546  alloc_(g%IdyBu(isdb:iedb,jsdb:jedb)) ; g%IdyBu(:,:) = 0.0
547 
548  alloc_(g%areaT(isd:ied,jsd:jed)) ; g%areaT(:,:) = 0.0
549  alloc_(g%IareaT(isd:ied,jsd:jed)) ; g%IareaT(:,:) = 0.0
550  alloc_(g%areaBu(isdb:iedb,jsdb:jedb)) ; g%areaBu(:,:) = 0.0
551  alloc_(g%IareaBu(isdb:iedb,jsdb:jedb)) ; g%IareaBu(:,:) = 0.0
552 
553  alloc_(g%mask2dT(isd:ied,jsd:jed)) ; g%mask2dT(:,:) = 0.0
554  alloc_(g%mask2dCu(isdb:iedb,jsd:jed)) ; g%mask2dCu(:,:) = 0.0
555  alloc_(g%mask2dCv(isd:ied,jsdb:jedb)) ; g%mask2dCv(:,:) = 0.0
556  alloc_(g%mask2dBu(isdb:iedb,jsdb:jedb)) ; g%mask2dBu(:,:) = 0.0
557  alloc_(g%geoLatT(isd:ied,jsd:jed)) ; g%geoLatT(:,:) = 0.0
558  alloc_(g%geoLatCu(isdb:iedb,jsd:jed)) ; g%geoLatCu(:,:) = 0.0
559  alloc_(g%geoLatCv(isd:ied,jsdb:jedb)) ; g%geoLatCv(:,:) = 0.0
560  alloc_(g%geoLatBu(isdb:iedb,jsdb:jedb)) ; g%geoLatBu(:,:) = 0.0
561  alloc_(g%geoLonT(isd:ied,jsd:jed)) ; g%geoLonT(:,:) = 0.0
562  alloc_(g%geoLonCu(isdb:iedb,jsd:jed)) ; g%geoLonCu(:,:) = 0.0
563  alloc_(g%geoLonCv(isd:ied,jsdb:jedb)) ; g%geoLonCv(:,:) = 0.0
564  alloc_(g%geoLonBu(isdb:iedb,jsdb:jedb)) ; g%geoLonBu(:,:) = 0.0
565 
566  alloc_(g%dx_Cv(isd:ied,jsdb:jedb)) ; g%dx_Cv(:,:) = 0.0
567  alloc_(g%dy_Cu(isdb:iedb,jsd:jed)) ; g%dy_Cu(:,:) = 0.0
568 
569  alloc_(g%areaCu(isdb:iedb,jsd:jed)) ; g%areaCu(:,:) = 0.0
570  alloc_(g%areaCv(isd:ied,jsdb:jedb)) ; g%areaCv(:,:) = 0.0
571  alloc_(g%IareaCu(isdb:iedb,jsd:jed)) ; g%IareaCu(:,:) = 0.0
572  alloc_(g%IareaCv(isd:ied,jsdb:jedb)) ; g%IareaCv(:,:) = 0.0
573 
574  alloc_(g%bathyT(isd:ied, jsd:jed)) ; g%bathyT(:,:) = 0.0
575  alloc_(g%CoriolisBu(isdb:iedb, jsdb:jedb)) ; g%CoriolisBu(:,:) = 0.0
576  alloc_(g%dF_dx(isd:ied, jsd:jed)) ; g%dF_dx(:,:) = 0.0
577  alloc_(g%dF_dy(isd:ied, jsd:jed)) ; g%dF_dy(:,:) = 0.0
578 
579  alloc_(g%sin_rot(isd:ied,jsd:jed)) ; g%sin_rot(:,:) = 0.0
580  alloc_(g%cos_rot(isd:ied,jsd:jed)) ; g%cos_rot(:,:) = 1.0
581 
582  allocate(g%gridLonT(isg:ieg)) ; g%gridLonT(:) = 0.0
583  allocate(g%gridLonB(g%IsgB:g%IegB)) ; g%gridLonB(:) = 0.0
584  allocate(g%gridLatT(jsg:jeg)) ; g%gridLatT(:) = 0.0
585  allocate(g%gridLatB(g%JsgB:g%JegB)) ; g%gridLatB(:) = 0.0
586 
587 end subroutine allocate_metrics
588 
589 !> Release memory used by the ocean_grid_type and related structures.
590 subroutine mom_grid_end(G)
591  type(ocean_grid_type), intent(inout) :: g !< The horizontal grid type
592 
593  dealloc_(g%dxT) ; dealloc_(g%dxCu) ; dealloc_(g%dxCv) ; dealloc_(g%dxBu)
594  dealloc_(g%IdxT) ; dealloc_(g%IdxCu) ; dealloc_(g%IdxCv) ; dealloc_(g%IdxBu)
595 
596  dealloc_(g%dyT) ; dealloc_(g%dyCu) ; dealloc_(g%dyCv) ; dealloc_(g%dyBu)
597  dealloc_(g%IdyT) ; dealloc_(g%IdyCu) ; dealloc_(g%IdyCv) ; dealloc_(g%IdyBu)
598 
599  dealloc_(g%areaT) ; dealloc_(g%IareaT)
600  dealloc_(g%areaBu) ; dealloc_(g%IareaBu)
601  dealloc_(g%areaCu) ; dealloc_(g%IareaCu)
602  dealloc_(g%areaCv) ; dealloc_(g%IareaCv)
603 
604  dealloc_(g%mask2dT) ; dealloc_(g%mask2dCu)
605  dealloc_(g%mask2dCv) ; dealloc_(g%mask2dBu)
606 
607  dealloc_(g%geoLatT) ; dealloc_(g%geoLatCu)
608  dealloc_(g%geoLatCv) ; dealloc_(g%geoLatBu)
609  dealloc_(g%geoLonT) ; dealloc_(g%geoLonCu)
610  dealloc_(g%geoLonCv) ; dealloc_(g%geoLonBu)
611 
612  dealloc_(g%dx_Cv) ; dealloc_(g%dy_Cu)
613 
614  dealloc_(g%bathyT) ; dealloc_(g%CoriolisBu)
615  dealloc_(g%dF_dx) ; dealloc_(g%dF_dy)
616  dealloc_(g%sin_rot) ; dealloc_(g%cos_rot)
617 
618  if (g%bathymetry_at_vel) then
619  dealloc_(g%Dblock_u) ; dealloc_(g%Dopen_u)
620  dealloc_(g%Dblock_v) ; dealloc_(g%Dopen_v)
621  endif
622 
623  deallocate(g%gridLonT) ; deallocate(g%gridLatT)
624  deallocate(g%gridLonB) ; deallocate(g%gridLatB)
625 
626  deallocate(g%Domain%mpp_domain)
627  deallocate(g%Domain)
628 
629 end subroutine mom_grid_end
630 
631 !> \namespace mom_grid
632 !!
633 !! Grid metrics and their inverses are labelled according to their staggered location on a Arakawa C (or B) grid.
634 !! - Metrics centered on h- or T-points are labelled T, e.g. dxT is the distance across the cell in the x-direction.
635 !! - Metrics centered on u-points are labelled Cu (C-grid u location). e.g. dyCu is the y-distance between
636 !! two corners of a T-cell.
637 !! - Metrics centered on v-points are labelled Cv (C-grid v location). e.g. dyCv is the y-distance between two -points.
638 !! - Metrics centered on q-points are labelled Bu (B-grid u,v location). e.g. areaBu is the area centered on a q-point.
639 !!
640 !! \image html Grid_metrics.png "The labelling of distances (grid metrics) at various staggered
641 !! location on an T-cell and around a q-point."
642 !!
643 !! Areas centered at T-, u-, v- and q- points are `areaT`, `areaCu`, `areaCv` and `areaBu` respectively.
644 !!
645 !! The reciprocal of metrics are pre-calculated and also stored in the ocean_grid_type with a I prepended to the name.
646 !! For example, `1./areaT` is called `IareaT`, and `1./dyCv` is `IdyCv`.
647 !!
648 !! Geographic latitude and longitude (or model coordinates if not on a sphere) are stored in
649 !! `geoLatT`, `geoLonT` for T-points.
650 !! u-, v- and q- point coordinates are follow same pattern of replacing T with Cu, Cv and Bu respectively.
651 !!
652 !! Each location also has a 2D mask indicating whether the entire column is land or ocean.
653 !! `mask2dT` is 1 if the column is wet or 0 if the T-cell is land.
654 !! `mask2dCu` is 1 if both neighboring column are ocean, and 0 if either is land.
655 
656 end module mom_grid
mom_domains::get_domain_extent
subroutine, public get_domain_extent(Domain, isc, iec, jsc, jec, isd, ied, jsd, jed, isg, ieg, jsg, jeg, idg_offset, jdg_offset, symmetric, local_indexing, index_offset)
Returns various data that has been stored in a MOM_domain_type.
Definition: MOM_domains.F90:1769
mom_grid::set_first_direction
subroutine, public set_first_direction(G, y_first)
Store an integer indicating which direction to work on first.
Definition: MOM_grid.F90:502
mom_grid::mom_grid_end
subroutine, public mom_grid_end(G)
Release memory used by the ocean_grid_type and related structures.
Definition: MOM_grid.F90:591
mom_grid::allocate_metrics
subroutine allocate_metrics(G)
Allocate memory used by the ocean_grid_type and related structures.
Definition: MOM_grid.F90:520
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_hor_index::hor_index_init
subroutine, public hor_index_init(Domain, HI, param_file, local_indexing, index_offset)
Sets various index values in a hor_index_type.
Definition: MOM_hor_index.F90:58
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_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_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_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_grid::mom_grid_init
subroutine, public mom_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_vel)
MOM_grid_init initializes the ocean grid array sizes and grid memory.
Definition: MOM_grid.F90:184
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_grid::get_global_grid_size
subroutine, public get_global_grid_size(G, niglobal, njglobal)
Return global shape of horizontal grid.
Definition: MOM_grid.F90:510
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:15
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
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_domains::get_global_shape
subroutine, public get_global_shape(domain, niglobal, njglobal)
Returns the global shape of h-point arrays.
Definition: MOM_domains.F90:1927
mom_grid::rescale_grid_bathymetry
subroutine, public rescale_grid_bathymetry(G, m_in_new_units)
rescale_grid_bathymetry permits a change in the internal units for the bathymetry on the grid,...
Definition: MOM_grid.F90:381
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::set_derived_metrics
subroutine, public set_derived_metrics(G, US)
set_derived_metrics calculates metric terms that are derived from other metrics.
Definition: MOM_grid.F90:413
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ispointincell
logical function, public ispointincell(G, i, j, x, y)
Returns true if the coordinates (x,y) are within the h-cell (i,j)
Definition: MOM_grid.F90:467
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_domains::get_domain_extent_dsamp2
subroutine, public get_domain_extent_dsamp2(Domain, isc_d2, iec_d2, jsc_d2, jec_d2, isd_d2, ied_d2, jsd_d2, jed_d2, isg_d2, ieg_d2, jsg_d2, jeg_d2)
Definition: MOM_domains.F90:1827
mom_grid::adcroft_reciprocal
real function adcroft_reciprocal(val)
Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0.
Definition: MOM_grid.F90:459