Ticket #23: ana_grid.h

File ana_grid.h, 33.7 KB (added by m.hadfield, 19 years ago)
Line 
1 SUBROUTINE ana_grid (ng, tile, model)
2!
3!! svn $Id: ana_grid.h 34 2007-04-27 04:40:21Z arango $
4!!======================================================================
5!! Copyright (c) 2002-2007 The ROMS/TOMS Group !
6!! Licensed under a MIT/X style license !
7!! See License_ROMS.txt !
8!! !
9!=======================================================================
10! !
11! This routine sets model grid using an analytical expressions. !
12! !
13! On Output: stored in common blocks: !
14! !
15! "grid" (file grid.h) !
16! "scalars" (file scalar.h) !
17! !
18! el Length (m) of domain box in the ETA-direction. !
19! f Coriolis parameter (1/seconds) at RHO-points. !
20! h Bathymetry (meters; positive) at RHO-points. !
21! hmin Minimum depth of bathymetry (m). !
22! hmax Maximum depth of bathymetry (m). !
23! pm Coordinate transformation metric "m" (1/meters) !
24! associated with the differential distances in XI !
25! at RHO-points. !
26! pn Coordinate transformation metric "n" (1/meters) !
27! associated with the differential distances in ETA. !
28! at RHO-points. !
29! xl Length (m) of domain box in the XI-direction. !
30! xp XI-coordinates (m) at PSI-points. !
31! xr XI-coordinates (m) at RHO-points. !
32! yp ETA-coordinates (m) at PSI-points. !
33! yr ETA-coordinates (m) at RHO-points. !
34! !
35!=======================================================================
36!
37 USE mod_param
38 USE mod_grid
39 USE mod_ncparam
40!
41! Imported variable declarations.
42!
43 integer, intent(in) :: ng, tile, model
44
45#include "tile.h"
46!
47 CALL ana_grid_tile (ng, model, Istr, Iend, Jstr, Jend, &
48 & LBi, UBi, LBj, UBj, &
49 & GRID(ng) % angler, &
50#if defined CURVGRID && defined UV_ADV
51 & GRID(ng) % dmde, &
52 & GRID(ng) % dndx, &
53#endif
54#ifdef ICESHELF
55 & GRID(ng) % zice, &
56#endif
57#ifdef SPHERICAL
58 & GRID(ng) % lonp, &
59 & GRID(ng) % lonr, &
60 & GRID(ng) % lonu, &
61 & GRID(ng) % lonv, &
62 & GRID(ng) % latp, &
63 & GRID(ng) % latr, &
64 & GRID(ng) % latu, &
65 & GRID(ng) % latv, &
66#else
67 & GRID(ng) % xp, &
68 & GRID(ng) % xr, &
69 & GRID(ng) % xu, &
70 & GRID(ng) % xv, &
71 & GRID(ng) % yp, &
72 & GRID(ng) % yr, &
73 & GRID(ng) % yu, &
74 & GRID(ng) % yv, &
75#endif
76 & GRID(ng) % pn, &
77 & GRID(ng) % pm, &
78 & GRID(ng) % f, &
79 & GRID(ng) % h)
80!
81! Set analytical header file name used.
82!
83 IF (Lanafile) THEN
84 ANANAME( 7)='ROMS/Functionals/ana_grid.h'
85 END IF
86
87 RETURN
88 END SUBROUTINE ana_grid
89!
90!***********************************************************************
91 SUBROUTINE ana_grid_tile (ng, model, Istr, Iend, Jstr, Jend, &
92 & LBi, UBi, LBj, UBj, &
93 & angler, &
94#if defined CURVGRID && defined UV_ADV
95 & dmde, dndx, &
96#endif
97#ifdef ICESHELF
98 & zice, &
99#endif
100#ifdef SPHERICAL
101 & lonp, lonr, lonu, lonv, &
102 & latp, latr, latu, latv, &
103#else
104 & xp, xr, xu, xv, &
105 & yp, yr, yu, yv, &
106#endif
107 & pn, pm, f, h)
108!***********************************************************************
109!
110 USE mod_param
111 USE mod_parallel
112 USE mod_scalars
113!
114#ifdef DISTRIBUTE
115 USE distribute_mod, ONLY : mp_reduce
116#endif
117#if defined EW_PERIODIC || defined NS_PERIODIC
118 USE exchange_2d_mod, ONLY : exchange_r2d_tile
119#endif
120#ifdef DISTRIBUTE
121 USE mp_exchange_mod, ONLY : mp_exchange2d
122#endif
123!
124! Imported variable declarations.
125!
126 integer, intent(in) :: ng, model, Iend, Istr, Jend, Jstr
127 integer, intent(in) :: LBi, UBi, LBj, UBj
128!
129#ifdef ASSUMED_SHAPE
130 real(r8), intent(out) :: angler(LBi:,LBj:)
131# if defined CURVGRID && defined UV_ADV
132 real(r8), intent(out) :: dmde(LBi:,LBj:)
133 real(r8), intent(out) :: dndx(LBi:,LBj:)
134# endif
135# ifdef ICESHELF
136 real(r8), intent(out) :: zice(LBi:,LBj:)
137# endif
138# ifdef SPHERICAL
139 real(r8), intent(out) :: lonp(LBi:,LBj:)
140 real(r8), intent(out) :: lonr(LBi:,LBj:)
141 real(r8), intent(out) :: lonu(LBi:,LBj:)
142 real(r8), intent(out) :: lonv(LBi:,LBj:)
143 real(r8), intent(out) :: latp(LBi:,LBj:)
144 real(r8), intent(out) :: latr(LBi:,LBj:)
145 real(r8), intent(out) :: latu(LBi:,LBj:)
146 real(r8), intent(out) :: latv(LBi:,LBj:)
147# else
148 real(r8), intent(out) :: xp(LBi:,LBj:)
149 real(r8), intent(out) :: xr(LBi:,LBj:)
150 real(r8), intent(out) :: xu(LBi:,LBj:)
151 real(r8), intent(out) :: xv(LBi:,LBj:)
152 real(r8), intent(out) :: yp(LBi:,LBj:)
153 real(r8), intent(out) :: yr(LBi:,LBj:)
154 real(r8), intent(out) :: yu(LBi:,LBj:)
155 real(r8), intent(out) :: yv(LBi:,LBj:)
156# endif
157 real(r8), intent(out) :: pn(LBi:,LBj:)
158 real(r8), intent(out) :: pm(LBi:,LBj:)
159 real(r8), intent(out) :: f(LBi:,LBj:)
160 real(r8), intent(out) :: h(LBi:,LBj:)
161#else
162 real(r8), intent(out) :: angler(LBi:UBi,LBj:UBj)
163# if defined CURVGRID && defined UV_ADV
164 real(r8), intent(out) :: dmde(LBi:UBi,LBj:UBj)
165 real(r8), intent(out) :: dndx(LBi:UBi,LBj:UBj)
166# endif
167# ifdef ICESHELF
168 real(r8), intent(out) :: zice(LBi:UBi,LBj:UBj)
169# endif
170# ifdef SPHERICAL
171 real(r8), intent(out) :: lonp(LBi:UBi,LBj:UBj)
172 real(r8), intent(out) :: lonr(LBi:UBi,LBj:UBj)
173 real(r8), intent(out) :: lonu(LBi:UBi,LBj:UBj)
174 real(r8), intent(out) :: lonv(LBi:UBi,LBj:UBj)
175 real(r8), intent(out) :: latp(LBi:UBi,LBj:UBj)
176 real(r8), intent(out) :: latr(LBi:UBi,LBj:UBj)
177 real(r8), intent(out) :: latu(LBi:UBi,LBj:UBj)
178 real(r8), intent(out) :: latv(LBi:UBi,LBj:UBj)
179# else
180 real(r8), intent(out) :: xp(LBi:UBi,LBj:UBj)
181 real(r8), intent(out) :: xr(LBi:UBi,LBj:UBj)
182 real(r8), intent(out) :: xu(LBi:UBi,LBj:UBj)
183 real(r8), intent(out) :: xv(LBi:UBi,LBj:UBj)
184 real(r8), intent(out) :: yp(LBi:UBi,LBj:UBj)
185 real(r8), intent(out) :: yr(LBi:UBi,LBj:UBj)
186 real(r8), intent(out) :: yu(LBi:UBi,LBj:UBj)
187 real(r8), intent(out) :: yv(LBi:UBi,LBj:UBj)
188# endif
189 real(r8), intent(out) :: pn(LBi:UBi,LBj:UBj)
190 real(r8), intent(out) :: pm(LBi:UBi,LBj:UBj)
191 real(r8), intent(out) :: f(LBi:UBi,LBj:UBj)
192 real(r8), intent(out) :: h(LBi:UBi,LBj:UBj)
193#endif
194!
195! Local variable declarations.
196!
197#ifdef DISTRIBUTE
198# ifdef EW_PERIODIC
199 logical :: EWperiodic=.TRUE.
200# else
201 logical :: EWperiodic=.FALSE.
202# endif
203# ifdef NS_PERIODIC
204 logical :: NSperiodic=.TRUE.
205# else
206 logical :: NSperiodic=.FALSE.
207# endif
208#endif
209 integer :: Imin, Imax, Jmin, Jmax
210 integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
211 integer :: NSUB, i, j, k
212
213 real(r8), parameter :: twopi = 2.0_r8*pi
214
215 real(r8) :: Esize, Xsize, beta, cff, depth, dth
216 real(r8) :: dx, dy, f0, my_min, my_max, r, theta, val1, val2
217
218#ifdef DISTRIBUTE
219 real(r8), dimension(2) :: buffer
220 character (len=3), dimension(2) :: op_handle
221#endif
222#ifdef WEDDELL
223 real(r8) :: hwrk(-1:235), xwrk(-1:235), zwrk
224#endif
225 real(r8) :: wrkX(PRIVATE_2D_SCRATCH_ARRAY)
226 real(r8) :: wrkY(PRIVATE_2D_SCRATCH_ARRAY)
227!
228#include "set_bounds.h"
229!
230!-----------------------------------------------------------------------
231! Set grid parameters:
232!
233! Xsize Length (m) of domain box in the XI-direction.
234! Esize Length (m) of domain box in the ETA-direction.
235! depth Maximum depth of bathymetry (m).
236! f0 Coriolis parameter, f-plane constant (1/s).
237! beta Coriolis parameter, beta-plane constant (1/s/m).
238!-----------------------------------------------------------------------
239!
240#if defined BASIN
241 Xsize=3600.0E+03_r8
242 Esize=2800.0E+03_r8
243 depth=5000.0_r8
244 f0=1.0E-04_r8
245 beta=2.0E-11_r8
246#elif defined BENCHMARK
247 Xsize=360.0_r8 ! degrees of longitude
248 Esize=20.0_r8 ! degrees of latitude
249 depth=4000.0_r8
250 f0=-1.0E-04_r8
251 beta=2.0E-11_r8
252#elif defined BL_TEST
253 Xsize=100.0E+03_r8
254 Esize=5.0E+03_r8
255 depth=47.5_r8
256 f0=9.25E-04_r8
257 beta=0.0_r8
258#elif defined CANYON
259 Xsize=128.0E+03_r8
260 Esize=96.0E+03_r8
261 depth=4000.0_r8
262 f0=1.0E-04_r8
263 beta=0.0_r8
264#elif defined COUPLING_TEST
265 Xsize=6000.0_r8*REAL(Lm(ng),r8)
266 Esize=6000.0_r8*REAL(Mm(ng),r8)
267 depth=1500.0_r8
268 f0=5.0E-05_r8
269 beta=0.0_r8
270#elif defined DOUBLE_GYRE
271 Xsize=1000.0E+03_r8
272 Esize=2000.0E+03_r8
273 depth=500.0_r8
274!! depth=5000.0_r8
275 f0=7.3E-05_r8
276 beta=2.0E-11_r8
277#elif defined ESTUARY_TEST
278 Xsize=100000.0_r8
279 Esize=300.0_r8
280 depth=10.0_r8
281 f0=0.0_r8
282 beta=0.0_r8
283#elif defined KELVIN
284 Xsize=20000.0_r8*REAL(Lm(ng),r8)
285 Esize=20000.0_r8*REAL(Mm(ng),r8)
286 depth=100.0_r8
287 f0=1.0E-04_r8
288 beta=0.0_r8
289#elif defined FLT_TEST
290 Xsize=1.0E+03_r8*REAL(Lm(ng),r8)
291 Esize=1.0E+03_r8*REAL(Mm(ng),r8)
292 depth=10.0_r8
293 f0=0.0_r8
294 beta=0.0_r8
295#elif defined GRAV_ADJ
296 Xsize=64.0E+03_r8
297 Esize=2.0E+03_r8
298 depth=20.0_r8
299 f0=0.0_r8
300 beta=0.0_r8
301#elif defined LAB_CANYON
302 Xsize=0.55_r8 ! width of annulus
303 Esize=2.0_r8*pi ! azimuthal length (radians)
304 f0=4.0_r8*pi/25.0_r8
305 beta=0.0_r8
306#elif defined LAKE_SIGNELL
307 Xsize=50.0e3_r8
308 Esize=10.0e3_r8
309 depth=18.0_r8
310 f0=0.0E-04_r8
311 beta=0.0_r8
312#elif defined LMD_TEST
313 Xsize=100.0E+03_r8
314 Esize=100.0E+03_r8
315 depth=50.0_r8
316 f0=1.09E-04_r8
317 beta=0.0_r8
318# elif defined MIXED_LAYER
319 Xsize=500.0_r8
320 Esize=400.0_r8
321 depth=50.0_r8
322 f0=0.0_r8
323 beta=0.0_r8
324#elif defined OVERFLOW
325 Xsize=4.0E+03_r8
326 Esize=200.0E+03_r8
327 depth=4000.0_r8
328 f0=0.0_r8
329 beta=0.0_r8
330#elif defined RIVERPLUME1
331 Xsize=58.5E+03_r8
332 Esize=201.0E+03_r8
333 depth=150.0_r8
334 f0=1.0E-04_r8
335 beta=0.0_r8
336#elif defined RIVERPLUME2
337 Xsize=100.0E+03_r8
338 Esize=210.0E+03_r8
339 depth=190.0_r8
340 f0=1.0E-04_r8
341 beta=0.0_r8
342#elif defined SEAMOUNT
343 Xsize=320.0E+03_r8
344 Esize=320.0E+03_r8
345 depth=5000.0_r8
346 f0=1.0E-04_r8
347 beta=0.0_r8
348#elif defined SOLITON
349!! Xsize=0.5_r8*REAL(Lm(ng),r8)
350!! Esize=0.5_r8*REAL(Mm(ng),r8)
351 Xsize=48.0_r8
352 Esize=16.0_r8
353 depth=1.0_r8
354 f0=0.0_r8
355 beta=1.0_r8
356 g=1.0_r8
357#elif defined SED_TEST1
358 Xsize=300.0_r8
359 Esize=36.0_r8
360 depth=10.0_r8
361 f0=0.0_r8
362 beta=0.0_r8
363#elif defined SED_TOY
364 Xsize=40.0_r8
365 Esize=30.0_r8
366 depth=0.5_r8
367 f0=0.0_r8
368 beta=0.0_r8
369# elif defined SHOREFACE
370 Xsize=1180.0_r8
371 Esize=140.0_r8
372 depth=15.0_r8
373 f0=0.0E-04_r8
374 beta=0.0_r8
375#elif defined TEST_CHAN
376 Xsize=10000.0_r8
377 Esize=1000.0_r8
378 depth=10.0_r8
379 f0=0.0_r8
380 beta=0.0_r8
381#elif defined UPWELLING
382 Xsize=1000.0_r8*REAL(Lm(ng),r8)
383 Esize=1000.0_r8*REAL(Mm(ng),r8)
384 depth=150.0_r8
385 f0=-8.26E-05_r8
386 beta=0.0_r8
387#elif defined WEDDELL
388 Xsize=4000.0_r8*REAL(Lm(ng),r8)
389 Esize=4000.0_r8*REAL(Mm(ng),r8)
390 depth=4500.0_r8
391 f0=0.0_r8
392 beta=0.0_r8
393#elif defined WINDBASIN
394 Xsize=2000.0_r8*REAL(Lm(ng),r8)
395 Esize=1000.0_r8*REAL(Mm(ng),r8)
396 depth=50.0_r8
397 f0=1.0E-04_r8
398 beta=0.0_r8
399#else
400 ana_grid.h: no values provided for Xsize, Esize, depth, f0, beta.
401#endif
402!
403! Load grid parameters to global storage.
404!
405 IF (SOUTH_WEST_TEST) THEN
406 xl(ng)=Xsize
407 el(ng)=Esize
408 END IF
409!
410!-----------------------------------------------------------------------
411! Compute the (XI,ETA) coordinates at PSI- and RHO-points.
412! Set grid spacing (m).
413!-----------------------------------------------------------------------
414!
415! Determine I- and J-ranges for computing grid data. This ranges
416! are special in periodic boundary conditons since periodicity cannot
417! be imposed in the grid coordinates.
418!
419 IF (WESTERN_EDGE) THEN
420 Imin=Istr-1
421 ELSE
422 Imin=Istr
423 END IF
424 IF (EASTERN_EDGE) THEN
425 Imax=Iend+1
426 ELSE
427 Imax=Iend
428 END IF
429 IF (SOUTHERN_EDGE) THEN
430 Jmin=Jstr-1
431 ELSE
432 Jmin=Jstr
433 END IF
434 IF (NORTHERN_EDGE) THEN
435 Jmax=Jend+1
436 ELSE
437 Jmax=Jend
438 END IF
439
440#if defined BENCHMARK
441!
442! Spherical coordinates set-up.
443!
444 dx=Xsize/REAL(Lm(ng),r8)
445 dy=Esize/REAL(Mm(ng),r8)
446 spherical=.TRUE.
447 DO j=Jmin,Jmax
448 val1=-70.0_r8+dy*(REAL(j,r8)-0.5_r8)
449 val2=-70.0_r8+dy*REAL(j,r8)
450 DO i=Imin,Imax
451 lonr(i,j)=dx*(REAL(i,r8)-0.5_r8)
452 latr(i,j)=val1
453 lonu(i,j)=dx*REAL(i,r8)
454 lonp(i,j)=lonu(i,j)
455 latu(i,j)=latr(i,j)
456 lonv(i,j)=lonr(i,j)
457 latv(i,j)=val2
458 latp(i,j)=latv(i,j)
459 END DO
460 END DO
461#elif defined LAB_CANYON
462!
463! Polar coordinates set-up.
464!
465 dx=Xsize/REAL(Lm(ng),r8)
466 dy=Esize/REAL(Mm(ng),r8)
467!! dth=twopi/REAL(Mm(ng),r8) ! equal azimultal spacing
468 dth=0.01_r8 ! azimultal spacing
469 cff=(4.0_r8*pi/(dth*REAL(Mm(ng),r8)))-1.0_r8 ! F
470 DO j=Jmin,Jmax
471 DO i=Imin,Imax
472 r=0.35_r8+dx*REAL(i-1,r8)
473 theta=-pi+ &
474 & 0.5_r8*dth*((cff+1.0_r8)*REAL(j-1,r8)+ &
475 & (cff-1.0_r8)*(REAL(Mm(ng),r8)/twopi)* &
476 & SIN(twopi*REAL(j-1,r8)/REAL(Mm(ng),r8)))
477 xp(i,j)=r*COS(theta)
478 yp(i,j)=r*SIN(theta)
479 r=0.35_r8+dx*(REAL(i-1,r8)+0.5_r8)
480 theta=-pi+ &
481 & 0.5_r8*dth*((cff+1.0_r8)*(REAL(j-1,r8)+0.5_r8)+ &
482 & (cff-1.0_r8)*(REAL(Mm(ng),r8)/twopi)* &
483 & SIN(twopi*(REAL(j-1,r8)+0.5_r8)/ &
484 & REAL(Mm(ng),r8)))
485 xr(i,j)=r*COS(theta)
486 yr(i,j)=r*SIN(theta)
487 xu(i,j)=xp(i,j)
488 yu(i,j)=yr(i,j)
489 xv(i,j)=xr(i,j)
490 yv(i,j)=yp(i,j)
491 END DO
492 END DO
493#else
494 dx=Xsize/REAL(Lm(ng),r8)
495 dy=Esize/REAL(Mm(ng),r8)
496 DO j=Jmin,Jmax
497 DO i=Imin,Imax
498# ifdef BL_TEST
499 dx=0.5_r8*(4000.0_r8/REAL(Lm(ng)+1,r8))*REAL(i,r8)+675.0_r8
500# endif
501 xp(i,j)=dx*REAL(i-1,r8)
502 xr(i,j)=dx*(REAL(i-1,r8)+0.5_r8)
503 xu(i,j)=xp(i,j)
504 xv(i,j)=xr(i,j)
505 yp(i,j)=dy*REAL(j-1,r8)
506 yr(i,j)=dy*(REAL(j-1,r8)+0.5_r8)
507 yu(i,j)=yr(i,j)
508 yv(i,j)=yp(i,j)
509 END DO
510 END DO
511#endif
512#ifdef DISTRIBUTE
513# ifdef SPHERICAL
514 CALL mp_exchange2d (ng, model, 4, Istr, Iend, Jstr, Jend, &
515 & LBi, UBi, LBj, UBj, &
516 & NghostPoints, .FALSE., .FALSE., &
517 & lonp, lonr, lonu, lonv)
518 CALL mp_exchange2d (ng, model, 4, Istr, Iend, Jstr, Jend, &
519 & LBi, UBi, LBj, UBj, &
520 & NghostPoints, .FALSE., .FALSE., &
521 & latp, latr, latu, latv)
522# else
523 CALL mp_exchange2d (ng, model, 4, Istr, Iend, Jstr, Jend, &
524 & LBi, UBi, LBj, UBj, &
525 & NghostPoints, .FALSE., .FALSE., &
526 & xp, xr, xu, xv)
527 CALL mp_exchange2d (ng, model, 4, Istr, Iend, Jstr, Jend, &
528 & LBi, UBi, LBj, UBj, &
529 & NghostPoints, .FALSE., .FALSE., &
530 & yp, yr, yu, yv)
531# endif
532#endif
533!
534!-----------------------------------------------------------------------
535! Compute coordinate transformation metrics at RHO-points "pm" and
536! "pn" (1/m) associated with the differential distances in XI and
537! ETA, respectively.
538!-----------------------------------------------------------------------
539!
540#define J_RANGE MIN(JstrR,Jstr-1),MAX(Jend+1,JendR)
541#define I_RANGE MIN(IstrR,Istr-1),MAX(Iend+1,IendR)
542
543#if defined BENCHMARK
544!
545! Spherical coordinates set-up.
546!
547 val1=REAL(Lm(ng),r8)/(2.0_r8*pi*Eradius)
548 val2=REAL(Mm(ng),r8)*360.0_r8/(2.0_r8*pi*Eradius*Esize)
549 DO j=J_RANGE
550 cff=1.0_r8/COS((-70.0_r8+dy*(REAL(j,r8)-0.5_r8))*deg2rad)
551 DO i=I_RANGE
552 wrkX(i,j)=val1*cff
553 wrkY(i,j)=val2
554 END DO
555 END DO
556#elif defined LAB_CANYON
557!
558! Polar coordinates set-up.
559!
560 DO j=J_RANGE
561 DO i=I_RANGE
562 r=0.35_r8+dx*(REAL(i-1,r8)+0.5_r8)
563 theta=0.5_r8*dth*((cff+1.0_r8)+ &
564 & (cff-1.0_r8)* &
565 & COS(twopi*REAL(j-1,r8)/REAL(Mm(ng),r8)))
566 wrkX(i,j)=1.0_r8/dx
567 wrkY(i,j)=1.0_r8/(r*theta)
568 END DO
569 END DO
570#else
571 DO j=J_RANGE
572 DO i=I_RANGE
573# ifdef BL_TEST
574 dx=0.5_r8*(4000.0_r8/REAL(Lm(ng)+1,r8))*REAL(i,r8)+675.0_r8
575# endif
576 wrkX(i,j)=1.0_r8/dx
577 wrkY(i,j)=1.0_r8/dy
578 END DO
579 END DO
580#endif
581#undef J_RANGE
582#undef I_RANGE
583 DO j=JstrR,JendR
584 DO i=IstrR,IendR
585 pm(i,j)=wrkX(i,j)
586 pn(i,j)=wrkY(i,j)
587 END DO
588 END DO
589#if defined EW_PERIODIC || defined NS_PERIODIC
590 CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend, &
591 & LBi, UBi, LBj, UBj, &
592 & pm)
593 CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend, &
594 & LBi, UBi, LBj, UBj, &
595 & pn)
596#endif
597#ifdef DISTRIBUTE
598 CALL mp_exchange2d (ng, model, 2, Istr, Iend, Jstr, Jend, &
599 & LBi, UBi, LBj, UBj, &
600 & NghostPoints, EWperiodic, NSperiodic, &
601 & pm, pn)
602#endif
603#if (defined CURVGRID && defined UV_ADV)
604!
605!-----------------------------------------------------------------------
606! Compute d(1/n)/d(xi) and d(1/m)/d(eta) at RHO-points.
607!-----------------------------------------------------------------------
608!
609 DO j=Jstr,Jend
610 DO i=Istr,Iend
611 dndx(i,j)=0.5_r8*((1.0_r8/wrkY(i+1,j ))- &
612 & (1.0_r8/wrkY(i-1,j )))
613 dmde(i,j)=0.5_r8*((1.0_r8/wrkX(i ,j+1))- &
614 & (1.0_r8/wrkX(i ,j-1)))
615 END DO
616 END DO
617# if defined EW_PERIODIC || defined NS_PERIODIC
618 CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend, &
619 & LBi, UBi, LBj, UBj, &
620 & dndx)
621 CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend, &
622 & LBi, UBi, LBj, UBj, &
623 & dmde)
624# endif
625# ifdef DISTRIBUTE
626 CALL mp_exchange2d (ng, model, 2, Istr, Iend, Jstr, Jend, &
627 & LBi, UBi, LBj, UBj, &
628 & NghostPoints, EWperiodic, NSperiodic, &
629 & dndx, dmde)
630# endif
631#endif
632!
633!-----------------------------------------------------------------------
634! Angle (radians) between XI-axis and true EAST at RHO-points.
635!-----------------------------------------------------------------------
636!
637#if defined LAB_CANYON
638 DO j=JstrR,JendR
639 DO i=IstrR,IendR
640 theta=-pi+ &
641 & 0.5_r8*dth*((cff+1.0_r8)*(REAL(j-1,r8)+0.5_r8)+ &
642 & (cff-1.0_r8)*(REAL(Mm(ng),r8)/twopi)* &
643 & SIN(twopi*(REAL(j-1,r8)+0.5_r8)/ &
644 & REAL(Mm(ng),r8)))
645 angler(i,j)=theta
646 END DO
647 END DO
648#elif defined WEDDELL
649 val1=90.0_r8*deg2rad
650 DO j=JstrR,JendR
651 DO i=IstrR,IendR
652 angler(i,j)=val1
653 END DO
654 END DO
655#else
656 DO j=JstrR,JendR
657 DO i=IstrR,IendR
658 angler(i,j)=0.0_r8
659 END DO
660 END DO
661#endif
662#if defined EW_PERIODIC || defined NS_PERIODIC
663 CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend, &
664 & LBi, UBi, LBj, UBj, &
665 & angler)
666#endif
667#ifdef DISTRIBUTE
668 CALL mp_exchange2d (ng, model, 1, Istr, Iend, Jstr, Jend, &
669 & LBi, UBi, LBj, UBj, &
670 & NghostPoints, EWperiodic, NSperiodic, &
671 & angler)
672#endif
673!
674!-----------------------------------------------------------------------
675! Compute Coriolis parameter (1/s) at RHO-points.
676!-----------------------------------------------------------------------
677!
678#if defined BENCHMARK
679 val1=2.0_r8*(2.0_r8*pi*366.25_r8/365.25_r8)/86400.0_r8
680 DO j=JstrR,JendR
681 DO i=IstrR,IendR
682 f(i,j)=val1*SIN(latr(i,j)*deg2rad)
683 END DO
684 END DO
685#elif defined WEDDELL
686 val1=10.4_r8/REAL(Lm(ng),r8)
687 DO j=JstrR,JendR
688 DO i=IstrR,IendR
689 f(i,j)=2.0_r8*7.2E-05_r8* &
690 & SIN((-79.0_r8+REAL(i-1,r8)*val1)*deg2rad)
691 END DO
692 END DO
693#else
694 val1=0.5_r8*Esize
695 DO j=JstrR,JendR
696 DO i=IstrR,IendR
697 f(i,j)=f0+beta*(yr(i,j)-val1)
698 END DO
699 END DO
700#endif
701#if defined EW_PERIODIC || defined NS_PERIODIC
702 CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend, &
703 & LBi, UBi, LBj, UBj, &
704 & f)
705#endif
706#ifdef DISTRIBUTE
707 CALL mp_exchange2d (ng, model, 1, Istr, Iend, Jstr, Jend, &
708 & LBi, UBi, LBj, UBj, &
709 & NghostPoints, EWperiodic, NSperiodic, &
710 & f)
711#endif
712!
713!-----------------------------------------------------------------------
714! Set bathymetry (meters; positive) at RHO-points.
715!-----------------------------------------------------------------------
716!
717#if defined BENCHMARK
718 DO j=JstrR,JendR
719 DO i=IstrR,IendR
720 h(i,j)=500.0_r8+1750.0_r8*(1.0+TANH((68.0_r8+latr(i,j))/dy))
721 END DO
722 END DO
723#elif defined BL_TEST
724 DO j=JstrR,JendR
725 DO i=IstrR,IendR
726 val1=(xr(i,j)+500.0_r8)/15000.0_r8
727 h(i,j)=14.0_r8+ &
728 & 25.0_r8*(1.0_r8-EXP(-pi*xr(i,j)*1.0E-05_r8))- &
729 & 8.0_r8*EXP(-val1*val1)
730 END DO
731 END DO
732#elif defined CANYON_A || defined CANYON_B
733 DO j=JstrR,JendR
734 DO i=IstrR,IendR
735 val1=32000.0_r8-16000.0_r8*(SIN(pi*xr(i,j)/Xsize))**24
736 h(i,j)=20.0_r8+0.5_r8*(depth-20.0_r8)* &
737 & (1.0_r8+TANH((yr(i,j)-val1)/10000.0_r8))
738 END DO
739 END DO
740#elif defined ESTUARY_TEST
741 DO j=JstrR,JendR
742 DO i=IstrR,IendR
743 h(i,j)=5.0_r8+(Xsize-xr(i,j))/Xsize*5.0_r8
744 END DO
745 END DO
746#elif defined LAB_CANYON
747 DO j=JstrR,JendR
748 DO i=IstrR,IendR
749 r=0.35_r8+dx*(REAL(i-1,r8)+0.5_r8)
750 theta=-pi+ &
751 & 0.5_r8*dth*((cff+1.0_r8)*(REAL(j-1,r8)+0.5_r8)+ &
752 & (cff-1.0_r8)*(REAL(Mm(ng),r8)/twopi)* &
753 & SIN(dth*(REAL(j-1,r8)+0.5_r8)/ &
754 & REAL(Mm(ng),r8)))
755 val1=0.55_r8-0.15_r8*(COS(pi*theta*0.55_r8/0.2_r8)**2) !r_small
756 val2=0.15_r8+0.15_r8*(COS(pi*theta*0.55_r8/0.2_r8)**2) !lambda
757 IF (ABS(theta).ge.0.181818181818_r8) THEN
758 IF (r.le.0.55_r8) THEN
759 h(i,j)=0.025_r8 ! shelf
760 ELSE IF (r.ge.0.7_r8) THEN
761 h(i,j)=0.125_r8 ! deep
762 ELSE
763 h(i,j)=0.125_r8-0.1_r8* &
764 & (COS(0.5_r8*pi*(r-0.55_r8)/0.15_r8)**2)
765 END IF
766 ELSE
767 IF (r.le.val1) THEN
768 h(i,j)=0.025_r8 ! shelf
769 ELSE IF (r.ge.0.7_r8) THEN
770 h(i,j)=0.125_r8 ! deep
771 ELSE
772 h(i,j)=0.125_r8-0.1_r8* &
773 & (COS(0.5_r8*pi*(r-val1)/val2)**2)
774 END IF
775 END IF
776 END DO
777 END DO
778#elif defined LAKE_SIGNELL
779 DO j=JstrR,JendR
780 DO i=IstrR,IendR
781 h(i,j)=18.0_r8-16.0_r8*FLOAT(Mm(ng)-j)/FLOAT(Mm(ng)-1)
782 END DO
783 END DO
784# elif defined MIXED_LAYER
785 DO j=JstrR,JendR
786 DO i=IstrR,IendR
787 h(i,j)=50.0_r8
788 END DO
789 END DO
790#elif defined OVERFLOW
791 val1=200.0_r8
792 DO j=JstrR,JendR
793 DO i=IstrR,IendR
794 h(i,j)=val1+0.5_r8*(depth-val1)* &
795 & (1.0_r8+TANH((yr(i,j)-100000.0_r8)/20000.0_r8))
796 END DO
797 END DO
798#elif defined RIVERPLUME1
799 DO j=JstrR,JendR
800 DO i=IstrR,MIN(5,IendR)
801 h(i,j)=15.0_r8
802 END DO
803 DO i=MAX(6,IstrR),IendR
804 h(i,j)=depth+REAL(Lm(ng)-i,r8)*(15.0_r8-depth)/ &
805 & REAL(Lm(ng)-6,r8)
806 END DO
807 END DO
808#elif defined RIVERPLUME2
809 DO j=JstrR,JendR
810 DO i=IstrR,MIN(5,IendR)
811 h(i,j)=15.0_r8
812 END DO
813 DO i=MAX(6,IstrR),IendR
814 h(i,j)=depth+REAL(Lm(ng)-i,r8)*(15.0_r8-depth)/ &
815 & REAL(Lm(ng)-6,r8)
816 END DO
817 END DO
818#elif defined SEAMOUNT
819 DO j=JstrR,JendR
820 DO i=IstrR,IendR
821 val1=(xr(i,j)-xr(Lm(ng)/2,Mm(ng)/2))/40000.0_r8
822 val2=(yr(i,j)-yr(Lm(ng)/2,Mm(ng)/2))/40000.0_r8
823 h(i,j)=depth-4500.0_r8*EXP(-(val1*val1+val2*val2))
824 END DO
825 END DO
826#elif defined SED_TOY
827 DO j=JstrR,JendR
828 DO i=IstrR,IendR
829 h(i,j)=20.0_r8
830 END DO
831 END DO
832#elif defined SHOREFACE
833 DO j=JstrR,JendR
834 DO i=IstrR,IendR
835 h(i,j)=11.75_r8-0.0125_r8*Xsize/FLOAT(Lm(ng)+1)*FLOAT(i)
836 END DO
837 END DO
838#elif defined TEST_CHAN
839 DO j=JstrR,JendR
840 DO i=IstrR,IendR
841 h(i,j)=10.0_r8+0.4040_r8*REAL(i,r8)/REAL(Lm(ng)+1,r8)
842 END DO
843 END DO
844#elif defined UPWELLING
845 DO j=JstrR,JendR
846 IF (j.le.Mm(ng)/2) THEN
847 val1=REAL(j,r8)
848 ELSE
849 val1=REAL(Mm(ng)+1-j,r8)
850 END IF
851 val2=MIN(depth,84.5_r8+66.526_r8*TANH((val1-10.0_r8)/7.0_r8))
852 DO i=IstrR,IendR
853 h(i,j)=val2
854 END DO
855 END DO
856#elif defined WEDDELL
857 val1=98.80_r8
858 val2=0.8270_r8
859 DO k=-1,26
860 xwrk(k)=REAL(k-1,r8)*15.0_r8*1000.0_r8
861 hwrk(k)=375.0_r8
862 END DO
863 DO k=27,232
864 zwrk=-2.0_r8+REAL(k-1,r8)*0.020_r8
865 xwrk(k)=(520.0_r8+val1+zwrk*val1+ &
866 & val1*val2*LOG(COSH(zwrk)))*1000.0_r8
867 hwrk(k)=-75.0_r8+2198.0_r8*(1.0_r8+val2*TANH(zwrk))
868 END DO
869 DO k=233,235
870 xwrk(k)=(850.0_r8+REAL(k-228,r8)*50.0_r8)*1000.0_r8
871 hwrk(k)=4000.0_r8
872 END DO
873 DO j=JstrR,JendR
874 DO i=IstrR,IendR
875 DO k=1,234
876 IF ((xwrk(k).le.xr(i,1)).and.(xr(i,1).lt.xwrk(k+1))) THEN
877 cff=1.0_r8/(xwrk(k+1)-xwrk(k))
878 h(i,j)=cff*(xwrk(k+1)-xr(i,j))*hwrk(k )+ &
879 & cff*(xr(i,j)-xwrk(k ))*hwrk(k+1)
880 END IF
881 END DO
882 END DO
883 END DO
884#elif defined WINDBASIN
885 DO i=IstrR,IendR
886 val1=1;
887 IF ((i-IstrR).lt.(INT(0.03_r8*REAL(IendR-IstrR,r8)))) THEN
888 val1=1.0_r8-(REAL((i-IstrR+1)- &
889 & INT(0.03_r8*REAL(IendR-IstrR,r8)),r8)/ &
890 & (0.03_r8*REAL(IendR-IstrR,r8)))**2
891 END IF
892 IF ((IendR-i).lt.(INT(0.03_r8*REAL(IendR-IstrR,r8)))) THEN
893 val1=1.0_r8-(REAL((IendR-i+1)- &
894 & INT(0.03_r8*REAL(IendR-IstrR,r8)),r8)/ &
895 & (0.03_r8*REAL(IendR-IstrR,r8)))**2
896 END IF
897 DO j=JstrR,JendR
898 val2=2.0_r8*REAL(j-(Mm(ng)+1)/2,r8)/REAL(Mm(ng)+1,r8)
899 h(i,j)=depth*(0.08_r8+0.92_r8*val1*(1.0_r8-val2*val2))
900 END DO
901 END DO
902#else
903 DO j=JstrR,JendR
904 DO i=IstrR,IendR
905 h(i,j)=depth
906 END DO
907 END DO
908#endif
909#if defined EW_PERIODIC || defined NS_PERIODIC
910 CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend, &
911 & LBi, UBi, LBj, UBj, &
912 & h)
913#endif
914#ifdef DISTRIBUTE
915 CALL mp_exchange2d (ng, model, 1, Istr, Iend, Jstr, Jend, &
916 & LBi, UBi, LBj, UBj, &
917 & NghostPoints, EWperiodic, NSperiodic, &
918 & h)
919#endif
920!
921! Determine minimum depth: first, determine minimum values of depth
922! within each subdomain (stored as private variable cff), then
923! determine global minimum by comparing these subdomain minima.
924!
925 my_min=h(IstrR,JstrR)
926 my_max=h(IstrR,JstrR)
927 DO j=JstrR,JendR
928 DO i=IstrR,IendR
929 my_min=MIN(my_min,h(i,j))
930 my_max=MAX(my_max,h(i,j))
931 END DO
932 END DO
933 IF (SOUTH_WEST_CORNER.and. &
934 & NORTH_EAST_CORNER) THEN
935 NSUB=1 ! non-tiled application
936 ELSE
937 NSUB=NtileX(ng)*NtileE(ng) ! tiled application
938 END IF
939!$OMP CRITICAL (H_RANGE)
940 IF (tile_count.eq.0) THEN
941 hmin(ng)=my_min
942 hmax(ng)=my_max
943 ELSE
944 hmin(ng)=MIN(hmin(ng),my_min)
945 hmax(ng)=MAX(hmax(ng),my_max)
946 END IF
947 tile_count=tile_count+1
948 IF (tile_count.eq.NSUB) THEN
949 tile_count=0
950#ifdef DISTRIBUTE
951 buffer(1)=hmin(ng)
952 buffer(2)=hmax(ng)
953 op_handle(1)='MIN'
954 op_handle(2)='MAX'
955 CALL mp_reduce (ng, model, 2, buffer, op_handle)
956 hmin(ng)=buffer(1)
957 hmax(ng)=buffer(2)
958#endif
959 END IF
960!$OMP END CRITICAL (H_RANGE)
961#ifdef ICESHELF
962!
963!-----------------------------------------------------------------------
964! Set depth of ice shelf (meters; negative) at RHO-points.
965!-----------------------------------------------------------------------
966!
967# ifdef WEDDELL
968 val1=340.0_r8/16.0_r8
969 DO j=JstrR,JendR
970 DO i=IstrR,IendR
971 IF (i.gt.20) THEN
972 zice(i,j)=0.0_r8
973 ELSE IF (i.gt.4) THEN
974 zice(i,j)=-340.0_r8+REAL(i-1,r8)*val1
975 ELSE
976 zice(i,j)=-340.0_r8
977 END IF
978 END DO
979 END DO
980# else
981 DO j=JstrR,JendR
982 DO i=IstrR,IendR
983 zice(i,j)=0.0_r8
984 END DO
985 END DO
986# endif
987# if defined EW_PERIODIC || defined NS_PERIODIC
988 CALL exchange_r2d_tile (ng, Istr, Iend, Jstr, Jend, &
989 & LBi, UBi, LBj, UBj, &
990 & zice)
991# endif
992# ifdef DISTRIBUTE
993 CALL mp_exchange2d (ng, model, 1, Istr, Iend, Jstr, Jend, &
994 & LBi, UBi, LBj, UBj, &
995 & NghostPoints, EWperiodic, NSperiodic, &
996 & zice)
997# endif
998#endif
999 RETURN
1000 END SUBROUTINE ana_grid_tile