PROGRAM SOUND C C *** This program accepts input from the sounding contained *** C *** in the file <> and calculates various *** C *** properties of samples of air raised or lowered *** C *** to different levels. *** C PARAMETER (NA=800,NK=20,NS=1000, NAT=900) C REAL T(NA), P(NA), R(NA), TLR(NK,NA), TLP(NK,NA) REAL PTEM(NAT), RTEM(NAT), TTEM(NAT),EVTEM(NAT),ESTEM(NAT) REAL EV(NA), ES(NA), TVRDIF(NK,NA), TVPDIF(NK,NA) REAL TLVR(NK,NA), TLVP(NK,NA), LW(NK,NA), TVD(NA) REAL CAPER(NS,NK), CAPEP(NS,NK), DCAPE(NS,NK), PAR(NS,NK) REAL CAPERAV(NK), CAPEPAV(NK), STDCAPEP(NK),STDCAPER(NK) REAL PAP(NS,NK), NAR(NS,NK), NAP(NS,NK), PL(NK) REAL ATVRD(NS,NK,NA),ATVPD(NS,NK,NA) REAL TRDBAR(NK,NA),TPDBAR(NK,NA),STDR(NK,NA),STDP(NK,NA) INTEGER JJIR(NK,NA),JJIP(NK,NA) CHARACTER*12 AB CHARACTER*16 AA C CHARACTER NAME*15, OBS*15, PREP*2, OTIME*3, MONTH*3, CNAM*4 CHARACTER*5 CNAM CHARACTER*15 DUMMY1, DUMMY2 CHARACTER*11 DATETIME C C *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *** C CPD=1005.7 CPV=1870.0 CL=2500.0 CPVMCL=2320.0 RV=461.5 RD=287.04 EPS=RD/RV ALV0=2.501E6 C P1BAR=0.0 NPC=0 NMAX=0 DO 8 I=1,NK DO 7 J=1,NA JJIR(I,J)=0 JJIP(I,J)=0 TRDBAR(I,J)=0.0 TPDBAR(I,J)=0.0 STDR(I,J)=0.0 STDP(I,J)=0.0 DO 6 K=1,NS ATVRD(K,I,J)=0.0 ATVPD(K,I,J)=0.0 6 CONTINUE 7 CONTINUE 8 CONTINUE NSTOT=0 C DO 900 JJ=1,1 C C C *** Read in the sounding from the file <> *** c OPEN(UNIT=11,FILE='sounding.txt',STATUS='OLD') READ(11,*) READ(11,*) READ(11,*) READ(11,*) READ(11,*) C KK=880 I=0 OPEN(UNIT=10,FILE='modsound.txt',STATUS='UNKNOWN') DO 20 KI=1,KK READ(11,15,ERR=21,END=21)PK,TK,RK 15 FORMAT(1X,F6.1,9X,F5.1,2X,F5.1) IF(PK.LT.78.0.OR.TK.GT.100.0.OR.RK.GT.100.0)GOTO 20 I=I+1 PTEM(I)=PK TTEM(I)=TK RTEM(I)=RK IF(TTEM(I).LT.-200.0)TTEM(I)=-200.0 EVTEM(I)=6.112*EXP(17.67*RTEM(I)/(243.5+RTEM(I))) ESTEM(I)=6.112*EXP(17.67*TTEM(I)/(243.5+TTEM(I))) EVTEM(I)=MIN(EVTEM(I),ESTEM(I)) RTEM(I)=0.622*EVTEM(I)/(PTEM(I)-EVTEM(I)) TTEM(I)=TTEM(I)+273.15 WRITE(10,*)PK,TK,EVTEM(I)/ESTEM(I) 20 CONTINUE 21 CONTINUE READ(11,*)DUMMY1,DUMMY2,CNAM if(dummy2.eq.'identifier: ')then READ(11,*)DUMMY1,DUMMY2,ID READ(11,*)DUMMY1,DUMMY2,idate,itime else READ(11,*)DUMMY1,DUMMY2,idate,itime end if N=I CLOSE(UNIT=11) CLOSE(10) C C *** Write header *** C IYEAR=IDATE/10000 IY1=IYEAR IF(IYEAR.LT.50)THEN IYEAR=IYEAR+2000 ELSE IYEAR=IYEAR+1900 END IF IMONTH=IDATE/100-100*IY1 IDAY=IDATE-10000*IY1-100*IMONTH OPEN(UNIT=12,FILE='header.txt',STATUS='UNKNOWN') WRITE(12,*)CNAM,' ',ITIME,' ',IDAY,' ',IMONTH,' ',IYEAR CLOSE(12) C c IF(PTEM(1).LT.995.0.OR.TTEM(1).LT.273.0.OR.TTEM(1).GT.306.0.OR. c 1 RTEM(1).LT.0.002.OR.PTEM(N).GE.500.0.OR.PTEM(1).GT.1050.0) c 2 GOTO 900 NSTOT=NSTOT+1 C C *** Interpolate soundings to 5 mb intervals C P(1)=PTEM(1) T(1)=TTEM(1) R(1)=RTEM(1) EV(1)=EVTEM(1) ES(1)=ESTEM(1) NL=(1005.-PTEM(N))/5.0+0.001 DO I=2,NL P(I)=1005.0-5.0*FLOAT(I) DO J=2,N IF(PTEM(J).LE.P(I))THEN T(I)=(TTEM(J)*(PTEM(J-1)-P(I))+TTEM(J-1)*(P(I)-PTEM(J)))/ 1 (PTEM(J-1)-PTEM(J)) R(I)=(RTEM(J)*(PTEM(J-1)-P(I))+RTEM(J-1)*(P(I)-PTEM(J)))/ 1 (PTEM(J-1)-PTEM(J)) EV(I)=(EVTEM(J)*(PTEM(J-1)-P(I))+EVTEM(J-1)*(P(I)-PTEM(J)))/ 1 (PTEM(J-1)-PTEM(J)) ES(I)=(ESTEM(J)*(PTEM(J-1)-P(I))+ESTEM(J-1)*(P(I)-PTEM(J)))/ 1 (PTEM(J-1)-PTEM(J)) GOTO 23 END IF END DO 23 CONTINUE END DO c print*, p(nl),t(nl),r(nl),p(1),t(1),r(1) N=NL NMAX=MAX(N,NMAX) C C *** Begin outer loop, which cycles through parcel origin levels I *** C DO 500 I=1,NK C C *** Define various conserved parcel quantities: reversible *** C *** entropy, S, pseudo-adiabatic entropy, SP, *** C *** and enthalpy, AH *** C PL(I)=P(I) RS=EPS*ES(I)/(P(I)-ES(I)) ALV=ALV0-CPVMCL*(T(I)-273.15) EM=MAX(EV(I),1.0E-6) S=(CPD+R(I)*CL)*LOG(T(I))-RD*LOG(P(I)-EV(I))+ 1 ALV*R(I)/T(I)-R(I)*RV*LOG(EM/ES(I)) SP=CPD*LOG(T(I))-RD*LOG(P(I)-EV(I))+ 1 ALV*R(I)/T(I)-R(I)*RV*LOG(EM/ES(I)) AH=(CPD+R(I)*CL)*T(I)+ALV*R(I) C C *** Find the temperature and mixing ratio of the parcel at *** C *** level I saturated by a wet bulb process *** C SLOPE=CPD+ALV*ALV*RS/(RV*T(I)*T(I)) TG=T(I) RG=RS DO 100 J=1,20 ALV1=ALV0-CPVMCL*(TG-273.15) AHG=(CPD+CL*RG)*TG+ALV1*RG TG=TG+(AH-AHG)/SLOPE TC=TG-273.15 ENEW=6.112*EXP(17.67*TC/(243.5+TC)) RG=EPS*ENEW/(P(I)-ENEW) 100 CONTINUE C C *** Calculate conserved variable at top of downdraft *** C EG=RG*P(I)/(EPS+RG) SPD=CPD*LOG(TG)-RD*LOG(P(I)-EG)+ 1 ALV1*RG/TG TVD(I)=TG*(1.+RG/EPS)/(1.+RG)-T(I)*(1.+R(I)/EPS)/ 1 (1.+R(I)) IF(P(I).LT.100.0)TVD(I)=0.0 RGD0=RG TGD0=TG C C *** Find lifted condensation pressure *** C RH=R(I)/RS RH=MIN(RH,1.0) CHI=T(I)/(1669.0-122.0*RH-T(I)) PLCL=1.0 IF(RH.GT.0.0)THEN PLCL=P(I)*(RH**CHI) END IF C C *** Begin updraft loop *** C SUM=0.0 RG0=R(I) TG0=T(I) DO 200 J=I,N C C *** Calculate estimates of the rates of change of the entropies *** C *** with temperature at constant pressure *** C RS=EPS*ES(J)/(P(J)-ES(J)) ALV=ALV0-CPVMCL*(T(J)-273.15) SL=(CPD+R(I)*CL+ALV*ALV*RS/(RV*T(J)*T(J)))/T(J) SLP=(CPD+RS*CL+ALV*ALV*RS/(RV*T(J)*T(J)))/T(J) C C *** Calculate lifted parcel temperature below its LCL *** C IF(P(J).GE.PLCL)THEN TLR(I,J)=T(I)*(P(J)/P(I))**(RD/CPD) TLP(I,J)=TLR(I,J) LW(I,J)=0.0 TLVR(I,J)=TLR(I,J)*(1.+R(I)/EPS)/(1.+R(I)) TLVP(I,J)=TLVR(I,J) TVRDIF(I,J)=TLVR(I,J)-T(J)*(1.+R(J)/EPS)/(1.+R(J)) TVPDIF(I,J)=TVRDIF(I,J) ELSE C C *** Iteratively calculate lifted parcel temperature and mixing *** C *** ratios for both reversible and pseudo-adiabatic ascent *** C TG=T(J) RG=RS DO 150 K=1,20 EM=RG*P(J)/(EPS+RG) ALV=ALV0-CPVMCL*(TG-273.15) SG=(CPD+R(I)*CL)*LOG(TG)-RD*LOG(P(J)-EM)+ 1 ALV*RG/TG TG=TG+(S-SG)/SL TC=TG-273.15 ENEW=6.112*EXP(17.67*TC/(243.5+TC)) RG=EPS*ENEW/(P(J)-ENEW) 150 CONTINUE TLR(I,J)=TG TLVR(I,J)=TG*(1.+RG/EPS)/(1.+R(I)) LW(I,J)=R(I)-RG LW(I,J)=MAX(0.0,LW(I,J)) TVRDIF(I,J)=TLVR(I,J)-T(J)*(1.+R(J)/EPS)/(1.+R(J)) c IF(TVRDIF(I,J).GE.7.0)GOTO 900 C C *** Now do pseudo-adiabatic ascent *** C TG=T(J) RG=RS DO 180 K=1,20 CPW=0.0 IF(J.GT.1)THEN CPW=SUM+CL*0.5*(RG0+RG)*(LOG(TG)-LOG(TG0)) END IF EM=RG*P(J)/(EPS+RG) ALV=ALV0-CPVMCL*(TG-273.15) SPG=CPD*LOG(TG)-RD*LOG(P(J)-EM)+CPW+ 1 ALV*RG/TG TG=TG+(SP-SPG)/SLP TC=TG-273.15 ENEW=6.112*EXP(17.67*TC/(243.5+TC)) RG=EPS*ENEW/(P(J)-ENEW) 180 CONTINUE TLP(I,J)=TG TLVP(I,J)=TG*(1.+RG/EPS)/(1.+RG) TVPDIF(I,J)=TLVP(I,J)-T(J)*(1.+R(J)/EPS)/(1.+R(J)) RG0=RG TG0=TG SUM=CPW END IF 200 CONTINUE IF(I.EQ.1)GOTO 500 C C *** Begin downdraft loop *** C SUM2=0.0 DO 300 J=I-1,1,-1 C C *** Calculate estimate of the rate of change of entropy *** C *** with temperature at constant pressure *** C RS=EPS*ES(J)/(P(J)-ES(J)) ALV=ALV0-CPVMCL*(T(J)-273.15) SLP=(CPD+RS*CL+ALV*ALV*RS/(RV*T(J)*T(J)))/T(J) TG=T(J) RG=RS C C *** Do iteration to find downdraft temperature *** C DO 250 K=1,20 CPW=SUM2+CL*0.5*(RGD0+RG)*(LOG(TG)-LOG(TGD0)) EM=RG*P(J)/(EPS+RG) ALV=ALV0-CPVMCL*(TG-273.15) SPG=CPD*LOG(TG)-RD*LOG(P(J)-EM)+CPW+ 1 ALV*RG/TG TG=TG+(SPD-SPG)/SLP TC=TG-273.15 ENEW=6.112*EXP(17.67*TC/(243.5+TC)) RG=EPS*ENEW/(P(J)-ENEW) 250 CONTINUE SUM2=CPW TGD0=TG RGD0=RG TLP(I,J)=TG TLVP(I,J)=TG*(1.+RG/EPS)/(1.+RG) TVPDIF(I,J)=TLVP(I,J)-T(J)*(1.+R(J)/EPS)/(1.+R(J)) IF(P(I).LT.100.0)TVPDIF(I,J)=0.0 TVPDIF(I,J)=MIN(TVPDIF(I,J),0.0) TLR(I,J)=T(J) TLVR(I,J)=T(J) TVRDIF(I,J)=0.0 LW(I,J)=0.0 300 CONTINUE 500 CONTINUE C C *** Begin loop to find CAPE, PA, and NA from reversible and *** C *** pseudo-adiabatic ascent, and DCAPE *** C DO 800 I=1,NK CAPER(JJ,I)=0.0 CAPEP(JJ,I)=0.0 DCAPE(JJ,I)=0.0 PAP(JJ,I)=0.0 PAR(JJ,I)=0.0 NAP(JJ,I)=0.0 NAR(JJ,I)=0.0 C C *** Find lifted condensation pressure *** C RS=EPS*ES(I)/(P(I)-ES(I)) RH=R(I)/RS RH=MIN(RH,1.0) CHI=T(I)/(1669.0-122.0*RH-T(I)) PLCL=1.0 IF(RH.GT.0.0)THEN PLCL=P(I)*(RH**CHI) END IF C C *** Find lifted condensation level and maximum level *** C *** of positive buoyancy *** C ICB=N INBR=1 INBP=1 DO 550 J=N,I,-1 IF(P(J).LT.PLCL)ICB=MIN(ICB,J) IF(TVRDIF(I,J).GT.0.0)INBR=MAX(INBR,J) IF(TVPDIF(I,J).GT.0.0)INBP=MAX(INBP,J) 550 CONTINUE C C *** Do updraft loops *** C IF(INBR.GT.I)THEN DO 600 J=I+1,INBR TVM=0.5*(TVRDIF(I,J)+TVRDIF(I,J-1)) PM=0.5*(P(J)+P(J-1)) IF(TVM.LE.0.0)THEN NAR(JJ,I)=NAR(JJ,I)-RD*TVM*(P(J-1)-P(J))/PM ELSE PAR(JJ,I)=PAR(JJ,I)+RD*TVM*(P(J-1)-P(J))/PM END IF 600 CONTINUE CAPER(JJ,I)=PAR(JJ,I)-NAR(JJ,I) END IF IF(INBP.GT.I)THEN DO 650 J=I+1,INBP TVM=0.5*(TVPDIF(I,J)+TVPDIF(I,J-1)) PM=0.5*(P(J)+P(J-1)) IF(TVM.LE.0.0)THEN NAP(JJ,I)=NAP(JJ,I)-RD*TVM*(P(J-1)-P(J))/PM ELSE PAP(JJ,I)=PAP(JJ,I)+RD*TVM*(P(J-1)-P(J))/PM END IF 650 CONTINUE CAPEP(JJ,I)=PAP(JJ,I)-NAP(JJ,I) END IF C C *** Find DCAPE *** C IF(I.EQ.1)GOTO 800 DO 700 J=I-1,1,-1 TVDIFM=TVPDIF(I,J+1) IF(I.EQ.(J+1))TVDIFM=TVD(I) TVM=0.5*(TVPDIF(I,J)+TVDIFM) PM=0.5*(P(J)+P(J+1)) IF(TVM.LT.0.0)THEN DCAPE(JJ,I)=DCAPE(JJ,I)-RD*TVM*(P(J)-P(J+1))/PM END IF 700 CONTINUE 800 CONTINUE C DO 520 I=1,NK INBR=1 INBP=1 DO 505 J=N,I,-1 IF(TVRDIF(I,J).GT.0.0)INBR=MAX(INBR,J) IF(TVPDIF(I,J).GT.0.0)INBP=MAX(INBP,J) 505 CONTINUE IMAX=MAX(INBR,I) IMAX=N IF(IMAX.GT.I)THEN DO 508 J=I,IMAX IF(TVRDIF(I,J).NE.0.0)THEN JJIR(I,J)=JJIR(I,J)+1 ATVRD(JJIR(I,J),I,J)=TVRDIF(I,J) END IF 508 CONTINUE END IF IMAX=MAX(I,INBP) IMAX=N DO 510 J=I,IMAX IF(TVPDIF(I,J).NE.0.0)THEN JJIP(I,J)=JJIP(I,J)+1 ATVPD(JJIP(I,J),I,J)=TVPDIF(I,J) END IF 510 CONTINUE 520 CONTINUE NPC=NPC+1 P1BAR=P1BAR+P(1) 900 CONTINUE 901 CONTINUE C CLOSE(UNIT=15) C P1BAR=P1BAR/FLOAT(NPC) P(1)=P1BAR PL(1)=P1BAR DO 920 I=1,NK PARS=0.0 PAPS=0.0 ANAPS=0.0 ANARS=0.0 CAPERS=0.0 CAPEPS=0.0 DCAPES=0.0 STDCAPEP(I)=0.0 STDCAPER(I)=0.0 JPAR=0 JPAP=0 JNAR=0 JNAP=0 JCAPER=0 JCAPEP=0 JDCAPE=0 DO 903 JJ=1,NSTOT IF(PAP(JJ,I).NE.0.0)THEN PAPS=PAPS+PAP(JJ,I) JPAP=JPAP+1 END IF IF(PAR(JJ,I).NE.0.0)THEN PARS=PARS+PAR(JJ,I) JPAR=JPAR+1 END IF IF(NAP(JJ,I).NE.0.0)THEN ANAPS=ANAPS+NAP(JJ,I) JNAP=JNAP+1 END IF IF(NAR(JJ,I).NE.0.0)THEN ANARS=ANARS+NAR(JJ,I) JNAR=JNAR+1 END IF IF(CAPER(JJ,I).NE.0.0)THEN CAPERS=CAPERS+CAPER(JJ,I) JCAPER=JCAPER+1 END IF IF(CAPEP(JJ,I).NE.0.0)THEN CAPEPS=CAPEPS+CAPEP(JJ,I) JCAPEP=JCAPEP+1 END IF IF(DCAPE(JJ,I).NE.0.0)THEN DCAPES=DCAPES+DCAPE(JJ,I) JDCAPE=JDCAPE+1 END IF 903 CONTINUE JNAP=MAX(JNAP,1) JNAR=MAX(JNAR,1) JPAP=MAX(JPAP,1) JPAR=MAX(JPAR,1) JCAPER=MAX(JCAPER,1) JCAPEP=MAX(JCAPEP,1) JDCAPE=MAX(JDCAPE,1) PAP(1,I)=PAPS/FLOAT(JPAP) PAR(1,I)=PARS/FLOAT(JPAR) NAP(1,I)=ANAPS/FLOAT(JNAP) NAR(1,I)=ANARS/FLOAT(JNAR) CAPEPAV(I)=CAPEPS/FLOAT(JCAPEP) CAPERAV(I)=CAPERS/FLOAT(JCAPER) DCAPE(1,I)=DCAPES/FLOAT(JDCAPE) DO 1010 JJ=1,JCAPEP IF(CAPEP(JJ,I).NE.0.0)THEN STDCAPEP(I)=STDCAPEP(I)+(CAPEP(JJ,I)-CAPEPAV(I))**2 END IF 1010 CONTINUE STDCAPEP(I)=SQRT(STDCAPEP(I)/FLOAT(JCAPEP)) DO 1020 JJ=1,JCAPER IF(CAPER(JJ,I).NE.0.0)THEN STDCAPER(I)=STDCAPER(I)+(CAPER(JJ,I)-CAPERAV(I))**2 END IF 1020 CONTINUE STDCAPER(I)=SQRT(STDCAPER(I)/FLOAT(JCAPER)) DO 915 J=1,NMAX IF(JJIR(I,J).GT.0)THEN DO 905 JJ=1,JJIR(I,J) TRDBAR(I,J)=TRDBAR(I,J)+ATVRD(JJ,I,J) 905 CONTINUE TRDBAR(I,J)=TRDBAR(I,J)/FLOAT(JJIR(I,J)) TRDBAR(I,J)=MAX(TRDBAR(I,J),-4.0) DO 907 JJ=1,JJIR(I,J) STDR(I,J)=STDR(I,J)+(ATVRD(JJ,I,J)-TRDBAR(I,J))**2 907 CONTINUE STDR(I,J)=SQRT(STDR(I,J)/FLOAT(JJIR(I,J))) END IF IF(JJIP(I,J).GT.0)THEN DO 908 JJ=1,JJIP(I,J) TPDBAR(I,J)=TPDBAR(I,J)+ATVPD(JJ,I,J) 908 CONTINUE TPDBAR(I,J)=TPDBAR(I,J)/FLOAT(JJIP(I,J)) TPDBAR(I,J)=MAX(TPDBAR(I,J),-4.0) DO 909 JJ=1,JJIP(I,J) STDP(I,J)=STDP(I,J)+(ATVPD(JJ,I,J)-TPDBAR(I,J))**2 909 CONTINUE STDP(I,J)=SQRT(STDP(I,J)/FLOAT(JJIP(I,J))) END IF 915 CONTINUE 920 CONTINUE C C *** Write values of PA, NA, CAPE and DCAPE *** C OPEN(UNIT=12, FILE='cape.out', STATUS='unknown') C WRITE(12,921)NPC 921 FORMAT(5X,'Number of soundings = ',I4) WRITE(12,801) 801 FORMAT(20X,'ALL AREAS IN UNITS OF J/kg',//) WRITE(12,802) 802 FORMAT(1X,'Origin',4X,'Rev.',4X,'P.A.',4X,'Rev.',4X,'P.A.', 1 4X,'Rev.',4X,'P.A.',4x,'Rev.',4x,'P.A.') WRITE(12,804) 804 FORMAT(1X,'p (mb)',4X,' PA ',4X,' PA ',4X,' NA ',4X,' NA ', 1 4X,'CAPE',4X,'CAPE',4X,'DCAPE',1X,'STDCAPE',4X,'STDCAPE') WRITE(12,806) 806 FORMAT(1X,'------',4X,'----',4X,'----',4X,'----',4X,'----', 1 4X,'----',4X,'----',4X,'-----',4X,'-------', 2 4X,'-------') DO 820 I=1,NK WRITE(12,810)PL(I),PAR(1,I),PAP(1,I),NAR(1,I), 1 NAP(1,I),CAPERAV(I),CAPEPAV(I),DCAPE(1,I),STDCAPEP(I), 2 STDCAPER(I) 810 FORMAT(1X,F6.1,9(F8.1)) 820 CONTINUE CLOSE(12) C C *** Output files for contour plots (e.g. for MATLAB) *** C OPEN(UNIT=12,FILE='p.out',STATUS='UNKNOWN') DO I=1,NMAX WRITE(12,*)P(I) END DO CLOSE(12) C OPEN(UNIT=12,FILE='porig.out',STATUS='UNKNOWN') DO I=1,NK WRITE(12,*)PL(I) END DO CLOSE(12) C OPEN(UNIT=12,FILE='tdifrev.out',STATUS='UNKNOWN') DO I=1,NMAX WRITE(12,1000)(TRDBAR(J,I),J=1,NK) END DO CLOSE(12) C OPEN(UNIT=12,FILE='tstdrev.out',STATUS='UNKNOWN') DO I=1,NMAX WRITE(12,1000)(STDR(J,I),J=1,NK) END DO CLOSE(12) C OPEN(UNIT=12,FILE='tdifpseudo.out',STATUS='UNKNOWN') DO I=1,NMAX WRITE(12,1000)(TPDBAR(J,I),J=1,NK) END DO CLOSE(12) C OPEN(UNIT=12,FILE='tstdpseudo.out',STATUS='UNKNOWN') DO I=1,NMAX WRITE(12,1000)(STDP(J,I),J=1,NK) END DO CLOSE(12) C 1000 FORMAT(21(1X,F7.3)) C STOP END