C**** SCC.FOR 2004/10/04 C**** Spatial Correlation Coefficients from IxJ DataFile records C**** C**** Compile in Executable module: FCE SCC.FOR C**** C**** SCC reads two IxJ DataFile records, calculates the global, C**** ocean, continent, ground, land ice, and lake Spatial Correlation C**** Coefficients, and writes the results to standard output. C**** DataFiles must have the format: TITLE(80),REAL4(IM,JM) C**** PARAMETER (DATMIS=-999999.) CHARACTER FILEIN*80, ZFILE*80, TITLEX*80,TITLEY*80 REAL*4 X(720*360),Y(720*360), * FOCEAN(720,360),FLAKE(720,360),FGRND(720,360),FGICE(720,360) C**** NARGS = IARGC() IF(NARGS.lt.3) GO TO 800 C**** C**** Determine resolution C**** IARG = 1 CALL RESOLU (IARG, IM,JM,ZFILE) IF(IM.eq.0) GO TO 810 IF(ZFILE(1:4).eq.'None') GO TO 811 IF(NARGS.lt.IARG+2) GO TO 800 C WRITE (6,*) 'Derived resolution: IM,JM=',IM,JM CALL INPUT (IM,JM, ZFILE, FOCEAN,FLAKE,FGRND,FGICE) C**** C**** Open DataFile1 C**** CALL GETARG (IARG,FILEIN) OPEN (1,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=820) C**** Determine record number for DataFile1 CALL GETARG (IARG+1,FILEIN) READ (FILEIN,*,ERR=821) NREC1 DO 210 N=1,NREC1-1 210 READ (1,ERR=822,END=823) 220 READ (1,ERR=822,END=823) TITLEX,(X(I),I=1,IM*JM) CLOSE (1) WRITE (6,924) TRIM(TITLEX) C**** C**** Open DataFile2 C**** CALL GETARG (IARG+2,FILEIN) OPEN (2,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=820) C**** Determine record number for DataFile2 NREC2 = 1 IF(NARGS.lt.IARG+3) GO TO 320 CALL GETARG (IARG+3,FILEIN) READ (FILEIN,*,ERR=821) NREC2 DO 310 N=1,NREC2-1 310 READ (2,ERR=822,END=823) 320 READ (2,ERR=822,END=823) TITLEY,(Y(I),I=1,IM*JM) CLOSE (2) WRITE (6,924) TRIM(TITLEY) C**** C**** Calculate Spatial Correlation Coefficients C**** 500 CALL SCC (IM,JM, X,Y, FOCEAN,FLAKE,FGRND,FGICE) GO TO 999 C**** 800 WRITE (6,*) 'Usage: SCC [IMxJM] datafile1 rec1 datafile2 [rec2] ' * // ' 2004/10/04' WRITE (6,*) 'Spatial Correlation Coefficients for each hemisphere' * // ' and each surface type' GO TO 999 810 CALL GETARG (1,FILEIN) WRITE (0,*) 'Error opening DataFile: ',TRIM(FILEIN) STOP 810 811 WRITE (0,*) *'Ocean fraction file not available for this resolution:',IM,JM STOP 811 820 WRITE (0,*) 'Error opening DataFile: ',TRIM(FILEIN) STOP 820 821 WRITE (0,*) 'Unable to decipher record number: ',TRIM(FILEIN) STOP 821 822 WRITE (0,922) N,TRIM(FILEIN) STOP 822 823 WRITE (0,923) N,TRIM(FILEIN) STOP 823 C**** 922 FORMAT ('Error reading record',I5,' of DataFile: ',A) 923 FORMAT ('End-of-File reading record',I5,' of DataFile: ',A) 924 FORMAT (A) 999 END SUBROUTINE SCC (IM,JM, X,Y, FOCEN,FLAKE,FGRND,FGICE) C**** C**** Write Spatial Correlation Coefficients for each hemisphere and C**** each surface type to standard output C**** IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (DATMIS=-999999.) CHARACTER STYPE(6)*7, CM(0:2)*10 REAL*4 X(IM,JM),Y(IM,JM), * FOCEN(IM,JM),FLAKE(IM,JM),FGRND(IM,JM),FGICE(IM,JM) REAL*8 S(0:2,6,6), XM(0:2,6) COMMON /GEOMCB/ DXYP(360) COMMON /WORKCB/ SW(0:2,6), SXW(0:2,6), SYW(0:2,6), * SXXW(0:2,6),SYYW(0:2,6),SXYW(0:2,6) EQUIVALENCE (S,SW) DATA STYPE /'Global','Ocean ','Contint','Ground','LandIce','Lake'/ C**** C**** Accumulate weighted values C**** DO 130 JHEMI=1,2 C**** Zero out weighted values over a hemisphere DO 110 K=1,6 SW (JHEMI,K) = 0. SXW (JHEMI,K) = 0. SYW (JHEMI,K) = 0. SXXW(JHEMI,K) = 0. SXYW(JHEMI,K) = 0. 110 SYYW(JHEMI,K) = 0. C**** Accumulate weighted values over Ocean, Lakes, Ground and Land Ice DO 120 JH=1,JM/2 J=JH+(2-JHEMI)*JM/2 DO 120 I=1,IM IF(X(I,J).le.DATMIS .or. Y(I,J).le.DATMIS) GO TO 120 IF(FOCEN(I,J).gt.0.) then SW (JHEMI,2) = SW (JHEMI,2) + FOCEN(I,J)*DXYP(J) SXW (JHEMI,2) = SXW (JHEMI,2) + FOCEN(I,J)*DXYP(J)*X(I,J) SYW (JHEMI,2) = SYW (JHEMI,2) + FOCEN(I,J)*DXYP(J)*Y(I,J) SXXW(JHEMI,2) = SXXW(JHEMI,2) + FOCEN(I,J)*DXYP(J)*X(I,J)**2 SYYW(JHEMI,2) = SYYW(JHEMI,2) + FOCEN(I,J)*DXYP(J)*Y(I,J)**2 SXYW(JHEMI,2) = SXYW(JHEMI,2) + FOCEN(I,J)*DXYP(J)*X(I,J)*Y(I,J) endif IF(FGRND(I,J).gt.0.) then SW (JHEMI,4) = SW (JHEMI,4) + FGRND(I,J)*DXYP(J) SXW (JHEMI,4) = SXW (JHEMI,4) + FGRND(I,J)*DXYP(J)*X(I,J) SYW (JHEMI,4) = SYW (JHEMI,4) + FGRND(I,J)*DXYP(J)*Y(I,J) SXXW(JHEMI,4) = SXXW(JHEMI,4) + FGRND(I,J)*DXYP(J)*X(I,J)**2 SYYW(JHEMI,4) = SYYW(JHEMI,4) + FGRND(I,J)*DXYP(J)*Y(I,J)**2 SXYW(JHEMI,4) = SXYW(JHEMI,4) + FGRND(I,J)*DXYP(J)*X(I,J)*Y(I,J) endif IF(FGICE(I,J).gt.0.) then SW (JHEMI,5) = SW (JHEMI,5) + FGICE(I,J)*DXYP(J) SXW (JHEMI,5) = SXW (JHEMI,5) + FGICE(I,J)*DXYP(J)*X(I,J) SYW (JHEMI,5) = SYW (JHEMI,5) + FGICE(I,J)*DXYP(J)*Y(I,J) SXXW(JHEMI,5) = SXXW(JHEMI,5) + FGICE(I,J)*DXYP(J)*X(I,J)**2 SYYW(JHEMI,5) = SYYW(JHEMI,5) + FGICE(I,J)*DXYP(J)*Y(I,J)**2 SXYW(JHEMI,5) = SXYW(JHEMI,5) + FGICE(I,J)*DXYP(J)*X(I,J)*Y(I,J) endif IF(FLAKE(I,J).gt.0.) then SW (JHEMI,6) = SW (JHEMI,6) + FLAKE(I,J)*DXYP(J) SXW (JHEMI,6) = SXW (JHEMI,6) + FLAKE(I,J)*DXYP(J)*X(I,J) SYW (JHEMI,6) = SYW (JHEMI,6) + FLAKE(I,J)*DXYP(J)*Y(I,J) SXXW(JHEMI,6) = SXXW(JHEMI,6) + FLAKE(I,J)*DXYP(J)*X(I,J)**2 SYYW(JHEMI,6) = SYYW(JHEMI,6) + FLAKE(I,J)*DXYP(J)*Y(I,J)**2 SXYW(JHEMI,6) = SXYW(JHEMI,6) + FLAKE(I,J)*DXYP(J)*X(I,J)*Y(I,J) endif 120 continue C**** Calculate weighted values over Global and Continents DO 130 L=1,6 S(JHEMI,3,L) = S(JHEMI,4,L) + S(JHEMI,5,L) + S(JHEMI,6,L) 130 S(JHEMI,1,L) = S(JHEMI,2,L) + S(JHEMI,3,L) C**** Calculate global weighted values from hemispheric values DO 140 K=1,6*6 140 S(0,K,1) = S(1,K,1) + S(2,K,1) C**** C**** Calculate Spatial Correlation Coefficients C**** DO 210 K=1,6 DO 210 N=0,2 XM(N,K) = DATMIS 210 IF(SXXW(N,K).gt.0. .and. SYYW(N,K).gt.0.) * XM(N,K) = (SXYW(N,K)-SXW(N,K)*SYW(N,K)/SW(N,K)) / / SQRT ((SXXW(N,K)-SXW(N,K)**2/SW(N,K))* * (SYYW(N,K)-SYW(N,K)**2/SW(N,K))) C**** C**** Write Spatial Correlation Coefficients to standard output C**** DO 472 K=1,6 DO 471 N=0,2 IF(XM(N,K).gt.DATMIS) WRITE (CM(N),947) XM(N,K) 471 IF(XM(N,K).le.DATMIS) CM(N) = ' * ' 472 WRITE (6,948) STYPE(K),CM RETURN C**** 947 FORMAT (1X,F9.7) 948 FORMAT (A8,' SpaCorCoef Both:',A10,' NH:',A10,' SH:',A10) END SUBROUTINE RESOLU (IARG, IM,JM,ZFILE) C**** C**** RESOLU obtains the IARG-th argument from the command line. C**** If that argument is a resolution character string (e.g. 72x46), C**** then IM and JM are filled. C**** If that argument is a datafile, then the first four byte integer C**** is read from the file, and IM and JM are determined from it. C**** The ZFILE is filled for standard resolutions. C**** C**** Input: IARG = location among the command line aruments C**** Output: IARG = incremented if ARG is a resolution character string C**** IM,JM = horizontal resolution C**** ZFILE = name of the ocean fraction file C**** CHARACTER ARG*80, ZFILE*(*) LOGICAL*4 QEXIST ZFILE = 'None' CALL GETARG (IARG,ARG) C**** C**** Is IARG-th command line argument a resolution character string? C**** LOCAX = SCAN(ARG,'Xx') IF(LOCAX.le.1) GO TO 100 READ (ARG(1:LOCAX-1) ,*,IOSTAT=IOSQ) IM IF(IOSQ.ne.0) GO TO 100 READ (ARG(LOCAX+1:80),*,IOSTAT=IOSQ) JM IF(IOSQ.ne.0) GO TO 100 IARG = IARG + 1 C**** Fill ocean fraction file IF(IM.eq. 72 .and. JM.eq. 46) ZFILE = 'Z72X46N' IF(IM.eq. 90 .and. JM.eq. 60) ZFILE = 'Z4X3N' IF(IM.eq.120 .and. JM.eq. 90) ZFILE = 'Z120X90N' IF(IM.eq.144 .and. JM.eq. 72) ZFILE = 'Z144X72' IF(IM.eq.144 .and. JM.eq. 90) ZFILE = 'Z144X90N' IF(IM.eq.360 .and. JM.eq.180) ZFILE = 'Z360X180' IF(IM.eq.720 .and. JM.eq.360) ZFILE = 'Z720X360' IF(ZFILE.eq.'None') RETURN GO TO 200 C**** C**** Derive resolution from DataFile C**** 100 OPEN (1,FILE=ARG,ACCESS='DIRECT',RECL=4,STATUS='OLD',ERR=800) READ (1,REC=1) IJT CLOSE (1) IJRES = (IJT-80)/4 C**** IF(IJRES.eq.46*9) then IM = 46 JM = 9 RETURN endif C**** IF(IJRES.eq.46*12) then IM = 46 JM = 12 RETURN endif C**** IF(IJRES.eq.72*46) then IM = 72 JM = 46 ZFILE = 'Z72X46N' GO TO 200 endif C**** IF(IJRES.eq.72*46*2) then IM = 72 JM = 46 ZFILE = 'Z72X46N' GO TO 200 endif C**** IF(IJRES.eq.90*60) then IM = 90 JM = 60 ZFILE = 'Z4X3N' GO TO 200 endif C**** IF(IJRES.eq.120*90) then IM = 120 JM = 90 ZFILE = 'Z120X90N' GO TO 200 endif C**** IF(IJRES.eq.144*72) then IM = 144 JM = 72 ZFILE = 'Z144X72' GO TO 200 endif C**** IF(IJRES.eq.144*73) then IM = 144 JM = 73 RETURN endif C**** IF(IJRES.eq.144*90) then IM = 144 JM = 90 ZFILE = 'Z144X90N' GO TO 200 endif C**** IF(IJRES.eq.360*180) then IM = 360 JM = 180 ZFILE = 'Z360X180' GO TO 200 endif C**** IF(IJRES.eq.720*360) then IM = 720 JM = 360 ZFILE = 'Z720X360' GO TO 200 endif C**** Unable to derive resolution from record length IM = IJRES JM = 1 RETURN C**** C**** Determine path to ZFILE C**** 200 INQUIRE (FILE=ZFILE,EXIST=QEXIST) IF(QEXIST) RETURN CALL GETENV ('IFDIR',ARG) ARG = TRIM(ARG) // '/' // ZFILE INQUIRE (FILE=ARG,EXIST=QEXIST) IF(.not.QEXIST) GO TO 210 ZFILE = ARG RETURN 210 ZFILE = '/u/cmrun/' // ZFILE RETURN C**** Unable to open datafile 800 IM = 0 GO TO 200 END SUBROUTINE INPUT (IM,JM, ZFILE, FOCEAN,FLAKE,FGRND,FGICE) IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (TWOPI=6.283185307179586477d0) CHARACTER ZFILE*80, TITLE*80 REAL*4 FOCEAN(IM,JM),FLAKE(IM,JM),FGRND(IM,JM),FGICE(IM,JM) COMMON /GEOMCB/ DXYP(360) C**** C**** Calculate the spherical geometry for the C grid C**** DLON = TWOPI/IM DLAT = TWOPI*NINT(360./(JM-1))/720. FJEQ = (1+JM)/2. C**** Geometric parameters centered at primary latitudes VSINS = -1. DO 10 J=1,JM VSINN = SIN(DLAT*(J+.5-FJEQ)) IF(J.eq.JM) VSINN = 1. DXYP(J) = DLON*(VSINN-VSINS) 10 VSINS = VSINN C**** C**** Read in FOCEAN, FLAKE, FGRND and FGICE from ZFILE C**** OPEN (11,FILE=ZFILE,FORM='UNFORMATTED',STATUS='OLD',ERR=810) READ (11) TITLE,FOCEAN READ (11) TITLE,FLAKE READ (11) TITLE,FGRND READ (11) TITLE,FGICE CLOSE (11) RETURN C**** 810 WRITE (0,*) 'Error opening: ',TRIM(ZFILE) STOP 810 END