C**** SUM.FOR Print mean and values of selected boxes 2004/10/18 C**** C**** Compile: FCE SUM.FOR C**** C**** For each record of an IxJ datafile, SUM prints the year and C**** the sum of the values of the specified grid cells. C**** PARAMETER (IMxJM=720*360) CHARACTER FILEIN*80,ZFILE*80,TITLE*80 REAL*4 X(IMxJM) COMMON /FIXDCB/ FOCEAN(IMxJM) C**** NARGS = IARGC() IF(NARGS.lt.3) GO TO 800 CALL GETARG (1,FILEIN) C**** C**** Determine whether data is for ocean only or land only C**** IARG = 1 KOorL = 0 IF(FILEIN(1:2).eq.'-O' .or. FILEIN(1:2).eq.'-o') KOorL = 1 IF(FILEIN(1:2).eq.'-L' .or. FILEIN(1:2).eq.'-l') KOorL = 2 IF(KOorL.gt.0) IARG = 2 IF(KOorL.gt.0) CALL GETARG (2,FILEIN) CALL RESOLU (FILEIN, IM,JM,ZFILE, *801) C**** C**** Read in Ocean Fraction File if requested C**** IF(KOorL.eq.0) GO TO 100 IF(ZFILE(1:4).eq.'None') GO TO 802 CALL GETENV ('IFDIR',TITLE) ZFILE = TRIM(TITLE) // '/' // ZFILE OPEN (11,FILE=ZFILE,FORM='UNFORMATTED',STATUS='OLD',ERR=803) READ (11) TITLE,(FOCEAN(I),I=1,IM*JM) CLOSE (11) C**** C**** Open datafile C**** 100 OPEN (1,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=801) CALL XMAIN (IM,JM,X,FOCEAN,KOorL,IARG) GO TO 999 C**** 800 WRITE (6,*) 'Usage: SUM datafile I1 J1 [I2min-I2max J2 ...]', * ' 2004/10/18' WRITE (6,*) ' or: SUM -O datafile I1 J1 [I2min-I2max J2 ...]', * ' Ocean data only' WRITE (6,*) ' or: SUM -L datafile I1 J1 [I2min-I2max J2 ...]', * ' Land data only' WRITE (6,*) 'Print year, sum, and individual values from each ', * 'record of DataFile' GO TO 999 801 WRITE (0,*) 'Error accessing ',FILEIN STOP 801 802 WRITE (0,*) 'ZFILE was not defined for this resolution: IM,JM=', * IM,JM STOP 802 803 WRITE (0,*) 'Error accessing ZFILE: ',ZFILE WRITE (0,*) 'Make sure that the Input File DIRectory is set.' WRITE (0,*) 'In your .profile add: "export IFDIR=/work/gcm"' STOP 803 999 END SUBROUTINE XMAIN (IM,JM, X,FOCEAN, KOorL,IARG) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (KMAX=12000, IMxJM=720*360) CHARACTER TITLEI*80, TITLE*80,XAXIS*80,YAXIS*80,LABEL*80, * CYEAR*4, MONTH(34)*3, CY(0:5)*10 INTEGER*4 IS(IMxJM),JS(IMxJM) REAL*4 X(IM,JM),FOCEAN(IM,JM) REAL*8 Y(KMAX,-1:5), TMONTH(34) DATA MONTH /'Ann','Jan','Feb','Mar','Apr','May','Jun', * 'Jul','Aug','Sep','Oct','Nov','Dec', * 'ANN','JAN','FEB','MAR','APR','MAY','JUN', * 'JUL','AUG','SEP','OCT','NOV','DEC', * 'DJF','MAM','JJA','SON','JFM','AMJ','JAS','OND'/ DATA TMONTH/.000 ,.042 ,.125 ,.208 ,.292 ,.375 ,.453 , * .542 ,.625 ,.708 ,.792 ,.875 ,.953 , * .000 ,.042 ,.125 ,.208 ,.292 ,.375 ,.453 , * .542 ,.625 ,.708 ,.792 ,.875 ,.953 , * .042 ,.292 ,.542 ,.792 ,.125 ,.375 ,.625 ,.875 / C**** C**** Determine grid points from arguments C**** NARGS = IARGC() NUM = 0 C**** Decode longitudes 100 CALL GETARG (IARG+1,TITLE) MINUS = INDEX (TITLE,'-') IF(MINUS.gt.0) GO TO 110 C**** Single longitude received READ (TITLE,*) IMIN IF(IMIN.lt.1 .or. IMIN.gt.IM) GO TO 820 IMAX = IMIN GO TO 120 C**** Minimum and maximum longitudes received, separated by minus sign 110 READ (TITLE(1:MINUS-1),*) IMIN IF(IMIN.lt.1 .or. IMIN.gt.IM) GO TO 820 READ (TITLE(MINUS+1:80),*) IMAX IF(IMAX.lt.IMIN .or. IMAX.gt.IM) GO TO 820 C**** Decode latitudes 120 CALL GETARG (IARG+2,TITLE) MINUS = INDEX (TITLE,'-') IF(MINUS.gt.0) GO TO 130 C**** Single latitude received READ (TITLE,*) JMIN IF(JMIN.lt.1 .or. JMIN.gt.JM) GO TO 830 JMAX = JMIN GO TO 140 C**** Minimum and maximum latitudes received, separated by minus sign 130 READ (TITLE(1:MINUS-1),*) JMIN IF(JMIN.lt.1 .or. JMIN.gt.JM) GO TO 830 READ (TITLE(MINUS+1:80),*) JMAX IF(JMAX.lt.JMIN .or. JMAX.gt.JM) GO TO 830 C**** 140 DO 150 J=JMIN,JMAX DO 150 I=IMIN,IMAX IF(KOorL.eq.1 .and. FOCEAN(I,J).le..5) GO TO 150 IF(KOorL.eq.2 .and. FOCEAN(I,J).gt..5) GO TO 150 NUM = NUM + 1 IS(NUM) = I JS(NUM) = J 150 continue IARG = IARG + 2 IF(IARG+2.le.NARGS) GO TO 100 NUMX = MIN(NUM,5) IF(NUM.eq.1) NUMX = 0 C**** C**** Loop over records in the file C**** DO 430 K=1,KMAX READ (1,IOSTAT=IOSQ) TITLEI,X IF(IOSQ.ne.0) GO TO 440 IF(K.gt.1) GO TO 300 C**** C**** Calculate 4 header lines of LinePlot file C**** TITLE = TITLEI READ (TITLE(52:54),920,IOSTAT=IOSQ) NUMRUN IF(IOSQ.eq.0) TITLE = TITLEI(1:59) IF(TITLEI(68:68).ne.'X') XAXIS = 'Years' IF(TITLEI(68:68).eq.'X') XAXIS = 'Decades' NPAREN = INDEX (TITLEI,')') IF(NPAREN.le.5 .or. NPAREN.gt.50) NPAREN = 50 YAXIS = TITLEI(1:NPAREN) IF(NUM.eq.1 .and. TITLEI(68:68).ne.'X') * WRITE (LABEL,921) IS(1),JS(1) IF(NUM.eq.1 .and. TITLEI(68:68).eq.'X') * WRITE (LABEL,922) IS(1),JS(1) IF(NUM.gt.1 .and. TITLEI(68:68).ne.'X') * WRITE (LABEL,923) (IS(N),JS(N),N=1,NUMX) IF(NUM.gt.1 .and. TITLEI(68:68).eq.'X') * WRITE (LABEL,924) (IS(N),JS(N),N=1,NUMX) C**** C**** Decode year and month or decade and month, C**** calculate absisca = Y( ,-1) C**** 300 Y(K,-1) = K CYEAR = TITLEI(65:68) IF(CYEAR(4:4).eq.'X') CYEAR = ' ' // TITLEI(65:67) READ (CYEAR,*,IOSTAT=IOSQ) Y(K,-1) DO 310 M=1,34 IF(TITLEI(70:72).ne.MONTH(M)) GO TO 310 Y(K,-1) = Y(K,-1) + TMONTH(M) GO TO 400 310 continue C**** C**** Calculate the sum of the values = Y( ,0) C**** 400 SUM = 0. NV = 0 DO 410 N=1,NUM IF(X(IS(N),JS(N)).le.-999999.) GO TO 410 SUM = SUM + X(IS(N),JS(N)) NV = 1 410 continue Y(K,0) = -999999. IF(NV.gt.0) Y(K,0) = SUM C**** Fill in first five grid point values = Y( ,N) DO 420 N=1,NUMX 420 Y(K,N) = X(IS(N),JS(N)) 430 continue 440 KM = K-1 C**** C**** Write LinePlot output file C**** CALL WRITLP (6, KMAX,KM,1+NUMX, TITLE,XAXIS,YAXIS,LABEL, Y,Y(1,0)) RETURN C**** 820 WRITE (0,*) 'Unable to decipher I argument: ',TITLE(1:10) STOP 820 830 WRITE (0,*) 'Unable to decipher J argument: ',TITLE(1:10) STOP 830 C**** 920 FORMAT (I3) 921 FORMAT ('Year.Mon',I6,',',I3.3) 922 FORMAT (' Dec.Mon',I6,',',I3.3) 923 FORMAT ('Year.Mon Mean',5(I6,',',I3.3)) 924 FORMAT (' Dec.Mon Mean',5(I6,',',I3.3)) END SUBROUTINE RESOLU (FILEIN, IM,JM,ZFILE, *) C**** C**** RESOLU reads in the first four byte integer of a datafile and C**** from this information it determines the horizontal resolution. C**** C**** Input: FILEIN = character file name of a DATAFILE C**** Output: IM,JM = horizontal resolution C**** ZFILE = character file name of the ocean fraction file C**** * = use RETURN 1 if resolution is not determined C**** CHARACTER FILEIN*(*), ZFILE*(*) C**** OPEN (1,FILE=FILEIN,ACCESS='DIRECT',RECL=4,STATUS='OLD',ERR=801) READ (1,REC=1) IJT CLOSE (1) IJRES = (IJT-80)/4 C**** IF(IJRES.eq.46*9) then IM = 46 JM = 9 ZFILE = 'None' RETURN endif C**** IF(IJRES.eq.46*12) then IM = 46 JM = 12 ZFILE = 'None' RETURN endif C**** IF(IJRES.eq.72*46) then IM = 72 JM = 46 ZFILE = 'Z5X4N' RETURN endif C**** IF(IJRES.eq.72*46*2) then IM = 72 JM = 46 ZFILE = 'Z5X4N' RETURN endif C**** IF(IJRES.eq.90*60) then IM = 90 JM = 60 ZFILE = 'Z4X3N' RETURN endif C**** IF(IJRES.eq.144*72) then IM = 144 JM = 72 ZFILE = 'Z144X72' RETURN endif C**** IF(IJRES.eq.144*73) then IM = 144 JM = 73 ZFILE = 'None' RETURN endif C**** IF(IJRES.eq.120*90) then IM = 120 JM = 90 ZFILE = 'Z3X2N' RETURN endif C**** IF(IJRES.eq.144*90) then IM = 144 JM = 90 ZFILE = 'Z144X90N' RETURN endif C**** IF(IJRES.eq.360*180) then IM = 360 JM = 180 ZFILE = 'Z360X180' RETURN endif C**** IF(IJRES.eq.720*360) then IM = 720 JM = 360 ZFILE = 'Z720X360' RETURN endif C**** Unable to derive resolution from record length IM = IJRES JM = 1 ZFILE = 'None' RETURN C**** Unable to open datafile 801 RETURN 1 END SUBROUTINE WRITLP (IUNIT, KMAX,KM,LM, TITLE,XAXIS,YAXIS,LABEL,X,Y) C**** C**** WRITLP writes the data lines of a LinePlot file to unit IUNIT, C**** using a reasonable format. Y values that are -999999. are C**** replaced with * in the output file. C**** C**** Input: KMAX = row dimension of X and Y arrays C**** KM = number of used rows in X and Y arrays C**** LM = number of Y columns C**** TITLE = title character string C**** XAXIS = X axis label character string C**** YAXIS = Y axis label character string C**** LABEL = column labels character string C**** X = one dimensional array of X coordinates C**** Y = two dimensional array of Y values C**** C**** Output: IUNIT = unit to send LinePlot data to, exclude titles C**** CHARACTER TITLE*(*),XAXIS*(*),YAXIS*(*),LABEL*(*), C(0:16)*10 REAL*8 X(KMAX),Y(KMAX,LM) C**** C**** Write 4 header lines of LinePlot file C**** WRITE (IUNIT,900) TRIM(TITLE) WRITE (IUNIT,900) TRIM(XAXIS) WRITE (IUNIT,900) TRIM(YAXIS) WRITE (IUNIT,900) TRIM(LABEL) C**** C**** Calculate maximum absolute Y value C**** ABSMAX = 0. DO 1 L=1,LM DO 1 K=1,KM IF(Y(K,L).eq.-999999.) GO TO 1 IF(ABSMAX.lt.ABS(Y(K,L))) ABSMAX = ABS(Y(K,L)) 1 continue KPOWER = LOG10(ABSMAX) + 2 IF(KPOWER.lt.1) KPOWER = 1 IF(KPOWER.gt.8) KPOWER = 8 C**** C**** Write data lines C**** GO TO (700,600,500,400,300,200,100),KPOWER C**** 1000000 < ABSMAX: output numbers have no decimal places DO 20 K=1,KM DO 10 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),901) NINT(Y(K,L)) 10 IF(Y(K,L).eq.-999999.) C(L) = ' *' 20 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 100000 < ABSMAX < 1000000: output numbers have 1 decimal place 100 DO 120 K=1,KM DO 110 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),911) Y(K,L) 110 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 120 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 10000 < ABSMAX < 100000: output numbers have 2 decimal places 200 DO 220 K=1,KM DO 210 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),921) Y(K,L) 210 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 220 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 1000 < ABSMAX < 10000: output numbers have 3 decimal places 300 DO 320 K=1,KM DO 310 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),931) Y(K,L) 310 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 320 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 100 < ABSMAX < 1000: output numbers have 4 decimal places 400 DO 420 K=1,KM DO 410 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),941) Y(K,L) 410 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 420 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 10 < ABSMAX < 100: output numbers have 5 decimal places 500 DO 520 K=1,KM DO 510 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),951) Y(K,L) 510 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 520 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 1 < ABSMAX < 10: output numbers have 6 decimal places 600 DO 620 K=1,KM DO 610 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),961) Y(K,L) 610 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 620 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** ABSMAX < 1: output numbers have 7 decimal places 700 DO 720 K=1,KM DO 710 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),961) Y(K,L) 710 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 720 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 900 FORMAT (A) 901 FORMAT (1X,I9) 902 FORMAT (F8.3,8A10) 911 FORMAT (1X,F9.1) 921 FORMAT (1X,F9.2) 931 FORMAT (1X,F9.3) 941 FORMAT (1X,F9.4) 951 FORMAT (1X,F9.5) 961 FORMAT (1X,F9.6) 971 FORMAT (1X,F9.7) END