119
120
126
128#ifdef DISTRIBUTE
130#endif
132
133
134
135 integer, intent(in) :: ng, tile, model
136 integer, intent(in) :: LBi, UBi, LBj, UBj
137 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
138
139#ifdef ASSUMED_SHAPE
140 real(r8), intent(out) :: angler(LBi:,LBj:)
141# if defined CURVGRID && defined UV_ADV
142 real(r8), intent(out) :: dmde(LBi:,LBj:)
143 real(r8), intent(out) :: dndx(LBi:,LBj:)
144# endif
145# ifdef ICESHELF
146 real(r8), intent(out) :: zice(LBi:,LBj:)
147# endif
148# ifdef SPHERICAL
149 real(r8), intent(out) :: lonp(LBi:,LBj:)
150 real(r8), intent(out) :: lonr(LBi:,LBj:)
151 real(r8), intent(out) :: lonu(LBi:,LBj:)
152 real(r8), intent(out) :: lonv(LBi:,LBj:)
153 real(r8), intent(out) :: latp(LBi:,LBj:)
154 real(r8), intent(out) :: latr(LBi:,LBj:)
155 real(r8), intent(out) :: latu(LBi:,LBj:)
156 real(r8), intent(out) :: latv(LBi:,LBj:)
157# else
158 real(r8), intent(out) :: xp(LBi:,LBj:)
159 real(r8), intent(out) :: xr(LBi:,LBj:)
160 real(r8), intent(out) :: xu(LBi:,LBj:)
161 real(r8), intent(out) :: xv(LBi:,LBj:)
162 real(r8), intent(out) :: yp(LBi:,LBj:)
163 real(r8), intent(out) :: yr(LBi:,LBj:)
164 real(r8), intent(out) :: yu(LBi:,LBj:)
165 real(r8), intent(out) :: yv(LBi:,LBj:)
166# endif
167 real(r8), intent(out) :: pn(LBi:,LBj:)
168 real(r8), intent(out) :: pm(LBi:,LBj:)
169 real(r8), intent(out) :: f(LBi:,LBj:)
170 real(r8), intent(out) :: h(LBi:,LBj:)
171#else
172 real(r8), intent(out) :: angler(LBi:UBi,LBj:UBj)
173# if defined CURVGRID && defined UV_ADV
174 real(r8), intent(out) :: dmde(LBi:UBi,LBj:UBj)
175 real(r8), intent(out) :: dndx(LBi:UBi,LBj:UBj)
176# endif
177# ifdef ICESHELF
178 real(r8), intent(out) :: zice(LBi:UBi,LBj:UBj)
179# endif
180# ifdef SPHERICAL
181 real(r8), intent(out) :: lonp(LBi:UBi,LBj:UBj)
182 real(r8), intent(out) :: lonr(LBi:UBi,LBj:UBj)
183 real(r8), intent(out) :: lonu(LBi:UBi,LBj:UBj)
184 real(r8), intent(out) :: lonv(LBi:UBi,LBj:UBj)
185 real(r8), intent(out) :: latp(LBi:UBi,LBj:UBj)
186 real(r8), intent(out) :: latr(LBi:UBi,LBj:UBj)
187 real(r8), intent(out) :: latu(LBi:UBi,LBj:UBj)
188 real(r8), intent(out) :: latv(LBi:UBi,LBj:UBj)
189# else
190 real(r8), intent(out) :: xp(LBi:UBi,LBj:UBj)
191 real(r8), intent(out) :: xr(LBi:UBi,LBj:UBj)
192 real(r8), intent(out) :: xu(LBi:UBi,LBj:UBj)
193 real(r8), intent(out) :: xv(LBi:UBi,LBj:UBj)
194 real(r8), intent(out) :: yp(LBi:UBi,LBj:UBj)
195 real(r8), intent(out) :: yr(LBi:UBi,LBj:UBj)
196 real(r8), intent(out) :: yu(LBi:UBi,LBj:UBj)
197 real(r8), intent(out) :: yv(LBi:UBi,LBj:UBj)
198# endif
199 real(r8), intent(out) :: pn(LBi:UBi,LBj:UBj)
200 real(r8), intent(out) :: pm(LBi:UBi,LBj:UBj)
201 real(r8), intent(out) :: f(LBi:UBi,LBj:UBj)
202 real(r8), intent(out) :: h(LBi:UBi,LBj:UBj)
203#endif
204
205
206
207 logical, save :: first = .true.
208
209 integer :: Imin, Imax, Jmin, Jmax
210 integer :: i, ival, j, k
211
212 real(r8),
parameter :: twopi = 2.0_r8*
pi
213
214 real(r8) :: Esize, Xsize, beta, cff, depth, dth
215 real(r8) :: dx, dy, f0, r, theta, val1, val2
216
217#ifdef WEDDELL
218 real(r8) :: hwrk(-1:235), xwrk(-1:235), zwrk
219#endif
220 real(r8) :: wrkX(IminS:ImaxS,JminS:JmaxS)
221 real(r8) :: wrkY(IminS:ImaxS,JminS:JmaxS)
222
223 TYPE (T_STATS), save :: Stats(16)
224
225#include "set_bounds.h"
226
227
228
229
230
231
232
233
234
235
236
237#if defined BASIN
238 xsize=3600.0e+03_r8
239 esize=2800.0e+03_r8
240 depth=5000.0_r8
241 f0=1.0e-04_r8
242 beta=2.0e-11_r8
243#elif defined BENCHMARK
244 xsize=360.0_r8
245 esize=20.0_r8
246 depth=4000.0_r8
247 f0=-1.0e-04_r8
248 beta=2.0e-11_r8
249#elif defined BL_TEST
250 xsize=100.0e+03_r8
251 esize=5.0e+03_r8
252 depth=47.5_r8
253 f0=9.25e-04_r8
254 beta=0.0_r8
255#elif defined CHANNEL
256 xsize=600.0e+03_r8
257 esize=360.0e+03_r8
258 depth=500.0_r8
259 f0=1.0e-04_r8
260 beta=0.0_r8
261#elif defined CANYON
262 xsize=128.0e+03_r8
263 esize=96.0e+03_r8
264 depth=4000.0_r8
265 f0=1.0e-04_r8
266 beta=0.0_r8
267#elif defined COUPLING_TEST
268 xsize=6000.0_r8*real(
lm(ng),r8)
269 esize=6000.0_r8*real(
mm(ng),r8)
270 depth=1500.0_r8
271 f0=5.0e-05_r8
272 beta=0.0_r8
273#elif defined DOUBLE_GYRE
274 xsize=1000.0e+03_r8
275 esize=2000.0e+03_r8
276 depth=500.0_r8
277
278 f0=7.3e-05_r8
279 beta=2.0e-11_r8
280#elif defined ESTUARY_TEST
281 xsize=100000.0_r8
282 esize=300.0_r8
283 depth=10.0_r8
284 f0=0.0_r8
285 beta=0.0_r8
286#elif defined KELVIN
287 xsize=20000.0_r8*real(
lm(ng),r8)
288 esize=20000.0_r8*real(
mm(ng),r8)
289 depth=100.0_r8
290 f0=1.0e-04_r8
291 beta=0.0_r8
292#elif defined FLT_TEST
293 xsize=1.0e+03_r8*real(
lm(ng),r8)
294 esize=1.0e+03_r8*real(
mm(ng),r8)
295 depth=10.0_r8
296 f0=0.0_r8
297 beta=0.0_r8
298#elif defined GRAV_ADJ
299 xsize=64.0e+03_r8
300
301 esize=
mm(ng)*xsize/
lm(ng)
302 depth=20.0_r8
303 f0=0.0_r8
304 beta=0.0_r8
305#elif defined LAB_CANYON
306 xsize=0.55_r8
309 beta=0.0_r8
310#elif defined LAKE_SIGNELL
311 xsize=50.0e3_r8
312 esize=10.0e3_r8
313 depth=18.0_r8
314 f0=0.0e-04_r8
315 beta=0.0_r8
316#elif defined LMD_TEST
317 xsize=100.0e+03_r8
318 esize=100.0e+03_r8
319 depth=50.0_r8
320 f0=1.09e-04_r8
321 beta=0.0_r8
322# elif defined MIXED_LAYER
323 xsize=500.0_r8
324 esize=400.0_r8
325 depth=50.0_r8
326 f0=0.0_r8
327 beta=0.0_r8
328#elif defined OVERFLOW
329 xsize=4.0e+03_r8
330 esize=200.0e+03_r8
331 depth=4000.0_r8
332 f0=0.0_r8
333 beta=0.0_r8
334#elif defined RIVERPLUME1
335 xsize=58.5e+03_r8
336 esize=201.0e+03_r8
337 depth=150.0_r8
338 f0=1.0e-04_r8
339 beta=0.0_r8
340#elif defined RIVERPLUME2
341 xsize=100.0e+03_r8
342 esize=210.0e+03_r8
343 depth=190.0_r8
344 f0=1.0e-04_r8
345 beta=0.0_r8
346#elif defined SEAMOUNT
347 xsize=320.0e+03_r8
348 esize=320.0e+03_r8
349 depth=5000.0_r8
350
351 f0=0.0_r8
352 beta=0.0_r8
353#elif defined SOLITON
354
355
356 xsize=48.0_r8
357 esize=16.0_r8
358 depth=1.0_r8
359 f0=0.0_r8
360 beta=1.0_r8
362#elif defined SED_TEST1
363 xsize=300.0_r8
364 esize=36.0_r8
365 depth=10.0_r8
366 f0=0.0_r8
367 beta=0.0_r8
368#elif defined SED_TOY
369 xsize=40.0_r8
370 esize=30.0_r8
371 depth=0.5_r8
372 f0=0.0_r8
373 beta=0.0_r8
374# elif defined SHOREFACE
375 xsize=1180.0_r8
376 esize=140.0_r8
377 depth=15.0_r8
378 f0=0.0e-04_r8
379 beta=0.0_r8
380#elif defined TEST_CHAN
381 xsize=10000.0_r8
382 esize=1000.0_r8
383 depth=10.0_r8
384 f0=0.0_r8
385 beta=0.0_r8
386#elif defined UPWELLING
387 xsize=1000.0_r8*real(
lm(ng),r8)
388 esize=1000.0_r8*real(
mm(ng),r8)
389 depth=150.0_r8
390 f0=-8.26e-05_r8
391 beta=0.0_r8
392#elif defined WEDDELL
393 xsize=4000.0_r8*real(
lm(ng),r8)
394 esize=4000.0_r8*real(
mm(ng),r8)
395 depth=4500.0_r8
396 f0=0.0_r8
397 beta=0.0_r8
398#elif defined WINDBASIN
399 xsize=2000.0_r8*real(
lm(ng),r8)
400 esize=1000.0_r8*real(
mm(ng),r8)
401 depth=50.0_r8
402 f0=1.0e-04_r8
403 beta=0.0_r8
404#else
405 ana_grid.h: no values provided for xsize, esize, depth, f0, beta.
406#endif
407
408
409
410 IF (
domain(ng)%NorthEast_Test(tile))
THEN
413 END IF
414
415
416
417
418
419 IF (first) THEN
420 first=.false.
421 DO i=1,SIZE(stats,1)
422 stats(i) % checksum=0_i8b
423 stats(i) % count=0
425 stats(i) % max=-
large
426 stats(i) % avg=0.0_r8
427 stats(i) % rms=0.0_r8
428 END DO
429 END IF
430 IF (
domain(ng)%NorthEast_Corner(tile))
WRITE (
stdout,
'(1x)')
431
432
433
434
435
436
437
438
439
440
441 IF (
domain(ng)%Western_Edge(tile))
THEN
442 imin=istr-1
443 ELSE
444 imin=istr
445 END IF
446 IF (
domain(ng)%Eastern_Edge(tile))
THEN
447 imax=iend+1
448 ELSE
449 imax=iend
450 END IF
451 IF (
domain(ng)%Southern_Edge(tile))
THEN
452 jmin=jstr-1
453 ELSE
454 jmin=jstr
455 END IF
456 IF (
domain(ng)%Northern_Edge(tile))
THEN
457 jmax=jend+1
458 ELSE
459 jmax=jend
460 END IF
461
462#if defined BENCHMARK
463
464
465
466 dx=xsize/real(
lm(ng),r8)
467 dy=esize/real(
mm(ng),r8)
469 DO j=jmin,jmax
470 val1=-70.0_r8+dy*(real(j,r8)-0.5_r8)
471 val2=-70.0_r8+dy*real(j,r8)
472 DO i=imin,imax
473 lonr(i,j)=dx*(real(i,r8)-0.5_r8)
474 latr(i,j)=val1
475 lonu(i,j)=dx*real(i,r8)
476 lonp(i,j)=lonu(i,j)
477 latu(i,j)=latr(i,j)
478 lonv(i,j)=lonr(i,j)
479 latv(i,j)=val2
480 latp(i,j)=latv(i,j)
481 END DO
482 END DO
483#elif defined LAB_CANYON
484
485
486
487 dx=xsize/real(
lm(ng),r8)
488 dy=esize/real(
mm(ng),r8)
489
490 dth=0.01_r8
491 cff=(4.0_r8*
pi/(dth*real(
mm(ng),r8)))-1.0_r8
492 DO j=jmin,jmax
493 DO i=imin,imax
494 r=0.35_r8+dx*real(i-1,r8)
496 & 0.5_r8*dth*((cff+1.0_r8)*real(j-1,r8)+ &
497 & (cff-1.0_r8)*(real(
mm(ng),r8)/twopi)* &
498 & sin(twopi*real(j-1,r8)/real(
mm(ng),r8)))
499 xp(i,j)=r*cos(theta)
500 yp(i,j)=r*sin(theta)
501 r=0.35_r8+dx*(real(i-1,r8)+0.5_r8)
503 & 0.5_r8*dth*((cff+1.0_r8)*(real(j-1,r8)+0.5_r8)+ &
504 & (cff-1.0_r8)*(real(
mm(ng),r8)/twopi)* &
505 & sin(twopi*(real(j-1,r8)+0.5_r8)/ &
507 xr(i,j)=r*cos(theta)
508 yr(i,j)=r*sin(theta)
509 xu(i,j)=xp(i,j)
510 yu(i,j)=yr(i,j)
511 xv(i,j)=xr(i,j)
512 yv(i,j)=yp(i,j)
513 END DO
514 END DO
515#else
516 dx=xsize/real(
lm(ng),r8)
517 dy=esize/real(
mm(ng),r8)
518 DO j=jmin,jmax
519 DO i=imin,imax
520# ifdef BL_TEST
521 dx=0.5_r8*(4000.0_r8/real(
lm(ng)+1,r8))*real(i,r8)+675.0_r8
522# endif
523 xp(i,j)=dx*real(i-1,r8)
524 xr(i,j)=dx*(real(i-1,r8)+0.5_r8)
525 xu(i,j)=xp(i,j)
526 xv(i,j)=xr(i,j)
527 yp(i,j)=dy*real(j-1,r8)
528 yr(i,j)=dy*(real(j-1,r8)+0.5_r8)
529 yu(i,j)=yr(i,j)
530 yv(i,j)=yp(i,j)
531 END DO
532 END DO
533#endif
534
535
536
537#ifdef SPHERICAL
539 & lbi, ubi, lbj, ubj, lonp)
540 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
541 WRITE (
stdout,10)
'longitude of PSI-points: lon_psi', &
542 & ng, stats(1)%min, stats(1)%max
543 END IF
545 & lbi, ubi, lbj, ubj, latp)
546 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
547 WRITE (
stdout,10)
'latitude of PSI-points: lat_psi', &
548 & ng, stats(2)%min, stats(2)%max
549 END IF
550
552 & lbi, ubi, lbj, ubj, lonr)
553 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
554 WRITE (
stdout,10)
'longitude of RHO-points: lon_rho', &
555 & ng, stats(3)%min, stats(3)%max
556 END IF
558 & lbi, ubi, lbj, ubj, latr)
559 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
560 WRITE (
stdout,10)
'latitude of RHO-points: lat_rho', &
561 & ng, stats(4)%min, stats(4)%max
562 END IF
563
565 & lbi, ubi, lbj, ubj, lonu)
566 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
567 WRITE (
stdout,10)
'longitude of U-points: lon_u', &
568 & ng, stats(5)%min, stats(5)%max
569 END IF
571 & lbi, ubi, lbj, ubj, latu)
572 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
573 WRITE (
stdout,10)
'latitude of U-points: lat_u', &
574 & ng, stats(6)%min, stats(6)%max
575 END IF
576
578 & lbi, ubi, lbj, ubj, lonv)
579 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
580 WRITE (
stdout,10)
'longitude of V-points: lon_v', &
581 & ng, stats(7)%min, stats(7)%max
582 END IF
584 & lbi, ubi, lbj, ubj, latv)
585 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
586 WRITE (
stdout,10)
'latitude of V-points: lat_v', &
587 & ng, stats(8)%min, stats(8)%max
588 END IF
589#else
591 & lbi, ubi, lbj, ubj, xp)
592 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
593 WRITE (
stdout,10)
'x-location of PSI-points: x_psi', &
594 & ng, stats(1)%min, stats(1)%max
595 END IF
597 & lbi, ubi, lbj, ubj, yp)
598 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
599 WRITE (
stdout,10)
'y-location of PSI-points: y_psi', &
600 & ng, stats(2)%min, stats(2)%max
601 END IF
602
604 & lbi, ubi, lbj, ubj, xr)
605 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
606 WRITE (
stdout,10)
'x-location of RHO-points: x_rho', &
607 & ng, stats(3)%min, stats(3)%max
608 END IF
610 & lbi, ubi, lbj, ubj, yr)
611 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
612 WRITE (
stdout,10)
'y-location of RHO-points: y_rho', &
613 & ng, stats(4)%min, stats(4)%max
614 END IF
615
617 & lbi, ubi, lbj, ubj, xu)
618 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
619 WRITE (
stdout,10)
'x-location of U-points: x_u', &
620 & ng, stats(5)%min, stats(5)%max
621 END IF
623 & lbi, ubi, lbj, ubj, yu)
624 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
625 WRITE (
stdout,10)
'y-location of U-points: y_u', &
626 & ng, stats(6)%min, stats(6)%max
627 END IF
628
630 & lbi, ubi, lbj, ubj, xv)
631 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
632 WRITE (
stdout,10)
'x-location of V-points: x_v', &
633 & ng, stats(7)%min, stats(7)%max
634 END IF
636 & lbi, ubi, lbj, ubj, yv)
637 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
638 WRITE (
stdout,10)
'y-location of V-points: y_v', &
639 & ng, stats(8)%min, stats(8)%max
640 END IF
641#endif
642
643#ifdef DISTRIBUTE
644
645
646
647# ifdef SPHERICAL
649 & lbi, ubi, lbj, ubj, &
651 & lonp, lonr, lonu, lonv)
653 & lbi, ubi, lbj, ubj, &
655 & latp, latr, latu, latv)
656# else
658 & lbi, ubi, lbj, ubj, &
660 & xp, xr, xu, xv)
662 & lbi, ubi, lbj, ubj, &
664 & yp, yr, yu, yv)
665# endif
666#endif
667
668
669
670
671
672
673
674#define J_RANGE MIN(JstrT,Jstr-1),MAX(Jend+1,JendT)
675#define I_RANGE MIN(IstrT,Istr-1),MAX(Iend+1,IendT)
676
677#if defined BENCHMARK
678
679
680
682 val2=real(
mm(ng),r8)*360.0_r8/(2.0_r8*
pi*
eradius*esize)
683 DO j=j_range
684 cff=1.0_r8/cos((-70.0_r8+dy*(real(j,r8)-0.5_r8))*
deg2rad)
685 DO i=i_range
686 wrkx(i,j)=val1*cff
687 wrky(i,j)=val2
688 END DO
689 END DO
690#elif defined LAB_CANYON
691
692
693
694 DO j=j_range
695 DO i=i_range
696 r=0.35_r8+dx*(real(i-1,r8)+0.5_r8)
697 theta=0.5_r8*dth*((cff+1.0_r8)+ &
698 & (cff-1.0_r8)* &
699 & cos(twopi*real(j-1,r8)/real(
mm(ng),r8)))
700 wrkx(i,j)=1.0_r8/dx
701 wrky(i,j)=1.0_r8/(r*theta)
702 END DO
703 END DO
704#else
705 DO j=j_range
706 DO i=i_range
707# ifdef BL_TEST
708 dx=0.5_r8*(4000.0_r8/real(
lm(ng)+1,r8))*real(i,r8)+675.0_r8
709# endif
710 wrkx(i,j)=1.0_r8/dx
711 wrky(i,j)=1.0_r8/dy
712 END DO
713 END DO
714#endif
715#undef J_RANGE
716#undef I_RANGE
717 DO j=jstrt,jendt
718 DO i=istrt,iendt
719 pm(i,j)=wrkx(i,j)
720 pn(i,j)=wrky(i,j)
721 END DO
722 END DO
723
724
725
727 & lbi, ubi, lbj, ubj, pm)
728 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
729 WRITE (
stdout,10)
'reciprocal XI-grid spacing: pm', &
730 & ng, stats(9)%min, stats(9)%max
731 END IF
733 & lbi, ubi, lbj, ubj, pn)
734 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
735 WRITE (
stdout,10)
'reciprocal ETA-grid spacing: pn', &
736 & ng, stats(10)%min, stats(10)%max
737 END IF
738
739
740
743 & lbi, ubi, lbj, ubj, &
744 & pm)
746 & lbi, ubi, lbj, ubj, &
747 & pn)
748 END IF
749
750#ifdef DISTRIBUTE
752 & lbi, ubi, lbj, ubj, &
755 & pm, pn)
756#endif
757
758#if (defined CURVGRID && defined UV_ADV)
759
760
761
762
763
764 DO j=jstr,jend
765 DO i=istr,iend
766 dndx(i,j)=0.5_r8*((1.0_r8/wrky(i+1,j ))- &
767 & (1.0_r8/wrky(i-1,j )))
768 dmde(i,j)=0.5_r8*((1.0_r8/wrkx(i ,j+1))- &
769 & (1.0_r8/wrkx(i ,j-1)))
770 END DO
771 END DO
772
773
774
776 & lbi, ubi, lbj, ubj, dmde)
777 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
778 WRITE (
stdout,10)
'ETA-derivative of inverse metric '// &
779 & 'factor pm: dmde', &
780 & ng, stats(11)%min, stats(11)%max
781 END IF
783 & lbi, ubi, lbj, ubj, dndx)
784 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
785 WRITE (
stdout,10)
'XI-derivative of inverse metric '// &
786 & 'factor pn: dndx', &
787 & ng, stats(12)%min, stats(12)%max
788 END IF
789
790
791
794 & lbi, ubi, lbj, ubj, &
795 & dndx)
797 & lbi, ubi, lbj, ubj, &
798 & dmde)
799 END IF
800
801# ifdef DISTRIBUTE
803 & lbi, ubi, lbj, ubj, &
806 & dndx, dmde)
807# endif
808#endif
809
810
811
812
813
814#if defined LAB_CANYON
815 DO j=jstrt,jendt
816 DO i=istrt,iendt
818 & 0.5_r8*dth*((cff+1.0_r8)*(real(j-1,r8)+0.5_r8)+ &
819 & (cff-1.0_r8)*(real(
mm(ng),r8)/twopi)* &
820 & sin(twopi*(real(j-1,r8)+0.5_r8)/ &
822 angler(i,j)=theta
823 END DO
824 END DO
825#elif defined WEDDELL
827 DO j=jstrt,jendt
828 DO i=istrt,iendt
829 angler(i,j)=val1
830 END DO
831 END DO
832#else
833 DO j=jstrt,jendt
834 DO i=istrt,iendt
835 angler(i,j)=0.0_r8
836 END DO
837 END DO
838#endif
839
840
841
843 & lbi, ubi, lbj, ubj, angler)
844 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
845 WRITE (
stdout,10)
'angle between XI-axis and EAST: '// &
846 & 'angler', &
847 & ng, stats(13)%min, stats(13)%max
848 END IF
849
850
851
854 & lbi, ubi, lbj, ubj, &
855 & angler)
856 END IF
857
858#ifdef DISTRIBUTE
860 & lbi, ubi, lbj, ubj, &
863 & angler)
864#endif
865
866
867
868
869
870#if defined BENCHMARK
871 val1=2.0_r8*(2.0_r8*
pi*366.25_r8/365.25_r8)/86400.0_r8
872 DO j=jstrt,jendt
873 DO i=istrt,iendt
874 f(i,j)=val1*sin(latr(i,j)*
deg2rad)
875 END DO
876 END DO
877#elif defined WEDDELL
878 val1=10.4_r8/real(
lm(ng),r8)
879 DO j=jstrt,jendt
880 DO i=istrt,iendt
881 f(i,j)=2.0_r8*7.2e-05_r8* &
882 & sin((-79.0_r8+real(i-1,r8)*val1)*
deg2rad)
883 END DO
884 END DO
885#else
886 IF (beta.eq.0.0_r8) THEN
887 DO j=jstrt,jendt
888 DO i=istrt,iendt
889 f(i,j)=f0
890 END DO
891 END DO
892 ELSE
893 val1=0.5_r8*esize
894 DO j=jstrt,jendt
895 DO i=istrt,iendt
896 f(i,j)=f0+beta*(yr(i,j)-val1)
897 END DO
898 END DO
899 END IF
900#endif
901
902
903
905 & lbi, ubi, lbj, ubj, f)
906 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
907 WRITE (
stdout,10)
'Coriolis parameter at RHO-points: f', &
908 & ng, stats(14)%min, stats(14)%max
909 END IF
910
911
912
915 & lbi, ubi, lbj, ubj, &
916 & f)
917 END IF
918
919#ifdef DISTRIBUTE
921 & lbi, ubi, lbj, ubj, &
924 & f)
925#endif
926
927
928
929
930
931#if defined BENCHMARK
932 DO j=jstrt,jendt
933 DO i=istrt,iendt
934 h(i,j)=500.0_r8+1750.0_r8*(1.0+tanh((68.0_r8+latr(i,j))/dy))
935 END DO
936 END DO
937#elif defined BL_TEST
938 DO j=jstrt,jendt
939 DO i=istrt,iendt
940 val1=(xr(i,j)+500.0_r8)/15000.0_r8
941 h(i,j)=14.0_r8+ &
942 & 25.0_r8*(1.0_r8-exp(-
pi*xr(i,j)*1.0e-05_r8))- &
943 & 8.0_r8*exp(-val1*val1)
944 END DO
945 END DO
946#elif defined CANYON
947 DO j=jstrt,jendt
948 DO i=istrt,iendt
949 val1=32000.0_r8-16000.0_r8*(sin(
pi*xr(i,j)/xsize))**24
950 h(i,j)=20.0_r8+0.5_r8*(depth-20.0_r8)* &
951 & (1.0_r8+tanh((yr(i,j)-val1)/10000.0_r8))
952 END DO
953 END DO
954#elif defined ESTUARY_TEST
955 DO j=jstrt,jendt
956 DO i=istrt,iendt
957 h(i,j)=5.0_r8+(xsize-xr(i,j))/xsize*5.0_r8
958 END DO
959 END DO
960#elif defined LAB_CANYON
961 DO j=jstrt,jendt
962 DO i=istrt,iendt
963 r=0.35_r8+dx*(real(i-1,r8)+0.5_r8)
965 & 0.5_r8*dth*((cff+1.0_r8)*(real(j-1,r8)+0.5_r8)+ &
966 & (cff-1.0_r8)*(real(
mm(ng),r8)/twopi)* &
967 & sin(dth*(real(j-1,r8)+0.5_r8)/ &
969 val1=0.55_r8-0.15_r8*(cos(
pi*theta*0.55_r8/0.2_r8)**2)
970 val2=0.15_r8+0.15_r8*(cos(
pi*theta*0.55_r8/0.2_r8)**2)
971 IF (abs(theta).ge.0.181818181818_r8) THEN
972 IF (r.le.0.55_r8) THEN
973 h(i,j)=0.025_r8
974 ELSE IF (r.ge.0.7_r8) THEN
975 h(i,j)=0.125_r8
976 ELSE
977 h(i,j)=0.125_r8-0.1_r8* &
978 & (cos(0.5_r8*
pi*(r-0.55_r8)/0.15_r8)**2)
979 END IF
980 ELSE
981 IF (r.le.val1) THEN
982 h(i,j)=0.025_r8
983 ELSE IF (r.ge.0.7_r8) THEN
984 h(i,j)=0.125_r8
985 ELSE
986 h(i,j)=0.125_r8-0.1_r8* &
987 & (cos(0.5_r8*
pi*(r-val1)/val2)**2)
988 END IF
989 END IF
990 END DO
991 END DO
992#elif defined LAKE_SIGNELL
993 DO j=jstrt,jendt
994 DO i=istrt,iendt
995 h(i,j)=18.0_r8-16.0_r8*real(
mm(ng)-j,r8)/real(
mm(ng)-1,r8)
996 END DO
997 END DO
998# elif defined MIXED_LAYER
999 DO j=jstrt,jendt
1000 DO i=istrt,iendt
1001 h(i,j)=50.0_r8
1002 END DO
1003 END DO
1004#elif defined OVERFLOW
1005 val1=200.0_r8
1006 DO j=jstrt,jendt
1007 DO i=istrt,iendt
1008 h(i,j)=val1+0.5_r8*(depth-val1)* &
1009 & (1.0_r8+tanh((yr(i,j)-100000.0_r8)/20000.0_r8))
1010 END DO
1011 END DO
1012#elif defined RIVERPLUME1
1013 DO j=jstrt,jendt
1014 DO i=istrt,min(5,iendt)
1015 h(i,j)=15.0_r8
1016 END DO
1017 DO i=max(6,istrt),iendt
1018 h(i,j)=depth+real(
lm(ng)-i,r8)*(15.0_r8-depth)/ &
1020 END DO
1021 END DO
1022#elif defined RIVERPLUME2
1023 DO j=jstrt,jendt
1024 DO i=istrt,min(5,iendt)
1025 h(i,j)=15.0_r8
1026 END DO
1027 DO i=max(6,istrt),iendt
1028 h(i,j)=depth+real(
lm(ng)-i,r8)*(15.0_r8-depth)/ &
1030 END DO
1031 END DO
1032#elif defined SEAMOUNT
1033 DO j=jstrt,jendt
1034 DO i=istrt,iendt
1035 val1=(xr(i,j)-0.5_r8*xsize)/40000.0_r8
1036 val2=(yr(i,j)-0.5_r8*esize)/40000.0_r8
1037 h(i,j)=depth-4500.0_r8*exp(-(val1*val1+val2*val2))
1038 END DO
1039 END DO
1040#elif defined SED_TOY
1041 DO j=jstrt,jendt
1042 DO i=istrt,iendt
1043 h(i,j)=20.0_r8
1044 END DO
1045 END DO
1046#elif defined SHOREFACE
1047 DO j=jstrt,jendt
1048 DO i=istrt,iendt
1049 h(i,j)=11.75_r8-0.0125_r8*xsize/real(
lm(ng)+1,r8)*real(i,r8)
1050 END DO
1051 END DO
1052#elif defined TEST_CHAN
1053 DO j=jstrt,jendt
1054 DO i=istrt,iendt
1055 h(i,j)=10.0_r8+0.4040_r8*real(i,r8)/real(
lm(ng)+1,r8)
1056 END DO
1057 END DO
1058#elif defined UPWELLING
1060 DO i=istrt,iendt
1061 IF (i.le.
lm(ng)/2)
THEN
1062 val1=real(i,r8)
1063 ELSE
1064 val1=real(
lm(ng)+1-i,r8)
1065 END IF
1066 val2=min(depth,84.5_r8+66.526_r8*tanh((val1-10.0_r8)/7.0_r8))
1067 DO j=jstrt,jendt
1068 h(i,j)=val2
1069 END DO
1070 END DO
1072 DO j=jstrt,jendt
1073 IF (j.le.
mm(ng)/2)
THEN
1074 val1=real(j,r8)
1075 ELSE
1076 val1=real(
mm(ng)+1-j,r8)
1077 END IF
1078 val2=min(depth,84.5_r8+66.526_r8*tanh((val1-10.0_r8)/7.0_r8))
1079 DO i=istrt,iendt
1080 h(i,j)=val2
1081 END DO
1082 END DO
1083 END IF
1084#elif defined WEDDELL
1085 val1=98.80_r8
1086 val2=0.8270_r8
1087 DO k=-1,26
1088 xwrk(k)=real(k-1,r8)*15.0_r8*1000.0_r8
1089 hwrk(k)=375.0_r8
1090 END DO
1091 DO k=27,232
1092 zwrk=-2.0_r8+real(k-1,r8)*0.020_r8
1093 xwrk(k)=(520.0_r8+val1+zwrk*val1+ &
1094 & val1*val2*log(cosh(zwrk)))*1000.0_r8
1095 hwrk(k)=-75.0_r8+2198.0_r8*(1.0_r8+val2*tanh(zwrk))
1096 END DO
1097 DO k=233,235
1098 xwrk(k)=(850.0_r8+real(k-228,r8)*50.0_r8)*1000.0_r8
1099 hwrk(k)=4000.0_r8
1100 END DO
1101 DO j=jstrt,jendt
1102 DO i=istrt,iendt
1103 h(i,j)=375.0_r8
1104 DO k=1,234
1105 IF ((xwrk(k).le.xr(i,1)).and.(xr(i,1).lt.xwrk(k+1))) THEN
1106 cff=1.0_r8/(xwrk(k+1)-xwrk(k))
1107 h(i,j)=cff*(xwrk(k+1)-xr(i,j))*hwrk(k )+ &
1108 & cff*(xr(i,j)-xwrk(k ))*hwrk(k+1)
1109 END IF
1110 END DO
1111 END DO
1112 END DO
1113#elif defined WINDBASIN
1114 DO i=istrt,iendt
1115 ival=int(0.03_r8*real(
lm(ng)+1,r8))
1116 IF (i.lt.ival) THEN
1117 val1=1.0_r8-(real((i+1)-ival,r8)/real(ival,r8))**2
1118 ELSE IF ((
lm(ng)+1-i).lt.ival)
THEN
1119 val1=1.0_r8-(real((
lm(ng)+1-i)-ival,r8)/real(ival,r8))**2
1120 ELSE
1121 val1=1.0_r8
1122 END IF
1123 DO j=jstrt,jendt
1124 val2=2.0_r8*real(j-(
mm(ng)+1)/2,r8)/real(
mm(ng)+1,r8)
1125 h(i,j)=depth*(0.08_r8+0.92_r8*val1*(1.0_r8-val2*val2))
1126 END DO
1127 END DO
1128#else
1129 DO j=jstrt,jendt
1130 DO i=istrt,iendt
1131 h(i,j)=depth
1132 END DO
1133 END DO
1134#endif
1135
1136
1137
1139 & lbi, ubi, lbj, ubj, h)
1140 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1141 WRITE (
stdout,10)
'bathymetry at RHO-points: h', &
1142 & ng, stats(15)%min, stats(15)%max
1143 END IF
1144 hmin(ng)=stats(15)%min
1145 hmax(ng)=stats(15)%max
1146
1147
1148
1151 & lbi, ubi, lbj, ubj, &
1152 & h)
1153 END IF
1154
1155#ifdef DISTRIBUTE
1157 & lbi, ubi, lbj, ubj, &
1160 & h)
1161#endif
1162#ifdef ICESHELF
1163
1164
1165
1166
1167
1168# ifdef WEDDELL
1169 val1=340.0_r8
1170 val2=val1/16.0_r8
1171 DO j=jstrt,jendt
1172 DO i=istrt,iendt
1173 IF (i.gt.20) THEN
1174 zice(i,j)=0.0_r8
1175 ELSE IF (i.gt.4) THEN
1176 zice(i,j)=-val1+real(i-1,r8)*val2
1177 ELSE
1178 zice(i,j)=-val1
1179 END IF
1180 END DO
1181 END DO
1182# else
1183 DO j=jstrt,jendt
1184 DO i=istrt,iendt
1185 zice(i,j)=0.0_r8
1186 END DO
1187 END DO
1188# endif
1189
1190
1191
1193 & lbi, ubi, lbj, ubj, zice)
1194 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1195 WRITE (
stdout,10)
'ice shelf thickness: zice', &
1196 & ng, stats(16)%min, stats(16)%max
1197 END IF
1198
1199
1200
1203 & lbi, ubi, lbj, ubj, &
1204 & zice)
1205 END IF
1206
1207# ifdef DISTRIBUTE
1209 & lbi, ubi, lbj, ubj, &
1212 & zice)
1213# endif
1214#endif
1215
1216 10 FORMAT (3x,' ANA_GRID - ',a,/,19x, &
1217 & '(Grid = ',i2.2,', Min = ',1p,e15.8,0p, &
1218 & ' Max = ',1p,e15.8,0p,')')
1219
1220 RETURN
subroutine ana_grid(ng, tile, model)
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
type(t_domain), dimension(:), allocatable domain
integer, dimension(:), allocatable lm
integer, parameter u2dvar
integer, parameter p2dvar
integer, dimension(:), allocatable mm
integer, parameter r2dvar
integer, parameter v2dvar
real(dp), dimension(:), allocatable hmin
real(r8), dimension(:), allocatable el
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
real(dp), parameter large
real(dp), parameter deg2rad
real(r8), dimension(:), allocatable xl
real(dp), dimension(:), allocatable hmax
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine, public stats_2dfld(ng, tile, model, gtype, s, extract_flag, lbi, ubi, lbj, ubj, f, fmask, debug)