cape.f90

计算几个复杂的天气诊断指数

!           ES =====> SATURATED VAPOR PRESSURE   (HPA)
	INTEGER LEV
	DIMENSION T(LEV),TD(LEV),RH(LEV)
      DIMENSION E(LEV),ES(LEV)
	CALL CAL_E(TD,LEV,E)
	CALL CAL_ES(T,LEV,ES)
      DO K=1,LEV
	    IF(E(K).LT.100.AND.ES(K).LT.100)THEN
	        RH(K)=E(K)/ES(K)
	    ELSE
	        RH(K)=99999.9
	    ENDIF
	ENDDO
	RETURN
	END

!	CAL_DCAPE.FOR BEGIN
!     DCAPE 下沉对流有效位能
!     从地面往上找，直到500HP，从最小湿球位温处下沉，若没有最小值，
!     则从此600处下沉
!      SUBROUTINE  CAL_DCAPE(FLOAT PP(M,N,40),FLOAT TT(M,N,40),FLOAT TD(M,N,40),INT NUM_PT,FLOAT *GETDCAPE)

!     参数：各点NUM_PT曾上的P、T、TD（三维数组），
!     所求的量DCAPEX（二维数组）
!     FLOAT TW_DEX(FLOAT,FLOAT,FLOAT)
!     //湿球温度
!     FLOAT PC_TC(FLOAT P0,FLOAT T0,FLOAT TD0,FLOAT *PC,FLOAT *TC)
!     //抬升凝结高度PC,TC
!     FLOAT EQU_T(FLOAT P,FLOAT T)//
!     //返回P,T处的假相当位温
!     FLOAT GET_T(FLOAT K,FLOAT P,FLOAT T) 
!     //参数：假相当位温，压强P及上一点的温度T，返回P处对应的T
!     FLOAT GET_T1(FLOAT K,FLOAT P,FLOAT T) 
!     //参数：假相当位温，压强P及上一点的温度T，返回P处对应的T

      SUBROUTINE  CAL_DCAPE(PP,TT,TD,NUM_PT,DCAPEX)
      REAL PP,TT,TD
      DIMENSION PP(NUM_PT),TT(NUM_PT),TD(NUM_PT)

      STEPP=10.0
      	 
      IFT0=0
      DCAPE=-999
      
      IF(PP(NUM_PT).GT.500) GOTO 799

      PMIN=PP(1)
      CALL TW_DEX(PP(1),TT(1),TD(1),TW600)						
      CALL EQU_T(PMIN,TW600,QMIN)
      P600=-999
      T600=-999
      TD600=-999
      DO 10,WHILE(PMIN.GT.500)
      	
        PMIN=PMIN-STEPP
!       FOR(I=0I.LE.NUM_PTI++)
        DO 2 I=1,NUM_PT-1
         IF(PP(I).GE.PMIN.AND.PP(I+1).LE.PMIN)THEN
       TMIN=TT(I)+LOG(PMIN/PP(I))/LOG(PP(I+1)/PP(I))* (TT(I+1)-TT(I))
       TDMIN=TD(I)+LOG(PMIN/PP(I))/LOG(PP(I+1)/PP(I))*(TD(I+1)-TD(I))
      	 CALL TW_DEX(PMIN,TMIN,TDMIN,TWMIN)
      	 CALL EQU_T(PMIN,TWMIN,QX) 
      	 GOTO 5
         END IF
2       CONTINUE
5       IF(QX.LT.QMIN)THEN
            QMIN=QX
            P600=PMIN
            T600=TMIN
            TD600=TDMIN
        END IF
10    CONTINUE

      IF(P600.LT.500.OR.ABS(T600).GT.800.OR.ABS(TD600).GT.800)THEN
!     FIRST P600=T600=TD600=-999
!       FOR(I=0I.LT.NUM_PT-1I++)
        DO 20 I=1,NUM_PT-1
           IF(PP(I).GE.600.AND.PP(I+1).LT.600)THEN
           P600=600
           T600=TT(I)+LOG(P600/PP(I))/LOG(PP(I+1)/PP(I))*(TT(I+1)-TT(I))
           TD600=TD(I)+LOG(P600/PP(I))/LOG(PP(I+1)/PP(I))*(TD(I+1)-TD(I))
           GOTO 30
           END IF
