subroutine hvgrid c c======================================================================= c === c This routine sets the auxiliary arrays for the horizontal grid. === c It also sets the auxiliary arrays for the "flat" vertical grid. === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #include #include #include #include #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,j,k FLOAT * cenlat,cenlon,dely0,fac,tng,x,y,yc c c----------------------------------------------------------------------- c Begin executable code. c----------------------------------------------------------------------- c c----------------------------------------------------------------------- c Set grid spacing in Y. c----------------------------------------------------------------------- c if(coord.eq.0) then fac=gridy else fac=gridy*deg2rad*re endif do 10 j=1,jmt dyt(j)=fac 10 continue c c----------------------------------------------------------------------- c Set grid spacing in X. c----------------------------------------------------------------------- c if(coord.eq.0) then fac=gridx else fac=gridx*deg2rad*re endif do 20 i=1,imt dxt(i)=fac 20 continue #ifdef cyclic c c Set Cyclic boundary conditions on DXT c dxt(1 )=dxt(imtm1) dxt(imt)=dxt(2 ) #endif #ifdef gridold c c----------------------------------------------------------------------- c Calculate topography at U-V points. c----------------------------------------------------------------------- c do 30 j=1,jmtm1 do 30 i=1,imtm1 hv(i,j)=p25*(hr(i,j)+hr(i+1,j)+hr(i,j+1)+hr(i+1,j+1)) 30 continue do 40 j=1,jmtm1 hv(imt,j)=p5*(hr(imt,j)+hr(imt,j+1)) 40 continue do 50 i=1,imtm1 hv(i,jmt)=p5*(hr(i,jmt)+hr(i+1,jmt)) 50 continue hv(imt,jmt)=p5*(hv(imtm1,jmt)+hv(imt,jmtm1)) #endif #ifdef cyclic c c Set Cyclic conditions on HV. c do 60 j=1,jmt hv(imt,j)=hv(2,j) hv( 1,j)=hv(imt-1,j) 60 continue #endif c c----------------------------------------------------------------------- c Compute reciprocal depths HV and build bottom depth HDV at uv points. c----------------------------------------------------------------------- c do 70 j=1,jmt do 70 i=1,imt hdv(i,j)=hv(i,j) hv(i,j)=c1/hv(i,j) 70 continue c c----------------------------------------------------------------------- c Set the vertical "flat" grid box thicknesses. c----------------------------------------------------------------------- c have=c0 hc=c0 kc=iflag(8) do 90 k=1,kc hc=hc+hz(k) 90 continue do 100 k=kc+1,km have=have+hz(k) 100 continue do 110 k=1,kc dz(k)=hz(k)/hc 110 continue do 120 k=kc+1,km dz(k)=hz(k)/have 120 continue c c----------------------------------------------------------------------- c Set vertical levels to enable vectorization. c----------------------------------------------------------------------- c do 130 k=1,km kar(k)=k 130 continue c c----------------------------------------------------------------------- c Set the depths at the center of the vertical box for the "flat" grid. c----------------------------------------------------------------------- c refz(1)=p5*hz(1) #ifndef barotropic do 140 k=2,km refz(k)=refz(k-1)+p5*(hz(k-1)+hz(k)) 140 continue #endif c c----------------------------------------------------------------------- c Compute vertical grid at uv-points, and reset bottom depth HD c at tracer points. c----------------------------------------------------------------------- c call set_depths c c----------------------------------------------------------------------- c Compute reciprocal depths HR at tracer points. c----------------------------------------------------------------------- c do 80 j=1,jmt do 80 i=1,imt hr(i,j)=c1/hd(i,j) 80 continue c c======================================================================= c Compute auxiliary arrays based upon the spacing specified above: c======================================================================= c c----------------------------------------------------------------------- c Set horizontal grid spacing reciprocal factors. c----------------------------------------------------------------------- c do 150 j=1,jmtm1 dyu(j)=p5*(dyt(j)+dyt(j+1)) 150 continue dyu(jmt)=dyt(jmt) c do 160 j=1,jmt dytr(j)=c1/dyt(j) dyt2r(j)=p5/dyt(j) dyt4r(j)=p25/dyt(j) dyur(j)=c1/dyu(j) dyu2r(j)=p5/dyu(j) dyu4r(j)=p25/dyu(j) 160 continue c do 170 i=1,imtm1 dxu(i)=p5*(dxt(i)+dxt(i+1)) 170 continue #ifdef cyclic dxu(imt)=dxu(2) #else dxu(imt)=dxt(imt) #endif c do 180 i=1,imt dxtr(i)=c1/dxt(i) dxt2r(i)=p5/dxt(i) dxt4r(i)=p25/dxt(i) dxur(i)=c1/dxu(i) dxu2r(i)=p5/dxu(i) dxu4r(i)=p25/dxu(i) 180 continue c c----------------------------------------------------------------------- c Compute longitude/latitude at the tracers and velocity points. c Compute latitude in radians. c----------------------------------------------------------------------- c do 190 j=1,jmt do 190 i=1,imt x=FLoaT(i) y=FLoaT(j) call xy2ll (x,y,coord,imt,jmt,gridx,gridy,rlngd,rlatd,delx,dely, * thetad,tlon(i,j),tlat(i,j)) phit(i,j)=tlat(i,j)*deg2rad x=x+p5 y=y+p5 call xy2ll (x,y,coord,imt,jmt,gridx,gridy,rlngd,rlatd,delx,dely, * thetad,vlon(i,j),vlat(i,j)) phi(i,j)=vlat(i,j)*deg2rad 190 continue c c----------------------------------------------------------------------- c Set cosine and tangent metric factors. Also compute special metric c factors for the laplacian of velocity. c----------------------------------------------------------------------- c yc=p5*FLoaT(jmt+1) fac=gridy*deg2rad if ((coord.eq.1).or.(coord.eq.2)) dely0 = dely*deg2rad do 200 j=1,jmt if(coord.eq.0) then cst(j)=c1 cs(j)=c1 tng=c0 lpmtvl(j)=c0 elseif ((coord.eq.1).or.(coord.eq.2)) then y=dely0+(FLoaT(j)-yc)*fac cst(j)=cos(y) y=dely0+(FLoaT(j)-yc+p5)*fac cs(j)=cos(y) tng=tan(y)/re lpmtvl(j)=cm1/(re*re)-tng*tng endif lpmtgd(j)=c2*tng/cs(j) cstr(j)=c1/cst(j) csr(j)=c1/cs(j) 200 continue c c----------------------------------------------------------------------- c Set sine and tangent of U,V point latitudes for Coriolis terms: c IOPT(7)=0 full curvature c IOPT(7)=1 f-plane approximation c IOPT(7)=2 beta-plane approximation c----------------------------------------------------------------------- c x=FLoaT(imt+2)*p5 y=FLoaT(jmt+2)*p5 call xy2ll (x,y,coord,imt,jmt,gridx,gridy,rlngd,rlatd,delx,dely, * thetad,cenlon,cenlat) do 210 j=1,jmt do 210 i=1,imt if(iopt(7).eq.0) sine(i,j)=sin(phi(i,j)) if(iopt(7).eq.1) sine(i,j)=sin(cenlat) if(iopt(7).eq.2) sine(i,j)=sin(cenlat)+ * (phi(i,j)-cenlat)*cos(cenlat) 210 continue c c----------------------------------------------------------------------- c Write out grid geometry arrays c----------------------------------------------------------------------- c write(stdout,901) 901 format(/' Flat Grid Box Thickness (nondimensional): DZ'/) write(stdout,910) (dz(k),k=1,km) write(stdout,902) 902 format(/' Flat Grid Depths (cm) at center of vertical boxes: ', * 'REFZ'/) write(stdout,910) (refz(k),k=1,km) c if(ifprnt.eq.0) return c write(stdout,905) 905 format(/' Latitude of T,S Points (radians): PHIT'/) write(stdout,910) ((phit(i,j), i=1,imt), j=1,jmt) write(stdout,906) 906 format(/' Latitude of U,V Points (radians): PHI'/) write(stdout,910) ((phi(i,j), i=1,imt), j=1,jmt) write(stdout,907) 907 format(/' Cosine metric factor at T,S points: CST'/) write(stdout,910) (cst(j), j=1,jmt) write(stdout,908) 908 format(/' Cosine metric factor at U,V points: CS'/) write(stdout,910) (cs(j), j=1,jmt) write(stdout,909) 909 format(/' Coriolis sine of U,V Latitude: SINE'/) write(stdout,910) ((sine(i,j), i=1,imt), j=1,jmt) 910 format(5(1x,1pe13.5)) return end