C CALC_SOUND.F C C DESCRIPTION: SUBROUTINE COMPUTES CAPE, CIN, DCAPE AND RETURNS C VALUES C C INPUT: C NA - number of vertical levels C T - temperature [degrees C] C Q - mixing ratio [g/kg] C P - pressure [hPa] C C *** PROGRAM ADAPTED FROM KERRY EMANUEL'S PROGRAM calcsound.f *** C *** AVAILABLE FROM ftp://texmex.mit.edu/pub/emanuel/BOOK/ *** C *** ADAPTED TO WORK WITH NCL BY ANDREW PENNY - OCT. 27, 2011 *** C======================================= C *** COMMENTS FROM ORIGINAL PROGRAM (RELATES TO ORIGINAL UNMODIFIED VERSION) *** C This program reads the input file "sounding", which can be created using C the program "readsound.f", and calculates the temperatures and virtual C temperatures of parcels lifted reversibly and also pseudo-adiabatically C from each level (but without ice processes). It also calculates the C temperatures and virtual temperatures of parcels that are first saturated by C a wet bulb process and then descend pseudo-adiabatically to the surface; C this is done for each level in the sounding. The Convective Available C Potential Energy (CAPE), Positive Area (PA) and Negative Area (NA) are C calculated for parcel ascent from each level, for both reversible and C pseudo-adiabatic processes, and the Downdraft Convective Available C Potential Energy (DCAPE) is calculated for the downdraft process described C above. C Output is directed to the output file "calcsound.out" and to an NCAR graphics C file. The output file "calcsound.out" contains the values of CAPE, PA, NA C and DCAPE, while the NCAR graphics file contains two contour plots showing, C the difference between the virtual temperatures (C) of reversibly and C pseudo-adiabatically lifted parcels, respectively, as a function of the C pressure level to which the parcel is lifted (ordinate) and the parcel C origin pressure (abscissa). When the pressure is greater than the origin C pressure (on the second graph only), the quantity graphed is the buoyancy C (degrees C) of air that has first been saturated by the wet bulb process at C its origin level, and then descends pseudo-adiabatically. C If the user does not have NCAR graphics, other contouring routines can be C substituted for the NCAR graphics calls toward the end of the program. C This program needs only to be compiled and run; no input values other than C the input sounding need be supplied. C======================================= C NCLFORTSTART SUBROUTINE CALCSOUND(NA,T,P,R,PAR,PAP,NAR,NAP,CAPER,CAPEP,DCAPE) IMPLICIT NONE INTEGER NA REAL T(NA), P(NA), R(NA), TLR(NA,NA), TLP(NA,NA) REAL EV(NA), ES(NA), TVRDIF(NA,NA), TVPDIF(NA,NA) REAL TLVR(NA,NA), TLVP(NA,NA), LW(NA,NA), TVD(NA) REAL CAPER(NA), CAPEP(NA), DCAPE(NA), PAR(NA) REAL PAP(NA), NAR(NA), NAP(NA) C NCLEND C C *** DECLARE LOCAL VARIABLES C INTEGER I, J, K, ICB, INBR, INBP, N, IMAX REAL AH, AHG, ALV, ALV0, ALV1, CHI, CL, CPD REAL CPV, CPVMCL, CPW, EG, EM, ENEW, EPS, PLCL, PM REAL RD, RGD0, RG, RG0, RH, RS, RV, S, SG, SL, SLOPE REAL SLP, SP, SPD, SPG, SUM, SUM2, TC, TG, TG0, TGD0 REAL TVDIFM, TVM N = NA C C *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *** C CPD=1005.7 CPV=1870.0 CL=4190.0 CPVMCL=2320.0 RV=461.5 RD=287.04 EPS=RD/RV ALV0=2.501E6 DO 20 I=1,N R(I)=R(I)*0.001 EV(I)=R(I)*P(I)/(EPS+R(I)) ES(I)=6.112*EXP(17.67*T(I)/(243.5+T(I))) T(I)=T(I)+273.15 20 CONTINUE C C *** Begin outer loop, which cycles through parcel origin levels I *** C DO 500 I=1,N C C *** Define various conserved parcel quantities: reversible *** C *** entropy, S, pseudo-adiabatic entropy, SP, *** C *** and enthalpy, AH *** C 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 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,N CAPER(I)=0.0 CAPEP(I)=0.0 DCAPE(I)=0.0 PAP(I)=0.0 PAR(I)=0.0 NAP(I)=0.0 NAR(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 IMAX=MAX(INBR+1,I) DO 555 J=IMAX,N TVRDIF(I,J)=0.0 555 CONTINUE IMAX=MAX(INBP+1,I) DO 565 J=IMAX,N TVPDIF(I,J)=0.0 565 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(I)=NAR(I)-RD*TVM*(P(J-1)-P(J))/PM ELSE PAR(I)=PAR(I)+RD*TVM*(P(J-1)-P(J))/PM END IF 600 CONTINUE CAPER(I)=PAR(I)-NAR(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(I)=NAP(I)-RD*TVM*(P(J-1)-P(J))/PM ELSE PAP(I)=PAP(I)+RD*TVM*(P(J-1)-P(J))/PM END IF 650 CONTINUE CAPEP(I)=PAP(I)-NAP(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(I)=DCAPE(I)-RD*TVM*(P(J)-P(J+1))/PM END IF 700 CONTINUE 800 CONTINUE RETURN END