20      CONTINUE
30	END IF

      IF(P600.LT.500.OR.ABS(T600).GT.800.OR.ABS(TD600).GT.800)GOTO 799

      DCAPE=0.0

      CALL TW_DEX(P600,T600,TD600,TW600)
      P1=P600
      T1=TW600
      CALL EQU_T(P600,TW600,QC)

      P2=P600
      P1=P600

      DO 50,WHILE(P2.LE.PP(1)) 
        CALL EQU_T(P1,T1,Q1)
        CALL EQU_T(P1,T1-0.1,Q2)
        T2=T1-(Q1-QC)/(Q1-Q2)*0.1	 
        
!       FOR(I=NUM_PT,I.GE.1,I--) //计算DCAPE
        DO 40 I=NUM_PT,2,-1
           IF(P1.GE.PP(I).AND.P1.LT.PP(I-1)) THEN
            T11=TT(I)+LOG(P1/PP(I))/LOG(PP(I-1)/PP(I))*(TT(I-1)-TT(I))
!                 T11,T22环境温度
            T22=TT(I)+LOG(P2/PP(I))/LOG(PP(I-1)/PP(I))*(TT(I-1)-TT(I))
              DCAPE=DCAPE+287.0/2.*(T11+T22-T1-T2)*LOG(P2/P1)
              GOTO 41
           END IF
40      CONTINUE
41      P1=P2
        T1=T2
        P2=P2+STEPP	  
50    CONTINUE
799   DCAPEX=DCAPE
800   CONTINUE
	if(dcapex<0) dcapex=0
      END 
!	CAL_DCAPE.FOR END

      SUBROUTINE EQU_T(P,T,K)
      REAL K,LV
      A=7.5*T/(237.3+T)
      LV=2500000.-2368*T
      E=6.11*10**A
      B=1000.0/(P-E)
      K=(273.16+T)*(B**0.286)*EXP(LV*0.622*E/(1004.*(273.16+T)*(P-E)))
      END

      SUBROUTINE	TW_DEX(P0,T0,TD0,TW)
      REAL L
      L=597.4
      CPD=0.2403
      TW=0 
      IF(P0.LT.10.OR.ABS(T0).GT.85.OR.ABS(TD0).GT.85)RETURN
      P=P0
      T1=TD0
      T2=TD0-0.1
      CALL COM_QW(P,TD0,QS)
!   	环境比湿
      I=1

10	T1=T1+0.1*ABS(T1-T2)
      CALL COM_QW(P,T1,QWS) 
!     饱和比湿
      T2=T0-L*(QWS-QS)/CPD
      IF(ABS(T1-T2).LT.0.05) THEN
      	TW=T2
      	RETURN
      END IF
      I=I+1
      IF(T2.GT.TD0)GOTO 10
      END
      SUBROUTINE COM_QW(P,T,QS)
      IF(T.GE.0)THEN
         A=7.5*T/(237.3+T)
      ELSE
         A=9.5*T/(265.5+T)
      END IF
      E=6.11*(10**A)
      QS=0.622*E/P
      END

	SUBROUTINE CAL_CCL(TD,T,P,LEV,CCL,TCON)
