c$Source: /usr/our/senya/work/model/MC_PG/senq/RCS/dyn_subs.f,v $ c$Author: senya $ c$Revision: 0.4 $ c$Date: 94/01/24 11:04:57 $ c$State: Exp $ c ------------------------------------------------------------------ subroutine aarea(npt,nz,lxxk,lyxk,emx,emy,area,basin,isk,nzi,h) c ------------------------------------------------------------------ c compute the stretched grid area factors for computing integrals c over the grid in the stretched coordinate system. c dx*dy = (d(x)/dpsi1(x) * d(y)/dpsi2(y)) * dpsi1 * dpsi2 c for transformation (stretching) functions dpsi1 and dpsi2. c c from common block grid: c nxp,nyp = (input) grid dimensions in the x and y directions. c nxyc = (input) number of ocean grid points. c land = (input) land storage flag. c nbx,nby = (input) number of x and y boundaries. c ncs = (input) number of interior corner boundary grid points. c c lxxk = (input) nbx indices of the x-boundaries for an x c (compressed or regular) sort. c lyxk = (input) nby indices of the y-boundaries for an x c (compressed or regular) sort. c emx = (input) factor for x-differencing = d(psi1)/d(x)*(1./delx) c emy = (input) factor for y-differencing = d(psi2)/d(y)*(1./dely) c isk = (input) index from compressed k to compressed k=1 points c area = (output) .5*dx*dy. c basin = (output) .5*(total basin area). c include 'comm_para.h' implicit real(a-h,o-z),integer(i-n) dimension lxxk(MXBDY,nz),lyxk(MXBDY,nz) dimension emx(1),emy(1),area(npt,nz),basin(nz) dimension nzi(npt),h(npt,nz) dimension isk(npt,nz) common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) c r2 = 0.5 do k = 1, nz c c find 1/2 the area (1/2 factor used on energy calculations) c as if all grid squares were completely ocean. c do i=1,nptk(k) j = isk(i,k) area(j,k) = r2/(emx(j)*emy(j)) enddo c c correct the x-boundary grid point areas. c do i=1,nbxk(k) area(lxxk(i,k),k) = area(lxxk(i,k),k)*.5 enddo c c correct the y-boundary grid point areas. c do i=1,nbyk(k) area(lyxk(i,k),k) = area(lyxk(i,k),k)*.5 enddo c c correct the interior corner point areas. c do i=nbxk(k)+1,nbxk(k)+ncsk(k) area(lxxk(i,k),k) = area(lxxk(i,k),k)*.75 enddo basin(k) = 0. do i=1,nptk(k) j = isk(i,k) basin(k) = basin(k) + area(j,k) enddo enddo c end of aarea. end c c ------------------------------------------------------------------ subroutine aarea_new(npt,nxp,nyp,nz,lxxk,lyyk,snxk,snyk,mask,iox, * emx,emy,area,basin,isk,nzi,h,isyk,lpbcwk,lpbcek, * saxmk,saxpk,saymk,saypk,tp,sxp,syp,sxm,sym,mbox) c ------------------------------------------------------------------ c compute the stretched grid area factors for computing integrals c over the grid in the stretched coordinate system. c dx*dy = (d(x)/dpsi1(x) * d(y)/dpsi2(y)) * dpsi1 * dpsi2 c for transformation (stretching) functions dpsi1 and dpsi2. c c from common block grid: c nxp,nyp = (input) grid dimensions in the x and y directions. c nxyc = (input) number of ocean grid points. c land = (input) land storage flag. c nbx,nby = (input) number of x and y boundaries. c ncs = (input) number of interior corner boundary grid points. c c lxxk = (input) nbx indices of the x-boundaries for an x c (compressed or regular) sort. c lyxk = (input) nby indices of the y-boundaries for an x c (compressed or regular) sort. c emx = (input) factor for x-differencing = d(psi1)/d(x)*(1./delx) c emy = (input) factor for y-differencing = d(psi2)/d(y)*(1./dely) c isk = (input) index from compressed k to compressed k=1 points c area = (output) .5*dx*dy. c basin = (output) .5*(total basin area). c include 'comm_para.h' implicit real(a-h,o-z),integer(i-n) dimension lxxk(MXBDY,nz),lyyk(MXBDY,nz),snxk(MAXNB,nz),snyk(MAXNB,nz) dimension lpbcwk(MAXSID,nz),lpbcek(MAXSID,nz) dimension emx(1),emy(1),area(npt,nz),basin(nz) dimension nzi(npt),h(npt,nz) dimension saxmk(npt,nz),saxpk(npt,nz),saymk(npt,nz),saypk(npt,nz), * tp(npt,4),sxp(npt,nz),syp(npt,nz),sxm(npt,nz),sym(npt,nz) dimension isk(npt,nz),isyk(npt,nz),mask(nxp,nyp,nz),iox(npt) common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ) + ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ) dimension ig(4) c r2 = 0.5 do k = 1, nz c c find 1/2 the area (1/2 factor used on energy calculations) c as if all grid squares were completely ocean. c npk = nptk(k) npbk = npbck(k) do i=1,npk j = isk(i,k) area(j,k) = r2/(emx(j)*emy(j)) saxmk(j,k) = emx(j) saxpk(j,k) = emx(j) saymk(j,k) = emy(j) saypk(j,k) = emy(j) enddo if (mbox.eq.2) then do j = 2, nyp- 1 do i = 2, nxp- 1 ipt = mask(i,j,k) if (ipt.gt.0) then jpt = isk(ipt,k) do m=1,4 ig(m) = 0 enddo c find out what type of point if (mask(i-1,j-1,k).eq.0) ig(1) = 1 if (mask(i+1,j-1,k).eq.0) ig(2) = 1 if (mask(i-1,j+1,k).eq.0) ig(3) = 1 if (mask(i+1,j+1,k).eq.0) ig(4) = 1 if (mask(i-1,j,k)*mask(i,j-1,k).eq.0) ig(1) = 1 if (mask(i+1,j,k)*mask(i,j-1,k).eq.0) ig(2) = 1 if (mask(i-1,j,k)*mask(i,j+1,k).eq.0) ig(3) = 1 if (mask(i+1,j,k)*mask(i,j+1,k).eq.0) ig(4) = 1 la = ig(1) + ig(2) + ig(3) + ig(4) c set metrics if (la.eq.3) then area(jpt,k) = 0.25*area(jpt,k) saxmk(jpt,k) = 2.0* emx(jpt) saxpk(jpt,k) = 2.0* emx(jpt) saymk(jpt,k) = 2.0* emy(jpt) saypk(jpt,k) = 2.0* emy(jpt) elseif (la.eq.2) then area(jpt,k) = 0.50*area(jpt,k) if (ig(1).ne.ig(2)) then saxmk(jpt,k) = 2.0* emx(jpt) saxpk(jpt,k) = 2.0* emx(jpt) if (ig(1).ne.ig(3)) then saxmk(jpt,k) = emx(jpt) saxpk(jpt,k) = emx(jpt) endif else saymk(jpt,k) = 2.0* emy(jpt) saypk(jpt,k) = 2.0* emy(jpt) endif elseif (la.eq.1) then area(jpt,k) = 0.75*area(jpt,k) if (ig(1).eq.1) then saxmk(jpt,k) = 2./3. * emx(jpt) saxpk(jpt,k) = 4./3. * emx(jpt) saymk(jpt,k) = 2./3. * emy(jpt) saypk(jpt,k) = 4./3. * emy(jpt) elseif (ig(2).eq.1) then saxmk(jpt,k) = 4./3. * emx(jpt) saxpk(jpt,k) = 2./3. * emx(jpt) saymk(jpt,k) = 2./3. * emy(jpt) saypk(jpt,k) = 4./3. * emy(jpt) elseif (ig(3).eq.1) then saxmk(jpt,k) = 2./3. * emx(jpt) saxpk(jpt,k) = 4./3. * emx(jpt) saymk(jpt,k) = 4./3. * emy(jpt) saypk(jpt,k) = 2./3. * emy(jpt) else saxmk(jpt,k) = 4./3. * emx(jpt) saxpk(jpt,k) = 2./3. * emx(jpt) saymk(jpt,k) = 4./3. * emy(jpt) saypk(jpt,k) = 2./3. * emy(jpt) endif endif endif enddo enddo c do north/south boundary points, as if not periodic c southern boundary j = 1 jp= 2 do i = 2 , nxp-1 ipt = mask(i,j,k) if (ipt.gt.0) then jpt = isk(ipt,k) area (jpt,k) = 0.50*area(jpt,k) saymk(jpt,k) = 2.0*emy(jpt) saypk(jpt,k) = 2.0*emy(jpt) la = 0 if (mask(i-1,j,k)*mask(i-1,jp,k)*mask(i,jp,k) * *mask(i+1,j,k)*mask(i+1,jp,k).eq.0) la = 1 if (la.eq.1) then ! a corner point area (jpt,k) = 0.50*area(jpt,k) saxmk(jpt,k) = 2.0*emx(jpt) saxpk(jpt,k) = 2.0*emx(jpt) endif endif enddo c northern boundary j = nyp jm= nyp-1 do i = 2 , nxp-1 ipt = mask(i,j,k) if (ipt.gt.0) then jpt = isk(ipt,k) area (jpt,k) = 0.50*area(jpt,k) saymk(jpt,k) = 2.0*emy(jpt) saypk(jpt,k) = 2.0*emy(jpt) la = 0 if (mask(i-1,j,k)*mask(i-1,jm,k)*mask(i,jm,k) * *mask(i+1,j,k)*mask(i+1,jm,k).eq.0) la = 1 if (la.eq.1) then ! a corner point area (jpt,k) = 0.50*area(jpt,k) saxmk(jpt,k) = 2.0*emx(jpt) saxpk(jpt,k) = 2.0*emx(jpt) endif endif enddo c do east/west boundary points, as if not periodic c western boundary i = 1 i2= 2 do j = 2 , nyp-1 ipt = mask(i,j,k) if (ipt.gt.0) then jpt = isk(ipt,k) area (jpt,k) = 0.50*area(jpt,k) saxmk(jpt,k) = 2.0*emx(jpt) saxpk(jpt,k) = 2.0*emx(jpt) la = 0 if (mask(i2,j-1,k)*mask(i,j-1,k)*mask(i2,j+1,k) * *mask(i2,j ,k)*mask(i,j+1,k).eq.0) la = 1 if (la.eq.1) then ! a corner point area (jpt,k) = 0.50*area(jpt,k) saymk(jpt,k) = 2.0*emy(jpt) saypk(jpt,k) = 2.0*emy(jpt) endif endif enddo c eastern boundary i = nxp i2= nxp-1 do j = 2 , nyp-1 ipt = mask(i,j,k) if (ipt.gt.0) then jpt = isk(ipt,k) area (jpt,k) = 0.50*area(jpt,k) saxmk(jpt,k) = 2.0*emx(jpt) saxpk(jpt,k) = 2.0*emx(jpt) la = 0 if (mask(i2,j-1,k)*mask(i,j-1,k)*mask(i2,j+1,k) * *mask(i2,j ,k)*mask(i,j+1,k).eq.0) la = 1 if (la.eq.1) then ! a corner point area (jpt,k) = 0.50*area(jpt,k) saymk(jpt,k) = 2.0*emy(jpt) saypk(jpt,k) = 2.0*emy(jpt) endif endif enddo ipt = mask(1,1,k) if (ipt.gt.0) then jpt = isk(ipt,k) area (jpt,k) = 0.25*area(jpt,k) saxmk(jpt,k) = 2.0*emx(jpt) saxpk(jpt,k) = 2.0*emx(jpt) saymk(jpt,k) = 2.0*emy(jpt) saypk(jpt,k) = 2.0*emy(jpt) endif ipt = mask(nxp,1,k) if (ipt.gt.0) then jpt = isk(ipt,k) area (jpt,k) = 0.25*area(jpt,k) saxmk(jpt,k) = 2.0*emx(jpt) saxpk(jpt,k) = 2.0*emx(jpt) saymk(jpt,k) = 2.0*emy(jpt) saypk(jpt,k) = 2.0*emy(jpt) endif ipt = mask(1,nyp,k) if (ipt.gt.0) then jpt = isk(ipt,k) area (jpt,k) = 0.25*area(jpt,k) saxmk(jpt,k) = 2.0*emx(jpt) saxpk(jpt,k) = 2.0*emx(jpt) saymk(jpt,k) = 2.0*emy(jpt) saypk(jpt,k) = 2.0*emy(jpt) endif ipt = mask(nxp,nyp,k) if (ipt.gt.0) then jpt = isk(ipt,k) area (jpt,k) = 0.25*area(jpt,k) saxmk(jpt,k) = 2.0*emx(jpt) saxpk(jpt,k) = 2.0*emx(jpt) saymk(jpt,k) = 2.0*emy(jpt) saypk(jpt,k) = 2.0*emy(jpt) endif c fix periodic boundary points do iper = 1, npbk jpt = lpbcwk(iper,k) area(jpt,k) = r2/(emx(jpt)*emy(jpt)) saxmk(jpt,k) = emx(jpt) saxpk(jpt,k) = emx(jpt) saymk(jpt,k) = emy(jpt) saypk(jpt,k) = emy(jpt) ixy = iox(jpt) i = 1 j = (ixy-1)/nxp + 1 jp = j+1 jm = j-1 im = nxp ip = 2 do m=1,4 ig(m) = 0 enddo if (jm.gt.0) then if (mask(im,j,k)*mask(i,jm,k)*mask(im,jm,k).eq.0) ig(1) = 1 if (mask(ip,j,k)*mask(i,jm,k)*mask(ip,jm,k).eq.0) ig(2) = 1 else ig(1) = 1 ig(2) = 1 endif if (jp.lt.nyp+1) then if (mask(im,j,k)*mask(i,jp,k)*mask(im,jp,k).eq.0) ig(3) = 1 if (mask(ip,j,k)*mask(i,jp,k)*mask(ip,jp,k).eq.0) ig(4) = 1 else ig(3) = 1 ig(4) = 1 endif la = ig(1) + ig(2) + ig(3) + ig(4) if (la.eq.3) then area(jpt,k) = 0.25*area(jpt,k) saxmk(jpt,k) = 2.0* emx(jpt) saxpk(jpt,k) = 2.0* emx(jpt) saymk(jpt,k) = 2.0* emy(jpt) saypk(jpt,k) = 2.0* emy(jpt) elseif (la.eq.2) then area(jpt,k) = 0.50*area(jpt,k) if (ig(1).ne.ig(2)) then saxmk(jpt,k) = 2.0* emx(jpt) saxpk(jpt,k) = 2.0* emx(jpt) if (ig(1).ne.ig(3)) then saxmk(jpt,k) = emx(jpt) saxpk(jpt,k) = emx(jpt) endif else saymk(jpt,k) = 2.0* emy(jpt) saypk(jpt,k) = 2.0* emy(jpt) endif elseif (la.eq.1) then area(jpt,k) = 0.75*area(jpt,k) saxmk(jpt,k) = 2./3. * emx(jpt) saxpk(jpt,k) = 2./3. * emx(jpt) saymk(jpt,k) = 2./3. * emy(jpt) saypk(jpt,k) = 2./3. * emy(jpt) if (ig(1).eq.1) then saxpk(jpt,k) = 2.* saxpk(jpt,k) saypk(jpt,k) = 2.* saypk(jpt,k) elseif (ig(2).eq.1) then saxmk(jpt,k) = 2.* saxmk(jpt,k) saypk(jpt,k) = 2.* saypk(jpt,k) elseif (ig(3).eq.1) then saxpk(jpt,k) = 2.* saxpk(jpt,k) saymk(jpt,k) = 2.* saymk(jpt,k) else saxmk(jpt,k) = 2.* saxmk(jpt,k) saymk(jpt,k) = 2.* saymk(jpt,k) endif endif jpt = lpbcek(iper,k) area(jpt,k) = r2/(emx(jpt)*emy(jpt)) saxmk(jpt,k) = emx(jpt) saxpk(jpt,k) = emx(jpt) saymk(jpt,k) = emy(jpt) saypk(jpt,k) = emy(jpt) ixy = iox(jpt) i = nxp j = (ixy-1)/nxp + 1 jp = j+1 jm = j-1 im = nxp - 1 ip = 1 do m=1,4 ig(m) = 0 enddo if (jm.gt.0) then if (mask(im,j,k)*mask(i,jm,k)*mask(im,jm,k).eq.0) ig(1) = 1 if (mask(ip,j,k)*mask(i,jm,k)*mask(ip,jm,k).eq.0) ig(2) = 1 else ig(1) = 1 ig(2) = 1 endif if (jp.lt.nyp+1) then if (mask(im,j,k)*mask(i,jp,k)*mask(im,jp,k).eq.0) ig(3) = 1 if (mask(ip,j,k)*mask(i,jp,k)*mask(ip,jp,k).eq.0) ig(4) = 1 else ig(3) = 1 ig(4) = 1 endif la = ig(1) + ig(2) + ig(3) + ig(4) if (la.eq.3) then area(jpt,k) = 0.25*area(jpt,k) saxmk(jpt,k) = 2.0* emx(jpt) saxpk(jpt,k) = 2.0* emx(jpt) saymk(jpt,k) = 2.0* emy(jpt) saypk(jpt,k) = 2.0* emy(jpt) elseif (la.eq.2) then area(jpt,k) = 0.50*area(jpt,k) if (ig(1).ne.ig(2)) then saxmk(jpt,k) = 2.0* emx(jpt) saxpk(jpt,k) = 2.0* emx(jpt) if (ig(1).ne.ig(3)) then saxmk(jpt,k) = emx(jpt) saxpk(jpt,k) = emx(jpt) endif else saymk(jpt,k) = 2.0* emy(jpt) saypk(jpt,k) = 2.0* emy(jpt) endif elseif (la.eq.1) then area(jpt,k) = 0.75*area(jpt,k) saxmk(jpt,k) = 2./3. * emx(jpt) saxpk(jpt,k) = 2./3. * emx(jpt) saymk(jpt,k) = 2./3. * emy(jpt) saypk(jpt,k) = 2./3. * emy(jpt) if (ig(1).eq.1) then saxpk(jpt,k) = 2.* saxpk(jpt,k) saypk(jpt,k) = 2.* saypk(jpt,k) elseif (ig(2).eq.1) then saxmk(jpt,k) = 2.* saxmk(jpt,k) saypk(jpt,k) = 2.* saypk(jpt,k) elseif (ig(3).eq.1) then saxpk(jpt,k) = 2.* saxpk(jpt,k) saymk(jpt,k) = 2.* saymk(jpt,k) else saxmk(jpt,k) = 2.* saxmk(jpt,k) saymk(jpt,k) = 2.* saymk(jpt,k) endif endif enddo endif basin(k) = 0. do i=1,npk j = isk(i,k) basin(k) = basin(k) + area(j,k) sxm(j,k) = 1./emx(j) sym(j,k) = 1./emy(j) enddo if (mbox.eq.2) then do i=1,nbxk(k) j = lxxk(i,k) sxm(j,k) = sxm(j,k)/2. enddo do i=1,nbyk(k) j = lyyk(i,k) sym(j,k) = sym(j,k)/2. enddo endif do j = 1, nyp do i = 1, nxp ipt = mask(i,j,k) if (ipt.gt.0) then jpt = isk(ipt,k) ip = min(nxp,i+1) jp = min(nyp,j+1) jn = mask(i,jp,k) je = mask(ip,j,k) if (jn.eq.0) jn = jpt if (je.eq.0) je = jpt sxp(jpt,k) = 0.25/max(saxpk(jpt,k),saxmk(je,k)) syp(jpt,k) = 0.25/max(saypk(jpt,k),saymk(jn,k)) endif enddo enddo if (mbox.eq.2) then do i=1,nbxk(k) j = lxxk(i,k) if (snxk(i,k).gt.0) sxp(j,k) = 0.50/saxpk(j,k) if (snxk(i,k).lt.0) sxp(j-1,k) = 0.50/saxmk(j,k) enddo nyoff = nbyk(k) + ncsk(k) do i=1,nbyk(k) j = lyyk(i,k) jm = lyyk(nyoff+i,k) if (snyk(i,k).gt.0) syp(j,k) = 0.50/saypk(j,k) if (snyk(i,k).lt.0) syp(jm,k) = 0.50/saymk(j,k) enddo endif enddo return c end of aarea. end c subroutine scaset(iox,x,y,xp,yp,f,emx,emy,cosl,cosli,emx2,emy2,tp, * dx2i,dy2i,ycos) c ------------------------------------------------------------------ c subroutine called by wdrivn that computes the coordinates of the c grid points and length scale variables with and without stretching, c e.g. arrays emx, emy, and f. c c from common block grid: c c nxp,nyp = (input) grid dimensions in the x and y directions. c nxyc = (input) number of ocean grid points. c c from common block coords: c c alon,blon = (input) min,max x grid coordinates in degrees. c alat,blat = (input) min,max y grid coordinates in degrees. c rlat = (input) reference latitude for beta plane and f plane. c mgrid = (input) determines coordinate system (stretching c determined by nsx, nsy). c = 1; beta plane with delx = (blon-alon)/((nxp-1)*rearth). c = 2; spherical coords with c delx = (blon-alon)((nxp-1)*rearth*cos(y(i)) c (for spherical coords, delx is a function of the c convergence of meridians away from the equator.) c = 3; f-plane with c delx = (blon-alon)((nxp-1)*rearth). c nsx,nsy = (input) see routine stretch for a description. c = 0; no stretching. c c iox = (input) nxyc indices of the x-sorted ocean grid points. c x,y = (output) nxp x and nyp y grid point coordinates (degrees). c xp = (output) nxp derivatives of the x-transformation function c if coordinate stretching is used (see gridxy). c yp = (output) nyp derivatives of the y-transformation function c if coordniate stretching is used (see gridxy). c emx = (output) factor for x-differencing= d(psi1)/d(x)*(1./delx) c emy = (output) factor for y-differencing= d(psi2)/d(y)*(1./dely) c f = (output) coriolis factor for routine dhoriz. c modified to account for the Pole's shift (Senya Basin, 1996) c include 'comm_para.h' implicit real(a-h,o-z),integer(i-n) parameter (REARTH = 6378000., * DTOSEC = 86400., * RAD = 3.14159265/180., * TOMEGA = 4.*3.14159265/DTOSEC) common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch common/strech/xs(maxxs),alpha(maxxs),beta(maxxs) common /new_geom/ ixs_type,iys_type,ipole,pole_alp,pole_bet,pole_gam c c dimension iox(nxyc),x(nxp),y(nyp),xp(nxp),yp(nyp),emx(nxyc),emy(nxyc) * ,f(nxyc),emx2(1), emy2(1), tp(nxyc), cosl(nxyc),cosli(nxyc) * ,dx2i(0:nxp+1,nyp),dy2i(nxp,nyp),ycos(0:nyp+1) c c compute the grid point x and y coordinates in degrees. c call gridxy(nxp,nyp,alon,blon,alat,blat,nsx,nsy,nystrch,xs,alpha,beta, + x,y,xp,yp,tp,tp(nxp+1)) delx = (blon-alon)*RAD/float(nxp-1) dely = (blat-alat)*RAD/float(nyp-1) c c convert to inverse delta x and delta y scale in m. c xfac = 1./(delx*REARTH) yfac = 1./(dely*REARTH) rlat = 0.5*(alat + blat) if ( (mgrid.eq.1) .or. (mgrid.eq.3) ) then c........beta plane or f plane. fz = sin(rlat*rad) betap = cos(rlat*rad) cosy = 1./betap c.............................for on f-plane. if (mgrid .eq.3 ) betap = 0. do k = 1, nxyc j = (iox(k)-1)/nxp + 1 i = iox(k) - (j-1)*nxp f(k) = (fz + betap*RAD*(y(j)-rlat))*TOMEGA emx(k) = xfac*xp(i)*cosy emy(k) = yfac*yp(j) cosl(k) = 1. cosli(k) = 1. emx2(k) = xfac*xfac*tp(i)*cosy*cosy emy2(k) = yfac*yfac*tp(nxp+j) enddo do j = 1, nyp do i = 1, nxp dx2i(i,j) = 0.5*xfac*xp(i)*cosy dy2i(i,j) = 0.5*yfac*yp(j) ycos(j) = 1. enddo dx2i(0 ,j) = dx2i(1,j) dx2i(nxp+1,j) = dx2i(nxp,j) enddo elseif (mgrid .eq. 2) then c.......................for spherical coordinates. do k = 1, nxyc c convert the x-sort index to i,j. j = (iox(k)-1)/nxp + 1 i = iox(k) - (j-1)*nxp yrad = y(j)*RAD fj = sin(yrad)*TOMEGA cosy = 1./cos(yrad) xyfac = -tan(yrad)/REARTH cos2 = cosy*cosy f(k) = fj emx(k) = xfac*xp(i)*cosy emy(k) = yfac*yp(j) emx2(k) = xfac*xfac*tp(i)*cos2 emy2(k) = yfac*yfac*tp(nxp+j) + xyfac cosli(k) = cosy cosl(k) = 1./cosy enddo do j = 1, nyp yrad = y(j)*RAD cosy = 1./cos(yrad) do i = 1, nxp dx2i(i,j) = 0.5*xfac*xp(i)*cosy dy2i(i,j) = 0.5*yfac*yp(j) ycos(j) = 1./cosy enddo dx2i(0 ,j) = dx2i(nxp,j) dx2i(nxp+1,j) = dx2i(1,j) enddo endif ycos(0) = ycos(1) ycos(nyp+1) = ycos(nyp) if (ipole .eq. 1) then c...........................re-compute Coriolis for a rotated Pole: do k = 1, nxyc j = (iox(k)-1)/nxp + 1 i = iox(k) - (j-1)*nxp f(k) = TOMEGA * rot_fcr2g(x(i), y(j)) enddo endif return c end of scaset. end c c ------------------------------------------------------------------ subroutine stretch(nx,xmin,xmax,ns,xs,alpha,beta,x,xp,xpp) c ------------------------------------------------------------------ c compute grid point locations for a stretched coordinate system. c c nx = (input) # of grid points in a coordinate direction. c xmin,xmax = (input) coordinate range of the nx grid points. c ns = (input) # of atan's defining the transformation from c the stretched coordinates to the regular coordinates. c xs = (input) ns locations of the atan's. c alpha = (input) ns scaling parameters for the atan's in degrees. c beta = (input) ns scaling parameters for the atan's in degrees. c x = (output) nx stretched grid point locations. c xp = (output) nx derivatives of the transformation function c d(psi(x))/d(x). c xpp = (output) nx second derivatives of the transformation c function: d^2(psi(x))/d(x)^2. (by Senya Basin) c c the transformation from the stretched coordinates x to the regular c coordinates xx is c xx = psi(x) = a + b*(x + sum of alpha(i)*atan((x-xs(i))/beta(i)) c + c*x**2 c c the extra degree of freedom provided by c*x**2 is used to locate c a grid point right at the equator. note the definition of psi c differs slightly from cane & gent & ncar code. alpha can be c modified to reproduce ncar stretching. c c for nx equally spaced values of xx from xmin to xmax, this routine c finds the corresponding values of the coord x. the stretched grid c point spacing will be the smallest near x=xs(i). increasing the c parameter alpha(i) increases the # of stretched grid points in the c vicinity of xs(i). increasing beta(i) widens the region around c xs(i) in which grid points are concentrated. if xmin*xmax.ge.0. c (basin does not span the equator) then the parameter c is set to c 0., and a and b are chosen by this routine to force the ends of c the computational and physical grids to be the same: c psi(xmin)=xmin, psi(xmax)=xmax. if xmin*xmax .lt. 0., then the c parameters a, b, and c are chosen to force psi(xmin)=xmin, c psi(0.)=0., and psi(xmax)=xmax. c c nbegns : the function psi(x) is initially computed at nbegns c equally spaced points between xmin and xmax to provide c a table of starting values for the iterative solution. c c maxit = the maximum allowed number of iterations. c c eps = the iterative solution for x is continued until c abs( (x_this_iter. - x_last_iter.)/x_last_iter. ).le.eps c implicit real(a-h,o-z),integer(i-n) parameter (nbegns=100, maxit=1000, eps=1.e-6) dimension xs(1),alpha(1),beta(1),x(1),xp(1),xpp(1) c c solve for a, b, (and c) by construing the physical and c computational grids at the end points (and equator). c if(xmin*xmax.ge.0.) then c c get the scaling factors a and b, which will force c psi(xmin)=xmin and psi(xmax)=xmax. c c evaluate psi at the end points as if a=0 and b=1. c ymin = psi(xmin,0.,1.,0.,ns,xs,alpha,beta) ymax = psi(xmax,0.,1.,0.,ns,xs,alpha,beta) c solving for two equations and two unknowns yields: a = (xmin*ymax-ymin*xmax)/(ymax-ymin) b = (xmax-xmin)/(ymax-ymin) c = 0. else c c for a basin that includes the equator. c c get the scaling factors a,b and c, which will force c psi(xmin)=xmin, psi(xmax)=xmax, and psi(0.) = 0. c c evaluate psi at the end points and equator as if c a=c=0 and b=1. c ymin = psi(xmin,0.,1.,0.,ns,xs,alpha,beta) ymax = psi(xmax,0.,1.,0.,ns,xs,alpha,beta) y0 = psi( 0.,0.,1.,0.,ns,xs,alpha,beta) x1 = xmin*xmin x2 = xmax*xmax c solving three equations for three unknowns: b = (xmin*x2-xmax*x1)/(x2*(ymin-y0)-x1*(ymax-y0)) c = (xmin - b*(ymin-y0))/x1 a = -b*y0 endif c kstart = 1 delx = (xmax-xmin)/float(nbegns-1) psi1 = psi(xmin,a,b,c,ns,xs,alpha,beta) c c loop over nx values of xx evenly spaced from xmin to xmax c and find the corresponding values of x such that xx = psi(x). c dx = (xmax-xmin)/float(nx-1) do 50 i = 1,nx xx = xmin + (i-1)*dx c c use Newton's method to find x for f(x) = psi(x)-xx = 0: c f(x+delx) ~= f(x) + f'(x)*delx + f"(x)/2*delx**2 + ... c for f(x+delx)=0, delx = -f(x)/f'(x) c xj1 = xj - f(xj))/fp(xj); fp(x) = dpsi = d(psi(x))/d(x) c c first find psi1 and psi2 straddling xy such that c psi1.le.xx .and. xx.le.psi2 and iterate on x from there. c 20 psi2 = psi(xmin+kstart*delx,a,b,c,ns,xs,alpha,beta) if (psi2 .ge. xx) goto 30 if (kstart .lt. nbegns-1) then c if not, increase x. kstart = kstart + 1 psi1 = psi2 goto 20 endif c c interpolate between two values of psi to get a starting xj. c 30 xj1 = xmin + (kstart-1)*delx + (xx-psi1)*delx/(psi2-psi1) c c loop until abs((xj1-xj)/xj).le.eps or iter.eq.maxit c iter = 0 40 xj = xj1 f = psi(xj,a,b,c,ns,xs,alpha,beta) - xx fp = dpsi(xj,a,b,c,ns,xs,alpha,beta) xj1 = xj - f/fp iter = iter + 1 if(xj .ne. 0.) then if((abs((xj1-xj)/xj).gt.eps).and.(iter.lt.maxit)) goto 40 else if(abs(xj1) .gt. eps) goto 40 endif c store x location. x(i)=xj c store psi derivative: xp = d(psi)/d(x) xp(i) = fp c xpp = d^2(psi)/dx^2 xpp(i) = d2psi(xj,a,b,c,ns,xs,alpha,beta) 50 continue c ccc fix the ends: x(1) = xmin x(nx) = xmax ccc return c end of stretch. end c ------------------------------------------------------------------ real function d2psi(x,a,b,c,ns,xs,alpha,beta) c ------------------------------------------------------------------ c Second derivative of the stretching function. (Senya Basin) c d2x = b * sum { -alpha/beta^2 * (2e/(1+e^2)^2 } + 2*c; e = (x-xs)/beta implicit real(a-h,o-z),integer(i-n) dimension xs(1),alpha(1),beta(1) c sum = 0. do i = 1, ns binv = 1./beta(i) e = binv*(x - xs(i)) e2 = 1. + e*e sum = sum + alpha(i)*binv*binv*(e + e)/(e2 * e2) enddo d2psi = -b*sum + 2.*c return c end of function d2psi. end c ------------------------------------------------------------------ real function dpsi(x,a,b,c,ns,xs,alpha,beta) c ------------------------------------------------------------------ c derivative of the stretching function. c dimension xs(1),alpha(1),beta(1) c sum = 0. do 10 i=1,ns e = (x - xs(i))/beta(i) 10 sum = sum + (alpha(i)/beta(i)) * 1./(1.+ e*e) dpsi = b*(1. + sum) + 2.*c*x return c end of function dpsi. end c c ------------------------------------------------------------------ real function psi(x,a,b,c,ns,xs,alpha,beta) c ------------------------------------------------------------------ c coordinate stretching function. c dimension xs(1),alpha(1),beta(1) c sum = x do 10 i=1,ns 10 sum = sum + alpha(i)*atan((x-xs(i))/beta(i)) psi = a + b*sum + c*x*x return c end of function psi. end c subroutine wspace(parm,need) c ------------------------------------------------------------------ c write an error message about the required array space and exit. c c parm = (input) six-character name of the dimension parameter. c need = (input) needed value for the parameter parm, or zero. c = 0, then the needed value is not included in the message. c implicit real(a-h,o-z),integer(i-n) character*6 parm character*72 msg if(need.ne.0) then write(msg,1) parm,need 1 format('insufficient space, increase dimension parameter ', + a6,' to ',i6,'$') else write(msg,2) parm 2 format('insufficient space, increase dimension parameter ',a6, + '$') endif call perror1(msg,1) return c end of wspace. end c ------------------------------------------------------------------ subroutine perror1(msg,istop) c ------------------------------------------------------------------ c print the character string msg and exit the program if istop.ne.0 c c msg = (input) character string containing a message to be c printed. c istop = (input) stop flag c = 0, continue execution c = otherwise, exit the program c ioerr = (input) unit number for error messages. c character*(*) msg character*72 err common/errors/ioerr,nstep c c check if ioerr is reasonable. if not, write to unit 6. c if(ioerr.ge.1 .and. ioerr.le.99) then io = ioerr else io = 6 endif c if(nstep .gt. 0) then if(istop.eq.0) then write(io,1) nstep 1 format(1x,'warning on step ',i10) else write(io,2) nstep 2 format(1x,'error exit on step ',i10) endif endif if(len(msg) .gt. 0) then l = len(msg) iend = 0 do 10 i=1,72 if(msg(i:i) .eq. '$') iend = 1 if(i .gt. l .or. iend .eq. 1) then err(i:i) = ' ' else err(i:i) = msg(i:i) endif 10 continue write(io,11) err 11 format(1x,a72) endif c c close all open output data files and stop execution. c if(istop.ne.0) then call cstop stop endif return c end of perror. end c