ssi.f90

计算修正的全总指数、K指数、A指数

      SUBROUTINE CAL_SSI(H0,H,U,V,lev,CAPE,shr,SSI)
!***********************************************************************
!    Compute INDEX OF STORM  INTENCITY                      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)
      
	REAL cape,u0,v0,H0,SSI,shr,h01,h1,u1,v1,vs0,u2,v2,vs1,vs2
      PI=3.1415926
      H01=H0
	 h1=500
	 vs1=sqrt(u(1)**2+v(1)**2)
	CALL CAL_H0_UV(H0,H,U,V,lev,U0,V0)
	vs0=sqrt(u0**2+v0**2)
	CALL CAL_H0_UV(H1,H,U,V,lev,U1,V1)
	u2=u1+u(1)
	v2=v1+v(1)
	vs2=sqrt(u2**2+v2**2)
!	print*, vs0,vs1,vs2
	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
		shr=(vs0-vs1-0.5*vs2)/h01
	end if
!	WRITE(*,*)'SHR=',SHR
!	PAUSE
	if(shr.gt.9000) then
	   ssi=9999.9
	else
	   SSI=100*(2+0.276*ALOG(abs(SHR))+2.011*CAPE*0.0001)
	end if

	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=9
	v0=9	
	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,v0
!	pause	
	RETURN
	END