!     THIS SUBROUTINE CALCULATE CONVECTIVE CONDENSATION LEVEL
!    INPUT TD =========> DEW POINT  (CELSIUS)
!         P ==========> PRESSURE   (HPA)
!           T ==========> TEMPERATURE  (CELSIUS)
!           LEV ========> LEVEL
!     OUTPUT CCL =======> CONVECTIVE CONDENSATION LEVEL (HPA)
!            TCON ======> TEMPERATURE OF CONVECTIVE (CELSIUS)
!     TEMP ARRAY
!          TE ==========> ENVIRONMENT TEMPERATURE (CELSIUS)
!         TDS =========> SATURATED DEW POINT OF PARCEL AT EACH LEVEL
	DIMENSION TD(LEV),T(LEV),P(LEV)
	DIMENSION TE(110),TDS(110)
	INTEGER LEV
	REAL CCL
      CALL CAL_TE(P,T,TE,LEV)
	TD0=TD(1)
	CALL CAL_TDS(TD0,P,LEV,TDS)
	KLFC=0
	IF(TDS(1).GT.TE(1).OR.(TDS(1).EQ.TE(1).AND.TDS(2).GT.TE(2)))THEN
	    CCL=P(1)
	    TCON=T(1)
	ELSEIF((TDS(1).EQ.TE(1).AND.TDS(2).LE.TE(2)).OR.TDS(2).LE.TE(2))THEN
          P_L=INT(P(1)/10.)*10
	    P_T=P_L-10
          DO K=2,90
	        IF(TDS(K).GT.TE(K))THEN
	            KLFC=K
                  P_B=P_T+10
	            DELT1=TE(K-1)-TDS(K-1)
	            DELT2=TDS(K)-TE(K)
	            CCL=(P_B*DELT2+P_T*DELT1)/(DELT1+DELT2)
	            CCT=TE(K-1)-DELT1/(DELT1+DELT2)*(TE(K-1)-TE(K))
	            THCL=(CCT+273.15)*(1000./CCL)**(287/1005.)
	            TCON=THCL*(1000/P(1))**(-287./1005.)-273.15
	            GOTO 290
	         ENDIF
	    P_T=P_T-10
	    ENDDO
	   IF(KLFC.EQ.0)CCL=99999.9
	ENDIF
290	CONTINUE
	RETURN
	END

	SUBROUTINE CAL_TDS(TD0,P,LEV,TDS)
!     THIS SUBROUTINE CALCULATE LAPSE RATE OF DEW POINT
!     INPUT TDO ======> DEW POINT OF THE FIRST PRESSURE LEVEL (CELSIUS)
!           P ========> PRESSURE  (HPA)
!           LEV ======> NUMBER OF LEVEL
!     OUTPUT TDS =====> DEW POINT ALONG THE SPECIFIC Q LINE
	DIMENSION P(LEV),TDS(110)
	INTEGER LEV
	REAL C
	C=461/2500000.
	TDS(1)=TD0
!	WRITE(*,*)' TD0 = ',TD0
	P1=INT(P(1)/10.)*10
	P2=P1-10.
	DO K=2,110
	   TDS(K)=1/(C*ALOG(P1/P2)+1/(TDS(K-1)+273.15))-273.15
	    P1=P1-10.
	    P2=P2-10.
          IF(P2.LT.300.)RETURN
	ENDDO
	RETURN
	END

	SUBROUTINE CAL_TW(T,TD,P,LEV,TW)
!     THIS SUBROUTINE CALCULATE WET-BULB TEMPERATURE
!     NOTICE THIS SUBROUTINE IS SUPPORTED BY CAL_DHZHF
!     INPUT   T ========> TEMPERATURE   (CELSIUS DEGREE)
!             TD =======> DEW POINT     (CELSIUS DEGREE)
!             P ========> PRESSURE      (HPA)
!            LEV ======> LEVEL
!    OUTPUT	TW =======> WET-BULB TEMPERATURE (CELSIUS DEGREE)
!     TEMP ARRAY
!             Q ========> SPECIFIC HUMIDITY
	DIMENSION T(LEV),TD(LEV),P(LEV)
	DIMENSION TW(LEV),Q(LEV)
	INTEGER LEV
	REAL CPD,TWW
	CPD=1005.
	DO K=1,LEV
	    Q(K)=0.622*6.112*10**(7.5*TD(K)/(237.3+TD(K)))/P(K)
	ENDDO
	DO K=1,LEV
          TOTAL0=T(K)+(2500000.-2368*TD(K))*Q(K)/CPD
	    PK=P(K)
	    T0=T(K)
	    TD0=TD(K)
	    CALL CAL_DHZHF(TOTAL0,T0,TD0,PK,TWW)
	    TW(K)=TWW
	ENDDO
      RETURN
	END

	SUBROUTINE CAL_DHZHF(TOTAL0,T0,TD0,PK,TW)
