C**** RMS.FOR Global, ocean, continent RMS of IxJ data 2004/06/17 C**** C**** Compile: FCE RMS.FOR C**** C**** RMS reads one or two IxJ datafile records, calculates the global, C**** ocean, continent, ground, land ice, and lake RMS values or C**** differences, 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, DIF*3 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.1) 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 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 NREC1 = 1 IF(NARGS.le.IARG) GO TO 220 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) DIF = ' ' IF(NARGS.le.IARG+1) GO TO 500 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.le.IARG+2) 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) DIF = 'Dif' C**** C**** Subtract record2 from record1 C**** 400 DO 410 I=1,IM*JM IF(X(I).le.DATMIS) GO TO 410 IF(Y(I).le.DATMIS) X(I) = DATMIS IF(Y(I).gt.DATMIS) X(I) = X(I) - Y(I) 410 continue C**** C**** Calculate RMSs C**** 500 CALL RMS (IM,JM, X, FOCEAN,FLAKE,FGRND,FGICE, DIF) GO TO 999 C**** 800 WRITE (6,*) 'Usage: RMS [IMxJM] datafile [recn] ' * // ' RMS value 2004/06/17' WRITE (6,*) ' or: RMS [IMxJM] datafile1 rec1 datafile2 [rec2]' * // ' RMS of rec1 minus rec2' WRITE (6,*) ' RMSs for each hemisphere and surface type a' * // 're sent to standard output' 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 RMS (IM,JM, X, FOCEAN,FLAKE,FGRND,FGICE, DIF) C**** C**** Write RMSs for each hemisphere and each surface type to C**** standard output C**** IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (DATMIS=-999999.) CHARACTER STYPE(6)*7, CM(0:2)*10, DIF*3 REAL*4 X(IM,JM), * FOCEAN(IM,JM),FLAKE(IM,JM),FGRND(IM,JM),FGICE(IM,JM) REAL*8 Q(0:2,6),W(0:2,6),XM(0:2,6) COMMON /GEOMCB/ DXYP(360) DATA STYPE /'Global','Ocean','Contint','Ground','LandIce','Lake'/ C**** C**** Accumulate weighted values and weights C**** DO 130 JHEMI=1,2 C**** Zero out weighted values and weights over a hemisphere DO 110 K=1,6 Q(JHEMI,K) = 0. 110 W(JHEMI,K) = 0. C**** Accumulate weighted values and weights over a hemisphere 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,2) = Q(JHEMI,2) + FOCEAN(I,J)*DXYP(J)*X(I,J)**2 W(JHEMI,2) = W(JHEMI,2) + FOCEAN(I,J)*DXYP(J) Q(JHEMI,4) = Q(JHEMI,4) + FGRND(I,J)*DXYP(J)*X(I,J)**2 W(JHEMI,4) = W(JHEMI,4) + FGRND(I,J)*DXYP(J) Q(JHEMI,5) = Q(JHEMI,5) + FGICE(I,J)*DXYP(J)*X(I,J)**2 W(JHEMI,5) = W(JHEMI,5) + FGICE(I,J)*DXYP(J) Q(JHEMI,6) = Q(JHEMI,6) + FLAKE(I,J)*DXYP(J)*X(I,J)**2 W(JHEMI,6) = W(JHEMI,6) + FLAKE(I,J)*DXYP(J) 120 continue Q(JHEMI,3) = Q(JHEMI,4) + Q(JHEMI,5) + Q(JHEMI,6) W(JHEMI,3) = W(JHEMI,4) + W(JHEMI,5) + W(JHEMI,6) Q(JHEMI,1) = Q(JHEMI,2) + Q(JHEMI,3) W(JHEMI,1) = W(JHEMI,2) + W(JHEMI,3) 130 continue C**** C**** Calculate RMS values C**** DO 210 K=1,6 XM(0,K) = DATMIS XM(1,K) = DATMIS XM(2,K) = DATMIS IF(W(1,K)+W(2,K).gt.0.) * XM(0,K) = SQRT((Q(1,K)+Q(2,K))/(W(1,K)+W(2,K))) IF(W(1,K).gt.0.) XM(1,K) = SQRT(Q(1,K)/W(1,K)) 210 IF(W(2,K).gt.0.) XM(2,K) = SQRT(Q(2,K)/W(2,K)) C**** C**** Calculate maximum absolute value of XM C**** ABSMAX = 0. DO 310 K=1,6 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 RMS 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,6 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,6 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,6 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,6 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,6 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,6 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,6 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,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),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,' RMS ',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 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 = 'Z5X4N' IF(IM.eq. 90 .and. JM.eq. 60) ZFILE = 'Z4X3N' IF(IM.eq.120 .and. JM.eq. 90) ZFILE = 'Z3X2N' 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 = 'Z5X4N' GO TO 200 ; endif C**** IF(IJRES.eq.72*46*2) then 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 = 'Z3X2N' 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 CALL GETENV ('IFDIR',ARG) ZFILE = TRIM(ARG) // '/' // 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