c subroutine ccm2gcm(imxtgcm,jmxtgcm,imxccm,jmxccm,nflds, + xccm,yccm,xtgcm,ytgcm,fccm,ftgcm,ivec) c c Transform from ccm2 to tgcm grid. c If ccm2 has more lat and lon points than tgcm, call fine2coarse c (e.g., when ccm2 is running in T42 resolution) c If ccm2 has fewer longitudes but more latitudes than tgcm c (e.g., ccm2 is running as maccm2), call fine2coarse, ignoring c the minor difference of only 8 longitudes. c If ccm2 has fewer lat and lon points than tgcm, call coarse2fine c (e.g., when ccm2 is running in T21 resolution) 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) real ccmlon(imxccm+1),ccmlat(jmxccm+2) ! for coarse2fine c c ccm -> gcm is fine -> coarse: c SUBROUTINE FINE2COURSE(IMXF,JMXF,IMXC,JMXC,XF,YF,XC,YC,FF,FC, c 1 NFLDS) c if ((imxccm.gt.imxtgcm.and.jmxccm.gt.jmxtgcm).or. ! T42 + (imxccm.lt.imxtgcm.and.jmxccm.gt.jmxtgcm)) then ! maccm2 call fine2coarse(imxccm,jmxccm,imxtgcm,jmxtgcm, + xccm,yccm,xtgcm,ytgcm,fccm,ftgcm,nflds) c c ccm -> gcm is coarse -> fine: c SUBROUTINE COARSE2FINE(FC,ICMX,JCMX,FF,IFMX,JFMX,XC,YC,XF,YF, c 1 ITYPE) c (note ccmlon(icmx+1) and ccmlat(jcmx+2) for wrap around lon and poles): c elseif (imxccm.lt.imxtgcm.and.jmxccm.lt.jmxtgcm) then ! T21 do i=1,imxccm ccmlon(i) = xccm(i) enddo ccmlon(imxccm+1) = 180. ccmlat(1) = -90. do j=1,jmxccm ccmlat(j+1) = yccm(j) enddo ccmlat(jmxccm+2) = 90. do ip=1,nflds call coarse2fine(fccm(1,1,ip),imxccm,jmxccm,ftgcm(1,1,ip), + imxtgcm,jmxtgcm,ccmlon,ccmlat,xtgcm,ytgcm,ivec(ip)) enddo else write(6,"('>>> WARNING ccm2gcm: do not know whether to ', + 'call fine2coarse or coarse2fine:')") write(6,"(' imxtgcm,jmxtgcm=',2i4)") imxtgcm,jmxtgcm write(6,"(' imxccm,jmxccm=',2i4)") imxccm,jmxccm stop 'ccm2gcm' endif 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 c gcm -> ccm is coarse -> fine c SUBROUTINE COARSE2FINE(FC,ICMX,JCMX,FF,IFMX,JFMX,XC,YC,XF,YF, c 1 ITYPE) c Must put pole points latitudes and wrap-around longitude in tgcm c when it is the coarse model: c if ((imxccm.gt.imxtgcm.and.jmxccm.gt.jmxtgcm).or. ! t42 + (imxccm.lt.imxtgcm.and.jmxccm.gt.jmxtgcm)) then ! maccm2 do i=1,imxtgcm gcmlon(i) = xtgcm(i) enddo gcmlon(imxtgcm+1) = 180. gcmlat(1) = -90. do j=1,jmxtgcm gcmlat(j+1) = ytgcm(j) enddo gcmlat(jmxtgcm+2) = 90. do ip=1,nflds call coarse2fine(ftgcm(1,1,ip),imxtgcm,jmxtgcm,fccm(1,1,ip), + imxccm,jmxccm,gcmlon,gcmlat,xccm,yccm,ivec(ip)) enddo c c gcm -> ccm is fine -> coarse c SUBROUTINE FINE2COARSE(IMXF,JMXF,IMXC,JMXC,XF,YF,XC,YC,FF,FC, c 1 NFLDS) c elseif (imxccm.lt.imxtgcm.and.jmxccm.lt.jmxtgcm) then ! t21 call fine2coarse(imxtgcm,jmxtgcm,imxccm,jmxccm,xtgcm,ytgcm, + xccm,yccm,ftgcm,fccm,nflds) else write(6,"('>>> WARNING gcm2ccm: do not know whether to ', + 'call fine2coarse or coarse2fine:')") write(6,"(' imxtgcm,jmxtgcm=',2i4)") imxtgcm,jmxtgcm write(6,"(' imxccm,jmxccm=',2i4)") imxccm,jmxccm stop 'gcm2ccm' endif return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE FINE2COARSE(IMXF,JMXF,IMXC,JMXC,XF,YF,XC,YC,FF,FC, 1 NFLDS) C **** C **** Performs transformation from a fine global grid to a C **** coarse one. C **** C **** Input: C **** C **** IMXF = no. of independent longitude values for fine grid C **** C **** JMXF = no. of latitude values for fine grid C **** C **** IMXC = no. of independent longitude values for coarse C **** grid C **** JMXC = no. of latitude values for coarse grid C **** C **** XF(IMXF) = longitude values for fine grid C **** C **** YF(JMXF) = latitude values for fine grid C **** C **** XC(IMXC) = longitude values for coarse grid C **** C **** YC(JMXC) = latitude values for coarse grid C **** C **** FF(IMXF,JMXF,NFLDS) = NFLDS fine grid arrays to be C **** transformed to coarse grid C **** Output: C **** C **** FC(IMXC,JMXC,NFLDS) = NFLDS transformed arrays on C **** coarse grid C **** C **** On first entry initialise interpolation process. C **** SAVE DATA IFIRST/1/ DIMENSION XF(IMXF), YF(JMXF), XC(IMXC), YC(JMXC), FF(IMXF,JMXF,1), 1 FC(IMXC,JMXC,1) C **** C **** At first entry, set up work arrays using heap allocation C **** C **** Define pointers: C **** REAL MSKF, MSKC POINTER (PXVF,XVF(IMXF+1)), (PYVF,YVF(IMXF+1)), 1 (PXVC,XVC(IMXF+1)), (PYVC,YVC(IMXF+1)), 2 (PMSKF,MSKF(IMXF,JMXF)), (PMSKC,MSKC(IMXC,JMXC)), 3 (PAREA,AREA(IMXC,JMXC)), 4 (PIBEG,IBEG(0:IMXC+1)), (PIEND,IEND(0:IMXC+1)), 5 (PJBEG,JBEG(0:JMXC+1)), (PJEND,JEND(0:JMXC+1)), 6 (PTMP,TMP(IMXF,JMXF)), (PWORK,WORK(IMXF,JMXF)) C **** C **** Where: C **** C **** XVF(IMXF+1) = longitudinal vertex grid points for fine C **** grid C **** YVF(JMXF+1) = latitudinal vertex grid points for fine C **** grid C **** XVC(IMXC+1) = longitudinal vertex grid points for coarse C **** grid C **** YVC(JMXC+1) = latitudinal vertex grid points for coarse C **** grid C **** MSKF(IMXF,JMXF) = mask for fines grid ( = 1.) C **** C **** MSKC(IMXC,JMXC) = same after transformation C **** ( should be 1.) C **** AREA(IMXC,JMXC) = area of each coarse grid cell allowing C **** for spherical geometry C **** C **** IBEG(0:IMXC+1), IEND(0:IMXC+1), C **** JBEG(0:JMXC+1), JEND(0:JMXC+1) C **** C **** These arrays are used to mark the beginning and end of C **** the sequence of fine grid vertex elements contributing C **** to each, larger, coarse grid element. The two C **** dimensions are trated separately. C **** C **** TMP(IMXF,JMXF) and WORK(IMXF,JMXF) are work arrays C **** C **** Logical variables: C **** LOGICAL SPHERE, PERIODX, PERIODY C **** C **** Where: C **** C **** SPHERE = .TRUE. for spherical geometry C **** C **** PERIODX = .TRUE. indicates functions are periodic C **** in the X-direction C **** C **** PERIODY = .FALSE. indicates functions are not C **** periodic in the Y-direction C **** DATA SPHERE, PERIODX, PERIODY / 2*.TRUE.,.FALSE./ C **** C **** Set MSKF to unity. Note that MSKC will be calculated C **** and should turn out to be unity also. C **** IF ( IFIRST.EQ.1) THEN IFIRST = 0 C **** C **** Initialisation: C **** C **** Allocate work space from heap: C **** CALL HPALLOC(PXVF,IMXF+1,IER1,0) CALL HPALLOC(PYVF,JMXF+1,IER2,0) CALL HPALLOC(PXVC,IMXC+1,IER3,0) CALL HPALLOC(PYVC,JMXC+1,IER4,0) CALL HPALLOC(PMSKF,IMXF*JMXF,IER5,0) CALL HPALLOC(PMSKC,IMXC*JMXC,IER6,0) CALL HPALLOC(PAREA,IMXC*JMXC,IER7,0) CALL HPALLOC(PIBEG,IMXC+2,IER8,0) CALL HPALLOC(PIEND,IMXC+2,IER9,0) CALL HPALLOC(PJBEG,JMXC+2,IER10,0) CALL HPALLOC(PJEND,JMXC+2,IER11,0) CALL HPALLOC(PTMP,IMXF*JMXF,IER12,0) CALL HPALLOC(PWORK,IMXF*JMXF,IER13,0) C **** C **** Check for successful allocation C **** IF(IABS(IER1)+IABS(IER2)+IABS(IER3)+IABS(IER4)+IABS(IER5)+ 1 IABS(IER6)+IABS(IER7)+IABS(IER8)+IABS(IER9)+IABS(IER10)+ 2 IABS(IER11)+IABS(IER12)+IABS(IER13).NE.0)THEN WRITE(6,*) 1 'FINE2COARSE: Error allocating space from heap, ', + 'IER1 - 13 = ', + IER1,IER2,IER3,IER4,IER5,IER6,IER7,IER8,IER9,IER10,IER11, + IER12,IER13 write(6,"('imxf=',i3,' jmxf=',i3,' imxc=',i3,' jmxc=',i3)") + imxf,jmxf,imxc,jmxc STOP ENDIF DO I = 1,IMXF*JMXF MSKF(I,1) = 1. ENDDO C **** C **** Now generate vertex grids falling half way between C **** central grid points. C **** DO I = 2,IMXC XVC(I) = .5*(XC(I-1) + XC(I)) ENDDO XVC(1) = .5*(XC(IMXC) - 360. + XC(1)) XVC(IMXC+1) = .5*(XC(IMXC) + XC(1) + 360.) DO J = 2,JMXC YVC(J) = .5*(YC(J-1) + YC(J)) ENDDO YVC(1) = -90. YVC(JMXC+1) = 90. DO I = 2,IMXF XVF(I) = .5*(XF(I-1) + XF(I)) ENDDO XVF(1) = .5*(XF(IMXF) - 360. + XF(1)) XVF(IMXF+1) = .5*(XF(IMXF) + XF(1) + 360.) DO J = 2,JMXF YVF(J) = .5*(YF(J-1) + YF(J)) ENDDO YVF(1) = -90. YVF(JMXF+1) = 90. C **** C **** Initialize arrays IBEG, IEND, JBEG, JEND C **** CALL SETTBL( 1 IMXC, IMXC, JMXC, IMXF, IMXF, JMXF, 2 XVC, YVC, XVF, YVF, PERIODX, PERIODY, 3 IBEG, IEND, JBEG, JEND) C **** C **** Transform MSKF to coarse coordinate system in MSKC. C **** Integrate the mask over each large cell in AREA. C **** CALL GEOGWT( 1 IMXC, IMXC, JMXC, IMXF, IMXF, JMXF, 2 XVC, YVC, XVF, YVF, MSKF, 3 IBEG, IEND, JBEG, JEND, 4 SPHERE, TMP, WORK, AREA, MSKC) ENDIF C **** C **** Transform field FF(IMXF,JMXF) on fine grid to FC(IMXC,JMXC) C **** on coarse grid. C **** CALL AVG2GRD( 1 IMXC, IMXC, JMXC, IMXF, IMXF, JMXF, 2 NFLDS, IBEG, IEND, JBEG, JEND, 3 XVC, YVC, XVF, YVF, MSKF, AREA, 4 SPHERE, TMP, WORK, FF, FC) RETURN END c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE COARSE2FINE(FC,ICMX,JCMX,FF,IFMX,JFMX,XC,YC,XF,YF, 1 ITYPE) C **** C **** This subroutine transforms a coarse grid array to a C **** finer grid using bilinear interpolation. C **** C **** Input: C **** C **** FC(ICMX,JCMX) = coarse grid array to be transformed C **** C **** ICMX, JCMX = dimensions of coarse grid array C **** C **** XC(ICMX+1) = array of coarse grid longitude values. C **** (XC(ICMX+1) = XC(1) for wrap around) C **** C **** YC(JCMX+2) = array of coarse grid latitude values. C **** (YC(1) and YC(JCMX+2) are polar values) C **** C **** IFMX, JFMX are dimensions of fine grid, FF, to which C **** FC is to be transformed C **** C **** XF(IFMX) = array of fine grid longitude values. C **** C **** YF(JFMX) = array of fine grid latitude values. C **** C **** ITYPE = 0 for scalar FC C **** = 1 for vector FC C **** C **** Output: C **** C **** FF(IFMX,JFMX) = fine grid version of array FC C **** DIMENSION FC(ICMX,JCMX), FF(IFMX,JFMX), XC(ICMX+1), YC(JCMX+2), 1 XF(IFMX), YF(JFMX) C **** C **** Set up pointers to work space needed for interpolation C **** procedure. C **** POINTER (PICM,ICM(IFMX)), (PJCM,JCM(JFMX)), 1 (PDICM,DICM(IFMX)), (PDJCM,DJCM(JFMX)), 2 (PWORK,WORK(ICMX+1,0:JCMX+1)) C **** C **** Where: C **** C **** ICM(IFMX), JCM(JFMX) will contain indices needed for C **** bilinear interpolation generated by TABULAT. C **** (ICM(I),JCM(J)) = indices of lower left hand corner of C **** coarse grid element containing fine grid point (I,J). C **** C **** DICM(IFMX), DJCN(JFMX) = weights needed for the bilinear C **** interpolation. These are generated by TABULAT C **** C **** WORK(ICMX+1,0:JCMX+1) = work space for FC after addition C **** of periodic and polar values. C **** C **** On first entry, generate arrays ICM, JCM, DICM, DJCM C **** needed for bilinear interpolation. Also set up work C **** array, WORK. C **** SAVE DATA IFIRST/1/ IF (IFIRST.EQ.1) THEN IFIRST = 0 C **** C **** Allocate needed space from heap C **** CALL HPALLOC(PICM,IFMX,IER1,0) CALL HPALLOC(PDICM,IFMX,IER2,0) CALL HPALLOC(PJCM,JFMX,IER3,0) CALL HPALLOC(PDJCM,JFMX,IER4,0) CALL HPALLOC(PWORK,(ICMX+1)*(JCMX+2),IER5,0) C **** C **** Check for successful allocation C **** IF(IABS(IER1)+IABS(IER2)+IABS(IER3)+IABS(IER4)+IABS(IER5).NE.0) 1 THEN WRITE(6,*) 1 'COARSE2FINE: Error allocating space from heap, IER1 - 5 = ', + IER1,IER2,IER3,IER4,IER5 write(6,"('imxf=',i3,' jmxf=',i3,' imxc=',i3,' jmxc=',i3)") + imxf,jmxf,imxc,jmxc STOP ENDIF C **** C **** Call TABULAT to tabulate indices and weights for C **** interpolation in x-direction. C **** CALL TABULAT(ICMX+1,XC,IFMX,XF,ICM,DICM) C **** C **** Now call TABULAT to tabulate indices and weights for C **** interpolation in y-direction C **** CALL TABULAT(JCMX+2,YC,JFMX,YF,JCM,DJCM) C **** C **** Subtract 1 from each value of JC to allow for initial C **** index of YC being 0 C **** DO J = 1,JFMX JCM(J) = JCM(J)-1 ENDDO C **** C **** Initialisation is now complete C **** ENDIF C **** C **** We are now ready to transform FC to FF C **** C **** Copy FC to WORK and insert polar and periodic values C **** DO J = 1,JCMX DO I = 1,ICMX WORK(I,J) = FC(I,J) ENDDO ENDDO C **** C **** Insert polar values C **** IF (ITYPE.EQ.0) THEN C **** C **** Field is scalar C **** C **** S and N poles C **** WORK(1,0) = SSUM(ICMX,FC(1,1),1)/FLOAT(ICMX) WORK(1,JCMX+1) = SSUM(ICMX,FC(1,JCMX),1)/FLOAT(ICMX) DO I = 2,ICMX WORK(I,0) = WORK(1,0) WORK(I,JCMX+1) = WORK(1,JCMX+1) ENDDO ELSE IF (ITYPE.EQ.1) THEN C **** C **** Field is vector (changes sign across pole) C **** C **** S and N poles C **** DO I = 1,ICMX WORK(I,0) = .5*(WORK(I,1) - WORK(1+MOD(I-1+ICMX/2,ICMX),1)) WORK(I,JCMX+1) = .5*(WORK(I,JCMX) - 1 WORK(1+MOD(I-1+ICMX/2,ICMX),JCMX)) ENDDO ELSE WRITE(6,*)'ITYPE must be 0 or 1' CALL EXIT ENDIF C **** C **** Wrap around point C **** DO J = 0,JCMX+1 WORK(ICMX+1,J) = WORK(1,J) ENDDO C **** C **** Now do interpolation C **** DO J = 1,JFMX DO I = 1,IFMX FF(I,J) = (1.-DICM(I))*(1.-DJCM(J))*WORK(ICM(I),JCM(J)) 1 +DICM(I)*(1.-DJCM(J))*WORK(ICM(I)+1,JCM(J)) 2 +DICM(I)*DJCM(J)*WORK(ICM(I)+1,JCM(J)+1) 3 +(1.-DICM(I))*DJCM(J)*WORK(ICM(I),JCM(J)+1) ENDDO ENDDO RETURN END C SUBROUTINE TABULAT(ICMX,XC,IFMX,XF,ILO,FRAC) C **** C **** Calculates indices and weights for transformation from C **** coarse to fine grid using linear interpolation. C **** C **** Input: C **** C **** ICMX = number of coarse grid points C **** C **** IFMX = number of fine grid points C **** C **** XC(ICMX) = values of X at each of ICMX coarse grid C **** points C **** C **** XF(IFMX) = values of X at each of IFMX fine grid C **** points C **** C **** Output: C **** C **** ILO(IFMX): C **** ILO(I) is defined by: C **** C **** XC(ILO(I)) .LE. XF(I) .LT. XC(ILO(I)+1) C **** C **** and gives the lower bound of the coarse index interval C **** within which XF(I) falls C **** C **** FRAC(I) = (XF(I)-XC(ILO(I)))/(XC(ILO(I)+1)-XC(ILO(I))) C **** DIMENSION XC(ICMX),XF(IFMX),ILO(IFMX),FRAC(IFMX) C **** C **** Get indices of coarse grid points surrounding each fine C **** grid point C **** DO I = 1,IFMX ILO(I) = ISRCHFGT(ICMX,XC,1,XF(I)) - 1 c write(6,"('tabulat: ifmx=',i2,' i=',i2,' icmx=',i2,' xf(i)=', c + f16.10,' xc(1)=',f16.10,' ilo(i)=',i2)") c + ifmx,i,icmx,xf(i),xc(1),ilo(i) C **** C **** Now calculate FRAC C **** FRAC(I) = (XF(I) - XC(ILO(I)))/(XC(ILO(I)+1) - XC(ILO(I))) ENDDO RETURN END subroutine avg2grd(icdim,icmax,jcmax, $ ifdim,ifmax,jfmax,nflds, $ ibeg,iend,jbeg,jend, $ xc,yc,xf,yf,mskfn,area,sphere, $ tmp,work,fldsf,fldsc) implicit none !---- Arguments ------------------------------------------------------- c Input integer ifdim, ! fine grid dimension $ ifmax,jfmax, ! fine grid extent to be averaged $ icdim, ! coarse grid dimension $ icmax,jcmax, ! coarse grid size $ nflds, ! number of fields $ ibeg(0:icmax+1), ! start - stop indices on fine grid $ iend(0:icmax+1), $ jbeg(0:jcmax+1), $ jend(0:jcmax+1) real fldsf(ifdim,jfmax,nflds), ! fine grid fields to average $ mskfn(ifdim,jfmax), ! fine grid surface type concentration $ xf(ifmax+1), ! fine vertex grid x-coordinates $ yf(jfmax+1), ! fine vertex grid y-coordinates $ xc(icmax+1), ! coarse grid x-coordinates $ yc(jcmax+1), ! coarse grid y-coordinates $ area(icdim,jcmax) ! area of coarse cell spanned by vertex points real $ tmp(ifdim,jfmax), $ work(ifdim,jfmax) ! work arrays for averaging logical sphere ! flag to use spherical coordinates ! (cosine factor taken into account) ! NB grids should be in radians on input c Output real fldsc(icdim,jcmax,nflds) ! coarse grid average field !---- Local Variables ------------------------------------------------- integer nf !---------------------------------------------------------------------- do nf=1,nflds call average(icdim,icmax,jcmax,ifdim,ifmax,jfmax, $ xc,yc,xf,yf, $ ibeg,iend,jbeg,jend,mskfn, $ area,sphere,tmp,work, $ fldsf(1,1,nf),fldsc(1,1,nf)) enddo return end c c======================================================================= c This module contains code required for averaging from a finer c grid to a coarser grid c c Two kinds of weights are required. The first kind, which we may c call "grid weights," are functions only of the two grids. For each c small cell, the weights are computed as follows: c weight = f*da/dA c where c f is the fraction of the small cell contained within the current large cell c da is the area of the small cell c dA is the area of the large cell c A normalization sum is computed for each larger cell c c Weights are computed based on the existence of "vertex" grids for each c mesh. These vertex grids are defined from the "central" grids such that c they are equidistant between two consecutive central points, in each c direction. If a vertex grid does not already exist, it must be constructed. c Areas are computed based on the vertex grids, but centered on the central c grids. c c The second kind of weight, the "geographical" weight, depends on the c type of surface represented by a given cell; for this code that may c be ocean, land, or sea ice. c c the values on the coarser grid. c This method is constructed such that the global average of the function c on the finer "ocean-only" grid is equal within roundoff to the global average c on the coarser "ocean-only" grid. c======================================================================= subroutine settbl (icd,icmax,jcmax,ifd,ifmax,jfmax, $ xc,yc,xf,yf,wrapx, wrapy, $ ibeg, iend,jbeg, jend ) c--------------------------------------------------------------------------- c Input c icd dimension of coarse grid c icmax, jcmax extent of coarse grid to be examined (may be subset) c xc, yc coarser grid, vertex values c ifd dimension of fine grid c ifmax, jfmax extent of fine grid to be examined (may be subset) c xf, yf finer grid, vertex values c wrapx periodic in longitude c wrapy periodic in latitude c Output c ibeg, iend c jbeg, jend c--------------------------------------------------------------------------- implicit none integer $ icd,icmax,jcmax,ifd,ifmax,jfmax real $ xc(icmax+1), $ yc(jcmax+1), $ xf(ifmax+1), $ yf(jfmax+1) logical $ wrapx, $ wrapy integer $ ibeg(0:icmax+1), $ iend(0:icmax+1), $ jbeg(0:jcmax+1), $ jend(0:jcmax+1), $ i,j integer $ isrchfge, $ isrchfgt external $ isrchfge, $ isrchfgt C***************************************** C DO I = 1,JCMAX+1 C WRITE(6,*)I,YF(I),YC(I) C ENDDO C DO I = JCMAX+2,JFMAX+1 C WRITE(6,*)I,YF(I) C ENDDO C***************************************** ibeg(0)=1 iend(0)=0 ibeg(icmax+1)=1 iend(icmax+1)=0 do i=1,icmax ibeg(i)= isrchfgt(ifmax+1,xf,1,xc(i)) - 1 enddo do i=1,icmax-1 iend(i) = ibeg(i+1) enddo iend(icmax)= isrchfge(ifmax+1,xf,1,xc(icmax+1)) - 1 if (wrapx) then if (xc(1).lt.xf(1)) then ibeg(0) = iend(icmax) iend(0) = ifmax ibeg(1) = 1 else if (xc(icmax+1).gt.xf(ifmax+1)) then iend(icmax) = ifmax ibeg(icmax+1) = 1 iend(icmax+1) = ibeg(1) endif endif jbeg(0)=1 jend(0)=0 jbeg(jcmax+1)=1 jend(jcmax+1)=0 do j=1,jcmax jbeg(j)= isrchfgt(jfmax+1,yf,1,yc(j)) - 1 c******************************* c write(6,*) jbeg(j) c******************************* enddo do j=1,jcmax-1 jend(j) = jbeg(j+1) enddo jend(jcmax)= isrchfge(jfmax+1,yf,1,yc(jcmax+1)) - 1 if (wrapy) then if (yc(1).lt.yf(1)) then jbeg(0) = jend(jcmax) jend(0) = jfmax jbeg(1) = 1 else if (yc(jcmax+1).gt.yf(jfmax+1)) then jend(jcmax) = jfmax jbeg(jcmax+1) = 1 jend(jcmax+1) = jbeg(1) endif endif return end c======================================================================= subroutine average(icd,icmax,jcmax,ifd,ifmax,jfmax, $ xc,yc,xf,yf, $ ibeg,iend,jbeg,jend,mskfn, $ wt,sphere,tmp,work,fldf,fldc) implicit none !---- Arguments ------------------------------------------------------- c Input integer ifd, ! fine grid dimension $ ifmax,jfmax, ! fine grid extent to be averaged $ icd, ! coarse grid dimension $ icmax,jcmax, ! coarse grid size $ ibeg(0:icmax+1), ! start - stop indices on fine grid $ iend(0:icmax+1), $ jbeg(0:jcmax+1), $ jend(0:jcmax+1) real fldf(ifd,jfmax), ! fine grid field to average $ mskfn(ifd,jfmax), ! fine grid surface type concentration $ xf(ifmax+1), ! fine vertex grid x-coordinates $ yf(jfmax+1), ! fine vertex grid y-coordinates $ xc(icmax+1), ! coarse grid x-coordinates $ yc(jcmax+1), ! coarse grid y-coordinates $ wt(icd,jcmax) ! area of coarse cell spanned by vertex points logical sphere ! flag to use spherical coordinates c Output real fldc(icd,jcmax) ! coarse grid average field !---- Local Variables ------------------------------------------------- integer if, jf, ic, jc real tmp(ifd,jfmax) real work(ifd,jfmax) !---------------------------------------------------------------------- !---- Mask fine grid field by surface percentage do jf=1,jfmax do if=1,ifmax tmp(if,jf) = fldf(if,jf)*mskfn(if,jf) enddo enddo !---- Integrate over each coarse grid cell call intgr2d(ifd,ifmax,jfmax,tmp,xf,yf, $ icd,icmax,jcmax,ibeg,iend,jbeg,jend,xc,yc, $ sphere,work,fldc) !---- Normalize by integral area occupied by this surface type do jc=1,jcmax do ic=1,icmax if ( wt(ic,jc) .ne. 0. ) then fldc(ic,jc) = fldc(ic,jc)/wt(ic,jc) else fldc(ic,jc)=0. endif enddo enddo return end c======================================================================= c subroutine geogwt(icd,icmax,jcmax,ifd,ifmax,jfmax, $ xc,yc,xf,yf,mskfn, $ ibeg,iend,jbeg,jend, $ sphere,tmp,work,area,wgt) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Note that when sea ice is included, the mask and hence the c geographical weights changes with time c c The ocean mask used should agree with CCM conventions c (boundaries are not tagged separately) c c index 1 = ocean (mask=0) c index 2 = land (mask=1) c index 3 = ice (mask=2) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none !---- Arguments c Input integer ifd, ! dimension of fine grid $ ifmax, jfmax, $ icd, ! dimension of coarse grid $ icmax, jcmax, $ ibeg(0:icmax+1), iend(0:icmax+1), $ jbeg(0:jcmax+1), jend(0:jcmax+1) real $ xf(ifmax+1), $ yf(jfmax+1), $ xc(icmax+1), $ yc(jcmax+1), $ mskfn(ifd,jfmax) logical sphere c Output real $ area(icd,jcmax), $ wgt(icd,jcmax) !---- Local variables integer if, jf, ic, jc real tmp(ifd,jfmax) real work(ifd,jfmax) !---------------------------------------------------------------------- !---- Compute area of each coarse cell covered by this surface type do jf=1,jfmax do if=1,ifmax tmp(if,jf) = 1. enddo enddo !---- determine area of each coarse grid cell (temp. store in f) call intgr2d(ifd,ifmax,jfmax,tmp,xf,yf, $ icd,icmax,jcmax,ibeg,iend,jbeg,jend,xc,yc, $ sphere,work,wgt) !---- Integrate concentration each coarse grid cell call intgr2d(ifd,ifmax,jfmax,mskfn,xf,yf, $ icd,icmax,jcmax,ibeg,iend,jbeg,jend,xc,yc, $ sphere,work,area) do jc=1,jcmax do ic=1,icmax if ( wgt(ic,jc) .ne. 0. ) then wgt(ic,jc) = area(ic,jc)/wgt(ic,jc) endif enddo enddo return end c======================================================================= c subroutine intgr2d(ifd,ifmax,jfmax,af,xf,yf, $ icd,icmax,jcmax, $ ibeg,iend,jbeg,jend,xc,yc, $ sphere,wrk,ac) implicit none !---- Arguments ------------------------------------------------------- integer ifd, ! fine grid size $ ifmax,jfmax ! fine grid size real af(ifd,jfmax), ! fine grid field to average $ xf(ifmax+1), ! fine grid x-coordinates $ yf(jfmax+1) ! fine grid y-coordinates integer icd, ! coarse grid dimension $ icmax,jcmax, ! coarse grid size $ ibeg(0:icmax+1), ! start - stop indices on fine grid $ iend(0:icmax+1), $ jbeg(0:jcmax+1), $ jend(0:jcmax+1) real xc(icmax+1), ! coarse grid x-coordinates $ yc(jcmax+1) ! coarse grid y-coordinates real wrk(ifd,jcmax) logical sphere ! flag to use spherical coordinates real ac(icd,jcmax) ! coarse grid average field !---- Local Variables ------------------------------------------------- integer i, j, if, jf, ic, jc real dx, dy real xclft, xcrgt, yclft, ycrgt real radian data radian/0.017453292/ real del !---------------------------------------------------------------------- !---- Initialize accumulators do jc=1,jcmax do ic=1,icmax ac(ic,jc) = 0. enddo do if=1,ifmax wrk(if,jc) = 0. enddo enddo !---- Integrate in Y-direction do j=0,jcmax+1 if ( j .eq. 0 ) then jc = 1 ! map into first coarse cell yclft = yc(jcmax+1) ycrgt = yc(jcmax+1) + (yc(2) - yc(1)) elseif ( j .eq. jcmax+1 ) then jc = jcmax ! map into last cell yclft = yc(1) - (yc(jcmax+1) - yc(jcmax)) ycrgt = yc(1) else jc = j yclft = yc(jc) ycrgt = yc(jc+1) endif do jf=jbeg(j),jend(j) if ( sphere ) then dy = sin(radian*min(ycrgt,yf(jf+1))) $ - sin(radian*max(yclft,yf(jf))) else dy = min(ycrgt,yf(jf+1)) - max(yclft,yf(jf)) endif do if=1,ifmax wrk(if,jc) = wrk(if,jc) + af(if,jf)*dy enddo enddo enddo !---- Integrate in X-direction do i=0,icmax if ( i .eq. 0 ) then ic = 1 del = 360. xclft = xc(icmax+1) xcrgt = xc(icmax+1) + (xc(2) - xc(1)) else ic = i del = 0. xclft = xc(ic) xcrgt = xc(ic+1) endif do if=ibeg(i),iend(i) if(sphere) then dx = radian*min(xc(ic+1) + del,xf(if+1)) $ - radian*max(xc(ic) + del,xf(if)) else dx = min(xc(ic+1) + del,xf(if+1)) $ - max(xc(ic) + del,xf(if)) endif do jc=1,jcmax ac(ic,jc) = ac(ic,jc) + wrk(if,jc)*dx enddo enddo enddo return end