!     THIS SUBROUTINE CALCULATE A PROPER TEMPERATUE BY USING
!     ITERATIVE METHOD
!     NOTICE IT IS SUPPORTED BY CAL_QL
!     INPUT   TOTAL0 ==> TARGET FOR ITERATION
!             T0 ======> TEMPERATURE    (CELSIUS DEGREE)
!             TD0 =====> DEW POINT      (CELSIUS DEGREE)
!            PK ======> PRESSURE       (HPA)
!     OUTPUT  TW ======> WET-BULB TEMPERATURE (CELSIUS DEGREE)
	REAL TOTAL0,T0,TD0,PK,TW,LQ
	CPD=1005.
	TM0=T0
	TM1=TD0
101	TM=(TM0+TM1)/2.
	LQ=(2500000.-2368*TM)*6.112*10**(7.5*TM/(237.3+TM))
	TOTALM=TM+0.622*LQ/PK/CPD
      IF(ABS(TOTALM-TOTAL0).GT.0.001)THEN
          IF(TOTALM.GT.TOTAL0)THEN
	        TM0=TM
	        TM1=TM1
	    ELSE
	        TM0=TM0
	        TM1=TM
	    ENDIF
	    GOTO 101
	ELSE
	    TW=TM
      ENDIF
	RETURN
	END
	

      SUBROUTINE CAL_SSI(p,H0,H,U,V,lev,CAPE,shr,SSI)
!***********************************************************************
!    Compute INDEX OF STORM  INTENSITY                      by liuhz 2003.9   *
!    INPUT   U =====> WIND U (M/S)
!            V =====> WIND V (M/S)
!           H ======> GEOPOTENTIAL HEIEGHT (gpm) 
!            H0 =====> GEOPOTENTIAL HEIEGHT ON GIVEN HIGHT  (gpm) 
!          CAPE =====> CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG)
!    OUTPUT  SSI  ====> INDEX OF STORM INTENCITY
!***********************************************************************
   
      real H(LEV),U(LEV),V(LEV)
!      DIMENSION P(K0),H1(K0),U1(K0),V1(K0)
      
	REAL cape,u0,v0,H0,SSI,shr,h01
!	REAL SHR2,H01
      PI=3.1415926
50	SSI=9999.
	
      H01=H0

	CALL CAL_H0_UV(H0,H,U,V,lev,U0,V0)
	
	if((u0.gt.9000).or.(u(1).gt.9000).or.(v0.gt.9000).or.(v(1).gt.9000)) then
          shr=9999.9
	else 
	   SHR=SQRT((U0-U(1))**2+(V0-V(1))**2)/H01
	end if
	WRITE(*,*)'SHR=',SHR
!	PAUSE
	if(shr.gt.9000) then
	   ssi=9999.9
	else
	   SSI=100*(2+0.276*ALOG(SHR)+2.011*CAPE*0.0001)
	end if
! 200  CONTINUE

	RETURN
	END
!-----------------
	SUBROUTINE CAL_H0_UV(H0,H,U,V,LEV,U0,V0)
!    THIS SUBROUTINE CALCULATE WIND SPEED(U0,V0) ON GIVEN HIGHT H0
	DIMENSION H(LEV),U(LEV),V(LEV)
	real hmin,umin,vmin,hmax,umax,vmax
	INTEGER LEV
	REAL H0,U0,V0
	u0=9999.0
	v0=9999.0	
	DO K=1,LEV
	if((h(k).lt.99990).and.(h(k).LE.H0)) then
	   hmin=h(k)
	   umin=u(k)
	   vmin=v(k)
	endif
	if((h(k).lt.99990).and.(h(k).GT.H0)) then
	   hmax=h(k)
	   umax=u(k)
	   vmax=v(k)
	   goto 123
	endif
	ENDDO
123	continue
	      CC=(Hmax-H0)/(Hmax-Hmin)
	      U0=Umax+CC*(Umin-Umax)
	      V0=Vmax+CC*(Vmin-Vmax)
!	print*,'u0=',u0,'v0=',v0
!	pause
	RETURN
	END