SUBROUTINE CALCAPE(ITYPE,P1D,T1D,Q1D,L1D,CAPE,CINS) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . C SUBPROGRAM: CALCAPE COMPUTES CAPE AND CINS C PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-10 C C ABSTRACT: C C THIS ROUTINE COMPUTES CAPE AND CINS GIVEN TEMPERATURE, C PRESSURE, AND SPECIFIC HUMIDTY. IN "STORM AND CLOUD C DYNAMICS" (1989, ACADEMIC PRESS) COTTON AND ANTHES DEFINE C CAPE (EQUATION 9.16, P501) AS C C EL C CAPE = SUM G * LN(THETAP/THETAA) DZ C LCL C C WHERE, C EL = EQUILIBRIUM LEVEL, C LCL = LIFTING CONDENSTATION LEVEL, C G = GRAVITATIONAL ACCELERATION, C THETAP = LIFTED PARCEL POTENTIAL TEMPERATURE, C THETAA = AMBIENT POTENTIAL TEMPERATURE. C C NOTE THAT THE INTEGRAND LN(THETAP/THETAA) APPROXIMATELY C EQUALS (THETAP-THETAA)/THETAA. THIS RATIO IS OFTEN USED C IN THE DEFINITION OF CAPE/CINS. C C TWO TYPES OF CAPE/CINS CAN BE COMPUTED BY THIS ROUTINE. THE C SUMMATION PROCESS IS THE SAME FOR BOTH CASES. WHAT DIFFERS C IS THE DEFINITION OF THE PARCEL TO LIFT. FOR ITYPE=1 THE C PARCEL WITH THE WARMEST THETA-E IN A DPBND PASCAL LAYER ABOVE C THE MODEL SURFACE IS LIFTED. THE ARRAYS P1D, T1D, Q1D, AND C L1D ARE NOT USED. FOR ITYPE=2 THE ARRAYS P1D, T1D, Q1D, AND C L1D DEFINE THE PARCEL TO LIFT IN EACH COLUMN. BOTH TYPES OF C CAPE/CINS MAY BE COMPUTED IN A SINGLE EXECUTION OF THE POST C PROCESSOR. C C THIS ALGORITHM PROCEEDS AS FOLLOWS. C FOR EACH COLUMN, C (1) INITIALIZE RUNNING CAPE AND CINS SUM TO 0.0 C (2) COMPUTE TEMPERATURE AND PRESSURE AT THE LCL USING C LOOK UP TABLE (PTBL). USE EITHER PARCEL THAT GIVES C MAX THETAE IN LOWEST DPBND ABOVE GROUND (ITYPE=1) C OR GIVEN PARCEL FROM T1D,Q1D,...(ITYPE=2). C (3) COMPUTE THE TEMP OF A PARCEL LIFTED FROM THE LCL. C WE KNOW THAT THE PARCEL'S C EQUIVALENT POTENTIAL TEMPERATURE (THESP) REMAINS C CONSTANT THROUGH THIS PROCESS. WE CAN C COMPUTE TPAR USING THIS KNOWLEDGE USING LOOK C UP TABLE (SUBROUTINE TTBLEX). C (4) FIND THE EQUILIBRIUM LEVEL. THIS IS DEFINED AS THE C HIGHEST POSITIVELY BUOYANT LAYER. C (IF THERE IS NO POSITIVELY BUOYANT LAYER, CAPE/CINS C WILL BE ZERO) C (5) COMPUTE CAPE/CINS. C (A) COMPUTE THETAP. WE KNOW TPAR AND P. C (B) COMPUTE THETAA. WE KNOW T AND P. C (6) ADD G*(THETAP-THETAA)*DZ TO THE RUNNING CAPE OR CINS SUM. C (A) IF THETAP > THETAA, ADD TO THE CAPE SUM. C (B) IF THETAP < THETAA, ADD TO THE CINS SUM. C (7) ARE WE AT EQUILIBRIUM LEVEL? C (A) IF YES, STOP THE SUMMATION. C (B) IF NO, CONTIUNUE THE SUMMATION. C (8) ENFORCE LIMITS ON CAPE AND CINS (I.E. NO NEGATIVE CAPE) C C PROGRAM HISTORY LOG: C 93-02-10 RUSS TREADON C 93-06-19 RUSS TREADON - GENERALIZED ROUTINE TO ALLOW FOR C TYPE 2 CAPE/CINS CALCULATIONS. C 94-09-23 MIKE BALDWIN - MODIFIED TO USE LOOK UP TABLES C INSTEAD OF COMPLICATED EQUATIONS. C 94-10-13 MIKE BALDWIN - MODIFIED TO CONTINUE CAPE/CINS CALC C UP TO AT HIGHEST BUOYANT LAYER. C 98-06-12 T BLACK - CONVERSION FROM 1-D TO 2-D C 98-08-18 T BLACK - COMPUTE APE INTERNALLY C 00-01-04 JIM TUCCILLO - MPI VERSION C C USAGE: CALL CALCAPE(ITYPE,P1D,T1D,Q1D,L1D,CAPE,CINS) C INPUT ARGUMENT LIST: C ITYPE - INTEGER FLAG SPECIFYING HOW PARCEL TO LIFT IS C IDENTIFIED. SEE COMMENTS ABOVE. C P1D - ARRAY OF PRESSURE OF PARCELS TO LIFT. C T1D - ARRAY OF TEMPERATURE OF PARCELS TO LIFT. C Q1D - ARRAY OF SPECIFIC HUMIDITY OF PARCELS TO LIFT. C L1D - ARRAY OF MODEL LEVEL OF PARCELS TO LIFT. C C OUTPUT ARGUMENT LIST: C CAPE - CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG) C CINS - CONVECTIVE INHIBITION (J/KG) C C OUTPUT FILES: C STDOUT - RUN TIME STANDARD OUT. C C SUBPROGRAMS CALLED: C UTILITIES: C P2FILT - SMOOTH DATA ON STAGGERED E-GRID. C BOUND - BOUND (CLIP) DATA BETWEEN UPPER AND LOWER LIMTS. C TTBLEX - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P C C LIBRARY: C COMMON - VRBLS C PHYS C EXTRA C IOUNIT C LOOPS C MASKS C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE : CRAY C-90 C$$$ C C C C INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980). INCLUDE "parmeta" INCLUDE "parm.tbl" INCLUDE "params" PARAMETER (DPBND=70.E2, & ELIVW=2.72E6,ELOCP=ELIVW/CP) PARAMETER (ISMTHP=2,ISMTHT=2,ISMTHQ=2) C C DECLARE VARIABLES. C INTEGER L1D(IM,JM),IEQL(IM,JM),IPTB(IM,JM),ITHTB(IM,JM) INTEGER KLRES(IM,JM),KHRES(IM,JM),LCL(IM,JM) C INTEGER KNUML(LM),KNUMH(LM) C REAL P1D(IM,JM),T1D(IM,JM),Q1D(IM,JM) REAL CAPE(IM,JM),CINS(IM,JM) REAL THESP(IM,JM),TPAR(IM,JM,LM),PSP(IM,JM) REAL QQ(IM,JM),PP(IM,JM) INTEGER IDX(IM,JM) C C INCLUDE COMMON BLOCKS. INCLUDE "VRBLS.comm" INCLUDE "EXTRA.comm" INCLUDE "IOUNIT.comm" INCLUDE "LOOPS.comm" INCLUDE "PHYS.comm" INCLUDE "MASKS.comm" INCLUDE "CTLBLK.comm" C C************************************************************** C START CALCAPE HERE. C C C COMPUTE CAPE/CINS C C WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF C G * (LN(THETAP) - LN(THETAA)) * DZ C C (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS) C C WHERE: C THETAP IS THE PARCEL THETA C THETAA IS THE AMBIENT THETA C DZ IS THE THICKNESS OF THE LAYER C C USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT C AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL. C C IEQL = EQ LEVEL C C INITIALIZE CAPE AND CINS ARRAYS C !$omp parallel do DO J=JSTA,JEND DO I=1,IM CAPE(I,J) = D00 CINS(I,J) = D00 LCL(I,J) = D00 THESP(I,J)= D00 IEQL(I,J) = LM+1 ENDDO ENDDO C !$omp parallel do DO L=1,LM DO J=JSTA,JEND DO I=1,IM TPAR(I,J,L)= D00 ENDDO ENDDO ENDDO C C FOR TYPE 2 CAPE/CINS SMOOTH DATA PRIOR TO CALCULATION. C NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D C ARE DUMMY ARRAYS. C IF (ITYPE.EQ.2) THEN CALL P2FILT(ISMTHP,HBM2,P1D) CALL P2FILT(ISMTHT,HBM2,T1D) CALL P2FILT(ISMTHQ,HBM2,Q1D) CALL BOUNDL(Q1D,H1M12,H99999,IM,JM) ENDIF C-------FOR ITYPE=1----------------------------------------------------- C---------FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- C-------FOR ITYPE=2----------------------------------------------------- C---------FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D--------------------- C--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES------------------- do20: DO KB=1,LM IF (ITYPE.EQ.2.AND.KB.GT.1) CYCLE do20 doout10: DO J=JSTA_M,JEND_M doin10: DO I=2,IM-1 LMHK =LMH(I,J) PSFCK =AETA(LMHK)*PDSL(I,J)+PT PKL = AETA(KB)*PDSL(I,J)+PT IF (ITYPE.EQ.1.AND.(PKL.LT.PSFCK-DPBND.OR.PKL.GT.PSFCK)) & CYCLE doin10 IF (ITYPE.EQ.1) THEN TBTK =T(I,J,KB) QBTK =Q(I,J,KB) APEBTK =(H10E5/(PDSL(I,J)*AETA(KB)+PT))**CAPA ELSE PKL =P1D(I,J) TBTK =T1D(I,J) QBTK =Q1D(I,J) APEBTK =(H10E5/PKL)**CAPA ENDIF C--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX-------------- TTHBTK =TBTK*APEBTK TTHK =(TTHBTK-THL)*RDTH QQ(I,J)=TTHK-AINT(TTHK) ITTBK =INT(TTHK)+1 C--------------KEEPING INDICES WITHIN THE TABLE------------------------- IF(ITTBK.LT.1) THEN ITTBK =1 QQ(I,J)=D00 ENDIF IF(ITTBK.GE.JTB) THEN ITTBK =JTB-1 QQ(I,J)=D00 ENDIF C--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY--------------- BQS00K=QS0(ITTBK) SQS00K=SQS(ITTBK) BQS10K=QS0(ITTBK+1) SQS10K=SQS(ITTBK+1) C--------------SCALING SPEC. HUMIDITY & TABLE INDEX--------------------- BQK =(BQS10K-BQS00K)*QQ(I,J)+BQS00K SQK =(SQS10K-SQS00K)*QQ(I,J)+SQS00K TQK =(QBTK-BQK)/SQK*RDQ PP(I,J)=TQK-AINT(TQK) IQ =INT(TQK)+1 C--------------KEEPING INDICES WITHIN THE TABLE------------------------- IF(IQ.LT.1) THEN IQ =1 PP(I,J)=D00 ENDIF IF(IQ.GE.ITB) THEN IQ =ITB-1 PP(I,J)=D00 ENDIF C--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.------- IT=ITTBK P00K=PTBL(IQ ,IT ) P10K=PTBL(IQ+1,IT ) P01K=PTBL(IQ ,IT+1) P11K=PTBL(IQ+1,IT+1) C--------------SATURATION POINT VARIABLES AT THE BOTTOM----------------- TPSPK=P00K+(P10K-P00K)*PP(I,J)+(P01K-P00K)*QQ(I,J) 2 +(P00K-P10K-P01K+P11K)*PP(I,J)*QQ(I,J) APESPK=(H10E5/TPSPK)**CAPA TTHESK=TTHBTK*EXP(ELOCP*QBTK*APESPK/TTHBTK) C--------------CHECK FOR MAXIMUM THETA E-------------------------------- IF(TTHESK.GT.THESP(I,J)) THEN PSP (I,J)=TPSPK THESP(I,J)=TTHESK ENDIF END DO doin10 END DO doout10 END DO do20 C-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------ C-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------ C-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)--------------------- DO L=1,LM !$omp parallel do !$omp& private(pkl) DO J=JSTA_M,JEND_M DO I=2,IM-1 PKL = AETA(L)*PDSL(I,J)+PT IF (PKL.LT.PSP(I,J)) LCL(I,J)=L+1 ENDDO ENDDO ENDDO !$omp parallel do DO J=JSTA_M,JEND_M DO I=2,IM-1 IF (LCL(I,J).GT.LMH(I,J)) LCL(I,J)=LMH(I,J) ENDDO ENDDO C----------------------------------------------------------------------- C---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)--------- C----------------------------------------------------------------------- !!$omp parallel do DO 30 L=LM,1,-1 C--------------SCALING PRESSURE & TT TABLE INDEX------------------------ KNUML=0 KNUMH=0 DO J=JSTA_M,JEND_M DO I=2,IM-1 KLRES(I,J)=0 KHRES(I,J)=0 PKL=AETA(L)*PDSL(I,J)+PT IF(L.LE.LCL(I,J)) THEN IF(PKL.LT.PLQ)THEN KNUML=KNUML+1 KLRES(I,J)=1 ELSE KNUMH=KNUMH+1 KHRES(I,J)=1 ENDIF ENDIF ENDDO ENDDO C*** C*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSUREPLQ C** IF(KNUMH.GT.0)THEN CALL TTBLEX(TPAR(1,1,L),TTBLQ,ITBQ,JTBQ,KHRES,PDSL,AETA(L) 1, HTM(1,1,L),PT,PLQ,QQ,PP,RDPQ,THE0Q,STHEQ,RDTHEQ 2, THESP,IPTB,ITHTB) ENDIF C------------SEARCH FOR EQ LEVEL---------------------------------------- DO J=JSTA_M,JEND_M DO I=2,IM-1 IF(KHRES(I,J).GT.0)THEN IF(TPAR(I,J,L).GT.T(I,J,L)) IEQL(I,J)=L ENDIF ENDDO ENDDO C DO J=JSTA_M,JEND_M DO I=2,IM-1 IF(KLRES(I,J).GT.0)THEN IF(TPAR(I,J,L).GT.T(I,J,L)) IEQL(I,J)=L ENDIF ENDDO ENDDO C----------------------------------------------------------------------- 30 CONTINUE C------------COMPUTE CAPE AND CINS-------------------------------------- c DO K=KHL01,KHH01 c LCLK=LCL(K) c IEQK=IEQL(K) cCDIR$ SHORTLOOP c DO L=IEQK,LCLK c PRESK=AETA(L)*PDSL(K)+PT c DZKL=T(K,L)*(Q(K,L)*D608+H1)*ROG*PDSL(K)*DETA(L)/PRESK c THETAP=TPAR(K,L)*(H10E5/PRESK)**CAPA c THETAA=T(K,L)*(H10E5/PRESK)**CAPA c IF (THETAP.LT.THETAA) c & CINS(K)=CINS(K)+G*(ALOG(THETAP)-ALOG(THETAA))*DZKL c IF (THETAP.GT.THETAA) c & CAPE(K)=CAPE(K)+G*(ALOG(THETAP)-ALOG(THETAA))*DZKL c ENDDO c ENDDO LBEG=100 LEND=0 C !$omp parallel do !$omp& private(lbeg,lend) DO J=JSTA_M,JEND_M DO I=2,IM-1 LBEG=MIN(IEQL(I,J),LBEG) LEND=MAX(LCL(I,J),LEND) ENDDO ENDDO C DO L=LBEG,LEND DO J=JSTA_M,JEND_M DO I=2,IM-1 IDX(I,J)=0 IF(L.GE.IEQL(I,J).AND.L.LE.LCL(I,J))THEN IDX(I,J)=1 ENDIF ENDDO C ENDDO C !$omp parallel do !$omp& private(dzkl,presk,thetaa,thetap) DO J=JSTA_M,JEND_M DO I=2,IM-1 IF(IDX(I,J).GT.0)THEN PRESK=AETA(L)*PDSL(I,J)+PT DZKL=T(I,J,L)*(Q(I,J,L)*D608+H1) * *ROG*PDSL(I,J)*DETA(L)/PRESK THETAP=TPAR(I,J,L)*(H10E5/PRESK)**CAPA THETAA=T(I,J,L)*(H10E5/PRESK)**CAPA IF(THETAP.LT.THETAA)THEN CINS(I,J)=CINS(I,J) cjjt * +G*(ALOG(THETAP/THETAA))*DZKL * +G*(ALOG(THETAP)-ALOG(THETAA))*DZKL ELSEIF(THETAP.GT.THETAA)THEN CAPE(I,J)=CAPE(I,J) cjjt * +G*(ALOG(THETAP/THETAA))*DZKL * +G*(ALOG(THETAP)-ALOG(THETAA))*DZKL ENDIF ENDIF ENDDO ENDDO ENDDO C C ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER C LIMIT OF 0.0 ON CINS. C !$omp parallel do doout40: DO J=JSTA,JEND doin40: DO I=1,IM c IF (CAPE(K).LT.D00) c X WRITE(STDOUT,*) c X 'CALCAPE: NEGATIVE CAPE AT K=',K,CAPE(K) c IF (CINS(K).GT.D00) c X WRITE(STDOUT,*) c X 'CALCAPE: POSTIVE CINS AT K = ',K,CINS(K) CAPE(I,J) = AMAX1(D00,CAPE(I,J)) CINS(I,J) = AMIN1(CINS(I,J),D00) END DO doin40 END DO doout40 C C END OF ROUTINE. C RETURN END