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
18 use mom_io,
only : corner, north_face, east_face
21 use mpp_domains_mod,
only : mpp_get_domain_extents, mpp_deallocate_domain
23 implicit none ;
private
25 public set_grid_metrics, initialize_masks, adcroft_reciprocal
33 type,
public ::
gps ;
private
41 real :: lat_enhance_factor
43 real :: lat_eq_enhance
51 logical :: equator_reference
62 subroutine set_grid_metrics(G, param_file, US)
68 #include "version_variable.h"
70 character(len=256) :: config
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.)
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)
91 case (
"cartesian");
call set_grid_metrics_cartesian(g, param_file)
92 case (
"spherical");
call set_grid_metrics_spherical(g, param_file)
93 case (
"mercator");
call set_grid_metrics_mercator(g, param_file)
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))
102 call calltree_enter(
"set_derived_metrics(), MOM_grid_initialize.F90")
103 call set_derived_dyn_horgrid(g)
104 call calltree_leave(
"set_derived_metrics()")
106 if (debug)
call grid_metrics_chksum(
'MOM_grid_init/set_grid_metrics',g)
108 call calltree_leave(
"set_grid_metrics()")
109 end subroutine set_grid_metrics
115 subroutine grid_metrics_chksum(parent, G)
116 character(len=*),
intent(in) :: parent
121 halo = min(g%ied-g%iec, g%jed-g%jec, 1)
124 g%dxT, g%dyT, g%HI, haloshift=halo)
126 call uvchksum(trim(parent)//
': dxC[uv]', g%dxCu, g%dyCv, g%HI, haloshift=halo)
128 call uvchksum(trim(parent)//
': dxC[uv]', &
129 g%dyCu, g%dxCv, g%HI, haloshift=halo)
132 g%dxBu, g%dyBu, g%HI, haloshift=halo)
135 g%IdxT, g%IdyT, g%HI, haloshift=halo)
137 call uvchksum(trim(parent)//
': Id[xy]C[uv]', &
138 g%IdxCu, g%IdyCv, g%HI, haloshift=halo)
140 call uvchksum(trim(parent)//
': Id[xy]C[uv]', &
141 g%IdyCu, g%IdxCv, g%HI, haloshift=halo)
144 g%IdxBu, g%IdyBu, g%HI, haloshift=halo)
146 call hchksum(g%areaT, trim(parent)//
': areaT',g%HI, haloshift=halo)
147 call bchksum(g%areaBu, trim(parent)//
': areaBu',g%HI, haloshift=halo)
149 call hchksum(g%IareaT, trim(parent)//
': IareaT',g%HI, haloshift=halo)
150 call bchksum(g%IareaBu, trim(parent)//
': IareaBu',g%HI, haloshift=halo)
152 call hchksum(g%geoLonT,trim(parent)//
': geoLonT',g%HI, haloshift=halo)
153 call hchksum(g%geoLatT,trim(parent)//
': geoLatT',g%HI, haloshift=halo)
155 call bchksum(g%geoLonBu, trim(parent)//
': geoLonBu',g%HI, haloshift=halo)
156 call bchksum(g%geoLatBu, trim(parent)//
': geoLatBu',g%HI, haloshift=halo)
158 call uvchksum(trim(parent)//
': geoLonC[uv]', &
159 g%geoLonCu, g%geoLonCv, g%HI, haloshift=halo)
161 call uvchksum(trim(parent)//
': geoLatC[uv]', &
162 g%geoLatCu, g%geoLatCv, g%HI, haloshift=halo)
164 end subroutine grid_metrics_chksum
169 subroutine set_grid_metrics_from_mosaic(G, param_file)
173 real,
dimension(G%isd :G%ied ,G%jsd :G%jed ) :: tempH1, tempH2, tempH3, tempH4
174 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: tempQ1, tempQ2, tempQ3, tempQ4
175 real,
dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: tempE1, tempE2
176 real,
dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: tempN1, tempN2
179 real,
dimension(G%isd :G%ied ,G%jsd :G%jed ) :: dxT, dyT, areaT
180 real,
dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: dxCu, dyCu
181 real,
dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: dxCv, dyCv
182 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: dxBu, dyBu, areaBu
184 real,
dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpT
185 real,
dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpU
186 real,
dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpV
187 real,
dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpZ
188 real,
dimension(:,:),
allocatable :: tmpGlbl
189 character(len=200) :: filename, grid_file, inputdir
190 character(len=64) :: mdl =
"MOM_grid_init set_grid_metrics_from_mosaic"
191 integer :: err=0, ni, nj, global_indices(4)
194 integer :: i, j, i2, j2
196 integer,
dimension(:),
allocatable :: exni,exnj
197 integer :: start(4), nread(4)
199 call calltree_enter(
"set_grid_metrics_from_mosaic(), MOM_grid_initialize.F90")
201 call get_param(param_file, mdl,
"GRID_FILE", grid_file, &
202 "Name of the file from which to read horizontal grid data.", &
203 fail_if_missing=.true.)
204 call get_param(param_file, mdl,
"USE_TRIPOLAR_GEOLONB_BUG", lon_bug, &
205 "If true, use older code that incorrectly sets the longitude "//&
206 "in some points along the tripolar fold to be off by 360 degrees.", &
208 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
209 inputdir = slasher(inputdir)
210 filename = trim(adjustl(inputdir)) // trim(adjustl(grid_file))
211 call log_param(param_file, mdl,
"INPUTDIR/GRID_FILE", filename)
213 call mom_error(fatal,
" set_grid_metrics_from_mosaic: Unable to open "//&
217 dxcu(:,:) = 0.0 ; dycu(:,:) = 0.0
218 dxcv(:,:) = 0.0 ; dycv(:,:) = 0.0
219 dxbu(:,:) = 0.0 ; dybu(:,:) = 0.0 ; areabu(:,:) = 0.0
222 ni = 2*(g%iec-g%isc+1)
223 nj = 2*(g%jec-g%jsc+1)
226 npei = g%domain%layout(1) ; npej = g%domain%layout(2)
227 allocate(exni(npei)) ;
allocate(exnj(npej))
228 call mpp_get_domain_extents(g%domain%mpp_domain, exni, exnj)
229 allocate(sgdom%mpp_domain)
230 sgdom%nihalo = 2*g%domain%nihalo+1
231 sgdom%njhalo = 2*g%domain%njhalo+1
232 sgdom%niglobal = 2*g%domain%niglobal
233 sgdom%njglobal = 2*g%domain%njglobal
234 sgdom%layout(:) = g%domain%layout(:)
235 sgdom%io_layout(:) = g%domain%io_layout(:)
236 global_indices(1) = 1+sgdom%nihalo
237 global_indices(2) = sgdom%niglobal+sgdom%nihalo
238 global_indices(3) = 1+sgdom%njhalo
239 global_indices(4) = sgdom%njglobal+sgdom%njhalo
240 exni(:) = 2*exni(:) ; exnj(:) = 2*exnj(:)
241 if (
associated(g%domain%maskmap))
then
242 call mom_define_domain(global_indices, sgdom%layout, sgdom%mpp_domain, &
243 xflags=g%domain%X_FLAGS, yflags=g%domain%Y_FLAGS, &
244 xhalo=sgdom%nihalo, yhalo=sgdom%njhalo, &
245 xextent=exni,yextent=exnj, &
246 symmetry=.true., name=
"MOM_MOSAIC", maskmap=g%domain%maskmap)
248 call mom_define_domain(global_indices, sgdom%layout, sgdom%mpp_domain, &
249 xflags=g%domain%X_FLAGS, yflags=g%domain%Y_FLAGS, &
250 xhalo=sgdom%nihalo, yhalo=sgdom%njhalo, &
251 xextent=exni,yextent=exnj, &
252 symmetry=.true., name=
"MOM_MOSAIC")
255 call mom_define_io_domain(sgdom%mpp_domain, sgdom%io_layout)
261 call mom_read_data(filename,
'x', tmpz, sgdom, position=corner)
264 call pass_var(tmpz, sgdom, position=corner)
266 call pass_var(tmpz, sgdom, position=corner, inner_halo=0)
268 call extrapolate_metric(tmpz, 2*(g%jsc-g%jsd)+2, missing=999.)
269 do j=g%jsd,g%jed ;
do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
270 g%geoLonT(i,j) = tmpz(i2-1,j2-1)
272 do j=g%JsdB,g%JedB ;
do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
273 g%geoLonBu(i,j) = tmpz(i2,j2)
275 do j=g%jsd,g%jed ;
do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
276 g%geoLonCu(i,j) = tmpz(i2,j2-1)
278 do j=g%JsdB,g%JedB ;
do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
279 g%geoLonCv(i,j) = tmpz(i2-1,j2)
286 call mom_read_data(filename,
'y', tmpz, sgdom, position=corner)
288 call pass_var(tmpz, sgdom, position=corner)
289 call extrapolate_metric(tmpz, 2*(g%jsc-g%jsd)+2, missing=999.)
290 do j=g%jsd,g%jed ;
do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
291 g%geoLatT(i,j) = tmpz(i2-1,j2-1)
293 do j=g%JsdB,g%JedB ;
do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
294 g%geoLatBu(i,j) = tmpz(i2,j2)
296 do j=g%jsd,g%jed ;
do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
297 g%geoLatCu(i,j) = tmpz(i2,j2-1)
299 do j=g%JsdB,g%JedB ;
do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
300 g%geoLatCv(i,j) = tmpz(i2-1,j2)
304 tmpu(:,:) = 0. ; tmpv(:,:) = 0.
305 call mom_read_data(filename,
'dx',tmpv,sgdom,position=north_face)
306 call mom_read_data(filename,
'dy',tmpu,sgdom,position=east_face)
307 call pass_vector(tmpu, tmpv, sgdom, to_all+scalar_pair, cgrid_ne)
308 call extrapolate_metric(tmpv, 2*(g%jsc-g%jsd)+2, missing=0.)
309 call extrapolate_metric(tmpu, 2*(g%jsc-g%jsd)+2, missing=0.)
311 do j=g%jsd,g%jed ;
do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
312 dxt(i,j) = tmpv(i2-1,j2-1) + tmpv(i2,j2-1)
313 dyt(i,j) = tmpu(i2-1,j2-1) + tmpu(i2-1,j2)
316 do j=g%jsd,g%jed ;
do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
317 dxcu(i,j) = tmpv(i2,j2-1) + tmpv(i2+1,j2-1)
318 dycu(i,j) = tmpu(i2,j2-1) + tmpu(i2,j2)
321 do j=g%JsdB,g%JedB ;
do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
322 dxcv(i,j) = tmpv(i2-1,j2) + tmpv(i2,j2)
323 dycv(i,j) = tmpu(i2-1,j2) + tmpu(i2-1,j2+1)
326 do j=g%JsdB,g%JedB ;
do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
327 dxbu(i,j) = tmpv(i2,j2) + tmpv(i2+1,j2)
328 dybu(i,j) = tmpu(i2,j2) + tmpu(i2,j2+1)
335 call extrapolate_metric(tmpt, 2*(g%jsc-g%jsd)+2, missing=0.)
337 do j=g%jsd,g%jed ;
do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
338 areat(i,j) = (tmpt(i2-1,j2-1) + tmpt(i2,j2)) + &
339 (tmpt(i2-1,j2) + tmpt(i2,j2-1))
341 do j=g%JsdB,g%JedB ;
do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
342 areabu(i,j) = (tmpt(i2,j2) + tmpt(i2+1,j2+1)) + &
343 (tmpt(i2,j2+1) + tmpt(i2+1,j2))
348 call mpp_deallocate_domain(sgdom%mpp_domain)
349 deallocate(sgdom%mpp_domain)
351 call pass_vector(dycu, dxcv, g%Domain, to_all+scalar_pair, cgrid_ne)
352 call pass_vector(dxcu, dycv, g%Domain, to_all+scalar_pair, cgrid_ne)
353 call pass_vector(dxbu, dybu, g%Domain, to_all+scalar_pair, bgrid_ne)
355 call pass_var(areabu, g%Domain, position=corner)
357 do i=g%isd,g%ied ;
do j=g%jsd,g%jed
358 g%dxT(i,j) = dxt(i,j) ; g%dyT(i,j) = dyt(i,j) ; g%areaT(i,j) = areat(i,j)
360 do i=g%IsdB,g%IedB ;
do j=g%jsd,g%jed
361 g%dxCu(i,j) = dxcu(i,j) ; g%dyCu(i,j) = dycu(i,j)
363 do i=g%isd,g%ied ;
do j=g%JsdB,g%JedB
364 g%dxCv(i,j) = dxcv(i,j) ; g%dyCv(i,j) = dycv(i,j)
366 do i=g%IsdB,g%IedB ;
do j=g%JsdB,g%JedB
367 g%dxBu(i,j) = dxbu(i,j) ; g%dyBu(i,j) = dybu(i,j) ; g%areaBu(i,j) = areabu(i,j)
372 start(:) = 1 ; nread(:) = 1
373 start(2) = 2 ; nread(1) = ni+1 ; nread(2) = 2
374 allocate( tmpglbl(ni+1,2) )
376 call read_data(filename,
"x", tmpglbl, start, nread, no_domain=.true.)
377 call broadcast(tmpglbl, 2*(ni+1), root_pe())
381 g%gridLonT(i) = tmpglbl(2*(i-g%isg)+2,2)
386 g%gridLonB(i) = tmpglbl(2*(i-g%isg)+3,1)
388 deallocate( tmpglbl )
390 allocate( tmpglbl(1, nj+1) )
391 start(:) = 1 ; nread(:) = 1
392 start(1) = int(ni/4)+1 ; nread(2) = nj+1
394 call read_data(filename,
"y", tmpglbl, start, nread, no_domain=.true.)
395 call broadcast(tmpglbl, nj+1, root_pe())
398 g%gridLatT(j) = tmpglbl(1,2*(j-g%jsg)+2)
401 g%gridLatB(j) = tmpglbl(1,2*(j-g%jsg)+3)
403 deallocate( tmpglbl )
405 call calltree_leave(
"set_grid_metrics_from_mosaic()")
406 end subroutine set_grid_metrics_from_mosaic
418 subroutine set_grid_metrics_cartesian(G, param_file)
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
429 character(len=80) :: units_temp
430 character(len=48) :: mdl =
"MOM_grid_init set_grid_metrics_cartesian"
432 niglobal = g%Domain%niglobal ; njglobal = g%Domain%njglobal
433 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
434 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
435 i1off = g%idg_offset ; j1off = g%jdg_offset
437 call calltree_enter(
"set_grid_metrics_cartesian(), MOM_grid_initialize.F90")
441 call get_param(param_file, mdl,
"AXIS_UNITS", units_temp, &
442 "The units for the Cartesian axes. Valid entries are: \n"//&
443 " \t degrees - degrees of latitude and longitude \n"//&
444 " \t m - meters \n \t k - kilometers", default=
"degrees")
445 call get_param(param_file, mdl,
"SOUTHLAT", g%south_lat, &
446 "The southern latitude of the domain or the equivalent "//&
447 "starting value for the y-axis.", units=units_temp, &
448 fail_if_missing=.true.)
449 call get_param(param_file, mdl,
"LENLAT", g%len_lat, &
450 "The latitudinal or y-direction length of the domain.", &
451 units=units_temp, fail_if_missing=.true.)
452 call get_param(param_file, mdl,
"WESTLON", g%west_lon, &
453 "The western longitude of the domain or the equivalent "//&
454 "starting value for the x-axis.", units=units_temp, &
456 call get_param(param_file, mdl,
"LENLON", g%len_lon, &
457 "The longitudinal or x-direction length of the domain.", &
458 units=units_temp, fail_if_missing=.true.)
459 call get_param(param_file, mdl,
"RAD_EARTH", g%Rad_Earth, &
460 "The radius of the Earth.", units=
"m", default=6.378e6)
462 if (units_temp(1:1) ==
'k')
then
463 g%x_axis_units =
"kilometers" ; g%y_axis_units =
"kilometers"
464 elseif (units_temp(1:1) ==
'm')
then
465 g%x_axis_units =
"meters" ; g%y_axis_units =
"meters"
467 call log_param(param_file, mdl,
"explicit AXIS_UNITS", g%x_axis_units)
472 g%gridLatB(j) = g%south_lat + g%len_lat*real(j-(g%jsg-1))/real(njglobal)
475 g%gridLatT(j) = g%south_lat + g%len_lat*(real(j-g%jsg)+0.5)/real(njglobal)
478 g%gridLonB(i) = g%west_lon + g%len_lon*real(i-(g%isg-1))/real(niglobal)
481 g%gridLonT(i) = g%west_lon + g%len_lon*(real(i-g%isg)+0.5)/real(niglobal)
485 grid_latb(j) = g%south_lat + g%len_lat*real(j+j1off-(g%jsg-1))/real(njglobal)
488 grid_latt(j) = g%south_lat + g%len_lat*(real(j+j1off-g%jsg)+0.5)/real(njglobal)
491 grid_lonb(i) = g%west_lon + g%len_lon*real(i+i1off-(g%isg-1))/real(niglobal)
494 grid_lont(i) = g%west_lon + g%len_lon*(real(i+i1off-g%isg)+0.5)/real(niglobal)
497 if (units_temp(1:1) ==
'k')
then
498 dx_everywhere = 1000.0 * g%len_lon / (real(niglobal))
499 dy_everywhere = 1000.0 * g%len_lat / (real(njglobal))
500 elseif (units_temp(1:1) ==
'm')
then
501 dx_everywhere = g%len_lon / (real(niglobal))
502 dy_everywhere = g%len_lat / (real(njglobal))
504 dx_everywhere = g%Rad_Earth * g%len_lon * pi / (180.0 * niglobal)
505 dy_everywhere = g%Rad_Earth * g%len_lat * pi / (180.0 * njglobal)
508 i_dx = 1.0 / dx_everywhere ; i_dy = 1.0 / dy_everywhere
510 do j=jsdb,jedb ;
do i=isdb,iedb
511 g%geoLonBu(i,j) = grid_lonb(i) ; g%geoLatBu(i,j) = grid_latb(j)
513 g%dxBu(i,j) = dx_everywhere ; g%IdxBu(i,j) = i_dx
514 g%dyBu(i,j) = dy_everywhere ; g%IdyBu(i,j) = i_dy
515 g%areaBu(i,j) = dx_everywhere * dy_everywhere ; g%IareaBu(i,j) = i_dx * i_dy
518 do j=jsd,jed ;
do i=isd,ied
519 g%geoLonT(i,j) = grid_lont(i) ; g%geoLatT(i,j) = grid_latt(j)
520 g%dxT(i,j) = dx_everywhere ; g%IdxT(i,j) = i_dx
521 g%dyT(i,j) = dy_everywhere ; g%IdyT(i,j) = i_dy
522 g%areaT(i,j) = dx_everywhere * dy_everywhere ; g%IareaT(i,j) = i_dx * i_dy
525 do j=jsd,jed ;
do i=isdb,iedb
526 g%geoLonCu(i,j) = grid_lonb(i) ; g%geoLatCu(i,j) = grid_latt(j)
528 g%dxCu(i,j) = dx_everywhere ; g%IdxCu(i,j) = i_dx
529 g%dyCu(i,j) = dy_everywhere ; g%IdyCu(i,j) = i_dy
532 do j=jsdb,jedb ;
do i=isd,ied
533 g%geoLonCv(i,j) = grid_lont(i) ; g%geoLatCv(i,j) = grid_latb(j)
535 g%dxCv(i,j) = dx_everywhere ; g%IdxCv(i,j) = i_dx
536 g%dyCv(i,j) = dy_everywhere ; g%IdyCv(i,j) = i_dy
539 call calltree_leave(
"set_grid_metrics_cartesian()")
540 end subroutine set_grid_metrics_cartesian
551 subroutine set_grid_metrics_spherical(G, param_file)
556 integer :: i, j, isd, ied, jsd, jed
557 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
558 integer :: i_offset, j_offset
559 real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB)
560 real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB)
561 real :: dLon,dLat,latitude,longitude,dL_di
562 character(len=48) :: mdl =
"MOM_grid_init set_grid_metrics_spherical"
564 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
565 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
566 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
567 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
568 i_offset = g%idg_offset ; j_offset = g%jdg_offset
570 call calltree_enter(
"set_grid_metrics_spherical(), MOM_grid_initialize.F90")
574 pi = 4.0*atan(1.0); pi_180 = atan(1.0)/45.
576 call get_param(param_file, mdl,
"SOUTHLAT", g%south_lat, &
577 "The southern latitude of the domain.", units=
"degrees", &
578 fail_if_missing=.true.)
579 call get_param(param_file, mdl,
"LENLAT", g%len_lat, &
580 "The latitudinal length of the domain.", units=
"degrees", &
581 fail_if_missing=.true.)
582 call get_param(param_file, mdl,
"WESTLON", g%west_lon, &
583 "The western longitude of the domain.", units=
"degrees", &
585 call get_param(param_file, mdl,
"LENLON", g%len_lon, &
586 "The longitudinal length of the domain.", units=
"degrees", &
587 fail_if_missing=.true.)
588 call get_param(param_file, mdl,
"RAD_EARTH", g%Rad_Earth, &
589 "The radius of the Earth.", units=
"m", default=6.378e6)
591 dlon = g%len_lon/g%Domain%niglobal
592 dlat = g%len_lat/g%Domain%njglobal
597 latitude = g%south_lat + dlat*(real(j-(g%jsg-1)))
598 g%gridLatB(j) = min(max(latitude,-90.),90.)
601 latitude = g%south_lat + dlat*(real(j-g%jsg)+0.5)
602 g%gridLatT(j) = min(max(latitude,-90.),90.)
605 g%gridLonB(i) = g%west_lon + dlon*(real(i-(g%isg-1)))
608 g%gridLonT(i) = g%west_lon + dlon*(real(i-g%isg)+0.5)
612 latitude = g%south_lat + dlat* real(j+j_offset-(g%jsg-1))
613 grid_latb(j) = min(max(latitude,-90.),90.)
616 latitude = g%south_lat + dlat*(real(j+j_offset-g%jsg)+0.5)
617 grid_latt(j) = min(max(latitude,-90.),90.)
620 grid_lonb(i) = g%west_lon + dlon*real(i+i_offset-(g%isg-1))
623 grid_lont(i) = g%west_lon + dlon*(real(i+i_offset-g%isg)+0.5)
626 dl_di = (g%len_lon * 4.0*atan(1.0)) / (180.0 * g%Domain%niglobal)
627 do j=jsdb,jedb ;
do i=isdb,iedb
628 g%geoLonBu(i,j) = grid_lonb(i)
629 g%geoLatBu(i,j) = grid_latb(j)
633 g%dxBu(i,j) = g%Rad_Earth * cos( g%geoLatBu(i,j)*pi_180 ) * dl_di
635 g%dyBu(i,j) = g%Rad_Earth * dlat*pi_180
636 g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
639 do j=jsdb,jedb ;
do i=isd,ied
640 g%geoLonCv(i,j) = grid_lont(i)
641 g%geoLatCv(i,j) = grid_latb(j)
645 g%dxCv(i,j) = g%Rad_Earth * cos( g%geoLatCv(i,j)*pi_180 ) * dl_di
647 g%dyCv(i,j) = g%Rad_Earth * dlat*pi_180
650 do j=jsd,jed ;
do i=isdb,iedb
651 g%geoLonCu(i,j) = grid_lonb(i)
652 g%geoLatCu(i,j) = grid_latt(j)
656 g%dxCu(i,j) = g%Rad_Earth * cos( g%geoLatCu(i,j)*pi_180 ) * dl_di
658 g%dyCu(i,j) = g%Rad_Earth * dlat*pi_180
661 do j=jsd,jed ;
do i=isd,ied
662 g%geoLonT(i,j) = grid_lont(i)
663 g%geoLatT(i,j) = grid_latt(j)
667 g%dxT(i,j) = g%Rad_Earth * cos( g%geoLatT(i,j)*pi_180 ) * dl_di
669 g%dyT(i,j) = g%Rad_Earth * dlat*pi_180
674 g%areaT(i,j) = g%dxT(i,j) * g%dyT(i,j)
677 call calltree_leave(
"set_grid_metrics_spherical()")
678 end subroutine set_grid_metrics_spherical
687 subroutine set_grid_metrics_mercator(G, param_file)
691 integer :: i, j, isd, ied, jsd, jed
692 integer :: I_off, J_off
694 character(len=128) :: warnmesg
695 character(len=48) :: mdl =
"MOM_grid_init set_grid_metrics_mercator"
697 real :: y_q, y_h, jd, x_q, x_h, id
698 real,
dimension(G%isd:G%ied,G%jsd:G%jed) :: &
700 real,
dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: &
702 real,
dimension(G%isd:G%ied,G%JsdB:G%JedB) :: &
704 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: &
709 integer :: itt1, itt2
710 logical :: debug = .false., simple_area = .true.
711 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
716 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
717 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
718 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
719 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
720 i_off = g%idg_offset ; j_off = g%jdg_offset
722 gp%niglobal = g%Domain%niglobal
723 gp%njglobal = g%Domain%njglobal
725 call calltree_enter(
"set_grid_metrics_mercator(), MOM_grid_initialize.F90")
729 pi = 4.0*atan(1.0) ; pi_2 = 0.5*pi
731 call get_param(param_file, mdl,
"SOUTHLAT", gp%south_lat, &
732 "The southern latitude of the domain.", units=
"degrees", &
733 fail_if_missing=.true.)
734 call get_param(param_file, mdl,
"LENLAT", gp%len_lat, &
735 "The latitudinal length of the domain.", units=
"degrees", &
736 fail_if_missing=.true.)
737 call get_param(param_file, mdl,
"WESTLON", gp%west_lon, &
738 "The western longitude of the domain.", units=
"degrees", &
740 call get_param(param_file, mdl,
"LENLON", gp%len_lon, &
741 "The longitudinal length of the domain.", units=
"degrees", &
742 fail_if_missing=.true.)
743 call get_param(param_file, mdl,
"RAD_EARTH", gp%Rad_Earth, &
744 "The radius of the Earth.", units=
"m", default=6.378e6)
745 g%south_lat = gp%south_lat ; g%len_lat = gp%len_lat
746 g%west_lon = gp%west_lon ; g%len_lon = gp%len_lon
747 g%Rad_Earth = gp%Rad_Earth
748 call get_param(param_file, mdl,
"ISOTROPIC", gp%isotropic, &
749 "If true, an isotropic grid on a sphere (also known as "//&
750 "a Mercator grid) is used. With an isotropic grid, the "//&
751 "meridional extent of the domain (LENLAT), the zonal "//&
752 "extent (LENLON), and the number of grid points in each "//&
753 "direction are _not_ independent. In MOM the meridional "//&
754 "extent is determined to fit the zonal extent and the "//&
755 "number of grid points, while grid is perfectly isotropic.", &
757 call get_param(param_file, mdl,
"EQUATOR_REFERENCE", gp%equator_reference, &
758 "If true, the grid is defined to have the equator at the "//&
759 "nearest q- or h- grid point to (-LOWLAT*NJGLOBAL/LENLAT).", &
761 call get_param(param_file, mdl,
"LAT_ENHANCE_FACTOR", gp%Lat_enhance_factor, &
762 "The amount by which the meridional resolution is "//&
763 "enhanced within LAT_EQ_ENHANCE of the equator.", &
764 units=
"nondim", default=1.0)
765 call get_param(param_file, mdl,
"LAT_EQ_ENHANCE", gp%Lat_eq_enhance, &
766 "The latitude range to the north and south of the equator "//&
767 "over which the resolution is enhanced.", units=
"degrees", &
775 if (gp%equator_reference)
then
779 jref = (g%jsg-1) + 0.5*floor(gp%njglobal*((-1.0*gp%south_lat*2.0)/gp%len_lat)+0.5)
780 fnref = int_dj_dy(0.0, gp)
784 fnref = int_dj_dy((gp%south_lat*pi/180.0), gp)
793 jd = fnref + (j - jref)
794 y_q = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt2)
795 g%gridLatB(j) = y_q*180.0/pi
800 jd = fnref + (j - jref) - 0.5
801 y_h = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt1)
802 g%gridLatT(j) = y_h*180.0/pi
806 do j=jsdb+j_off,jedb+j_off
807 jd = fnref + (j - jref)
808 y_q = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt2)
809 do i=isdb,iedb ; yq(i,j-j_off) = y_q ;
enddo
810 do i=isd,ied ; yv(i,j-j_off) = y_q ;
enddo
812 do j=jsd+j_off,jed+j_off
813 jd = fnref + (j - jref) - 0.5
814 y_h = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt1)
815 if ((j >= jsd+j_off) .and. (j <= jed+j_off))
then
816 do i=isd,ied ; yh(i,j-j_off) = y_h ;
enddo
817 do i=isdb,iedb ; yu(i,j-j_off) = y_h ;
enddo
824 iref = (g%isg-1) + gp%niglobal
825 fnref = int_di_dx(((gp%west_lon+gp%len_lon)*pi/180.0), gp)
831 id = fnref + (i - iref)
832 x_q = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt2)
833 g%gridLonB(i) = x_q*180.0/pi
836 id = fnref + (i - iref) - 0.5
837 x_h = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt1)
838 g%gridLonT(i) = x_h*180.0/pi
840 do i=isdb+i_off,iedb+i_off
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 do j=jsdb,jedb ; xq(i-i_off,j) = x_q ;
enddo
844 do j=jsd,jed ; xu(i-i_off,j) = x_q ;
enddo
846 do i=isd+i_off,ied+i_off
847 id = fnref + (i - iref) - 0.5
848 x_h = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt1)
849 do j=jsd,jed ; xh(i-i_off,j) = x_h ;
enddo
850 do j=jsdb,jedb ; xv(i-i_off,j) = x_h ;
enddo
853 do j=jsdb,jedb ;
do i=isdb,iedb
854 g%geoLonBu(i,j) = xq(i,j)*180.0/pi
855 g%geoLatBu(i,j) = yq(i,j)*180.0/pi
856 g%dxBu(i,j) = ds_di(xq(i,j), yq(i,j), gp)
857 g%dyBu(i,j) = ds_dj(xq(i,j), yq(i,j), gp)
859 g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
860 g%IareaBu(i,j) = 1.0 / g%areaBu(i,j)
863 do j=jsd,jed ;
do i=isd,ied
864 g%geoLonT(i,j) = xh(i,j)*180.0/pi
865 g%geoLatT(i,j) = yh(i,j)*180.0/pi
866 g%dxT(i,j) = ds_di(xh(i,j), yh(i,j), gp)
867 g%dyT(i,j) = ds_dj(xh(i,j), yh(i,j), gp)
869 g%areaT(i,j) = g%dxT(i,j)*g%dyT(i,j)
870 g%IareaT(i,j) = 1.0 / g%areaT(i,j)
873 do j=jsd,jed ;
do i=isdb,iedb
874 g%geoLonCu(i,j) = xu(i,j)*180.0/pi
875 g%geoLatCu(i,j) = yu(i,j)*180.0/pi
876 g%dxCu(i,j) = ds_di(xu(i,j), yu(i,j), gp)
877 g%dyCu(i,j) = ds_dj(xu(i,j), yu(i,j), gp)
880 do j=jsdb,jedb ;
do i=isd,ied
881 g%geoLonCv(i,j) = xv(i,j)*180.0/pi
882 g%geoLatCv(i,j) = yv(i,j)*180.0/pi
883 g%dxCv(i,j) = ds_di(xv(i,j), yv(i,j), gp)
884 g%dyCv(i,j) = ds_dj(xv(i,j), yv(i,j), gp)
887 if (.not.simple_area)
then
888 do j=jsdb+1,jed ;
do i=isdb+1,ied
889 g%areaT(i,j) = gp%Rad_Earth**2 * &
890 (dl(xq(i-1,j-1),xq(i-1,j),yq(i-1,j-1),yq(i-1,j)) + &
891 (dl(xq(i-1,j),xq(i,j),yq(i-1,j),yq(i,j)) + &
892 (dl(xq(i,j),xq(i,j-1),yq(i,j),yq(i,j-1)) + &
893 dl(xq(i,j-1),xq(i-1,j-1),yq(i,j-1),yq(i-1,j-1)))))
895 if ((isdb == isd) .or. (jsdb == jsq))
then
899 g%areaT(isd+1,jsd) = g%areaT(isd+1,jsd+1)
900 do j=jsd,jed ; g%areaT(isd,j) = g%areaT(isd+1,j) ;
enddo
901 do i=isd,ied ; g%areaT(i,jsd) = g%areaT(i,jsd+1) ;
enddo
905 do j=jsd,jed ;
do i=isd,ied
906 g%IareaT(i,j) = 1.0 / g%areaT(i,j)
910 call calltree_leave(
"set_grid_metrics_mercator()")
911 end subroutine set_grid_metrics_mercator
915 function ds_di(x, y, GP)
916 real,
intent(in) :: x
917 real,
intent(in) :: y
918 type(
gps),
intent(in) :: gp
922 ds_di = gp%Rad_Earth * cos(y) * dx_di(x,gp)
929 function ds_dj(x, y, GP)
930 real,
intent(in) :: x
931 real,
intent(in) :: y
932 type(
gps),
intent(in) :: gp
936 ds_dj = gp%Rad_Earth * dy_dj(y,gp)
945 function dl(x1, x2, y1, y2)
946 real,
intent(in) :: x1
947 real,
intent(in) :: x2
948 real,
intent(in) :: y1
949 real,
intent(in) :: y2
956 if (abs(dy) > 2.5e-8)
then
957 r = ((1.0 - cos(dy))*cos(y1) + sin(dy)*sin(y1)) / dy
959 r = (0.5*dy*cos(y1) + sin(y1))
968 function find_root( fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax)
971 real,
external :: dy_df
972 type(
gps),
intent(in) :: gp
973 real,
intent(in) :: fnval
974 real,
intent(in) :: y1
975 real,
intent(in) :: ymin
976 real,
intent(in) :: ymax
977 integer,
intent(out) :: ittmax
980 real :: ybot, ytop, fnbot, fntop
982 character(len=256) :: warnmesg
984 real :: dy_dfn, dy, fny
990 fnbot = fn(ybot,gp) - fnval ; itt = 0
991 do while (fnbot > 0.0)
992 if ((ybot - 2.0*dy_df(ybot,gp)) < (0.5*(ybot+ymin)))
then
994 ybot = ybot - 2.0*dy_df(ybot,gp)
996 ybot = 0.5*(ybot+ymin) ; itt = itt + 1
998 fnbot = fn(ybot,gp) - fnval
1000 if ((itt > 50) .and. (fnbot > 0.0))
then
1001 write(warnmesg,
'("PE ",I2," unable to find bottom bound for grid function. &
1002 &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4,&
1003 &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') &
1004 pe_here(),ybot,ymin,fn(ybot,gp),dy_df(ybot,gp),fnval, fnbot
1005 call mom_error(fatal,warnmesg)
1010 fntop = fn(ytop,gp) - fnval ; itt = 0
1011 do while (fntop < 0.0)
1012 if ((ytop + 2.0*dy_df(ytop,gp)) < (0.5*(ytop+ymax)))
then
1014 ytop = ytop + 2.0*dy_df(ytop,gp)
1016 ytop = 0.5*(ytop+ymax) ; itt = itt + 1
1018 fntop = fn(ytop,gp) - fnval
1020 if ((itt > 50) .and. (fntop < 0.0))
then
1021 write(warnmesg,
'("PE ",I2," unable to find top bound for grid function. &
1022 &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4, &
1023 &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') &
1024 pe_here(),ytop,ymax,fn(ytop,gp),dy_df(ytop,gp),fnval,fntop
1025 call mom_error(fatal,warnmesg)
1031 if ((fntop < 0.0) .or. (fnbot > 0.0) .or. (ytop < ybot))
then
1032 write(warnmesg,
'("PE ",I2," find_root failed to bracket function. y = ",&
1033 &2ES10.4,", fn = ",2ES10.4,".")') pe_here(),ybot,ytop,fnbot,fntop
1034 call mom_error(fatal, warnmesg)
1037 if (fntop == 0.0)
then ; y = ytop ; fny = fntop
1038 elseif (fnbot == 0.0)
then ; y = ybot ; fny = fnbot
1040 y = (ybot*fntop - ytop*fnbot) / (fntop - fnbot)
1041 fny = fn(y,gp) - fnval
1042 if (fny < 0.0)
then ; fnbot = fny ; ybot = y
1043 else ; fntop = fny ; ytop = y ;
endif
1047 dy_dfn = dy_df(y,gp)
1049 dy = -1.0* fny * dy_dfn
1051 if ((y_next >= ytop) .or. (y_next <= ybot))
then
1056 if (abs(fntop - fnbot) > epsilon(y) * (abs(fntop) + abs(fnbot))) &
1057 y_next = (ybot*fntop - ytop*fnbot) / (fntop - fnbot)
1061 if (abs(dy) < (2.0*epsilon(y)*(abs(y) + abs(y_next)) + 1.0e-20))
then
1066 fny = fn(y,gp) - fnval
1067 if (fny > 0.0)
then ; ytop = y ; fntop = fny
1068 elseif (fny < 0.0)
then ; ybot = y ; fnbot = fny
1072 if (abs(y) < 1e-12) y = 0.0
1076 end function find_root
1080 function dx_di(x, GP)
1081 real,
intent(in) :: x
1082 type(
gps),
intent(in) :: gp
1085 dx_di = (gp%len_lon * 4.0*atan(1.0)) / (180.0 * gp%niglobal)
1091 function int_di_dx(x, GP)
1092 real,
intent(in) :: x
1093 type(
gps),
intent(in) :: gp
1096 int_di_dx = x * ((180.0 * gp%niglobal) / (gp%len_lon * 4.0*atan(1.0)))
1098 end function int_di_dx
1102 function dy_dj(y, GP)
1103 real,
intent(in) :: y
1104 type(
gps),
intent(in) :: gp
1110 real :: y_eq_enhance
1113 if (gp%isotropic)
then
1114 c0 = (gp%len_lon * pi) / (180.0 * gp%niglobal)
1115 y_eq_enhance = pi*abs(gp%lat_eq_enhance)/180.0
1116 if (abs(y) < y_eq_enhance)
then
1117 dy_dj = c0 * (cos(y) / (1.0 + 0.5*cos(y) * (gp%lat_enhance_factor - 1.0) * &
1118 (1.0+cos(pi*y/y_eq_enhance)) ))
1123 c0 = (gp%len_lat * pi) / (180.0 * gp%njglobal)
1131 function int_dj_dy(y, GP)
1132 real,
intent(in) :: y
1133 type(
gps),
intent(in) :: gp
1140 real :: y_eq_enhance
1147 if (gp%isotropic)
then
1148 i_c0 = (180.0 * gp%niglobal) / (gp%len_lon * pi)
1149 y_eq_enhance = pi*abs(gp%lat_eq_enhance)/180.0
1152 r = i_c0 * log((1.0 + sin(y))/cos(y))
1154 r = -1.0 * i_c0 * log((1.0 - sin(y))/cos(y))
1157 if (y >= y_eq_enhance)
then
1158 r = r + i_c0*0.5*(gp%lat_enhance_factor - 1.0)*y_eq_enhance
1159 elseif (y <= -y_eq_enhance)
then
1160 r = r - i_c0*0.5*(gp%lat_enhance_factor - 1.0)*y_eq_enhance
1162 r = r + i_c0*0.5*(gp%lat_enhance_factor - 1.0) * &
1163 (y + (y_eq_enhance/pi)*sin(pi*y/y_eq_enhance))
1166 i_c0 = (180.0 * gp%njglobal) / (gp%len_lat * pi)
1171 end function int_dj_dy
1174 subroutine extrapolate_metric(var, jh, missing)
1175 real,
dimension(:,:),
intent(inout) :: var
1176 integer,
intent(in) :: jh
1177 real,
optional,
intent(in) :: missing
1182 badval = 0.0 ;
if (
present(missing)) badval = missing
1185 do j=lbound(var,2)+jh,lbound(var,2),-1 ;
do i=lbound(var,1),ubound(var,1)
1186 if (var(i,j)==badval) var(i,j) = 2.0*var(i,j+1)-var(i,j+2)
1190 do j=ubound(var,2)-jh,ubound(var,2) ;
do i=lbound(var,1),ubound(var,1)
1191 if (var(i,j)==badval) var(i,j) = 2.0*var(i,j-1)-var(i,j-2)
1195 do j=lbound(var,2),ubound(var,2) ;
do i=lbound(var,1)+jh,lbound(var,1),-1
1196 if (var(i,j)==badval) var(i,j) = 2.0*var(i+1,j)-var(i+2,j)
1200 do j=lbound(var,2),ubound(var,2) ;
do i=ubound(var,1)-jh,ubound(var,1)
1201 if (var(i,j)==badval) var(i,j) = 2.0*var(i-1,j)-var(i-2,j)
1204 end subroutine extrapolate_metric
1208 function adcroft_reciprocal(val)
result(I_val)
1209 real,
intent(in) :: val
1213 if (val /= 0.0) i_val = 1.0/val
1214 end function adcroft_reciprocal
1224 subroutine initialize_masks(G, PF, US)
1229 real :: m_to_z_scale
1233 character(len=40) :: mdl =
"MOM_grid_init initialize_masks"
1236 call calltree_enter(
"initialize_masks(), MOM_grid_initialize.F90")
1237 m_to_z_scale = 1.0 ;
if (
present(us)) m_to_z_scale = us%m_to_Z
1238 call get_param(pf, mdl,
"MINIMUM_DEPTH", min_depth, &
1239 "If MASKING_DEPTH is unspecified, then anything shallower than "//&
1240 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
1241 "If MASKING_DEPTH is specified, then all depths shallower than "//&
1242 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
1243 units=
"m", default=0.0, scale=m_to_z_scale)
1244 call get_param(pf, mdl,
"MASKING_DEPTH", mask_depth, &
1245 "The depth below which to mask points as land points, for which all "//&
1246 "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", &
1247 units=
"m", default=-9999.0, scale=m_to_z_scale)
1250 if (mask_depth>=0.) dmin = mask_depth
1252 g%mask2dCu(:,:) = 0.0 ; g%mask2dCv(:,:) = 0.0 ; g%mask2dBu(:,:) = 0.0
1255 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
1256 if (g%bathyT(i,j) <= dmin)
then
1257 g%mask2dT(i,j) = 0.0
1259 g%mask2dT(i,j) = 1.0
1263 do j=g%jsd,g%jed ;
do i=g%isd,g%ied-1
1264 if ((g%bathyT(i,j) <= dmin) .or. (g%bathyT(i+1,j) <= dmin))
then
1265 g%mask2dCu(i,j) = 0.0
1267 g%mask2dCu(i,j) = 1.0
1271 do j=g%jsd,g%jed-1 ;
do i=g%isd,g%ied
1272 if ((g%bathyT(i,j) <= dmin) .or. (g%bathyT(i,j+1) <= dmin))
then
1273 g%mask2dCv(i,j) = 0.0
1275 g%mask2dCv(i,j) = 1.0
1279 do j=g%jsd,g%jed-1 ;
do i=g%isd,g%ied-1
1280 if ((g%bathyT(i+1,j) <= dmin) .or. (g%bathyT(i+1,j+1) <= dmin) .or. &
1281 (g%bathyT(i,j) <= dmin) .or. (g%bathyT(i,j+1) <= dmin))
then
1282 g%mask2dBu(i,j) = 0.0
1284 g%mask2dBu(i,j) = 1.0
1288 call pass_var(g%mask2dBu, g%Domain, position=corner)
1289 call pass_vector(g%mask2dCu, g%mask2dCv, g%Domain, to_all+scalar_pair, cgrid_ne)
1291 do j=g%jsd,g%jed ;
do i=g%IsdB,g%IedB
1292 g%dy_Cu(i,j) = g%mask2dCu(i,j) * g%dyCu(i,j)
1293 g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
1294 g%IareaCu(i,j) = g%mask2dCu(i,j) * adcroft_reciprocal(g%areaCu(i,j))
1297 do j=g%JsdB,g%JedB ;
do i=g%isd,g%ied
1298 g%dx_Cv(i,j) = g%mask2dCv(i,j) * g%dxCv(i,j)
1299 g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
1300 g%IareaCv(i,j) = g%mask2dCv(i,j) * adcroft_reciprocal(g%areaCv(i,j))
1303 call calltree_leave(
"initialize_masks()")
1304 end subroutine initialize_masks