102 & LBi, UBi, LBj, UBj, &
103 & IminS, ImaxS, JminS, JmaxS, &
105#if defined CURVGRID && defined UV_ADV
112 & lonp, lonr, lonu, lonv, &
113 & latp, latr, latu, latv, &
135 integer,
intent(in) :: ng, tile, model
136 integer,
intent(in) :: LBi, UBi, LBj, UBj
137 integer,
intent(in) :: IminS, ImaxS, JminS, JmaxS
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:)
146 real(r8),
intent(out) :: zice(LBi:,LBj:)
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:)
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:)
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:)
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)
178 real(r8),
intent(out) :: zice(LBi:UBi,LBj:UBj)
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)
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)
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)
207 logical,
save :: first = .true.
209 integer :: Imin, Imax, Jmin, Jmax
210 integer :: i, ival, j, k
212 real(r8),
parameter :: twopi = 2.0_r8*
pi
214 real(r8) :: Esize, Xsize, beta, cff, depth, dth
215 real(r8) :: dx, dy, f0, r, theta, val1, val2
218 real(r8) :: hwrk(-1:235), xwrk(-1:235), zwrk
220 real(r8) :: wrkX(IminS:ImaxS,JminS:JmaxS)
221 real(r8) :: wrkY(IminS:ImaxS,JminS:JmaxS)
223 TYPE (T_STATS),
save :: Stats(16)
225#include "set_bounds.h"
243#elif defined BENCHMARK
267#elif defined COUPLING_TEST
268 xsize=6000.0_r8*real(
lm(ng),r8)
269 esize=6000.0_r8*real(
mm(ng),r8)
273#elif defined DOUBLE_GYRE
280#elif defined ESTUARY_TEST
287 xsize=20000.0_r8*real(
lm(ng),r8)
288 esize=20000.0_r8*real(
mm(ng),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)
298#elif defined GRAV_ADJ
301 esize=
mm(ng)*xsize/
lm(ng)
305#elif defined LAB_CANYON
310#elif defined LAKE_SIGNELL
316#elif defined LMD_TEST
322# elif defined MIXED_LAYER
328#elif defined OVERFLOW
334#elif defined RIVERPLUME1
340#elif defined RIVERPLUME2
346#elif defined SEAMOUNT
362#elif defined SED_TEST1
374# elif defined SHOREFACE
380#elif defined TEST_CHAN
386#elif defined UPWELLING
387 xsize=1000.0_r8*real(
lm(ng),r8)
388 esize=1000.0_r8*real(
mm(ng),r8)
393 xsize=4000.0_r8*real(
lm(ng),r8)
394 esize=4000.0_r8*real(
mm(ng),r8)
398#elif defined WINDBASIN
399 xsize=2000.0_r8*real(
lm(ng),r8)
400 esize=1000.0_r8*real(
mm(ng),r8)
405 ana_grid.h: no values provided for xsize, esize, depth, f0, beta.
410 IF (
domain(ng)%NorthEast_Test(tile))
THEN
422 stats(i) % checksum=0_i8b
425 stats(i) % max=-
large
426 stats(i) % avg=0.0_r8
427 stats(i) % rms=0.0_r8
430 IF (
domain(ng)%NorthEast_Corner(tile))
WRITE (
stdout,
'(1x)')
441 IF (
domain(ng)%Western_Edge(tile))
THEN
446 IF (
domain(ng)%Eastern_Edge(tile))
THEN
451 IF (
domain(ng)%Southern_Edge(tile))
THEN
456 IF (
domain(ng)%Northern_Edge(tile))
THEN
466 dx=xsize/real(
lm(ng),r8)
467 dy=esize/real(
mm(ng),r8)
470 val1=-70.0_r8+dy*(real(j,r8)-0.5_r8)
471 val2=-70.0_r8+dy*real(j,r8)
473 lonr(i,j)=dx*(real(i,r8)-0.5_r8)
475 lonu(i,j)=dx*real(i,r8)
483#elif defined LAB_CANYON
487 dx=xsize/real(
lm(ng),r8)
488 dy=esize/real(
mm(ng),r8)
491 cff=(4.0_r8*
pi/(dth*real(
mm(ng),r8)))-1.0_r8
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)))
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)/ &
516 dx=xsize/real(
lm(ng),r8)
517 dy=esize/real(
mm(ng),r8)
521 dx=0.5_r8*(4000.0_r8/real(
lm(ng)+1,r8))*real(i,r8)+675.0_r8
523 xp(i,j)=dx*real(i-1,r8)
524 xr(i,j)=dx*(real(i-1,r8)+0.5_r8)
527 yp(i,j)=dy*real(j-1,r8)
528 yr(i,j)=dy*(real(j-1,r8)+0.5_r8)
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
649 & lbi, ubi, lbj, ubj, &
651 & lonp, lonr, lonu, lonv)
653 & lbi, ubi, lbj, ubj, &
655 & latp, latr, latu, latv)
658 & lbi, ubi, lbj, ubj, &
662 & lbi, ubi, lbj, ubj, &
674#define J_RANGE MIN(JstrT,Jstr-1),MAX(Jend+1,JendT)
675#define I_RANGE MIN(IstrT,Istr-1),MAX(Iend+1,IendT)
682 val2=real(
mm(ng),r8)*360.0_r8/(2.0_r8*
pi*
eradius*esize)
684 cff=1.0_r8/cos((-70.0_r8+dy*(real(j,r8)-0.5_r8))*
deg2rad)
690#elif defined LAB_CANYON
696 r=0.35_r8+dx*(real(i-1,r8)+0.5_r8)
697 theta=0.5_r8*dth*((cff+1.0_r8)+ &
699 & cos(twopi*real(j-1,r8)/real(
mm(ng),r8)))
701 wrky(i,j)=1.0_r8/(r*theta)
708 dx=0.5_r8*(4000.0_r8/real(
lm(ng)+1,r8))*real(i,r8)+675.0_r8
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
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
743 & lbi, ubi, lbj, ubj, &
746 & lbi, ubi, lbj, ubj, &
752 & lbi, ubi, lbj, ubj, &
758#if (defined CURVGRID && defined UV_ADV)
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)))
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
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
794 & lbi, ubi, lbj, ubj, &
797 & lbi, ubi, lbj, ubj, &
803 & lbi, ubi, lbj, ubj, &
814#if defined LAB_CANYON
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)/ &
843 & lbi, ubi, lbj, ubj, angler)
844 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
845 WRITE (
stdout,10)
'angle between XI-axis and EAST: '// &
847 & ng, stats(13)%min, stats(13)%max
854 & lbi, ubi, lbj, ubj, &
860 & lbi, ubi, lbj, ubj, &
871 val1=2.0_r8*(2.0_r8*
pi*366.25_r8/365.25_r8)/86400.0_r8
874 f(i,j)=val1*sin(latr(i,j)*
deg2rad)
878 val1=10.4_r8/real(
lm(ng),r8)
881 f(i,j)=2.0_r8*7.2e-05_r8* &
882 & sin((-79.0_r8+real(i-1,r8)*val1)*
deg2rad)
886 IF (beta.eq.0.0_r8)
THEN
896 f(i,j)=f0+beta*(yr(i,j)-val1)
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
915 & lbi, ubi, lbj, ubj, &
921 & lbi, ubi, lbj, ubj, &
934 h(i,j)=500.0_r8+1750.0_r8*(1.0+tanh((68.0_r8+latr(i,j))/dy))
940 val1=(xr(i,j)+500.0_r8)/15000.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)
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))
954#elif defined ESTUARY_TEST
957 h(i,j)=5.0_r8+(xsize-xr(i,j))/xsize*5.0_r8
960#elif defined LAB_CANYON
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
974 ELSE IF (r.ge.0.7_r8)
THEN
977 h(i,j)=0.125_r8-0.1_r8* &
978 & (cos(0.5_r8*
pi*(r-0.55_r8)/0.15_r8)**2)
983 ELSE IF (r.ge.0.7_r8)
THEN
986 h(i,j)=0.125_r8-0.1_r8* &
987 & (cos(0.5_r8*
pi*(r-val1)/val2)**2)
992#elif defined LAKE_SIGNELL
995 h(i,j)=18.0_r8-16.0_r8*real(
mm(ng)-j,r8)/real(
mm(ng)-1,r8)
998# elif defined MIXED_LAYER
1004#elif defined OVERFLOW
1008 h(i,j)=val1+0.5_r8*(depth-val1)* &
1009 & (1.0_r8+tanh((yr(i,j)-100000.0_r8)/20000.0_r8))
1012#elif defined RIVERPLUME1
1014 DO i=istrt,min(5,iendt)
1017 DO i=max(6,istrt),iendt
1018 h(i,j)=depth+real(
lm(ng)-i,r8)*(15.0_r8-depth)/ &
1022#elif defined RIVERPLUME2
1024 DO i=istrt,min(5,iendt)
1027 DO i=max(6,istrt),iendt
1028 h(i,j)=depth+real(
lm(ng)-i,r8)*(15.0_r8-depth)/ &
1032#elif defined SEAMOUNT
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))
1040#elif defined SED_TOY
1046#elif defined SHOREFACE
1049 h(i,j)=11.75_r8-0.0125_r8*xsize/real(
lm(ng)+1,r8)*real(i,r8)
1052#elif defined TEST_CHAN
1055 h(i,j)=10.0_r8+0.4040_r8*real(i,r8)/real(
lm(ng)+1,r8)
1058#elif defined UPWELLING
1061 IF (i.le.
lm(ng)/2)
THEN
1064 val1=real(
lm(ng)+1-i,r8)
1066 val2=min(depth,84.5_r8+66.526_r8*tanh((val1-10.0_r8)/7.0_r8))
1073 IF (j.le.
mm(ng)/2)
THEN
1076 val1=real(
mm(ng)+1-j,r8)
1078 val2=min(depth,84.5_r8+66.526_r8*tanh((val1-10.0_r8)/7.0_r8))
1084#elif defined WEDDELL
1088 xwrk(k)=real(k-1,r8)*15.0_r8*1000.0_r8
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))
1098 xwrk(k)=(850.0_r8+real(k-228,r8)*50.0_r8)*1000.0_r8
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)
1113#elif defined WINDBASIN
1115 ival=int(0.03_r8*real(
lm(ng)+1,r8))
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
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))
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
1144 hmin(ng)=stats(15)%min
1145 hmax(ng)=stats(15)%max
1151 & lbi, ubi, lbj, ubj, &
1157 & lbi, ubi, lbj, ubj, &
1175 ELSE IF (i.gt.4)
THEN
1176 zice(i,j)=-val1+real(i-1,r8)*val2
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
1203 & lbi, ubi, lbj, ubj, &
1209 & lbi, ubi, lbj, ubj, &
1216 10
FORMAT (3x,
' ANA_GRID - ',a,/,19x, &
1217 &
'(Grid = ',i2.2,
', Min = ',1p,e15.8,0p, &
1218 &
' Max = ',1p,e15.8,0p,
')')