📄 cape.f90
字号:
CALL CAL_TDE(P,TD,TDE,LEV)
! WRITE(*,'(10F8.2)')(TDE(K),K=1,90)
NP=1
P_L=INT(P(1)/10.)*10.
1 CONTINUE
IF(P_L.GT.100.)THEN
PDS(NP)=P_L
P_L=P_L-10.
NP=NP+1
GOTO 1
ENDIF
NP=NP-1
CALL CAL_Q(TDE,PDS,NP,QE)
CALL CAL_QS(TA,PDS,NP,QSE)
! WRITE(*,*)' NP = ',NP
! WRITE(*,'(10F8.4)')(QE(I),I=80,110)
! WRITE(*,'(10F8.4)')(QSE(I),I=80,110)
! PAUSE 110
CAPEV=0.
R=287.04
P_L=INT(P(1)/10.)*10
P_T=P_L-10.
DO K=1,NP
IF(TA(K).GT.TE(K).AND.TA(K+1).GT.TE(K+1)) then
CAPEV=CAPEV+R*((1+0.61*QSE(K))*(TA(K)+273.15)+(1+0.61*QSE(K+1))*(TA(K+1)+273.15)-(1+0.61*QE(K))*(TE(K)+273.15)-(1+0.61*QE(K+1))*(TE(K+1)+273.15))/2.*ALOG(P_L/P_T)
endif
! 1 CAPEV=CAPEV+R*(TA(K)+TA(K+1)-TE(K)-TE(K+1))/2.*ALOG(P_L/P_T)
P_L=P_T
P_T=P_T-10.
IF(P_T.LT.P(LEV))GOTO 10
ENDDO
10 CONTINUE
!print*,'pro capev',capev
! RETURN
END
SUBROUTINE CAL_CIN(PP0,TE,TA,CIN)
! CALCULATE CONVECTIVE INHIBITION
! INPUT PP0 ======> BOTTOM PRESSURE
! TE =======> ENVIRONMENTAL TEMPERATURE
! TA =======> PARCEL TEMPERATURE
! OUTPUT
! CIN ======> CONVECTIVE INHIBITION (J/KG)
REAL PP0,R,P_L,CIN
DIMENSION TE(110),TA(110)
CIN=0.
R=287.04
P_L=INT(PP0/10.)*10
P_T=P_L-10.
IF(TA(1).GT.TE(1).OR.TA(2).GT.TE(2))THEN
CIN=0.0
RETURN
ELSE
DO K=1,90
IF(TA(K).LE.TE(K).AND.TA(K+1).LT.TE(K+1).AND.P_L.GE.300) THEN
CIN=CIN-R*(TA(K)+TA(K+1)-TE(K)-TE(K+1))/2.*ALOG(P_L/P_T)
! print*,cin
! pause
ELSEIF(TA(K+2).GT.TE(K+2))THEN
RETURN
ELSEIF(P_L.LT.300.)THEN
CIN=999.9
ENDIF
P_L=P_T
P_T=P_T-10.
ENDDO
ENDIF
if(cin>999) cin=1000
if(cin<1) cin=0
RETURN
END
SUBROUTINE CAL_LI(PP0,TE,TA,LI)
! CALCULATE LIFTING INDEX
! INPUT PP0 ========> BOTTON PRESSURE
! TE =========> ENVIRONMENTAL TEMPERATURE
! TA =========> PARCEL TEMPERATURE
! OUTPUT
! LI =========> LIFTING INDEX
REAL PP0,LI,P_L
DIMENSION TE(110),TA(110)
NP=1
P_L=INT(PP0/10.)*10.
1 CONTINUE
IF(P_L.GT.500.)THEN
P_L=P_L-10.
NP=NP+1
GOTO 1
ENDIF
LI=TE(NP)-TA(NP)
! WRITE(*,*)' NP === ',NP
! WRITE(*,*)' LI === ',LI
RETURN
END
SUBROUTINE CAL_Q(TD,P,LEV,Q)
! THIS SUBROUTINE CALCULATE SPECIFIC HUMIDITY AND SATURATE SPECFIC HUMIDITY
! NOTICE THIS SUBROUTINE IS SUPPORTED BY SUBROUTINE CAL_E
! INPUT TD ======> DEW POINT TEMPERATURE (CELSIUS DEGREE)
! P =====> PRESSURE (HPA)
! LEV ====> LEVEL
! OUTPUT Q ====> SPECIFIC HUMIDITY
! TEMP ARRAY
! E =====> VAPOR PRESSURE (HPA)
INTEGER LEV
DIMENSION TD(LEV),P(LEV),Q(LEV)
DIMENSION E(LEV)
CALL CAL_E(TD,LEV,E)
DO K=1,LEV
IF(E(K).LT.100.)THEN
Q(K)=0.622*E(K)/(P(K)-0.378*E(K))
ELSE
Q(K)=99999.9
ENDIF
ENDDO
RETURN
END
SUBROUTINE CAL_QS(T,P,LEV,QS)
! THIS SUBROUTINE CALCULATE SPECIFIC HUMIDITY AND SATURATE SPECFIC HUMIDITY
! NOTICE THIS SUBROUTINE IS SUPPORTED BY SUBROUTINE CAL_ES
! INPUT TD =====> DEW POINT (CELSIUS DEGREE)
! P =====> PRESSURE (HPA)
! LEV ====> LEVEL
! QS ====> SATURATE SPECIFIC HUMIDITY
! TEMP ARRAY
! ES =====> SATURATED VAPOR PRESSURE (HPA)
INTEGER LEV
DIMENSION T(LEV),P(LEV),QS(LEV)
DIMENSION ES(LEV)
CALL CAL_ES(T,LEV,ES)
DO K=1,LEV
IF(ES(K).LT.100.)THEN
QS(K)=0.622*ES(K)/(P(K)-0.378*ES(K))
ELSE
QS(K)=99999.9
ENDIF
ENDDO
RETURN
END
SUBROUTINE CAL_ES(T,LEV,ES)
! THIS SUBROUTINE CALCULATE SATURATE VAPOR PRESSURE
! NOTICE THE FORMULA ARE DIFFERENT WHEN TEMPERATURE
! IS ABOVE OR BELOW ZERO DEGREE
! INPUT T ======> TEMPERATURE (CELSIUS DEGREE)
! LEV ====> LEVEL
! OUTPUT ES ====> SATURATE VAPOR PRESSURE (HPA)
INTEGER LEV
DIMENSION T(LEV),ES(LEV)
DO K=1,LEV
IF(T(K).GT.-150.AND.T(K).LT.60.)THEN
IF(T(K).GE.0)ES(K)=6.112*10**(7.5*T(K)/(237.3+T(K)))
IF(T(K).LT.0)ES(K)=6.112*10**(9.5*T(K)/(265.5+T(K)))
ELSE
ES(K)=99999.9
ENDIF
ENDDO
RETURN
END
SUBROUTINE CAL_E(TD,LEV,E)
! THIS SUBROUTINE CALCULATE VAPOR PRESURE
! NOTICE THE FORMULA ARE DIFFERENT WHEN DEW POINT TEMPERATURE
! IS ABOVE OR BELOW ZERO DEGREE
! INPUT TD =====> DEW POINT (CELSIUS DEGREE)
! LEV ====> LEVEL
! OUTPUT E ====> VAPOR PRESSURE (HPA)
INTEGER LEV
DIMENSION TD(LEV),E(LEV)
DO K=1,LEV
IF(TD(K).GT.-150.AND.TD(K).LT.60.)THEN
IF(TD(K).GE.0)E(K)=6.112*10**(7.5*TD(K)/(237.3+TD(K)))
IF(TD(K).LT.0)E(K)=6.112*10**(9.5*TD(K)/(265.5+TD(K)))
ELSE
E(K)=99999.9
ENDIF
ENDDO
RETURN
END
SUBROUTINE CAL_TDE(PP,TP,TE,KMAX)
! ********** CALCULATE TEMPERTURE AT EACH LEVEL *********
! PARAMETER (KMAX=30)
DIMENSION PP(KMAX),TP(KMAX),TE(110)
REAL P_L
INTEGER I_NL,IN
DO K=1,110
TE(K)=-120.
ENDDO
I_NL=KMAX
I=1
IN=1
101 P_L=INT(PP(I)/10.)*10.0
P_T=INT(PP(I+1)/10.)*10.+10.
10 TE(IN)=TP(I)+(TP(I)-TP(I+1))/ALOG(PP(I)/PP(I+1))*ALOG(P_L/PP(I))
IN=IN+1
P_L=P_L-10.
IF(P_T.GT.PP(I))IN=IN-1
IF(P_L.GE.P_T)GOTO 10
I=I+1
! IF(PP(I+1).GT.50.OR.(I+1).LT.(I_NL))THEN
IF((I+1).LE.(I_NL))THEN
GOTO 101
ELSE
GOTO 20
ENDIF
20 CONTINUE
IF((PP(I_NL).GE.100.).AND.(INT(PP(I_NL)/10.)*10.).EQ.(PP(I_NL))) TE(IN)=TP(I_NL)
RETURN
END
SUBROUTINE CAL_BCAPE(P,T,TD,LEV,BCAPE)
! THIS SUBROUTINE CALCULATE THE BEST COVECTIVE AVAILABLE POTENTIAL ENERGY
! THIS PROGRAM IS USED WITH SUB(CAL_TC_PC,CAL_TA,CAL_TE,CAL_ENGY,
! IN WHICH SUBROUTINE CAL_TA AND CAL_TC_PC ARE USED WITH CAL_THSE
! INPUT T ======> TEMPERATURE (CELSIUS DEGREE)
! TD =====> DEW POINT (CELSIUS DEGREE)
! P ======> PRESSURE (HPA)
! LEV ====> LEVEL
! OUTPUT
! BCAPE ==> BEST CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG)
! TEMP ARRAY
! INPUT TP =====> TEMPERATURE FROM ABOVE MOIST UNSTABLE LAYER
! (CELSIUS DEGREE)
! TDP ====> DEW POINT FROM ABOVE MOIST UNSTABLE LAYER
! (CELSIUS DEGREE)
! PP =====> PRESSURE FROM ABOVE MOIST UNSTABLE LAYER (HPA)
DIMENSION P(LEV),T(LEV),TD(LEV),PP(LEV),TP(LEV),TDP(LEV)
DIMENSION TE(120),TA(120)
REAL BCAPE,TC,PC,THSE,PC0,TC0,TD0,THSE0
INTEGER KP
KP=0
THSE=220.
DO K=1,LEV
P0=P(K)
T0=T(K)
TD0=TD(K)
CALL CAL_TC_PC(P0,T0,TD0,PC0,TC0,THSE0)
IF(P(1)-P0.LT.300.AND.THSE0.GT.THSE)THEN
THSE=THSE0
PC=PC0
TC=TC0
KC=K
ENDIF
ENDDO
KL=LEV-KC+1
DO K=1,KL
PP(K)=P(K+KC-1)
TP(K)=T(K+KC-1)
TDP(K)=TD(K+KC-1)
ENDDO
P0=PP(1)
T0=TP(1)
CALL CAL_TE(PP,TP,TE,KL)
CALL CAL_TA(P0,T0,PC,TC,THSE,TA)
CALL CAL_VIR_ENGY(PP,TP,TDP,KL,TE,TA,BCAPE)
! CALL CAL_ENGY(P0,TE,TA,BCAPE)
RETURN
END
SUBROUTINE CAL_BIC(P,T,TD,LEV,BIC)
! THIS SUBROUTINE CALCULATE BEST STABILITY INDEX OF CONVECTION
! NOTICE THIS PROGRAM IS USED WITH CAL_TE AND CAL_TDE
! INPUT T ======> TEMPERATURE (CELSIUS DEGREE)
! TD =====> DEW POINT (CELSIUS DEGREE)
! P ======> PRESSURE (HPA)
! LEV ====> LEVEL
! OUTPUT
! BIC ====> BEST INDEX OF CONVECTIVE INDEX (KELVIN)
! TEMP ARRAY
! THE =====> EQUIVLENT POTENTIAL TEMPERATURE (KELVIN)
DIMENSION P(LEV),T(LEV),TD(LEV)
DIMENSION THE(LEV)
REAL BIC,THMAX,THMIN
THMAX=220
THMIN=500
IF(P(1).LT.850.)THEN
BIC=99999.9
RETURN
ELSE
CALL CAL_THE(T,TD,P,LEV,THE)
DO K=1,LEV
IF(P(1)-P(K).LT.150.AND.THE(K).GT.THMAX)THEN
THMAX=THE(K)
ENDIF
IF(P(K).GE.500.AND.P(K).LE.650.AND.THE(K).LT.THMIN)THEN
THMIN=THE(K)
ENDIF
IF(P(K).LE.500)THEN
BIC=THMIN-THMAX
RETURN
ENDIF
ENDDO
ENDIF
RETURN
END
SUBROUTINE CAL_THE(T,TD,P,LEV,THE)
! THIS SUBROUTINE CALCULATE EQUIVALENT POTENTIAL TEMPERATURE
! NOTICE THIS SUBROUTINE IS SUPPORTED BY SUBROUTINE CAL_E, CAL_R AND CAL_RH
! INPUT T ======> TEMPERATURE (CELSIUS DEGREE)
! TD =====> DEW POINT (CELSIUS DEGREE)
! P ======> PRESSURE (HPA)
! LEV ====> LEVEL
! OUTPUT THE ====> EQUIVALENT POTENTIAL TEMPERATURE (KELVIN)
! TEMP ARRAY
! R ======> MIXING RATIO
! RH =====> RELATIVE HUMIDITY
! E ======> VAPOR PRESSURE
INTEGER LEV
DIMENSION T(LEV),TD(LEV),P(LEV),THE(LEV)
DIMENSION R(LEV),E(LEV),RH(LEV)
REAL RD,RV,CPD
RD=287.04
RV=461.50
CPD=1005.
CALL CAL_E(TD,LEV,E)
CALL CAL_R(TD,P,LEV,R)
CALL CAL_RH(T,TD,LEV,RH)
DO K=1,LEV
! P0=P(K)
! T0=T(K)
! TD0=TD(K)
! CALL CAL_TC_PC(P0,T0,TD0,PC0,TC0,THSE0)
! THE(K)=THSE0
IF(T(K).LT.70.AND.R(K).LT.1)THEN
THE(K)=(T(K)+273.15)*(1000./(P(K)-E(K)))**(RD/CPD)*RH(K)**(R(K)*RV/CPD)*EXP((2500000.-2368*T(K))*R(K)/(CPD*(T(K)+273.15)))
ELSE
THE(K)=99999.9
ENDIF
ENDDO
RETURN
END
SUBROUTINE CAL_R(TD,P,LEV,R)
! THIS SUBROUTINE CALCULATE MIXING RATIO
! NOTICE THIS SUBROUTINE IS SUPPORTED BY SUBROUTINE CAL_E
! INPUT TD =====> DEW POINT (CELSIUS DEGREE)
! P =====> PRESSURE (HPA)
! LEV ====> LEVEL
! OUTPUT R ====> MIXING RATION
! TEMP ARRAY
! E =====> VAPOR PRESSURE (HPA)
INTEGER LEV
DIMENSION TD(LEV),P(LEV),R(LEV)
DIMENSION E(LEV)
CALL CAL_E(TD,LEV,E)
DO K=1,LEV
IF(E(K).LT.100.)THEN
R(K)=0.622*E(K)/(P(K)-E(K))
ELSE
R(K)=99999.9
ENDIF
ENDDO
RETURN
END
SUBROUTINE CAL_RH(T,TD,LEV,RH)
! THIS SUBROUTINE CALCULATE RELATIVE HUMIDITY
! NOTICE THIS SUBROUTINE IS SUPPORTED BY SUBROUTINE CAL_E AND CAL_ES
! INPUT T ======> TEMPERATURE (CELSIUS DEGREE)
! TD =====> DEW POINT (CELSIUS DEGREE)
! LEV ====> LEVEL
! OUTPUT RH ====> RELATIVE HUMIDITY
! TEMP ARRAY
! E =====> VAPOR PRESSURE (HPA)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -