C**** AIJK070.S make Atmospheric IxJ data file 98/11/03 C**** C**** AIJK070 reads several Climate Model III D files, calculates C**** specified quantities from the AIJ and AIJL accumulating arrays, C**** and writes the scaled data to a disk file. Each record of the C**** output file contains an 80 byte title and a REAL*4 IMxJM array. 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$),(IDIAG,IDACC) C**** PARAMETER (PPVSF=610.571d0, TKF=273.16d0) INTEGER*4 KU(32), KV(32),KMW(32), KT(32),KTC(32),KTV(32),KTW(32), * KQ(32),KQC(32),KQV(32),KQW(32), KZ(32),KRH(32) INTEGER*4 IDACCS(15) CHARACTER FILEIN*80, TITLE*80, NAME*50, OUTYR*4,OUTMON*4 COMMON XK(IM,JM), X0(LMA,IM,JM),XE(0:LMA,IM,JM), MS(IM,JM), * XL(LMA-1), DXD2M(LMA-1), DS2(LMA-1), * byDS1(LMA),byDS2(LMA-1),byDS3(2:LMA-1),byDS4(2:LMA-2) COMMON /NAMECB/ NAME(14) C**** NAMELIST /INPUTZ/ KU,KV,KMW, KT,KTC,KTV,KTW, KQ,KQC,KQV,KQW,KZ,KRH DATA KU/32*0/, KV/32*0/, KMW/32*0/, KT/32*0/, KTC/32*0/, * KTV/32*0/, KTW/32*0/, KQ/32*0/, KQC/32*0/, KQV/32*0/, * KQW/32*0/, KZ/32*0/, KRH/32*0/ C**** C**** MS = dynamical column mass (kg/m^2) C**** XK = quantity interpolated to requested pressure level C**** X0 = average quantity over GCM layer C**** XE = quantity interpolated to GCM layer edges C**** XL = quantity linearly interpolated to GCM layer edges C**** C PPVSAT(TK,EL) = PPVSF*EXP(EL*(1./TKF-1./TK)/RVAP) NARGS = IARGC() IF(IARGC().le.0) GO TO 800 OPEN (6,FILE='AIJK070.PRT') C**** Read first input file to determine KAIJ0, KAIJL0, KDACCx CALL GETARG (1,FILEIN) OPEN (1,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) READ (1) IHOURX,LABEL,IPARM CLOSE (1) KAIJ0 = JM*KAJ$*6 + JM*LMA$*KAJL$ + JM*3*4 KAIJL0 = KAIJ0 + 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 C**** C**** Determine output quantities from Namelist input file C**** FILEIN = 'AIJK070.I' OPEN (5,FILE=FILEIN,STATUS='OLD',ERR=810) READ (5,INPUTZ) CLOSE (5) CALL SGEOM C**** C**** Zero out summing arrays C**** DO 10 N=1,15 10 IDACCS(N) = 0 DO 20 I=1,IM*JM AIJ(I,1,14) = 0. 20 AIJ(I,1,35) = 0. DO 30 I=1,IM*JM*LMA*12 30 AIJL(I,1,1,1) = 0. C**** C**** Read in fixed arrays C**** FILEIN = '/u/cmrun/Z72X46N' OPEN (11,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=810) CALL READR4 (11,IM*JM,FOCEAN,FOCEAN) READ (11) CALL READR4 (11,IM*JM,FGRND,FGRND) CALL READR4 (11,IM*JM,FGICE,FGICE) CALL READR4 (11,IM*JM,ZATMO,ZATMO) CLOSE (11) C**** C**** Loop over the input disk files C**** KFILE=0 100 KFILE=KFILE+1 IF(KFILE.gt.NARGS) GO TO 200 CALL GETARG (KFILE,FILEIN) 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 160 N=1,15 160 IDACCS(N) = IDACCS(N) + IDACC(N) DO 170 I=1,IM*JM AIJ(I,1,14) = AIJ(I,1,14) + DIAGR4(KAIJ0+IM*JM*13+I) 170 AIJ(I,1,35) = AIJ(I,1,35) + DIAGR4(KAIJ0+IM*JM*34+I) DO 180 I=1,IM*JM*LMA*12 180 AIJL(I,1,1,1) = AIJL(I,1,1,1) + DIAGR4(KAIJL0+I) GO TO 100 C**** C**** End of input files: print table of input years and months C**** 200 CALL PRNTYM (OUTYR,OUTMON) OPEN (2,FILE='AIJK070.O',FORM='UNFORMATTED') C**** C**** Calculate dynamical column mass MS (kg/m^2) C**** DO 201 I=1,IM AIJ(I, 1,14) = AIJ(IM,1,14) 201 AIJ(I,JM,14) = AIJ(1,JM,14) DO 202 I=1,IM*JM 202 MS(I,1) = (AIJ(I,1,14)/IDACCS(1) + 101325.)/GRAV - MSTRAT C**** Calculate DSIGA sums and their reciprocals DO 203 L=1,LMA DSIGA(L) = SIGEA(L-1) - SIGEA(L) 203 byDS1(L) = 1. / DSIGA(L) DO 204 L=1,LMA-1 DS2(L) = DSIGA(L) + DSIGA(L+1) 204 byDS2(L) = 1. / DS2(L) DO 205 L=2,LMA-1 205 byDS3(L) = 1. / (DSIGA(L-1)+DS2(L)) DO 206 L=2,LMA-2 206 byDS4(L) = 1. / (DS2(L-1)+DS2(L+1)) C**** C**** Eastward velocity (m/s) C**** IF(KU(1).le.0) GO TO 230 I=IM DO 224 J=2,JM-1 DO 224 Ip1=1,IM C**** Calculate X0 = UA (m/s) for GCM layers DO 221 L=1,LMA 221 X0(L,I,J) = 4.*AIJL(I,J,L,2)*byDS1(L) / / (NDYNA*IDACCS(1)*(MS(I,J)+MS(Ip1,J))*DYP(J)) C**** Calculate XL, the linearly interpolated value, C**** and DXD2M, the derivative of X0 divided by 2 DO 222 L=1,LMA-1 XL(L) = (X0(L,I,J)*DSIGA(L+1)+X0(L+1,I,J)*DSIGA(L))*byDS2(L) 222 DXD2M(L) = (X0(L,I,J)-X0(L+1,I,J))*byDS2(L) C**** Calculate edge values XE (m/s) using cubic, quadratic or linear C**** interpolation schemes XE(0 ,I,J) = X0(1 ,I,J) + DXD2M(1 )*DSIGA(1) XE(LMA,I,J) = X0(LMA,I,J) - DXD2M(LMA-1)*DSIGA(LMA) XE(1,I,J) = XL(1) - (DXD2M(1)-DXD2M(2))*DSIGA(1)*DSIGA(2)*byDS3(2) XE(LMA-1,I,J) = XL(LMA-1) + + (DXD2M(LMA-1)-DXD2M(LMA-2))*DSIGA(LMA)*DSIGA(LMA-1)*byDS3(LMA-1) DO 223 L=2,LMA-2 223 XE(L,I,J) = XL(L) + ((DXD2M(L)-DXD2M(L-1))*DS2(L+1)*byDS3(L) + + (DXD2M(L+1)-DXD2M(L))*DS2(L)*byDS3(L+1))* * DSIGA(L)*DSIGA(L+1)*byDS4(L) 224 I=Ip1 C**** Use quadratic vertical interpolation scheme to calulate quantity C**** at constant pressure coordinate DO 226 K=1,32 IF(KU(K).le.0) GO TO 230 I=IM DO 225 J=2,JM-1 DO 225 Ip1=1,IM SG = (KU(K)*100/GRAV - MSTRAT)*2. / (MS(I,J)+MS(Ip1,J)) XK(I,J) = * VNLtoK (LMA, SG, SIGEA,byDS1, X0(1,I,J),XE(0,I,J), -999999.d0) 225 I=Ip1 WRITE (NAME(2)(28:31),922) KU(K) 226 CALL WRITED (NAME(2),LABEL,OUTYR,OUTMON,'U') C**** C**** Northward velocity (m/s) C**** 230 IF(KV(1).le.0) GO TO 240 DO 234 J=1,JM-1 DO 234 I=1,IM C**** Calculate X0 = VA (m/s) for GCM layers DO 231 L=1,LMA 231 X0(L,I,J) = 4.*AIJL(I,J,L,3)*byDS1(L) / / (NDYNA*IDACCS(1)*(MS(I,J)+MS(I,J+1))*DXV(J)) C**** Calculate XL, the linearly interpolated value, C**** and DXD2M, the derivative of X0 divided by 2 DO 232 L=1,LMA-1 XL(L) = (X0(L,I,J)*DSIGA(L+1)+X0(L+1,I,J)*DSIGA(L))*byDS2(L) 232 DXD2M(L) = (X0(L,I,J)-X0(L+1,I,J))*byDS2(L) C**** Calculate edge values XE (m/s) using cubic, quadratic or linear C**** interpolation schemes XE(0 ,I,J) = X0(1 ,I,J) + DXD2M(1 )*DSIGA(1) XE(LMA,I,J) = X0(LMA,I,J) - DXD2M(LMA-1)*DSIGA(LMA) XE(1,I,J) = XL(1) - (DXD2M(1)-DXD2M(2))*DSIGA(1)*DSIGA(2)*byDS3(2) XE(LMA-1,I,J) = XL(LMA-1) + + (DXD2M(LMA-1)-DXD2M(LMA-2))*DSIGA(LMA)*DSIGA(LMA-1)*byDS3(LMA-1) DO 233 L=2,LMA-2 233 XE(L,I,J) = XL(L) + ((DXD2M(L)-DXD2M(L-1))*DS2(L+1)*byDS3(L) + + (DXD2M(L+1)-DXD2M(L))*DS2(L)*byDS3(L+1))* * DSIGA(L)*DSIGA(L+1)*byDS4(L) 234 continue C**** Use quadratic vertical interpolation scheme to calulate quantity C**** at constant pressure coordinate DO 236 K=1,32 IF(KV(K).le.0) GO TO 240 DO 235 J=1,JM-1 DO 235 I=1,IM SG = (KV(K)*100/GRAV - MSTRAT)*2. / (MS(I,J)+MS(I,J+1)) 235 XK(I,J) = * VNLtoK (LMA, SG, SIGEA,byDS1, X0(1,I,J),XE(0,I,J), -999999.d0) WRITE (NAME(3)(29:32),922) KV(K) 236 CALL WRITED (NAME(3),LABEL,OUTYR,OUTMON,'V') C**** 240 continue C**** C**** Temperature (K) C**** 250 IF(KT(1).le.0) GO TO 260 DO 254 J=1,JM IMAX=IM IF(J.eq.1 .or. J.eq.JM) IMAX=1 DO 254 I=1,IMAX C**** Calculate X0 = PSE (J/kg) for GCM layers DO 251 L=1,LMA 251 X0(L,I,J) = AIJL(I,J,L,5) / (AIJL(I,J,L,1)*DXYP(J)) C**** Calculate XL, the linearly interpolated value, C**** and DXD2M, the derivative of X0 divided by 2 DO 252 L=1,LMA-1 XL(L) = (X0(L,I,J)*DSIGA(L+1)+X0(L+1,I,J)*DSIGA(L))*byDS2(L) 252 DXD2M(L) = (X0(L,I,J)-X0(L+1,I,J))*byDS2(L) C**** Calculate edge values XE (J/kg) using cubic, quadratic or linear C**** interpolation schemes XE(0 ,I,J) = SHCD*(273.16 + AIJ(I,J,35)/(IDACCS(1)*NSURF)) / * (101325.+ AIJ(I,J,14)/ IDACCS(1))**RKAP C XE(0 ,I,J) = X0(1 ,I,J) + DXD2M(1 )*DSIGA(1) XE(LMA,I,J) = X0(LMA,I,J) - DXD2M(LMA-1)*DSIGA(LMA) XE(1,I,J) = XL(1) - (DXD2M(1)-DXD2M(2))*DSIGA(1)*DSIGA(2)*byDS3(2) XE(LMA-1,I,J) = XL(LMA-1) + + (DXD2M(LMA-1)-DXD2M(LMA-2))*DSIGA(LMA)*DSIGA(LMA-1)*byDS3(LMA-1) DO 253 L=2,LMA-2 253 XE(L,I,J) = XL(L) + ((DXD2M(L)-DXD2M(L-1))*DS2(L+1)*byDS3(L) + + (DXD2M(L+1)-DXD2M(L))*DS2(L)*byDS3(L+1))* * DSIGA(L)*DSIGA(L+1)*byDS4(L) 254 continue C**** Use quadratic vertical interpolation scheme to calulate quantity C**** at constant pressure coordinate DO 256 K=1,32 IF(KT(K).le.0) GO TO 260 DO 255 J=1,JM IMAX=IM IF(J.eq.1 .or. J.eq.JM) IMAX=1 DO 255 I=1,IMAX SG = (KT(K)*100/GRAV - MSTRAT) / MS(I,J) XK(I,J) = * VNLtoK (LMA, SG, SIGEA,byDS1, X0(1,I,J),XE(0,I,J), -999999.d0) 255 IF(XK(I,J).gt.-999999.) * XK(I,J) = ((KT(K)*100.d0)**RKAP / SHCD) * XK(I,J) WRITE (NAME(5)(20:23),922) KT(K) 256 CALL WRITED (NAME(5),LABEL,OUTYR,OUTMON,'A') C**** 260 continue C**** C**** Specific humidity (g/kg) C**** IF(KQ(1).le.0) GO TO 300 DO 294 J=1,JM IMAX=IM IF(J.eq.1 .or. J.eq.JM) IMAX=1 DO 294 I=1,IMAX C**** Calculate X0 = Q (m/s) for GCM layers DO 291 L=1,LMA 291 X0(L,I,J) = AIJL(I,J,L,9) / (AIJL(I,J,L,1)*DXYP(J)) C**** Calculate XL, the linearly interpolated value, C**** and DXD2M, the derivative of X0 divided by 2 DO 292 L=1,LMA-1 XL(L) = (X0(L,I,J)*DSIGA(L+1)+X0(L+1,I,J)*DSIGA(L))*byDS2(L) 292 DXD2M(L) = (X0(L,I,J)-X0(L+1,I,J))*byDS2(L) C**** Calculate edge values XE (m/s) using cubic, quadratic or linear C**** interpolation schemes XE(0 ,I,J) = AIJ(I,J,50) / IDACCS(3) C XE(0 ,I,J) = X0(1 ,I,J) + DXD2M(1 )*DSIGA(1) XE(LMA,I,J) = X0(LMA,I,J) - DXD2M(LMA-1)*DSIGA(LMA) XE(1,I,J) = XL(1) - (DXD2M(1)-DXD2M(2))*DSIGA(1)*DSIGA(2)*byDS3(2) XE(LMA-1,I,J) = XL(LMA-1) + + (DXD2M(LMA-1)-DXD2M(LMA-2))*DSIGA(LMA)*DSIGA(LMA-1)*byDS3(LMA-1) DO 293 L=2,LMA-2 293 XE(L,I,J) = XL(L) + ((DXD2M(L)-DXD2M(L-1))*DS2(L+1)*byDS3(L) + + (DXD2M(L+1)-DXD2M(L))*DS2(L)*byDS3(L+1))* * DSIGA(L)*DSIGA(L+1)*byDS4(L) 294 continue C**** Use quadratic vertical interpolation scheme to calulate quantity C**** at constant pressure coordinate DO 296 K=1,32 IF(KQ(K).le.0) GO TO 300 DO 295 J=1,JM IMAX=IM IF(J.eq.1 .or. J.eq.JM) IMAX=1 DO 295 I=1,IMAX SG = (KQ(K)*100/GRAV - MSTRAT) / MS(I,J) XK(I,J) = * VNLtoK (LMA, SG, SIGEA,byDS1, X0(1,I,J),XE(0,I,J), -999999.d0) 295 IF(XK(I,J).gt.-999999.) XK(I,J) = XK(I,J)*1000. WRITE (NAME(9)(29:32),922) KQ(K) 296 CALL WRITED (NAME(9),LABEL,OUTYR,OUTMON,'A') C**** 300 continue C**** C**** Geopotential Height (m) C**** 330 IF(KZ(1).le.0) GO TO 340 DO 334 J=1,JM IMAX=IM IF(J.eq.1 .or. J.eq.JM) IMAX=1 DO 334 I=1,IMAX C**** Calculate X0 = PSE (J/kg) for GCM layers DO 331 L=1,LMA 331 X0(L,I,J) = AIJL(I,J,L,5) / (AIJL(I,J,L,1)*DXYP(J)) C**** Calculate XL, the linearly interpolated value, C**** and DXD2M, the derivative of X0 divided by 2 DO 332 L=1,LMA-1 XL(L) = (X0(L,I,J)*DSIGA(L+1)+X0(L+1,I,J)*DSIGA(L))*byDS2(L) 332 DXD2M(L) = (X0(L,I,J)-X0(L+1,I,J))*byDS2(L) C**** Calculate edge values XE (J/kg) using cubic, quadratic or linear C**** interpolation schemes XE(0 ,I,J) = SHCD*(273.16 + AIJ(I,J,35)/(IDACCS(1)*NSURF)) / * (101325.+ AIJ(I,J,14)/ IDACCS(1))**RKAP C XE(0 ,I,J) = X0(1 ,I,J) + DXD2M(1 )*DSIGA(1) XE(LMA,I,J) = X0(LMA,I,J) - DXD2M(LMA-1)*DSIGA(LMA) XE(1,I,J) = XL(1) - (DXD2M(1)-DXD2M(2))*DSIGA(1)*DSIGA(2)*byDS3(2) XE(LMA-1,I,J) = XL(LMA-1) + + (DXD2M(LMA-1)-DXD2M(LMA-2))*DSIGA(LMA)*DSIGA(LMA-1)*byDS3(LMA-1) DO 333 L=2,LMA-2 333 XE(L,I,J) = XL(L) + ((DXD2M(L)-DXD2M(L-1))*DS2(L+1)*byDS3(L) + + (DXD2M(L+1)-DXD2M(L))*DS2(L)*byDS3(L+1))* * DSIGA(L)*DSIGA(L+1)*byDS4(L) 334 continue C**** Use quadratic vertical interpolation scheme to calulate quantity C**** at constant pressure coordinate DO 336 K=1,32 IF(KZ(K).le.0) GO TO 340 DO 335 J=1,JM IMAX=IM IF(J.eq.1 .or. J.eq.JM) IMAX=1 DO 335 I=1,IMAX PX = KZ(K)*100. 335 XK(I,J) = GHLtoK (LMA, PX, MS(I,J)*GRAV,MSTRAT*GRAV, SIGEA,byDS1, * ZATMO(I,J), X0(1,I,J),XE(0,I,J), -999999.d0) WRITE (NAME(13)(28:31),922) KZ(K) 336 CALL WRITED (NAME(13),LABEL,OUTYR,OUTMON,'A') C**** 340 continue C**** 500 CLOSE (2) CLOSE (6) GO TO 999 C**** 800 WRITE (0,*) 'Usage: AIJK070 /raid1/C070/DJa* ', * 'Make atmosphere IxJ data files' WRITE (0,*) 'Check the input file: AIJK070.I 98/11/03' 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) 922 FORMAT (I4) 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) 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 PLATS = -TWOPI/4. PSINS = SIN(PLATS) DO 10 J=1,JM-1 PLATN = DLAT*(J+1.-FJEQ) IF(J.eq.JM-1) PLATN = TWOPI/4. PSINN = SIN(PLATN) DXV(J) = DLON*RADIUS*(PSINN-PSINS) / (PLATN-PLATS) PLATS = PLATN 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) DXP(J) =.5*(DXV(J-1)+DXV(J)) DYP(J) = (VLATN-VLATS)*RADIUS 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 AIJK070: inconsistant year and months received' WRITE (0,*) 'From AIJK070: 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 FUNCTION VNLtoK (LM, MK, ME,byDME, X0,XE, BAD) C**** C**** VNLtoK locates MK in a GCM layer. VNLtoK then calculates the C**** quadratic distribution of X in the layer based on the values C**** XE at the layer edges and the mean values X0 of each layer. C**** VNLtoK then evaluates the quadratic polynomial at MK. C**** C**** Input: LM = number of GCM layers C**** MK = mass coordinate to which X is interpolated C**** ME = decreasing mass coordinate of GCM layer edges C**** byDME = reciprocal of mass coordinate layer thicknesses C**** X0 = mean value of X in each layer C**** XE = value of X at each layer edge C**** BAD = value returned by VNLtoK if MK is outside GCM layers C**** IMPLICIT REAL*8 (A-H,M,O-Z) REAL*8 ME(0:LM), byDME(LM), X0(LM), XE(0:LM) C**** IF(MK.gt.ME(0)) GO TO 20 DO 10 L=1,LM IF(MK.lt.ME(L)) GO TO 10 C**** ME(L-1) > MK > ME(L) A = 3.*(XE(L-1)-2.*X0(L)+XE(L))*byDME(L)*byDME(L) B = (XE(L-1)-XE(L))*byDME(L) - A*(ME(L-1)+ME(L)) C = X0(L) - A*(ME(L-1)*ME(L-1)+ME(L-1)*ME(L)+ME(L)*ME(L))/3. - - B*(ME(L-1)+ME(L))*.5 VNLtoK = A*MK*MK + B*MK + C RETURN 10 continue C**** Either MK > ME(0) or ME(LM) > MK: set VNLtoK to BAD value 20 VNLtoK = BAD RETURN END FUNCTION GHLtoK (LM, PX, DP,PT, SIGE,byDSIG, Z0, H0,HE, BAD) C**** C**** GHLtoK locates PX in a GCM layer. GHLtoK then calculates the C**** quadratic distribution of H in the layer based on the values HE C**** at the layer edges and the mean values H0 of each layer. GHLtoK C**** then integrates the quadratic polynomial from the surface to PX. C**** C**** Temp = PSE*P^RKAP/SHCD C**** P = Rho*RGAS*Temp => Rho = P^(1-RKAP)/RKAP*PSE C**** dP/dZ = - Rho*GRAV => dZ = - (RKAP/GRAV)*PSE*P^(RKAP-1)*dP C**** PSE = A*P^2 + B*P + C C**** dZ = - (RKAP/GRAV)*[A*P^(RKAP+1) + B*P^RKAP + C*P^(RKAP-1)]*dP C**** C**** Input: LM = number of GCM layers C**** PX = desired pressure level (Pa) C**** DP = GCM pressure difference (Pa) from surface to PT C**** PT = GCM pressure (Pa) at Sigma = 0 C**** SIGE = GCM Sigma coordinate layer edges C**** byDSIG = reciprocal of sigma coordinate layer thicknesses C**** Z0 = surface topography (m) C**** H0 = mean potential specific enthalpy in each layer C**** HE = potential specific enthalpy at each layer edge C**** BAD = value returned by GHLtoK if PX is outside GCM layers C**** IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (GRAV=9.81d0, RGAS=287., SHCD=1003.5, RKAP=RGAS/SHCD, * byK=1./RKAP, byKp1=1./(RKAP+1.), byKp2=1./(RKAP+2.)) REAL*8 SIGE(0:LM), byDSIG(LM), H0(LM), HE(0:LM) C**** IF(PX.gt.DP+PT .or. PX.lt.PT) GO TO 30 GHLtoK = Z0 byDP = 1./DP PDN = DP + PT PKDN = PDN**RKAP DO 10 L=1,LM PUP = SIGE(L)*DP + PT PKUP = PUP**RKAP byDPL = byDSIG(L)*byDP C**** Calculate quadratic distribution of potential specific enthalpy C**** PSE(Sigma) = A*P^2 + B*P + C A = 3.*(HE(L-1)-2.*H0(L)+HE(L))*byDPL*byDPL B = (HE(L-1)-HE(L))*byDPL - A*(PDN+PUP) C = H0(L) - A*(PDN*PDN+PDN*PUP+PUP*PUP)/3. - B*(PDN+PUP)*.5 IF(PX.ge.PUP) GO TO 20 C**** Integrate height from PDN to PUP GHLtoK = GHLtoK + (RKAP/GRAV)* * (PKDN*(A*PDN*PDN*byKp2 + B*PDN*byKp1 + C*byK) - - PKUP*(A*PUP*PUP*byKp2 + B*PUP*byKp1 + C*byK)) PDN = PUP 10 PKDN = PKUP STOP 'GHLtoK' C**** PDN > PX > PUP: integrate height from PDN to PX 20 GHLtoK = GHLtoK + (RKAP/GRAV)* * (PKDN*(A*PDN*PDN*byKp2 + B*PDN*bYKp1 + C*byK) - - PX**RKAP*(A*PX *PX *byKp2 + B*PX *bYKp1 + C*byK)) RETURN C**** Either PX > DP+PT or PX < PT: set GHLtoK to BAD value 30 GHLtoK = BAD RETURN END SUBROUTINE WRITED (NAME,LABEL,OUTYR,OUTMON,GRID) C**** C**** Write a record to a data file whose format is: TITLE(80),A(IM,JM) C**** Input: NAME = name of quantity C**** LABEL = run number C**** OUTYR = year C**** OUTMON = month C**** GRID = grid arrangement of variables C**** PARAMETER (IM=72,JM=46) REAL*8 Q CHARACTER NAME*50,LABEL*5,OUTYR*4,OUTMON*4,GRID*1, TITLE*80 COMMON Q(IM,JM) IF(GRID.eq.'A') GO TO 10 IF(GRID.eq.'B') GO TO 20 IF(GRID.eq.'U') GO TO 30 IF(GRID.eq.'V') GO TO 20 GO TO 50 C**** A-grid: replicate I=1 value to all longitudes at the poles 10 DO 11 I=2,IM Q(I, 1) = Q(1, 1) 11 Q(I,JM) = Q(1,JM) GO TO 50 C**** B-grid or V-grid: zero out values at North Pole 20 DO 21 I=1,IM 21 Q(I,JM) = 0. GO TO 50 C**** U-grid: zero out values at both poles 30 DO 31 I=1,IM Q(I, 1) = 0. 31 Q(I,JM) = 0. C**** Build TITLE 50 WRITE (TITLE,905) NAME,LABEL,OUTYR,OUTMON C**** Write record to disk WRITE (2) TITLE,(SNGL(Q(I,1)),I=1,IM*JM) WRITE (6,906) TITLE,(Q(IM/2+1,J),J=1,JM) RETURN C**** 905 FORMAT (A50,A5,9X,A4,1X,A3) 906 FORMAT ('0',A72/(12F11.3)) 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 BLOCK DATA C**** C**** Contents of latitude by longitude arrays C**** IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (LMA=9) CHARACTER*50 NAME COMMON /LAYACB/ DSIGA(LMA),SIGEA(0:LMA) COMMON /NAMECB/ NAME(14) C**** DATA SIGEA /1.,.948665d0,.866530d0,.728953d0,.554415d0, * .390144d0,.251540d0,.143737d0,.061602d0,0./ C**** DATA NAME / 1' ', 2'EASTWARD VELOCITY (m/s) at 1000 mb', 3'NORTHWARD VELOCITY (m/s) at 1000 mb', 4'UPWARD MASS FLUX (Pa/s) at 1000 mb', 5'TEMPERATURE (K) at 1000 mb', 6'CONVERGED HEAT TRANSPORT (K/s) at 1000 mb', 7'NORTHWARD HEAT FLUX (K*m/s) at 1000 mb', 8' ', 9'SPECIFIC HUMIDITY (g/kg) at 1000 mb', 1O'CONVERGED MOISTURE TRANSPORT (1/s) at 1000 mb', 11'NORTHWARD MOISTURE FLUX (m/s) at 1000 mb', 12'UPWARD MOISTURE FLUX (Pa/s) at 1000 mb', 13'GEOPOTENTIAL HEIGHT (m) at 1000 mb', 14'RELATIVE HUMIDITY (%) at 1000 mb'/ END