MODULE stability_indices use mpas_kind_types, only : RKIND ! STABILITY ANALYSIS USING PARCEL METHOD IMPLICIT NONE ! integer, parameter :: RKIND=8 PRIVATE PUBLIC :: INDX1 CONTAINS SUBROUTINE INDX1 (JNO, P, TS, TSD, RLI, RKI, RWI) !----------------------------------------------------------------------- ! ! This routine computes stability indices (Stone 1983, 1984) ! The following fields are computed by this routine: ! ! Lifted index ! K index ! Showalter index ! ! TERRAIN NOTE: levels 1 (~850 hPa) and 2 (~700 hPa) are adapted to the ! surface pressure when the standard level does not exist above ground ! (e.g. Andes, Tibetan Plateau). Level 3 (500 hPa) is extrapolated ! below the surface on the rare points where even it does not exist ! (extreme terrain, > ~5500 m). This keeps the routine as close as ! possible to the original fixed-level algorithm below, while still ! avoiding the out-of-bounds access that the fixed-level version hits ! on elevated terrain. ! ! Stone, H.M., 1983: Stability Analysis Program, NOAA Eastern Region ! Computer Programs and Problems NWS ERCP-No. 9, National Weather ! ·Service, Garden City, NY. ! ! Stone, H.M., 1984: Stability Analysis Plot Program, NOAA Eastern Region ! Computer Programs and Problems NWS ERCP-No. 16, National Weather ! ·Service, Garden City, NY. ! !----------------------------------------------------------------------- ! COMPUTES LIFTED INDEX, K INDEX, AND SHOWALTER INDEX IMPLICIT NONE INTEGER, INTENT(IN) :: JNO INTEGER :: I,J,K REAL (kind=RKIND), DIMENSION(JNO), INTENT(IN) :: P,TS,TSD REAL (kind=RKIND), DIMENSION(JNO) :: TW REAL (kind=RKIND) :: PL(3),TL(3),TDL(3) REAL (kind=RKIND), INTENT(OUT) :: RLI, RKI, RWI REAL (kind=RKIND) :: TC,TH,WTH,WTC,WSUM,THSUM,P1,P2,PFINISH REAL (kind=RKIND) :: PLOG1,PLOG2,PLOG3,PLOG4 REAL (kind=RKIND) :: DP1,E2,FACTOR,FACTORD,FACTORT,PPARCEL,TDE1,TDE2,TDPARCEL REAL (kind=RKIND) :: TE1,TE2,TH1,TH2,THAVG,THW,TP,TP500,TPARCEL,W,W1,W2,WAVEG REAL (kind=RKIND) :: WAVG,WFW,X REAL (kind=RKIND) :: DELTAZ REAL (kind=RKIND), PARAMETER :: RLAPSE = 0.0065_RKIND ! K/m, atm padrao REAL (kind=RKIND), PARAMETER :: RSTDEXP = 0.190297_RKIND ! 1/5.255 PL(1) = 850.0_RKIND PL(2) = 700.0_RKIND PL(3) = 500.0_RKIND ! ADAPTA OS NIVEIS 1 E 2 AO TERRENO QUANDO NECESSARIO (nunca deixam de ! ter um valor real: nivel 1 nunca fica abaixo da superficie, nivel 2 ! sempre fica entre a superficie e o nivel 3) IF (P(1) .LT. PL(1)) PL(1) = P(1) IF (P(1) .LT. PL(2)) & PL(2) = MIN( MAX(P(1)-150.0_RKIND, PL(3)+50.0_RKIND), P(1)-10.0_RKIND ) ! COMPUTE WET BULB POTENTIAL TEMP AT ALL SGFNT LEVELS ABV SFC DO K = 1, JNO TC=TCONOF(TS(K),TSD(K)) ! CONDENSATION TEMPERATURE TH=THETA(TS(K)+273.16,1000.0_RKIND,P(K))-273.16_RKIND ! POT TEMP DEG C WTH=WOBF(TH) WTC=WOBF(TC) TW(K)=TH-WTH+WTC ! WET BULB POT TEMP - DEG C ENDDO DP1=50. ! AVERAGES OVER FIRST -DPI- MBS. WSUM=0. THSUM=0. J=1 P1=P(1) TE1=TS(1) TDE1=TSD(1) TH1=THETA(TE1+273.16,1000.0_RKIND,P(1)) ! POT TEMP W1=WMROF(P(1),TDE1) PFINISH=P1-DP1 3 P2=PFINISH IF(P(J+1)-P2) 1,2,2 ! <0, 0, >0 2 P2=P(J+1) J=J+1 1 PLOG1=ALOG(P(J)/P(J+1)) FACTORT=(TS(J)-TS(J+1))/PLOG1 FACTORD=(TSD(J)-TSD(J+1))/PLOG1 PLOG2=ALOG(P2/P(J+1)) TE2=TS(J+1)+FACTORT*PLOG2 ! ENVIRONMENT TEMP AT P2 TDE2=TSD(J+1)+FACTORD*PLOG2 ! ENVIRONMENT DEWPT AT P2 TH2=THETA(TE2+273.16,1000.0_RKIND,P2) ! POT TEMP AT TE2,P2 W2=WMROF(P2,TDE2) ! MIXING RATIO AT P2 PLOG3=ALOG(P1/P2) TH=.5*(TH1+TH2)*PLOG3 ! AVG POT TEMP IN LYR P1-P2 W=.5*(W1+W2)*PLOG3 ! AVG MIX RATIO IN LYR P1-P2 THSUM=THSUM+TH WSUM=WSUM+W P1=P2 TH1=TH2 W1=W2 IF (P2.GT.PFINISH) GO TO 3 ! COMPUTE AVG VALUES FOR FIRST -DPI- MBS. PLOG4=ALOG(P(1)/PFINISH) THAVG=THSUM/PLOG4 WAVG=WSUM/PLOG4 PPARCEL=P(1)-.5*DP1 TPARCEL=THETA(THAVG,PPARCEL,1000.0_RKIND)-273.16 ! DEG C X=.0200*(TPARCEL-12.5+7500./PPARCEL) ! NON-IDEAL GAS CORRECTION WFW=1.+.0000045*PPARCEL+.00140*X*X ! NON-IDEAL GAS CORRECTION E2=.001*WAVG*PPARCEL/((WAVG*.001+.62197)*WFW) ! VAPOR PRES (MB) TDPARCEL=DPTOF(E2) TC=TCONOF(TPARCEL,TDPARCEL) TH=THAVG-273.16 ! POT TEMP DEG C WTH=WOBF(TH) WTC=WOBF(TC) THW=TH-WTH+WTC ! EQUIV WET BULB POT TEMP (DEG C) TP500=SATLFT(THW,500.0_RKIND) ! GET TEMP AND DEWPT AT LEVELS 1, 2 E 3 (POSSIVELMENTE ADAPTADOS) DO 5 J=1,3 DO 4 I=1,JNO IF (PL(J)-P(I)) 4,6,7 4 CONTINUE !PRINT*,"RAOB TERMINATES TOO SOON, P(JNO) = ",P(JNO) RETURN 6 TL(J)=TS(I) TDL(J)=TSD(I) GO TO 5 7 IF (I.NE.1) GO TO 8 ! So chega aqui para J=3 em terreno extremo (>500 hPa de ! altitude), ja que os niveis 1 e 2 foram adaptados acima e nunca ! ficam abaixo de P(1). Extrapola-se hidrostaticamente em vez de ! travar (mesmo principio do metodo de Shuell/NCEP para reducao ! de pressao ao nivel do mar). DELTAZ = 44330.0_RKIND*( (PL(J)/1013.25_RKIND)**RSTDEXP & - (P(1) /1013.25_RKIND)**RSTDEXP ) TL(J) = TS(I) + RLAPSE*DELTAZ TDL(J) = TSD(I) + RLAPSE*DELTAZ GO TO 5 8 FACTOR=ALOG(PL(J)/P(I))/ALOG(P(I-1)/P(I)) TL(J)=TS(I)+FACTOR*(TS(I-1)-TS(I)) TDL(J)=TSD(I)+FACTOR*(TSD(I-1)-TSD(I)) 5 CONTINUE ! COMPUTE LIFTED INDEX RLI=TL(3)-TP500 ! LIFTED INDEX ! COMPUTE K INDEX RKI=(TL(1)-TL(3))+TDL(1)-(TL(2)-TDL(2)) ! K INDEX ! COMPUTE SHOWALTER INDEX TC=TCONOF(TL(1),TDL(1)) TH=THETA(TL(1)+273.16,1000.0_RKIND,PL(1))-273.16_RKIND ! DEG C WTH=WOBF(TH) WTC=WOBF(TC) THW=TH-WTH+WTC ! EQUIV WET BULB POT TEMP TP=SATLFT(THW,500.0_RKIND) RWI=TL(3)-TP ! SHOWALTER INDEX RETURN END !-------------- REAL(kind=RKIND) FUNCTION WOBF(T) IMPLICIT NONE REAL (kind=RKIND), INTENT(IN) :: T REAL (kind=RKIND) :: X, POL ! COMPUTE BY DOUBLE ASYMPTOTIC APPROXIMATION ! CONSIDER SEPARATELY IF .GT. OR .LE. 20 DEG. ! CENT. FOR ALL TEMPS...THETU=THETA-UOBF(THETA)+UOBF(TEMPCON) ! CENT. FOR ALL TEMPS...THETM=THETA-U08F(THETA)+UOBF(TEMP) X=T-20.0 IF(X) 10,10,20 10 CONTINUE ! CURVE FIG FOR COOL TEMPERATURE RANGE POL=1.000+X*(-8.8416605E-3+X*(1.4714143E-4+X*(-9.6719890E-7 & +X*(-3.2607217E-8+X*(-3.8598073E-10))))) POL=POL*POL WOBF=15.130/(POL*POL) RETURN 20 CONTINUE ! CURVE FIT FOR WARMER TEMPERATURES POL=1.000+X*(3.6182989E-3+X*(-1.3603273E-5+X*(4.9618922E-7 & +X*(-6.1059365E-9+X*(3.9401551E-11+X*(-1.2588129E-13 & +X*(1.6688280E-16) )))))) POL=POL*POL WOBF=29.930/(POL*POL)+0.9600*X-14.800 RETURN END REAL(kind=RKIND) FUNCTION SATLFT (THM,P) ! COMPUTES TEMPERATURE (DEG C) WHERE THETA MOIST (DEG C) CROSSES P (MB) IMPLICIT NONE REAL (kind=RKIND), INTENT(IN) :: THM,P ! CONSIDER THE EXPONENTIAL FOR POTENTIAL TEMPERATURE AS ROCP REAL, PARAMETER:: ROCP=0.28571428 REAL (kind=RKIND) :: PWRP,TONE,TTWO,EONE,ETWO,RATE,EOR IF(ABS(P-1000.0)-0.0010) 100,100,200 100 SATLFT=THM RETURN 200 PWRP=(P/1000.0)**ROCP ! COMPUTE TEMPERATURE OF DRY ADIABATIC LIFT FOR FIRST GUESS TONE=(THM+273.16)*PWRP-273.16 ! CONSIDER PSEUDO-ADIABAT, EW1, THROUGH TONE AT P. ! COMPUTE EONE=EW1-THM EONE=WOBF(TONE)-WOBF(THM) RATE=1.0 GO TO 330 300 CONTINUE ! CONTRIBUTION TO ITERATION IS CHANGE IN T ! CORRESPONDING TO CHANGE IN E RATE=(TTWO-TONE)/(ETWO-EONE) TONE=TTWO EONE=ETWO 330 CONTINUE ! COMPUTE ESTIMATED SATLIFT, TTWO TTWO=TONE-EONE*RATE ! CONSIDER PSEUDO-ADIABAT, EW2, THROUGH TTWO AT P. ! COMPUTE ETW0=EW2-THM ETWO=(TTWO+273.16)/PWRP-273.16 ETWO=ETWO+WOBF(TTWO)-WOBF(ETWO )-THM ! CORRECTION TO TTWO IS EOR EOR=ETWO*RATE IF(ABS(EOR)-0.1000) 400,400,300 400 SATLFT=TTWO-EOR RETURN END REAL(kind=RKIND) FUNCTION TCONOF(TEMP,DEWPT) ! COMPUTES CONDENSATION TEMPERATURE (DEGREES CENT) BY LIFTING IMPLICIT NONE REAL (kind=RKIND), INTENT(IN) :: TEMP, DEWPT REAL (kind=RKIND) :: S, T, DLT S=TEMP-DEWPT ! CONSIDER TEMP AND DEWPT TO BE LIKE UNITS (C OR K) T=TEMP IF(100.-TEMP) 4,5,5 4 T=TEMP-273.16 ! COMPUTE CURVE FIT IN MOST EFFICIENT MANNER 5 DLT=S*(1.2185+0.001278*T+S*(-0.002190+11.73E-6*S-5.20E-6*T)) TCONOF=T-DLT RETURN END REAL(kind=RKIND) FUNCTION WMROF(P,TD) ! COMPUTE MIXING RATIO (G/KG)...DEWPOINT (DEGREES C OR K)... PRESSURE (MB) IMPLICIT NONE REAL (kind=RKIND), INTENT(IN) :: P, TD REAL (kind=RKIND) :: T, X, WFW, FWESW T=TD IF (100.-T) 3,4,4 3 T=T-273.16 ! CURVE FIT CORRECTION FOR NON-IDEAL GAS 4 X=0.0200*(T-12.5+7500.0/P) WFW=1.+0.0000045*P+0.00140*X*X ! COMPUTE ACCORDING TO STANDARD FORMULA FWESW=WFW*VAPFW(T) WMROF=621.97*(FWESW/(P-FWESW)) RETURN END REAL(kind=RKIND) FUNCTION THETA(T,P2,P1) REAL (kind=RKIND), INTENT(IN) :: T REAL (kind=RKIND), INTENT(IN) :: P2 REAL (kind=RKIND), INTENT(IN) :: P1 REAL (kind=RKIND) :: Rd, Cp, Rd_Cp Rd = 287.058_RKIND ! Specific gas constant for dry air (J/kg/K) Cp = 1004.0_RKIND ! Specific heat at constant pressure (J/kg/K) Rd_Cp = Rd/Cp ! DRY ADIABATIC (T,P1) TO (THETA,P2) THETA = T * (P2/P1)**Rd_Cp ! DRY ADIABATIC (T,P1) TO (THETA,P2) RETURN END REAL(kind=RKIND) FUNCTION DPTOF(EW) ! COMPUTE DEWPOINT, DPT, IN DEGREES C GIVEN WATER VAPOR PRESSURE (MB) IMPLICIT NONE REAL (kind=RKIND), INTENT(IN) :: EW REAL (kind=RKIND) :: TOL,X,BOT,DELTM, EDP,DTDE, DELT, DM ! CREATE TOLERANCE TO DEGREE DESIRED TOL=0.00010 IF (EW-0.21382876E-09) 20,20,30 20 DPTOF=-10000. RETURN 30 IF (1013.0-EW) 20, 100, 100 ! CREATE GUESS BY INVERTING TETEN-S FORMULA 100 X=ALOG(EW/6.1078) BOT=17.269388-X DPTOF=(237.3*X)/BOT BOT=BOT*EW DELTM=0. 200 EDP=VAPFW(DPTOF) ! CORRECT GUESS BY DERIVATIVE OF TEMPERATURE WITH RESPECT TO VAPOR PRES. ! CALCULATED FROM INVERSE OF TETEN-S FORMULA DTDE=(DPTOF+237.3)/BOT DELT=DTDE*(EW-EDP) DPTOF=DPTOF+DELT ! CHECK THAT ITERATION IS NOT IN AN ENDLESS CYCLE, A RARE SITUATION ! IF NEEDED. CHANGE -TOL- AND EXIT DM=DELT-DELTM IF(ABS(DM).GE.1.E-7) GO TO 10 ! IF DM VERY SMALL, ITERATION IS ENDLESS TOL=ABS(DELT) !PRINT*,"TOLERANCE (TOL) IN DPTOF CHANGED TO ",TOL," (NORMAL TOL = .00010)" 10 DELTM=-DELT ! CHECK TO SEE IF ANSWER CLOSE ENOUGH, IF NOT ITERATE OVER CORRECTION IF (ABS(DELT)-TOL) 300,300,200 ! CHANGE SO DEWPOINT IS ALWAYS LESS THAN THE TEMP. ! COMPATIBILITY WITH TOL IS FORCED 300 DPTOF=DPTOF-TOL RETURN END REAL(kind=RKIND) FUNCTION VAPFW(T) ! COMPUTE SATURATION VAPOR PRESSURE OVER WATER, VAPFW, IN MBS. IMPLICIT NONE REAL (kind=RKIND), INTENT(IN) :: T REAL (kind=RKIND) :: X, POL ! CONSIDER T(TEMPERATURE) IN DEGREES C OR DEGREES K. X=T IF (100.0-X) 3,4,4 3 X=X-273.16 ! CURVE FIT FOR RANGE -50 < T < 100 DEGREES C. 4 POL = 0.99999683E-00 + X *(-0.90826951E-02 + & X *(0.70736169E-04 + X *(-0.61117958E-06 + & X *(0.43884187E-08 + X *(-0.29883885E-10 + & X *(0.21874425E-12 + X *(-0.17892321E-14 + & X *(0.11112018E-16 + X *(-0.30994571E-19) )))))))) POL=POL*POL POL=POL*POL VAPFW=6.107800/ (POL*POL) RETURN END END MODULE stability_indices