www.pudn.com > cape.rar > cape.f90
PROGRAM MAIN
! parameter (ix=41,iy=31,ik=7,lev=7)
integer iii,m,len_file,len_name,lenst
INTEGER stanum1,stanum,lev1(800),lev11(800),LEV
INTEGER year,month,date,time,sta(800),station(800),number
character*128 dir_in,dir_out(5),vary_ch(5),time_ch
character*128 ch_date,ch3_micaps,file_in,file_out
INTEGER*2 LEN
REAL P(50),h(50),T(50),TD(50),ff(50),dd(50),TW(50)
real p2(800,50),t2(800,50),td2(800,50)
real h2(800,50),ff2(800,50),dd2(800,50)
real h22(800,50),ff22(800,50),dd22(800,50)
real p22(800,50),t22(800,50),td22(800,50)
real CAPE,PLFC,PE,CIN,LI,PC,TC,BCAPE,BIC,DCAPE,TCON,CCL
REAL CAPE1(800),TC1(800),PC1(800),bcape1(800),ncape1(800)
REAL PLFC1(800),PE1(800),LI1(800),CIN1(800),bic1(800)
REAL tcon1(800),ccl1(800)
REAL jd(800),wd(800),gd(800)
REAL jdd(800),wdd(800),gdd(800)
character head*100
character*100 filnam(100)
INTEGER*2 FNUM
! 读入文件列表
open(1,file='d:\wj\cape\filelist.txt',status='old')
read(1,*) FNUM
write(*,*) 'FNUM=',FNUM
do k=1,FNUM
read(1,*)filnam(k)
! filnam(k)='d:\wj\cape\dat\'//filnam(k)
! write(*,*)filnam(k)
enddo
close(1)
! file_in='E:\cape\dat\05090320.57494.text'
! 打开tlogp文件
! open(1,file='E:\cape\dat\05090320.57494.text',status='old')
do iii=1,FNUM
! head=filnam(iii)
! year=head(1:2)
! month=head(3:4)
! date=head(5:6)
! time=head(7:8)
head='d:\wj\cape\dat\'//filnam(iii)
write(*,*)head
open(1,file=head,status='old')
! read(1,*)head
! read(1,*)year,month,date,time,stanum1
stanum=1
! stanum=stanum1-2
! write(*,*)'stanum=',stanum
i=1
do 100 k=1,1
number=0
read(1,*)sta(k),jdd(k),wdd(k),gdd(k),number
write(*,*)sta(k),jdd(k),wdd(k),gdd(k),number
lev11(k)=number/6
do j=1,lev11(k)
read(1,*)p22(K,J),h22(k,J),t22(K,J),td22(K,J),ff22(k,J),dd22(k,J)
enddo
L=0
IF((FF22(K,1).LT.9990).AND.(DD22(K,1).LT.9990).and.(td22(k,1).lt.9990).and.(t22(k,1).lt.9990)) THEN
L=L+1
P2(I,L)=P22(K,1)
T2(I,L)=T22(K,1)
TD2(I,L)=TD22(K,1)
h2(i,L)=0.0
ff2(i,L)=ff22(k,1)
dd2(i,L)=dd22(k,1)
ENDIF
do j=2,lev11(k)
IF((H22(K,J).LT.9990).AND.(FF22(K,J).LT.9990).AND.(DD22(K,J).LT.9990).and.(td22(k,j).lt.9990).and.(t22(k,j).lt.9990)) THEN
L=L+1
P2(I,L)=P22(K,J)
T2(I,L)=T22(K,J)
TD2(I,L)=TD22(K,J)
h2(i,L)=h22(k,j)*10.0
ff2(i,L)=ff22(k,j)
dd2(i,L)=dd22(k,j)
ENDIF
enddo
LEV1(I)=L
if(l.gt.6)then
STATION(I)=STA(K)
JD(I)=JDd(K)
WD(I)=WDd(K)
GD(I)=GDd(K)
I=I+1
endif
GOTO 100
100 continue
!STANUM=I-2
CLOSE(1)
do 150 i=1,stanum
lev=lev1(i)
do j=1,lev
p(j)=p2(i,j)
t(j)=t2(i,j)
td(j)=td2(i,j)
h(j)=h2(i,j)
! u(j)=dd2(i,j)*sin(ff2(i,j)*3.1415926/180.)
! v(j)=dd2(i,j)*cos(ff2(i,j)*3.1415926/180.)
enddo
CALL CAL_CAPE(P,T,TD,lev,CAPE,PLFC,PE,CIN,LI,PC,TC)
cape1(i)=cape
cin1(i)=cin
li1(i)=li
write(*,*)'cape=',CAPE1(i),station(i),pc
if((plfc.lt.90000).and.(pe.lt.90000))then
ncape1(i)=cape/(plfc-pe)
else if(plfc.gt.90000) then
ncape1(i)=0.0
else if(pe.gt.90000) then
ncape1(i)=cape/(plfc-p(lev))
endif
!call CAL_BCAPE(P,T,TD,LEV,BCAPE)
!bcape1(i)=bcape
150 continue
!c 写入抬升指数
! open(2,file='d:\12\cape\抬升指数.txt')
! write(2,*)'diamond 3 ',ch3_micaps
! write(2,'(4i5,a36,i5)')year,month,date,time,' 1000 1 2000 1 2000 1 1.0 2.0 1 ',stanum
! do i=1,stanum
! write(2,'(i10,4f10.1)') station(i),jd(i),wd(i),gd(i),li1(i)
! enddo
! close(2)
!c 写入对流有效位能
! open(3,file='E:\cape\cape.txt')
! open(3,file='d:\wj\cape\cape.txt')
head='d:\wj\cape\result\'//filnam(iii)
write(*,*) head
open(3,file=head)
write(3,'(4i5,a36,i5)')year,month,date,time,' 1000 1 2000 1 2000 1 1.0 2.0 1 ',stanum
!200 format'(4i5,a36,i5)'
do i=1,stanum
write(3,'(i10,4f10.1)') station(i),jd(i),wd(i),gd(i),cape1(i)
enddo
close(3)
!c 写入最佳对流有效位能
! open(2,file='d:\12\cape\最佳对流有效位能.txt')
! write(2,*)'diamond 3 ',ch3_micaps
! write(2,'(4i5,a36,i5)')year,month,date,time,' 1000 1 2000 1 2000 1 1.0 2.0 1 ',stanum
! do i=1,stanum
! write(2,'(i10,4f10.1)') station(i),jd(i),wd(i),gd(i),bcape1(i)
! enddo
! close(2)
!c 写入归一化对流有效位能
! LEN=len_trim(DIR_OUT(4))
! file_out=dir_out(4)(1:LEN)//'\'//ch_date(1:len_name)
! ch3_micaps=ch_date(1:2)//'年'//ch_date(3:4)//'月'//ch_date(5:6)
! &//'日'//ch_date(7:8)//'时'//vary_ch(4)
! open(3,file=file_out)
! write(3,*)'diamond 3 ',ch3_micaps
!! write(3,200)year,month,date,time,' 1000 1 2000 1 2000 1
! &1.0 2.0 1 ',stanum
!200 format(4i5,a36,i5)
! do i=1,stanum
! write(3,300) station(i),jd(i),wd(i),gd(i),ncape1(i)
! write(3,300) ((station(i),jd(i),wd(i),gd(i),ncape1(i)),i=1,stanum)
! enddo
! 300 format(i10,4f10.1)
! close(3)
! 写入对流抑制能量
! LEN=len_trim(DIR_OUT(5))
! file_out=dir_out(5)(1:LEN)//'\'//ch_date(1:len_name)
! ch3_micaps=ch_date(1:2)//'年'//ch_date(3:4)//'月'//ch_date(5:6)
! &//'日'//ch_date(7:8)//'时'//vary_ch(5)
! open(4,file=file_out)
! write(4,*)'diamond 3 ',ch3_micaps
! write(4,200)year,month,date,time,' 1000 1 2000 1 2000 1
! &1.0 2.0 1 ',stanum
!200 format(4i5,a36,i5)
! do i=1,stanum
! write(4,300) station(i),jd(i),wd(i),gd(i),cin1(i)
! write(4,300) ((station(i),jd(i),wd(i),gd(i),cin1(i)),i=1,stanum)
! enddo
! 300 format(i10,4f10.1)
! close(4)
!110 continue
enddo
! close(998)
END
SUBROUTINE CAL_CAPE(P,T,TD,LEV,CAPE,PLFC,PE,CIN,LI,PC,TC)
! THIS SUBROUTINE CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY
! AND IT IS USED WITH 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 CAPE ===> EQUIVALENT POTENTIAL TEMPERATURE (J/KG)
! PLFC ===> FREE CONVECTION LEVEL
! PE =====> EMBALANCE LEVEL
! CIN ====> CONVECTIVE INHIBITION
! LI =====> LIFTING INDEX
! TEMP ARRAY
! TE =====> ENVIRONMENT TEMPERATURE (KELVIN)
! TA =====> PARCEL TEMPERATURE (KELVIN)
! TC
DIMENSION P(LEV),T(LEV),TD(LEV)
DIMENSION TE(120),TA(120)
REAL CAPE,TC,PC,THSE,PC0,TC0,TD0,THSE0
REAL PLFC,PE,LI,CIN
INTEGER LEV
! ASSUMING THE BOTTOM POINT AS THE LIFTING LAYER
THSE=220.
PC=P(1)
TC=T(1)
P0=P(1)
T0=T(1)
TD0=TD(1)
CALL CAL_TC_PC(P0,T0,TD0,PC0,TC0,THSE0)
PC=PC0
TC=TC0
CALL CAL_TE(P,T,TE,LEV)
CALL CAL_TA(P0,T0,PC,TC,THSE0,TA)
CALL CAL_PLFC_AND_PE(P,TE,TA,LEV,PLFC,PE)
open(1,file='pe_ta_te')
do k=1,90
ppp=int(p(1)/10.)*10-(k-1)*10.
write(1,'(3f10.2)')ppp,ta(k),te(k)
enddo
close(1)
CALL CAL_VIR_ENGY(P,T,TD,LEV,TE,TA,CAPE)
! print*,ta
! stop
CALL CAL_CIN(P0,TE,TA,CIN)
CALL CAL_LI(P0,TE,TA,LI)
! RETURN
END
SUBROUTINE CAL_TC_PC(P,T,TD,PC,TC,THSE)
! CALCULATE LIFTING CONDENSATION LAYER AND LIFTING CONDENSATION TEMPERATURE
! AND ALSO PSEUDO-POTENTIAL TEMPERATURE AT CONDESAATING POINT
! INPUT P ========> LIFTING PRESSURE (HPA)
! T ========> LIFTING TEMPERATURE (CELSIUS DEGREE)
! TD =======> DEW POINT (CELSIUS DEGREE)
! OUTPUT
! PC =======> CONDESATION PRESSURE
! TC =======> CONDENSATION TEMPERATURE
! THSE =====> PSEUSO EQUIVALENT POTENTION TEMPERATURE (KELVIN)
REAL P,T,TD,PC,TC,THSE
TC=T-(T-TD)*0.976/(0.976-0.000833*(237.3+TD)**2/(273.15+TD))
PC=P*((273.15+TC)/(273.15+T))**3.5
CALL CAL_THSE(PC,TC,THSE)
RETURN
END
SUBROUTINE CAL_THSE(P,T,THSE)
! 计算假相当位温，不考虑 0 度以下假相当位温计算公式的差别
REAL E,A,B,P,T,THSE
A=7.5*T/(237.3+T)
E=6.11*10.**A
B=1000.0/(P-E)
THSE=(273.16+T)*B**0.286*EXP((2.5*10.**6-2368.*T)*0.622*E/(1004.*(273.16+T)*(P-E)))
RETURN
END
SUBROUTINE CAL_THSE1(P,T,THSE)
! 本子程序考虑 0 度以下，假相当位温计算公式不同
REAL E,A,B,P,T,THSE
IF(T.GE.0)A=7.5*T/(237.3+T)
IF(T.LT.0)A=9.5*T/(265.5+T)
E=6.11*10.**A
B=1000.0/(P-E)
IF(T.GE.0) THSE=(273.16+T)*B**0.286*EXP((597.3+0.5*T)*0.622*E/(0.24*(273.16+T)*(P-E)))
IF(T.LT.0) THSE=(273.16+T)*B**0.286*EXP((667.0+0.5*T)*0.622*E/(0.24*(273.16+T)*(P-E)))
RETURN
END
SUBROUTINE CAL_TE(PP,TP,TE,KMAX)
! p ,T ,TE,LEV
! ********** CALCULATE TEMPERTURE AT EACH LEVEL *********
! 计算层结每隔 10HPA 的温度（对数气压线性内插）
! PARAMETER (KMAX=30)
DIMENSION PP(KMAX),TP(KMAX),TE(110)
REAL P_L
INTEGER I_NL,IN
DO K=1,110
TE(K)=-60.
ENDDO
! WRITE(*,*)(PP(K),K=1,KMAX)
! WRITE(*,*)(TP(K),K=1,KMAX)
! PAUSE 31
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_TA(PP0,TP0,PC,TC,THSE,TA)
! 计算状态曲线各整层高度（间隔 10HPA）上的温度
! INPUT PP0 ======> PRESSURE AT THE LIFTING LEVEL (HPA)
! TP0 ======> TEMPERATURE AT THE LIFTING LEVEL (CELSIUS DEGREE)
! PC =======> LIFTING CONDENSATION LEVEL (HPA)
! TC =======> LIFTING CONDESATION TEMPERATURE (CELSIUS DEGREE)
! THESE ====> PSEUDO EQUIVALENT POTENTIAL TEMPERATURE (KELVIN)
! OUTPUT ARRAY
! TA =======> PARCEL TEMPERATURE AT EACH LEVEL (CELSIUS DEGREE)
REAL PP0,TP0,PC,TC,THSE,THSE1,THSE2,P_L,T1,T2,RD,CPD
INTEGER KPC
DIMENSION TA(110)
RD=287.04
CPD=1005.
DO K=1,110
TA(K)=-130.
ENDDO
!C ************** DRY ADIABATIC PROCESS BEGIN
KPC=INT(PP0/10)-INT(PC/10)
P_L=INT(PP0/10.)*10.
THD=(TP0+273.15)*(1000./PP0)**(RD/CPD)
DO I=1,KPC
TA(I)=THD*(P_L/1000.)**(RD/CPD)-273.15
P_L=P_L-10.
ENDDO
! DO I=1,KPC
! TA(I)=TP0+(TC-TP0)/ALOG(PC/PP0)*ALOG(P_L/PP0)
! P_L=P_L-10.
! ENDDO
! ************** DRY ADIABATIC PROCESS END
IN=KPC
T1=TC
P1=PC
! WRITE(*,*)' IN = ',IN,' TC = ',TC,' PC = ',PC
KZERO=1
! KZERO IS A JUDGEMENT FOR TEMPERATURE
!*************** MOIST ADIABATIC PROCESS BEGIN
10 CALL CAL_THSE(P_L,T1,THSE1)
CALL CAL_THSE(P_L,T1-0.1,THSE2)
T2=T1-(THSE1-THSE)/(THSE1-THSE2)*0.1
! NOTICE DIFFERENCE OF THSE WHEN TEMPERERTURE OBOVE OR BELEOW ZERO
IF(T2.LT.0.AND.KZERO.EQ.1)THEN
CALL CAL_THSE(P_L,T2,THSE)
KZERO=0
ENDIF
IN=IN+1
TA(IN)=T2
T1=T2
P_L=P_L-10.
! WRITE(*,*)' T1 = ',T1
! PAUSE
IF(P_L.GE.90.)GOTO 10
! *************** MOIST ADIABATIC PROCESS END
RETURN
END
SUBROUTINE CAL_PLFC_AND_PE(P,TE,TA,LEV,PLFC,PE)
! THIS SUBROUTINE CALCULATE LAYER OF FREE CONVECTION AND EMBALANCE LAYER
! NOTICE THIS SUBROUTINE IS NOT INDEPENDENT!!!
! INPUT TE =====> ENVIRONMENT TEMPERATURE (KELVIN)
! TA =====> PARCEL TEMPERATURE (KELVIN)
! P =====> PRESSURE (HPA)
! LEV ====> NUMBER OF LAYER
! OUTPUT PLFC ===> FREE CONVECTION LAYER (HPA)
! PE =====> EMBALANCE PRESSURE (HPA)
REAL PLFC,PE
DIMENSION TE(110),TA(110),P(LEV)
INTEGER LEV,KLFC
REAL DELT1,DELT2
! WRITE(*,'(10F8.2)')(TA(K),K=1,10),(TE(K),K=1,10)
KLFC=0
! IF(TA(1).GE.TE(1).AND.TA(2).GT.TE(2))THEN
IF(TA(1).GT.TE(1).OR.(TA(1).EQ.TE(1).AND.TA(2).GT.TE(2)))THEN
PLFC=P(1)
ELSEIF((TA(1).EQ.TE(1).AND.TA(2).LE.TE(2)).OR.TA(2).LE.TE(2))THEN
P_L=INT(P(1)/10.)*10
P_T=P_L-10
DO K=2,90
IF(TA(K).GT.TE(K))THEN
KLFC=K
P_B=P_T+10
DELT1=TE(K-1)-TA(K-1)
DELT2=TA(K)-TE(K)
PLFC=(P_B*DELT2+P_T*DELT1)/(DELT1+DELT2)
GOTO 290
ENDIF
P_T=P_T-10
ENDDO
IF(KLFC.EQ.0)PLFC=99999.9
ENDIF
290 CONTINUE
NB=INT(P(1)/10.)
NT=INT((P(LEV)-0.001)/10.)+1
NE=NB-NT+1
KPE=NE
IF(INT(P(LEV))/10..EQ.P(LEV)/10.)THEN
P_T=P(LEV)
ELSE
P_T=INT(P(LEV)/10.)*10+10.
ENDIF
IF(TE(NE).LT.TA(NE))THEN
PE=99999.9
ELSE
DO K=NE-1,1,-1
IF(TE(K).LT.TA(K))THEN
KPE=K
P_B=P_T+10
DELT1=TE(K+1)-TA(K+1)
DELT2=TA(K)-TE(K)
PE=(P_B*DELT2+P_T*DELT1)/(DELT1+DELT2)
GOTO 390
ENDIF
P_T=P_T+10
ENDDO
IF(KPE.EQ.NE)PE=99999.9
ENDIF
390 CONTINUE
RETURN
END
SUBROUTINE CAL_VIR_ENGY(P,T,TD,LEV,TE,TA,CAPEV)
!
! NOTICE IT IS SUPPORTED BY CAL_TDE,CAL_Q AND CAL_QS
! INPUT T ======> TEMPERATURE (CELSIUS DEGREE)
! TD =====> DEW POINT (CELSIUS DEGREE)
! P ======> PRESSURE (HPA)
! LEV ====> LEVEL
! TE =====> ENVIRONMENT TEMPERATURE (CELSIUS DEGREE)
! TA =====> PARCEL TEMPERATURE (CELSIUS DEGREE)
! OUTPUT
! CAPEV ==> CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG)
! TEMP ARRAY
! PDS ====> PRESSURE (HPA)
! QE =====> SPECIFIC HUMICITY
! QSE ====> SATURATE SPECIFIC HUMIDITY
! TDE ====> DEW POINT (CELSIUS DEGREE)
REAL CAPEV,R,P_L
DIMENSION TE(110),TA(110),TDE(110),P(LEV),T(LEV),TD(LEV)
DIMENSION PDS(110),QE(110),QSE(110)
DO K=1,LEV
IF(TD(K).GT.90)TD(K)=T(K)-30.
ENDDO
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)
ELSEIF(TA(K+2).GT.TE(K+2))THEN
RETURN
ELSEIF(P_L.LT.300.)THEN
CIN=99999.9
RETURN
ENDIF
P_L=P_T
P_T=P_T-10.
ENDDO
ENDIF
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 RATION
! 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)
! 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
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 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)
! 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)
! SHR2=SQRT((U0-U(1))**2+(V0-V(1))**2)/(H0/1000)
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
! shr(i,j)=shr2
! 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)
RETURN
END