C**** MEAN.F0R Global, ocean, continent MEAN of IxJ data 2004/06/18 C**** C**** Compile: FCE MEAN.FOR C**** C**** MEAN reads one or two IxJ DataFile records, calculates the global, C**** ocean, continent, ground, land ice, and lake MEAN vlaues or C**** differences, and writes the results to standard output. C**** DataFiles must have the format: TITLE(80),REAL4(IM,JM) C**** IMPLICIT NONE CHARACTER FILEIN*80, ZFILE*80, TITLEX*80,TITLEY*80, DIF*3 INTEGER*4 IARGC,NARGS,IARG, IM,JM, NREC1,NREC2, N,I REAL*4 X(720*360),Y(720*360), FOCEN,FGRND,FGICE,FLAKE, DATMIS PARAMETER (DATMIS=-999999.) COMMON /FIXDCB/ FOCEN(720*360),FGRND(720*360), * FGICE(720*360),FLAKE(720*360) C**** NARGS = IARGC() IF(NARGS.lt.1) GO TO 800 C**** C**** Determine resolution C**** IARG = 1 CALL RESOLU (IARG, IM,JM,ZFILE) CALL INPUT (IM,JM, ZFILE, X) C**** C**** Open DataFile1 C**** CALL GETARG (IARG,FILEIN) OPEN (1,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) C**** Determine record number for DataFile1 NREC1 = 1 IF(NARGS.lt.IARG+1) GO TO 120 CALL GETARG (IARG+1,FILEIN) READ (FILEIN,*,ERR=811) NREC1 DO 110 N=1,NREC1-1 110 READ (1,ERR=812,END=813) 120 READ (1,ERR=812,END=813) TITLEX,(X(I),I=1,IM*JM) CLOSE (1) WRITE (6,914) TRIM(TITLEX) DIF = ' ' IF(NARGS.lt.IARG+2) GO TO 400 C**** C**** Open DataFile2 C**** CALL GETARG (IARG+2,FILEIN) OPEN (2,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) C**** Determine record number for DataFile2 NREC2 = 1 IF(NARGS.lt.IARG+3) GO TO 220 CALL GETARG (IARG+3,FILEIN) READ (FILEIN,*,ERR=811) NREC2 DO 210 N=1,NREC2-1 210 READ (2,ERR=812,END=813) 220 READ (2,ERR=812,END=813) TITLEY,(Y(I),I=1,IM*JM) CLOSE (2) DIF = 'Dif' WRITE (6,914) TRIM(TITLEY) C**** C**** Subtract record2 from record1 C**** DO 310 I=1,IM*JM IF(Y(I).le.DATMIS) X(I) = DATMIS 310 IF(X(I).gt.DATMIS) X(I) = X(I) - Y(I) C**** C**** Calculate MEANs C**** 400 CALL MEAN (IM,JM, X, FOCEN,FLAKE,FGRND,FGICE, DIF) GO TO 999 C**** 800 WRITE (6,*) 'Usage: MEAN [IMxJM] DataFile [recn] ' * // ' MEAN value 2004/06/18' WRITE (6,*) ' or: MEAN [IMxJM] DataFile1 rec1 DataFile2 [rec2]' * // ' MEAN of rec1 minus rec2' WRITE (6,*) ' MEANs for each hemisphere and surface type a' * // 're sent to standard output.' GO TO 999 810 WRITE (0,*) 'Error opening DataFile: ',FILEIN STOP 810 811 WRITE (0,*) 'Unable to decipher record number: ',FILEIN STOP 811 812 WRITE (0,912) N,TRIM(FILEIN) STOP 812 813 WRITE (0,913) N,TRIM(FILEIN) STOP 813 C**** 912 FORMAT ('Error reading record',I5,' of DataFile: ',A) 913 FORMAT ('End-of-File reading record',I5,' of DataFile: ',A) 914 FORMAT (A) 999 END SUBROUTINE MEAN (IM,JM, X, FOCEN,FLAKE,FGRND,FGICE, DIF) C**** C**** Write MEANs for each hemisphere and each surface type to C**** standard output C**** IMPLICIT NONE CHARACTER STYPE(6)*7, CM(0:2)*10, DIF*3 INTEGER*4 IM,JM, KMAX, I,J,K,N,JH,JHEMI,KPOWER REAL*4 X(IM,JM), * FOCEN(IM,JM),FLAKE(IM,JM),FGRND(IM,JM),FGICE(IM,JM) REAL*8 Q(2),W(2), XM(0:2,6), DATMIS, DXYP, ABSMAX PARAMETER (DATMIS=-999999.) COMMON /GEOMCB/ DXYP(360) DATA STYPE /'Global','Ocean','Contint','Ground','LandIce','Lake'/ C**** C**** Calculate Global mean values C**** KMAX = 1 DO 110 JHEMI=1,2 Q(JHEMI) = 0. W(JHEMI) = 0. DO 110 JH=1,JM/2 J=JH+(2-JHEMI)*JM/2 DO 110 I=1,IM IF(X(I,J).le.DATMIS) GO TO 110 Q(JHEMI) = Q(JHEMI) + DXYP(J)*X(I,J) W(JHEMI) = W(JHEMI) + DXYP(J) 110 continue XM(0,1) = DATMIS XM(1,1) = DATMIS XM(2,1) = DATMIS IF(W(1)+W(2).gt.0.) XM(0,1) = (Q(1)+Q(2))/(W(1)+W(2)) IF(W(1) .gt.0.) XM(1,1) = Q(1) / W(1) IF(W(2) .gt.0.) XM(2,1) = Q(2) / W(2) C**** C**** Calculate Ocean mean values C**** IF(FOCEN(1,1).eq.DATMIS) GO TO 300 KMAX = 3 DO 120 JHEMI=1,2 Q(JHEMI) = 0. W(JHEMI) = 0. DO 120 JH=1,JM/2 J=JH+(2-JHEMI)*JM/2 DO 120 I=1,IM IF(X(I,J).le.DATMIS) GO TO 120 Q(JHEMI) = Q(JHEMI) + FOCEN(I,J)*DXYP(J)*X(I,J) W(JHEMI) = W(JHEMI) + FOCEN(I,J)*DXYP(J) 120 continue XM(0,2) = DATMIS XM(1,2) = DATMIS XM(2,2) = DATMIS IF(W(1)+W(2).gt.0.) XM(0,2) = (Q(1)+Q(2))/(W(1)+W(2)) IF(W(1) .gt.0.) XM(1,2) = Q(1) / W(1) IF(W(2) .gt.0.) XM(2,2) = Q(2) / W(2) C**** C**** Calculate Continental mean values C**** DO 130 JHEMI=1,2 Q(JHEMI) = 0. W(JHEMI) = 0. DO 130 JH=1,JM/2 J=JH+(2-JHEMI)*JM/2 DO 130 I=1,IM IF(X(I,J).le.DATMIS) GO TO 130 Q(JHEMI) = Q(JHEMI) + (1.-FOCEN(I,J))*DXYP(J)*X(I,J) W(JHEMI) = W(JHEMI) + (1.-FOCEN(I,J))*DXYP(J) 130 continue XM(0,3) = DATMIS XM(1,3) = DATMIS XM(2,3) = DATMIS IF(W(1)+W(2).gt.0.) XM(0,3) = (Q(1)+Q(2))/(W(1)+W(2)) IF(W(1) .gt.0.) XM(1,3) = Q(1) / W(1) IF(W(2) .gt.0.) XM(2,3) = Q(2) / W(2) C**** C**** Calculate Ground mean values C**** IF(FGRND(1,1).eq.DATMIS) GO TO 300 KMAX = 4 DO 140 JHEMI=1,2 Q(JHEMI) = 0. W(JHEMI) = 0. DO 140 JH=1,JM/2 J=JH+(2-JHEMI)*JM/2 DO 140 I=1,IM IF(X(I,J).le.DATMIS) GO TO 140 Q(JHEMI) = Q(JHEMI) + FGRND(I,J)*DXYP(J)*X(I,J) W(JHEMI) = W(JHEMI) + FGRND(I,J)*DXYP(J) 140 continue XM(0,4) = DATMIS XM(1,4) = DATMIS XM(2,4) = DATMIS IF(W(1)+W(2).gt.0.) XM(0,4) = (Q(1)+Q(2))/(W(1)+W(2)) IF(W(1) .gt.0.) XM(1,4) = Q(1) / W(1) IF(W(2) .gt.0.) XM(2,4) = Q(2) / W(2) C**** C**** Calculate Glacial Ice mean values C**** IF(FGICE(1,1).eq.DATMIS) GO TO 300 KMAX = 5 DO 150 JHEMI=1,2 Q(JHEMI) = 0. W(JHEMI) = 0. DO 150 JH=1,JM/2 J=JH+(2-JHEMI)*JM/2 DO 150 I=1,IM IF(X(I,J).le.DATMIS) GO TO 150 Q(JHEMI) = Q(JHEMI) + FGICE(I,J)*DXYP(J)*X(I,J) W(JHEMI) = W(JHEMI) + FGICE(I,J)*DXYP(J) 150 continue XM(0,5) = DATMIS XM(1,5) = DATMIS XM(2,5) = DATMIS IF(W(1)+W(2).gt.0.) XM(0,5) = (Q(1)+Q(2))/(W(1)+W(2)) IF(W(1) .gt.0.) XM(1,5) = Q(1) / W(1) IF(W(2) .gt.0.) XM(2,5) = Q(2) / W(2) C**** C**** Calculate Lake mean values C**** IF(FLAKE(1,1).eq.DATMIS) GO TO 300 KMAX = 6 DO 160 JHEMI=1,2 Q(JHEMI) = 0. W(JHEMI) = 0. DO 160 JH=1,JM/2 J=JH+(2-JHEMI)*JM/2 DO 160 I=1,IM IF(X(I,J).le.DATMIS) GO TO 160 Q(JHEMI) = Q(JHEMI) + FLAKE(I,J)*DXYP(J)*X(I,J) W(JHEMI) = W(JHEMI) + FLAKE(I,J)*DXYP(J) 160 continue XM(0,6) = DATMIS XM(1,6) = DATMIS XM(2,6) = DATMIS IF(W(1)+W(2).gt.0.) XM(0,6) = (Q(1)+Q(2))/(W(1)+W(2)) IF(W(1) .gt.0.) XM(1,6) = Q(1) / W(1) IF(W(2) .gt.0.) XM(2,6) = Q(2) / W(2) C**** C**** Calculate maximum absolute value of XM C**** 300 ABSMAX = 0. DO 310 K=1,KMAX DO 310 N=0,2 IF(XM(N,K).le.DATMIS) GO TO 310 IF(ABS(XM(N,K)).gt.ABSMAX) ABSMAX = ABS(XM(N,K)) 310 continue IF(ABSMAX < 1) GO TO 470 IF(ABSMAX >= 1000000) GO TO 400 KPOWER = LOG10(ABSMAX) + 2 C**** C**** Write MEAN values to standard output C**** GO TO (470,460,450,440,430,420,410),KPOWER C**** 1000000 < ABSMAX: output numbers have no decimal places 400 DO 402 K=1,KMAX DO 401 N=0,2 IF(XM(N,K).gt.DATMIS) WRITE (CM(N),940) NINT(XM(N,K)) 401 IF(XM(N,K).le.DATMIS) CM(N) = ' *' 402 WRITE (6,948) STYPE(K),DIF,CM RETURN C**** 100000 < ABSMAX < 1000000: output numbers have 1 decimal place 410 DO 412 K=1,KMAX DO 411 N=0,2 IF(XM(N,K).gt.DATMIS) WRITE (CM(N),941) XM(N,K) 411 IF(XM(N,K).le.DATMIS) CM(N) = ' * ' 412 WRITE (6,948) STYPE(K),DIF,CM RETURN C**** 10000 < ABSMAX < 100000: output numbers have 2 decimal places 420 DO 422 K=1,KMAX DO 421 N=0,2 IF(XM(N,K).gt.DATMIS) WRITE (CM(N),942) XM(N,K) 421 IF(XM(N,K).le.DATMIS) CM(N) = ' * ' 422 WRITE (6,948) STYPE(K),DIF,CM RETURN C**** 1000 < ABSMAX < 10000: output numbers have 3 decimal places 430 DO 432 K=1,KMAX DO 431 N=0,2 IF(XM(N,K).gt.DATMIS) WRITE (CM(N),943) XM(N,K) 431 IF(XM(N,K).le.DATMIS) CM(N) = ' * ' 432 WRITE (6,948) STYPE(K),DIF,CM RETURN C**** 100 < ABSMAX < 1000: output numbers have 4 decimal places 440 DO 442 K=1,KMAX DO 441 N=0,2 IF(XM(N,K).gt.DATMIS) WRITE (CM(N),944) XM(N,K) 441 IF(XM(N,K).le.DATMIS) CM(N) = ' * ' 442 WRITE (6,948) STYPE(K),DIF,CM RETURN C**** 10 < ABSMAX < 100: output numbers have 5 decimal places 450 DO 452 K=1,KMAX DO 451 N=0,2 IF(XM(N,K).gt.DATMIS) WRITE (CM(N),945) XM(N,K) 451 IF(XM(N,K).le.DATMIS) CM(N) = ' * ' 452 WRITE (6,948) STYPE(K),DIF,CM RETURN C**** 1 < ABSMAX < 10: output numbers have 6 decimal places 460 DO 462 K=1,KMAX DO 461 N=0,2 IF(XM(N,K).gt.DATMIS) WRITE (CM(N),946) XM(N,K) 461 IF(XM(N,K).le.DATMIS) CM(N) = ' * ' 462 WRITE (6,948) STYPE(K),DIF,CM RETURN C**** ABSMAX < 1: output numbers have 7 decimal places 470 DO 472 K=1,KMAX 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),DIF,CM RETURN C**** 940 FORMAT (1X,I9) 941 FORMAT (1X,F9.1) 942 FORMAT (1X,F9.2) 943 FORMAT (1X,F9.3) 944 FORMAT (1X,F9.4) 945 FORMAT (1X,F9.5) 946 FORMAT (1X,F9.6) 947 FORMAT (1X,F9.7) 948 FORMAT (A8,' Mean ',A3,' 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 DataFile, 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 DataFile C**** IMPLICIT NONE CHARACTER ARG*80, ZFILE*(*) LOGICAL*4 QEXIST INTEGER*4 IARG, IM,JM, LOCAX,IOSQ,IJT,IJRES ZFILE = 'None' CALL GETARG (IARG,ARG) C**** C**** Is IARG-th command line argument a resolution character string? C**** IF(ARG(1:1).lt.'0' .or. ARG(1:1).gt.'9') GO TO 100 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 DataFile IF(IM.eq. 72 .and. JM.eq. 46) ZFILE = 'Z5X4N' 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=820) 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 .or. IJRES.eq.72*46*2 .or. * IJRES.eq.72*46+46+1) then ! I,J + LATave + GLOBave IM = 72 JM = 46 ZFILE = 'Z5X4N' 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 ('HOME',ARG) ARG = TRIM(ARG) // '/' // ZFILE INQUIRE (FILE=ARG,EXIST=QEXIST) IF(QEXIST) then ; ZFILE = ARG ; RETURN ; endif CALL GETENV ('IFDIR',ARG) ARG = TRIM(ARG) // '/' // ZFILE INQUIRE (FILE=ARG,EXIST=QEXIST) IF(QEXIST) then ; ZFILE = ARG ; RETURN ; endif ARG = '/u/cmrun/' // ZFILE INQUIRE (FILE=ARG,EXIST=QEXIST) IF(QEXIST) then ; ZFILE = ARG ; RETURN ; endif C**** Unable to locate Z DataFile WRITE (0,*) 'Unable to locate Z DataFile.' WRITE (0,*) 'Tried: ',TRIM(ZFILE) WRITE (0,*) 'Tried: $HOME/',TRIM(ZFILE) WRITE (0,*) 'Tried: $IFDIR/',TRIM(ZFILE) WRITE (0,*) 'Tried: /u/cmrun/',TRIM(ZFILE) STOP 821 C**** Unable to open DataFile 820 WRITE (0,*) 'Error opening DataFile: ',ARG STOP 820 END SUBROUTINE INPUT (IM,JM, ZFILE, X) IMPLICIT NONE CHARACTER ZFILE*80, TITLE*80 INTEGER*4 IM,JM, I,J,KF,KQ REAL*4 X(IM,JM), FOCEN REAL*8 TWOPI,DATMIS, DXYP, DLON,DLAT,FJEQ,VSINS,VSINN PARAMETER (TWOPI=6.283185307179586477d0, DATMIS=-999999.) COMMON /FIXDCB/ FOCEN(720*360,4) 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 FOCEN, FGRND, FGICE and FLAKE from ZFILE C**** FOCEN(1,1) = DATMIS FOCEN(1,2) = DATMIS FOCEN(1,3) = DATMIS FOCEN(1,4) = DATMIS IF(ZFILE(1:4).eq.'None') RETURN OPEN (11,FILE=ZFILE,FORM='UNFORMATTED',STATUS='OLD',ERR=830) DO 110 KQ=1,4 READ (11,END=120) TITLE,X KF=0 IF(TITLE(1:5).eq.'Ocean'.or.TITLE(1:5).eq.'FOCEA') KF=1 IF(TITLE(1:5).eq.'Groun'.or.TITLE(1:5).eq.'FGRND') KF=2 IF(TITLE(1:5).eq.'Glaci'.or.TITLE(1:5).eq.'FGICE') KF=3 IF(TITLE(1:5).eq.'Lake '.or.TITLE(1:5).eq.'FLAKE') KF=4 IF(KF.eq.0) GO TO 120 DO 110 I=1,IM*JM 110 FOCEN(I,KF) = X(I,1) 120 CLOSE (11) RETURN C**** 830 WRITE (0,*) 'Error opening Z DataFile: ',TRIM(ZFILE) STOP 830 END