C**** ATJ070.S 98/11/23 C**** Line Plot file of Atmospheric Energy Transport C**** C**** ATJ070 reads several Climate Model III D files, calculates C**** northward transport of dry static energy and latent energy, C**** sums them vertically and longitudinally, and C**** writes the scaled data of J as a line plot file. 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**** REAL*8 TOT(JM),DSE(JM),ENT(JM),GPE(JM),ELH(JM) CHARACTER FILEIN*80, TITLE*80, OUTYR*4,OUTMON*4, * NAME*50,XLAB*8,YLAB*21,LINE*72 COMMON /WORK01/ PE(0:LMA,IM,JM), PUP(LMA,IM,JM),PDN(LMA,IM,JM), * TUP(LMA,IM,JM),TDN(LMA,IM,JM), AUP(LMA,IM,JM),ADN(LMA,IM,JM), * EZ(0:LMA,IM,JM), Z(LMA,IM,JM), PK(LMA,IM,JM) C**** DATA NAME /'NORTHWARD TRANSPORT of STATIC ENERGY'/ DATA XLAB /'Latitude'/ DATA YLAB /'Energy Transport (PW)'/ DATA LINE / *' Lat Total Dry_Stat_E Enthalpy Geopo_Enrg Latent_Enrg'/ C**** NARGS = IARGC() IF(NARGS.le.0) GO TO 800 OPEN (6,FILE='ATJ070.PRT') C**** Read first input file to determine KAIJL0 and 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$ KDACCX = KAIJL0 + IM*JM*LMA$*KAIJL$ + IM*JM*LMO$*KOIJL$ + * LMO$*NMST$*KOLNS$ + 24*50*4 + JM*KCON$ WRITE (6,*) 'KDINT,KDACCX=',KDINT,KDACCX CALL SGEOM C**** C**** Read in atmospheric topography C**** FILEIN = '/u/cmrun/Z72X46N' OPEN (11,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) READ (11) READ (11) READ (11) READ (11) CALL READR4 (11,IM*JM,ZATMO,ZATMO) CLOSE (11) C**** C**** Zero out summing arrays C**** DO 10 N=1,15 10 IDACC(N) = 0 DO 20 K=1,11,2 DO 20 I=1,IM*JM*LMA 20 AIJL(I,1,1,K) = 0. C**** C**** Loop over input D files C**** KFILE = 0 100 KFILE = KFILE + 1 IF(KFILE.gt.NARGS) GO TO 200 CALL GETARG (KFILE,FILEIN) C**** Open and read input D 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 K=1,11,2 DO 120 I=1,IM*JM*LMA 120 AIJL(I,1,1,K) = AIJL(I,1,1,K) + DIAGR4(KAIJL0+IM*JM*LMA*(K-1)+I) GO TO 100 C**** C**** End of input D files: print table of input years and months C**** 200 CALL PRNTYM (OUTYR,OUTMON) C**** C**** Calculate mass MA (kg/m^2); pressures PE, PUP and PDN (Pa); C**** temperatures TUP and TDN (K); and C**** specific volumes AUP and ADN (m^3/kg) C**** DO 210 J=1,JM DO 210 I=1,IM PE(LMA,I,J) = MSTRAT DO 210 L=LMA,1,-1 MA(I,J,L) = AIJL(I,J,L,1)/IDACC(4) H0M(I,J,L) = AIJL(I,J,L,5)/IDACC(4) Q0M(I,J,L) = AIJL(I,J,L,9)/IDACC(4) PUP(L,I,J) = PE(L,I,J) + GRAV*(.5-.5*RRT3)*MA(I,J,L) PDN(L,I,J) = PE(L,I,J) + GRAV*(.5+.5*RRT3)*MA(I,J,L) PK(L,I,J) = .5*(PDN(L,I,J)**RKAP+PUP(L,I,J)**RKAP) TUP(L,I,J) = H0M(I,J,L)*PUP(L,I,J)**RKAP/(SHCD*DXYP(J)*MA(I,J,L)) TDN(L,I,J) = H0M(I,J,L)*PDN(L,I,J)**RKAP/(SHCD*DXYP(J)*MA(I,J,L)) AUP(L,I,J) = RDRY*TUP(L,I,J)/PUP(L,I,J) ADN(L,I,J) = RDRY*TDN(L,I,J)/PDN(L,I,J) 210 PE(L-1,I,J) = PE(L,I,J) + GRAV*MA(I,J,L) C**** C**** Calculate altitudes EZ and Z (m) C**** DO 220 J=1,JM DO 220 I=1,IM EZ(0,I,J) = ZATMO(I,J) DO 220 L=1,LMA Z(L,I,J) = EZ(L-1,I,J) + .25*MA(I,J,L)* * ((1.-RRT3)*ADN(L,I,J)+(1.+RRT3)*AUP(L,I,J)) 220 EZ(L,I,J) = EZ(L-1,I,J) + .5*MA(I,J,L)*(ADN(L,I,J)+AUP(L,I,J)) C**** C**** Calculate energy transports C**** DO 320 J=1,JM-1 EV = 0. GV = 0. QV = 0. DO 310 I=1,IM DO 310 L=1,LMA EV = EV + AIJL(I,J,L,7)*(PK(L,I,J)+PK(L,I,J+1)) GV = GV + AIJL(I,J,L,3)*( Z(L,I,J)+ Z(L,I,J+1)) 310 QV = QV + AIJL(I,J,L,11) ENT(J) = .5d-15* EV/(IDACC(1)*DTS) GPE(J) = 1.d-15*GRAV*GV/(IDACC(1)*NDYNA) ELH(J) = 1.d-15*ELHE*QV/(IDACC(1)*DTS) DSE(J) = ENT(J) + GPE(J) 320 TOT(J) = DSE(J) + ELH(J) C**** C**** Write line plot output file C**** WRITE (TITLE,940) NAME,LABEL(1:5),OUTYR,OUTMON OPEN (2,FILE='ATJ070.LP') WRITE (2,*) TITLE(1:72) WRITE (2,*) XLAB WRITE (2,*) YLAB WRITE (2,*) LINE DO 410 J=1,JM-1 WRITE (2,941) RLAT(J),TOT(J),DSE(J),ENT(J),GPE(J),ELH(J) 410 WRITE (6,941) RLAT(J),TOT(J),DSE(J),ENT(J),GPE(J),ELH(J) C**** CLOSE (2) CLOSE (6) GO TO 999 C**** 800 WRITE (0,*) 'Usage: ATJ070 /raid2/C070/DAnn199* ', * 'Northward Transports 98/11/23' GO TO 999 810 WRITE (0,*) 'Error accessing ',FILEIN STOP 810 820 WRITE (0,*) 'IHOURX and IHOURY do not match:',IHOURX,IHOURY STOP 820 C**** 911 FORMAT (' From',I6,A5,' to',I6,A5,' FILEIN=',A) 940 FORMAT (A50,A5,9X,A4,1X,A3) 941 FORMAT (F6.1,5F12.4) 999 END SUBROUTINE SGEOM IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (IM=72,JM=46,LMO=9, * 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) 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 DO 10 J=1,JM-1 10 RLAT(J) = 360.*DLAT*(J+.5-FJEQ)/TWOPI 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) VLATS = VLATN 20 VSINS = VSINN 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 ATJ070: inconsistant year and months received' WRITE (0,*) 'From ATJ070: 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