C**** CONS070.S print GCM CONServation and budget pages 1999/07/02 C**** C**** Compile: FCE CONS070.S $RUN/C071D.o C**** C**** CONS070 reads several Climate Model III Diagnostic Files, C**** averages the data, calls GCM subroutines DIAGCP and DIAGJ, and C**** may create a line plot file of the conservation diagnostics. C**** C**** The instantaneous conservation quantities are wrong unless a C**** single input D file is used. C**** INCLUDE '/u/cmrun/C070.COM' PARAMETER (KDINT = 15 + 1 + 4 + 2*4 + 42*50, * KDACC = JM*KAJ*6 + JM*LMA*KAJL + JM*3*4 + IM*JM*KAIJ + * IM*JM*LMA*KAIJL + IM*JM*LMO*KOIJL + * LMO*NMST*KOLNST + 24*50*4 + JM*KCON) INTEGER*4 IPARM(300),IDIAG(KDINT) REAL*4 DIAGR4(KDACC) EQUIVALENCE (IPARM,IM$) C**** CHARACTER FILEIN*80, TITLE*80, OUTYR*4,OUTMON*4, LPLAB(16)*11 INTEGER*4 INSTAN(21), LPINDX(16) COMMON /WORK01/ SCALE(KCON), FLAT(-2:JM,KCON+5) DATA INSTAN /1,10,20,26,32,38, 43,47, 51,55, 59,65, * 71,79,87,95,104, 109,117, 125,131/ C**** NARGS = IARGC() IF(NARGS.le.0) GO TO 800 OPEN (6,FILE='CONS070.PRT') C**** Read first input file to determine KAIJL0, KCON0, KDACCX CALL GETARG (1,FILEIN) OPEN (1,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) READ (1) IHOURX,LABEL,IPARM CLOSE (1) KAIJL0 = JM*KAJ$*6 + JM*LMA$*KAJL$ + JM*3*4 + IM*JM*KAIJ$ KCON0 = KAIJL0 + IM*JM*LMA$*KAIJL$ + IM*JM*LMO$*KOIJL$ + + LMO$*NMST$*KOLNS$ + 24*50*4 KDACCX = KCON0 + JM*KCON$ WRITE (6,*) 'KDINT,KDACCX =',KDINT,KDACCX C**** C**** Read in fixed arrays C**** CALL SGEOM FILEIN = '/u/cmrun/Z72X46N' OPEN (11,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) CALL READR4 (11,IM*JM,FOCEAN,FOCEAN) CALL READR4 (11,IM*JM,FLAKE, FLAKE) CALL READR4 (11,IM*JM,FGRND, FGRND) CALL READR4 (11,IM*JM,FGICE, FGICE) CLOSE (11) C**** C**** Zero out summing arrays C**** DO 10 N=1,15 10 IDACC(N) = 0 DO 20 J=1,JM*KAJ$*6 20 AJ(J,1) = 0. DO 30 I=1,IM*JM*LMA$ 30 AIJL(I,1,1,9) = 0. DO 40 J=1,JM*KCON$ 40 CONSRV(J,1) = 0. C**** C**** Loop over input disk files C**** IHOURL = -123456789 KFILE=0 100 KFILE=KFILE+1 IF(KFILE.gt.NARGS) GO TO 200 CALL GETARG (KFILE,FILEIN) IF(ICHAR(FILEIN).ge.Z'30' .and. ICHAR(FILEIN).le.Z'39') GO TO 200 IF(FILEIN(1:2).eq.'AT' .or. FILEIN(1:2).eq.'AG' .or. * FILEIN(1:2).eq.'AI' .or. FILEIN(1:2).eq.'AL' .or. * FILEIN(1:2).eq.'AO') GO TO 200 C**** Open and read input file OPEN (1,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) READ (1) IHOURX,LABEL,IPARM,IDIAG,(DIAGR4(K),K=1,KDACCX),IHOURY CLOSE (1) IF(IHOURX.ne.IHOURY) GO TO 820 CALL STORYM (JYEAR0,JMON0, JYEAR,JMON) WRITE (6,911) JYEAR0,JMON0, JYEAR,JMON, FILEIN(1:60) C**** Accumulate diagnostics for time averaging DO 110 N=1,15 110 IDACC(N) = IDACC(N) + IDIAG(N) DO 120 J=1,JM*KAJ$*6 120 AJ(J,1) = AJ(J,1) + DIAGR4(J) DO 130 I=1,IM*JM*LMA$ 130 AIJL(I,1,1,9) = AIJL(I,1,1,9) + DIAGR4(KAIJL0+IM*JM*LMA$*8+I) DO 160 K=1,KCON$ DO 140 I=1,21 140 IF(K.eq.INSTAN(I)) GO TO 160 DO 150 J=1,JM 150 CONSRV(J,K) = CONSRV(J,K) + DIAGR4(KCON0+JM*(K-1)+J) 160 continue IF(IHOURX.lt.IHOURL) GO TO 180 IHOURL = IHOURX DO 170 I=1,21 K = INSTAN(I) DO 170 J=1,JM 170 CONSRV(J,K) = DIAGR4(KCON0+JM*(K-1)+J) 180 GO TO 100 C**** C**** End of input files: print table of input months and years 200 CALL PRNTYM (OUTYR,OUTMON) C**** C**** Call diagnostic subroutines C**** CALL DIAGCP CALL DIAGJ CLOSE (6) IF(KFILE.gt.NARGS) GO TO 999 C**** C**** Interpret last few arguments from command line as conservation C**** diagnostics or as areas C**** LPMAX = NARGS - KFILE +1 IF(LPMAX.gt.16) LPMAX = 16 ! arrays are dimensioned by 16 DO 310 LP=1,LPMAX LPINDX(LP) = 0 CALL GETARG (KFILE,LPLAB(LP)) IF(ICHAR(LPLAB(LP)).ge.Z'30' .and. ICHAR(LPLAB(LP)).le.Z'39') * READ (LPLAB(LP),*) LPINDX(LP) IF(LPLAB(LP)(1:2).eq.'AT') LPINDX(LP) = KCON$ + 1 IF(LPLAB(LP)(1:2).eq.'AG') LPINDX(LP) = KCON$ + 2 IF(LPLAB(LP)(1:2).eq.'AI') LPINDX(LP) = KCON$ + 3 IF(LPLAB(LP)(1:2).eq.'AL') LPINDX(LP) = KCON$ + 4 IF(LPLAB(LP)(1:2).eq.'AO') LPINDX(LP) = KCON$ + 5 IF(LPINDX(LP).le.0 .or. LPINDX(LP).gt.KCON$+5) GO TO 830 310 KFILE = KFILE + 1 C**** Convert the areas from (m^2) to (10^10 m^2) DO 320 N=KCON+1,KCON+5 DO 320 J=1,JM 320 FLAT(J,N) = FLAT(J,N)*1.d-10 C**** C**** Write line plot file of conservation diagnostics C**** OPEN (2,FILE='CONS070.LP') WRITE (2,940) LABEL,OUTYR,OUTMON WRITE (2,*) 'Latitude' WRITE (2,*) 'Conservation Diagnostic or Area' WRITE (2,941) (LPLAB(LP),LP=1,LPMAX) DO 420 J=1,JM 420 WRITE (2,942) RLAT(J)*360./TWOPI, (FLAT(J,LPINDX(LP)),LP=1,LPMAX) WRITE (2,943) WRITE (2,943) 'NH', (FLAT(-1,LPINDX(LP)),LP=1,LPMAX) WRITE (2,943) 'SH', (FLAT( 0,LPINDX(LP)),LP=1,LPMAX) WRITE (2,943) 'GL', (FLAT(-2,LPINDX(LP)),LP=1,LPMAX) CLOSE (2) GO TO 999 C**** 800 WRITE (0,*) 'Usage: CONS070 /raid1/C070/DAnn* ', * '[39 40 ... AT AG AI AL AO] 1999/07/02' WRITE (0,*) 'print DIAGCP and DIAGJ Make op', * 'tional line plot file of conservation diags' GO TO 999 810 WRITE (0,*) 'Error accessing ',FILEIN STOP 810 820 WRITE (0,*) 'IHOURX and IHOURY do not match:',IHOURX,IHOURY STOP 820 830 WRITE (0,*) 'Unable to decipher conservation diagnostic: ', * LPLAB(LP) STOP 830 C**** 911 FORMAT (' From',I6,A5,' to',I6,A5,' FILEIN=',A) 940 FORMAT ('Conservation Diagnostics',5X,A5,A9,1X,A3) 941 FORMAT (' Lat',8X,16A11) 942 FORMAT (F5.1,16(1X,F11.3)) 943 FORMAT (A3,2X,16(1X,F11.3)) 999 END SUBROUTINE SGEOM IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (IM=72,JM=46, * TWOPI=6.283185307179586477d0, RADIUS=6375000.) COMMON /GEOMCB/ DXYP(JM),DXP(JM),DYP(JM),DXV(JM),DYV(JM), * COSV(0:JM),COSP(JM),SINP(JM),DXYV(JM),DXYS(JM),DXYN(JM), * RAMVS(JM),RAMVN(JM),RLAT(JM),SINI(IM),COSI(IM),BYDXYP(JM) C**** C**** Calculates the spherical geometry for the C grid C**** DLON = TWOPI/IM DLAT = TWOPI*NINT(360./(JM-1))/720. FJEQ = .5*(1+JM) C**** Geometric parameters centered at secondary latitudes RLAT(1)= -TWOPI/4. PSINS = SIN(RLAT(1)) DO 10 J=1,JM-1 RLAT(J+1) = DLAT*(J+1.-FJEQ) IF(J.eq.JM-1) RLAT(JM) = TWOPI/4. PSINN = SIN(RLAT(J+1)) DXV(J) = DLON*RADIUS*(PSINN-PSINS) / (RLAT(J+1)-RLAT(J)) 10 PSINS = PSINN C**** Geometric parameters centered at primary latitudes VLATS = -TWOPI/4. VSINS = SIN(VLATS) DO 20 J=1,JM VLATN = DLAT*(J+.5-FJEQ) IF(J.eq.JM) VLATN = TWOPI/4. VSINN = SIN(VLATN) DXYP(J)= DLON*RADIUS*RADIUS*(VSINN-VSINS) BYDXYP(J) = 1./DXYP(J) DXP(J) =.5*(DXV(J-1)+DXV(J)) DYP(J) = (VLATN-VLATS)*RADIUS VLATS = VLATN 20 VSINS = VSINN DXP(1) = .5*DXV(1) DXP(JM)= .5*DXV(JM-1) RETURN END SUBROUTINE STORYM (JYEAR0,JMON0, JYEAR,JMON) C**** C**** STORYM receives years and months and stores them in an array. C**** C**** Start of accumulation period: JYEAR0, JMON0, day 1, hour 0 C**** End of accumulation period: JYEAR, JMON, day 1, hour 0 C**** INTEGER*4 NYofM(12) CHARACTER*1 QMY(12,0:2500) CHARACTER*4 JMON0,JMON, MON(12), OUTYR,OUTMON DATA QMY /30012*' '/, NYofM/12*0/, MINYR/2501/, MAXYR/-1/ DATA MON /'Jan','Feb','Mar','Apr','May','Jun', * 'Jul','Aug','Sep','Oct','Nov','Dec'/ C**** IF(JMON0.eq.'Ann ') JMON0 = 'Jan ' ! compatible with prior vers IF(JYEAR0.gt.2500) JYEAR0 = 1900 + JYEAR0/100 DO 10 I=1,12 10 IF(JMON0.eq.MON(I)) GO TO 20 WRITE (6,901) JYEAR0,JMON0, JYEAR,JMON, 'JMON0 is not a month.' WRITE (0,901) JYEAR0,JMON0, JYEAR,JMON, 'JMON0 is not a month.' STOP 10 20 IMON0 = I DO 30 I=1,12 30 IF(JMON.eq.MON(I)) GO TO 40 WRITE (6,901) JYEAR0,JMON0, JYEAR,JMON, 'JMON is not a month.' WRITE (0,901) JYEAR0,JMON0, JYEAR,JMON, 'JMON is not a month.' STOP 30 40 IMON = I IF(IMON.gt.IMON0) then IMON = IMON-1 IYEAR = JYEAR else IMON = IMON+11 IYEAR = JYEAR-1 endif C**** DO 50 IY=JYEAR0,IYEAR DO 50 IM=IMON0 ,IMON IF(QMY(IM,IY).eq.'*') GO TO 805 QMY(IM,IY) = '*' JM = 1 + MOD(IM-1,12) NYofM(JM) = NYofM(JM) + 1 IF(JYEAR0.lt.MINYR) MINYR = JYEAR0 50 IF(JYEAR .gt.MAXYR) MAXYR = JYEAR RETURN C**** C**** C**** ENTRY PRNTYM (OUTYR,OUTMON) C**** C**** PRNTYM prints a table of input years and month, and C**** calculates OUTYR and OUTMON C**** C**** Output: OUTYR = character describing input years C**** OUTMON = character describing input months C**** C**** Produce printer table of months and years received C**** WRITE (6,910) DO 110 IY=MINYR,MAXYR 110 WRITE (6,911) IY,(QMY(IM,IY),IM=1,12) C**** C**** Determine output month: a single month, three months, or an annual C**** DO 210 M=1,12 210 IF(NYofM(M).gt.0) GO TO 220 GO TO 800 C**** Is output month a single month ? 220 M1 = M DO 230 M=M1+1,12 230 IF(NYofM(M).gt.0) GO TO 240 OUTMON = MON(M1) GO TO 400 C**** Is output twelve months ? 240 IF(M1.ne.1) GO TO 260 DO 250 M=2,12 250 IF(NYofM(M).ne.NYofM(M1)) GO TO 260 OUTMON = 'Ann ' GO TO 400 C**** Is output month three consecutive months ? 260 IF(M1.ge.11) GO TO 800 IF(NYofM(M1+1).ne.NYofM(M1)) GO TO 280 IF(NYofM(M1+2).ne.NYofM(M1)) GO TO 300 DO 270 M=M1+3,12 270 IF(NYofM(M).ne.0) GO TO 800 OUTMON = MON(M1)(1:1) // MON(M1+1)(1:1) // MON(M1+2)(1:1) // ' ' GO TO 400 C**** Is output 3 months November, December and January ? 280 IF(NYofM(11).ne.NYofM(1) .or. NYofM(12).ne.NYofM(1)) GO TO 800 DO 290 M=2,10 290 IF(NYofM(M).ne.0) GO TO 800 OUTMON = 'NDJ ' GO TO 400 C**** Is output 3 months December, January and February ? 300 IF(NYofM(12).ne.NYofM(1) .or. NYofM(2).ne.NYofM(1)) GO TO 800 DO 310 M=3,11 310 IF(NYofM(M).ne.0) GO TO 800 OUTMON = 'DJF ' C**** C**** Calculate the first and last years received, C**** and determine a character string describing those years C**** 400 MAXYR = MINYR + NYofM(M1) - 1 WRITE (OUTYR,940) MINYR IF(MINYR.lt.MAXYR) * WRITE (OUTYR,941) MOD(MINYR,100),MOD(MAXYR,100) IF(MINYR+8.le.MAXYR .and. MINYR/10.eq.MAXYR/10) * WRITE (OUTYR,942) MINYR/10 WRITE (6,943) OUTYR,OUTMON RETURN C**** 800 WRITE (6,*) 'From CONSC070: inconsistant year and months received' WRITE (0,*) 'From CONSC070: inconsistant year and months received' STOP 800 805 WRITE (6,901) JYEAR0,JMON0, JYEAR,JMON, 'accumulated already' WRITE (0,901) JYEAR0,JMON0, JYEAR,JMON, 'accumulated already' STOP 805 C**** 901 FORMAT (' From',I6,A5,' to',I6,A5,1X,A) 910 FORMAT ('0The following months and years were received:' / '0', * 9X,'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec') 911 FORMAT (1X,I6,21A5) 940 FORMAT (I4) 941 FORMAT (2I2.2) 942 FORMAT (I3,'X') 943 FORMAT ('0Description of the years and months above:',2A5) END SUBROUTINE READR4 (IUNIT,NM,DATAR4,DATAR8) C**** C**** READR4 reads a record from unit IUNIT containing TITLE,DATAR4, C**** converts the REAL*4 array DATAR4 to the REAL*8 array DATAR8, and C**** writes a line to unit 6 containing the TITLE just read. C**** REAL*4 DATAR4(NM) REAL*8 DATAR8(NM) CHARACTER*80 TITLE C**** READ (IUNIT) TITLE,DATAR4 DO 10 I=NM,1,-1 10 DATAR8(I) = DATAR4(I) WRITE (6,901) IUNIT,TITLE RETURN C**** 901 FORMAT (' Read on unit',I3,': ',A80) END SUBROUTINE TIMER (MNOW,MSUM) C**** C**** Output: MNOW (.01 s) = current CPU time C**** MSUM (.01 s) = MSUM + time since last call to TIMER C**** DATA MLAST /0/ MNOW = MCLOCK () MSUM = MNOW - MLAST + MSUM MLAST = MNOW RETURN C**** C**** C**** TIMER0 and TIMER1 are used to time an individual subroutine and C**** to attribute its cost to a short process (e.g. diagnostics) while C**** a longer process is being timed both before and after the short C**** subroutine was called. C**** ENTRY TIMER0 MNOW0 = MCLOCK () RETURN C**** C**** ENTRY TIMER1 (MSUM) MNOW1 = MCLOCK () MSUM = MNOW1 - MNOW0 + MSUM MLAST = MNOW1 - MNOW0 + MLAST RETURN END