c subroutine ccm2gcm(imxtgcm,jmxtgcm,imxccm,jmxccm,nflds, + xccm,yccm,xtgcm,ytgcm,fccm,ftgcm,ivec) c c Transform from ccm2 to tgcm grid. c c tgcm : imxtgcm = 72, jmxtgcm = 36 c c T21 : imxccm = 64, jmxccm = 32 c maccm2: imxccm = 64, jmxccm = 64 c T42 : imxccm = 128, jmxccm = 64 c real xccm(imxccm),yccm(jmxccm),xtgcm(imxtgcm),ytgcm(jmxtgcm) real fccm(imxccm,jmxccm,nflds),ftgcm(imxtgcm,jmxtgcm,nflds) integer ivec(nflds) c do ip=1,nflds call transgrd(fccm(1,1,ip),imxccm,jmxccm,ftgcm(1,1,ip),imxtgcm, + jmxtgcm,xccm,yccm,xtgcm,ytgcm,ivec(ip)) enddo return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine gcm2ccm(imxtgcm,jmxtgcm,imxccm,jmxccm,nflds, + xccm,yccm,xtgcm,ytgcm,fccm,ftgcm,ivec) c c Transform fields ftgcm on the tgcm grid to fccm on the ccm2 grid: c c tgcm : imxtgcm = 72, jmxtgcm = 36 c T21 : imxccm = 64, jmxccm = 32 c maccm2: imxccm = 64, jmxccm = 64 c T42 : imxccm = 128, jmxccm = 64 c real xccm(imxccm),yccm(jmxccm),xtgcm(imxtgcm),ytgcm(jmxtgcm) real fccm(imxccm,jmxccm,nflds),ftgcm(imxtgcm,jmxtgcm,nflds) integer ivec(nflds) real gcmlat(jmxtgcm+2),gcmlon(imxtgcm+1) ! for coarse2fine c do ip=1,nflds call transgrd(ftgcm(1,1,ip),imxtgcm,jmxtgcm,fccm(1,1,ip), + imxccm,jmxccm,xtgcm,ytgcm,xccm,yccm,ivec(ip)) enddo return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE TRANSGRD(F1,NX1,NY1,F2,NX2,NY2,X1,Y1,X2,Y2,ITYPE) SAVE C **** C **** This subroutine transforms a function defined on global C **** grid(1) to global grid(2) using John Adam's method C **** (bicubic but not bicubic spline) C **** C **** Input: C **** C **** F1(NX1,NY1) = array on grid(1) C **** C **** NX1, NY1 = dimensions of array F1 C **** C **** X1(NX1) = array of longitude values for grid(1) C **** in degrees starting at -180. (no wrap around) C **** C **** Y1(NY1) = array of latitude values for grid(1) in C **** degrees. (no polar values) C **** C **** NX2,NY2 = dimensions of array F2 on grid(2) to which C **** F1 is to be transformed C **** C **** X2(NX2) = array of longitude values for grid(2) starting C **** at -180. (no wrap around) C **** C **** Y2(NY2) = array of latitude values for grid(2) C **** (no polar values) C **** C **** ITYPE = 0 for scalar F1 C **** = 1 for vector F1 (e.g. horizontal velocities) C **** C **** Output: C **** C **** F2(NX2,NY2) = transformed version of F1 C **** DIMENSION F1(NX1,NY1), F2(NX2,NY2), X1(NX1), Y1(NY1), X2(NX2), 1 Y2(NY2), INTPOL(2) C **** C **** Where: IBD defines bicubic spline boundary conditions. C **** (periodic in this case) C **** C **** Set up pointers to work space needed for interpolation C **** procedure. C **** POINTER (PFF1,FF1(0:NX1+2,2*NY1+1)), (PXX1,XX1(0:NX1+2)), 1 (PYY1,YY1(2*NY1+1)), 4 (PWK,WK(5*NY2+9*NX2)) C **** C **** Where: C **** C **** FF1 = space for extended version of F1. This is C **** periodic in both latitude and longitude, C **** wrap around points are added. C **** C **** XX1 = array of longitude values with wrap around points. C **** (in degrees) C **** C **** YY1 = array of latitude values running from -180 degrees C **** to +180 degrees. Periodic points are C **** added. C **** WK = work space for tlc package C **** DATA INTPOL/2*3/ C **** C **** Initialization C **** C **** Allocate work space from heap C **** CALL HPALLOC(PFF1,(NX1+3)*(2*NY1+1),IER1,0) CALL HPALLOC(PXX1,NX1+3,IER2,0) CALL HPALLOC(PYY1,2*NY1+1,IER3,0) CALL HPALLOC(PWK,5*NY2+9*NX2,IER4,0) C **** C **** Check for successful allocation C **** IF(IABS(IER1)+IABS(IER2)+IABS(IER3)+IABS(IER4) 1 .NE.0)THEN WRITE(6,*) 1 'TRANSGRD: Error allocating space from heap, IER1 thru 4 = ', 2 IER1,IER2,IER3,IER4 STOP ENDIF C **** C **** Set up extended longitude and latitude arrays, XX1 & YY1 C **** DO I = 1,NX1 XX1(I) = X1(I) ENDDO DO I = 1,2 XX1(NX1+I) = XX1(I) + 360. ENDDO XX1(0) = XX1(NX1) - 360. DO J = 1,NY1 YY1(NY1/2+J) = Y1(J) ENDDO DO J =1,NY1/2 YY1(3*NY1/2+J) = 180. - Y1(NY1+1-J) YY1(J) = -180. + Y1(NY1/2+J) ENDDO C **** C **** Insert values for (180. + DY/2) C **** YY1(2*NY1+1) = 180.+ Y1(NY1/2+1) C **** C **** Now transfer array, F1, to extended array , FF1, taking C **** account of scalar or vector properties. C **** C **** SIGN = 1. for scalar field, -1. for vector field C **** SIGN = 1. IF(ITYPE.EQ.1)SIGN = -1. C **** C **** Copy array F1 to central part of array FF1 C **** DO J = 1,NY1 DO I = 1,NX1 FF1(I,NY1/2+J) = F1(I,J) ENDDO ENDDO C **** C **** Now fill remainder of FF1 C **** DO J = 1,NY1/2 DO I = 1,NX1 FF1(I,NY1/2+1-J) = SIGN*F1(1+MOD(I+NX1/2-1,NX1),J) FF1(I,2*NY1+1-J) = SIGN*F1(1+MOD(I+NX1/2-1,NX1),NY1/2+J) ENDDO ENDDO C **** C **** Insert periodic points C **** DO I = 1,NX1 FF1(I,2*NY1+1) = FF1(I,1) ENDDO DO J = 1,2*NY1+1 FF1(NX1+1,J) = FF1(1,J) FF1(NX1+2,J) = FF1(2,J) FF1(0,J) = FF1(NX1,J) ENDDO C **** C **** Perform bicubic interpolation C **** CALL TLC2(NX1+3,2*NY1+1,XX1(0),YY1,FF1(0,1),NX2,NY2,X2,Y2,F2, 1 INTPOL,5*NY2+9*NX2,WK,LWMIN,IER) C **** C **** Release heap arrays C **** CALL HPDEALLC(PFF1,IER1,0) CALL HPDEALLC(PXX1,IER2,0) CALL HPDEALLC(PYY1,IER3,0) CALL HPDEALLC(PWK,IER4,0) C **** C **** Check for successful deallocation C **** IF(IABS(IER1)+IABS(IER2)+IABS(IER3)+IABS(IER4) 1 .NE.0)THEN WRITE(6,*) 1 'TRANSGRD: Error deallocating heap space, IER1 thru 7 = ', 2 IER1,IER2,IER3,IER4 STOP ENDIF RETURN END c subroutine tlc2(nx,ny,x,y,p,mx,my,xx,yy,q,intpol,lw,w,lwmin,ier) dimension x(nx),y(ny),p(nx,ny),xx(mx),yy(my),q(mx,my),w(lw) dimension intpol(2) c c check input arguments c ier = 1 c c check (xx,yy) grid resolution c if (min0(mx,my) .lt. 1) return c c check intpol c ier = 6 if (intpol(1).ne.1 .and. intpol(1).ne.3) return if (intpol(2).ne.1 .and. intpol(2).ne.3) return c c check (x,y) grid resolution c ier = 2 if (intpol(1).eq.1 .and. nx.lt.2) return if (intpol(1).eq.3 .and. nx.lt.4) return if (intpol(2).eq.1 .and. ny.lt.2) return if (intpol(2).eq.3 .and. ny.lt.4) return c c check work space length input and set minimum c if (intpol(1).eq.1) then lwx = 2*mx else lwx = 5*mx end if if (intpol(2).eq.1) then lwy = 2*(my+mx) else lwy = 5*my+4*mx end if lwmin = lwx+lwy if (lw .lt. lwmin) then ier = 5 return end if c c check (xx,yy) grid contained in (x,y) grid c ier = 3 if (xx(1).lt.x(1) .or. xx(mx).gt.x(nx)) return if (yy(1).lt.y(1) .or. yy(my).gt.y(ny)) return c c check montonicity of grids c ier = 4 do 1 i=2,nx if (x(i-1).ge.x(i)) return 1 continue do 2 j=2,ny if (y(j-1).ge.y(j)) return 2 continue do 4 ii=2,mx if (xx(ii-1).gt.xx(ii)) return 4 continue do 5 jj=2,my if (yy(jj-1).gt.yy(jj)) return 5 continue c c arguments o.k. c ier = 0 if (intpol(2) .eq.1) then c c linearly interpolate in y c j1 = 1 j2 = j1+my j3 = j2 j4 = j3+my j5 = j4 j6 = j5 j7 = j6 j8 = j7+mx j9 = j8+mx c c set y interpolation indices and scales and linearly interpolate c call linmx(ny,y,my,yy,w(j1),w(j3)) i1 = j9 c c set work space portion and indices which depend on x interpolation c if (intpol(1) .eq. 1) then i2 = i1+mx i3 = i2 i4 = i3 i5 = i4 call linmx(nx,x,mx,xx,w(i1),w(i3)) else i2 = i1+mx i3 = i2+mx i4 = i3+mx i5 = i4+mx call cubnmx(nx,x,mx,xx,w(i1),w(i2),w(i3),w(i4),w(i5)) end if call lint2(nx,ny,p,mx,my,q,intpol(1),w(j1),w(j3), + w(j7),w(j8),w(i1),w(i2),w(i3),w(i4),w(i5)) return else c c cubically interpolate in y, set indice pointers c j1 = 1 j2 = j1+my j3 = j2+my j4 = j3+my j5 = j4+my j6 = j5+my j7 = j6+mx j8 = j7+mx j9 = j8+mx call cubnmx(ny,y,my,yy,w(j1),w(j2),w(j3),w(j4),w(j5)) i1 = j9+mx c c set work space portion and indices which depend on x interpolation c if (intpol(1) .eq. 1) then i2 = i1+mx i3 = i2 i4 = i3 i5 = i4 call linmx(nx,x,mx,xx,w(i1),w(i3)) else i2 = i1+mx i3 = i2+mx i4 = i3+mx i5 = i4+mx call cubnmx(nx,x,mx,xx,w(i1),w(i2),w(i3),w(i4),w(i5)) end if call cubt2(nx,ny,p,mx,my,q,intpol,w(j1),w(j2),w(j3), +w(j4),w(j5),w(j6),w(j7),w(j8),w(j9),w(i1),w(i2),w(i3),w(i4),w(i5)) return end if end subroutine lint1(nx,p,mx,q,ix,dx) dimension p(nx),q(mx),ix(mx),dx(mx) c c linearly interpolate p on x onto q on xx c do 1 ii=1,mx i = ix(ii) q(ii) = p(i)+dx(ii)*(p(i+1)-p(i)) 1 continue return end subroutine cubt1(nx,p,mx,q,ix,dxm,dx,dxp,dxpp) dimension p(nx),q(mx) dimension ix(mx),dxm(mx),dx(mx),dxp(mx),dxpp(mx) c c cubically interpolate p on x to q on xx c do 2 ii=1,mx i = ix(ii) q(ii) = dxm(ii)*p(i-1)+dx(ii)*p(i)+dxp(ii)*p(i+1)+dxpp(ii)*p(i+2) 2 continue return end subroutine linmx(nx,x,mx,xx,ix,dx) c c set x grid pointers for xx grid and interpolation scale terms c dimension x(nx),xx(mx) dimension ix(mx),dx(mx) isrt = 1 do 1 ii=1,mx c c find x(i) s.t. x(i) < xx(ii) <= x(i+1) c do 2 i=isrt,nx-1 if (x(i+1) .ge. xx(ii)) then isrt = i ix(ii) = i go to 3 end if 2 continue 3 continue 1 continue c c set linear scale term c do 4 ii=1,mx i = ix(ii) dx(ii) = (xx(ii)-x(i))/(x(i+1)-x(i)) 4 continue return end subroutine cubnmx(nx,x,mx,xx,ix,dxm,dx,dxp,dxpp) dimension x(nx),xx(mx),ix(mx),dxm(mx),dx(mx),dxp(mx),dxpp(mx) isrt = 1 do 1 ii=1,mx c c set i in [2,nx-2] closest s.t. c x(i-1),x(i),x(i+1),x(i+2) can interpolate xx(ii) c do 2 i=isrt,nx-1 if (x(i+1) .ge. xx(ii)) then ix(ii) = min0(nx-2,max0(2,i)) isrt = ix(ii) go to 3 end if 2 continue 3 continue 1 continue c c set cubic scale terms c do 4 ii=1,mx i = ix(ii) dxm(ii) = (xx(ii)-x(i))*(xx(ii)-x(i+1))*(xx(ii)-x(i+2))/ + ((x(i-1)-x(i))*(x(i-1)-x(i+1))*(x(i-1)-x(i+2))) dx(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i+1))*(xx(ii)-x(i+2))/ + ((x(i)-x(i-1))*(x(i)-x(i+1))*(x(i)-x(i+2))) dxp(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i))*(xx(ii)-x(i+2))/ + ((x(i+1)-x(i-1))*(x(i+1)-x(i))*(x(i+1)-x(i+2))) dxpp(ii) = (xx(ii)-x(i-1))*(xx(ii)-x(i))*(xx(ii)-x(i+1))/ + ((x(i+2)-x(i-1))*(x(i+2)-x(i))*(x(i+2)-x(i+1))) 4 continue return end subroutine cubt2(nx,ny,p,mx,my,q,intpol,jy,dym,dy,dyp, +dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp) dimension p(nx,ny) dimension q(mx,my) dimension pjm(mx),pj(mx),pjp(mx),pjpp(mx) dimension jy(my),dym(my),dy(my),dyp(my),dypp(my) dimension ix(mx),dxm(mx),dx(mx),dxp(mx),dxpp(mx) if (intpol.eq.1) then c c linear in x c jsave = -3 do 2 jj=1,my c c load closest four j lines containing interpolate on xx mesh c for j-1,j,j+1,j+2 in pjm,pj,pjp,pjpp c j = jy(jj) if (j.eq.jsave) then c c j pointer has not moved since last pass (no updates or interpolation) c else if (j.eq.jsave+1) then c c update j-1,j,j+1 and interpolate j+2 c do 3 ii=1,mx pjm(ii) = pj(ii) pj(ii) = pjp(ii) pjp(ii) = pjpp(ii) 3 continue call lint1(nx,p(1,j+2),mx,pjpp,ix,dx) else if (j.eq.jsave+2) then c c update j-1,j and interpolate j+1,j+2 c do 4 ii=1,mx pjm(ii) = pjp(ii) pj(ii) = pjpp(ii) 4 continue call lint1(nx,p(1,j+1),mx,pjp,ix,dx) call lint1(nx,p(1,j+2),mx,pjpp,ix,dx) else if (j.eq.jsave+3) then c c update j-1 and interpolate j,j+1,j+2 c do 5 ii=1,mx pjm(ii) = pjpp(ii) 5 continue call lint1(nx,p(1,j),mx,pj,ix,dx) call lint1(nx,p(1,j+1),mx,pjp,ix,dx) call lint1(nx,p(1,j+2),mx,pjpp,ix,dx) else c c interpolate all four j-1,j,j+1,j+2 c call lint1(nx,p(1,j-1),mx,pjm,ix,dx) call lint1(nx,p(1,j),mx,pj,ix,dx) call lint1(nx,p(1,j+1),mx,pjp,ix,dx) call lint1(nx,p(1,j+2),mx,pjpp,ix,dx) end if c c save j pointer for next pass c jsave = j c c cubically interpolate q(ii,jj) from pjm,pj,pjp,pjpp in y direction c do 6 ii=1,mx q(ii,jj) = dym(jj)*pjm(ii)+dy(jj)*pj(ii)+dyp(jj)*pjp(ii)+ + dypp(jj)*pjpp(ii) 6 continue 2 continue return else c c cubic in x c jsave = -3 do 12 jj=1,my c c load closest four j lines containing interpolate on xx mesh c for j-1,j,j+1,j+2 in pjm,pj,pjp,pjpp c j = jy(jj) if (j.eq.jsave) then c c j pointer has not moved since last pass (no updates or interpolation) c else if (j.eq.jsave+1) then c c update j-1,j,j+1 and interpolate j+2 c do 13 ii=1,mx pjm(ii) = pj(ii) pj(ii) = pjp(ii) pjp(ii) = pjpp(ii) 13 continue call cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp) else if (j.eq.jsave+2) then c c update j-1,j and interpolate j+1,j+2 c do 14 ii=1,mx pjm(ii) = pjp(ii) pj(ii) = pjpp(ii) 14 continue call cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp) call cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp) else if (j.eq.jsave+3) then c c update j-1 and interpolate j,j+1,j+2 c do 15 ii=1,mx pjm(ii) = pjpp(ii) 15 continue call cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp) call cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp) call cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp) else c c interpolate all four j-1,j,j+1,j+2 c call cubt1(nx,p(1,j-1),mx,pjm,ix,dxm,dx,dxp,dxpp) call cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp) call cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp) call cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp) end if c c save j pointer for next pass c jsave = j c c cubically interpolate q(ii,jj) from pjm,pj,pjp,pjpp in y direction c do 16 ii=1,mx q(ii,jj) = dym(jj)*pjm(ii)+dy(jj)*pj(ii)+dyp(jj)*pjp(ii)+ + dypp(jj)*pjpp(ii) 16 continue 12 continue return end if end subroutine lint2(nx,ny,p,mx,my,q,intpol,jy,dy,pj,pjp, + ix,dxm,dx,dxp,dxpp) dimension p(nx,ny) dimension q(mx,my) dimension pj(mx),pjp(mx),jy(my),dy(my) dimension ix(mx),dxm(mx),dx(mx),dxp(mx),dxpp(mx) c c linearly interpolate in y c if (intpol.eq.1) then c c linear in x c jsave = -1 do 1 jj=1,my j = jy(jj) if (j.eq.jsave) then c c j pointer has not moved since last pass (no updates or interpolation) c else if (j.eq.jsave+1) then c c update j and interpolate j+1 c do 2 ii=1,mx pj(ii) = pjp(ii) 2 continue call lint1(nx,p(1,j+1),mx,pjp,ix,dx) else c c interpolate j,j+1in pj,pjp on xx mesh c call lint1(nx,p(1,j),mx,pj,ix,dx) call lint1(nx,p(1,j+1),mx,pjp,ix,dx) end if c c save j pointer for next pass c jsave = j c c linearly interpolate q(ii,jj) from pjp,pj in y direction c do 3 ii=1,mx q(ii,jj) = pj(ii)+dy(jj)*(pjp(ii)-pj(ii)) 3 continue 1 continue else c c cubic in x c jsave = -1 do 11 jj=1,my j = jy(jj) if (j.eq.jsave) then c c j pointer has not moved since last pass (no updates or interpolation) c else if (j.eq.jsave+1) then c c update j and interpolate j+1 c do 12 ii=1,mx pj(ii) = pjp(ii) 12 continue call cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp) else c c interpolate j,j+1 in pj,pjp on xx mesh c call cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp) call cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp) end if c c save j pointer for next pass c jsave = j c c linearly interpolate q(ii,jj) from pjp,pj in y direction c do 13 ii=1,mx q(ii,jj) = pj(ii)+dy(jj)*(pjp(ii)-pj(ii)) 13 continue 11 continue return end if end