indexes.f90

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

!---------------主程序开始
program main
implicit none
integer,parameter :: x=61,y=34,z=21
real,parameter :: T0=273.16
integer i,j,k
real tk(x,y,z),f(x,y,z),h(x,y,z),u(x,y,z),v(x,y,z)    
real citase(x,y,z),q(x,y,z),tdk(x,y,z)     
real cape1(x,y),tti1(x,y),li1(x,y),ki1(x,y),ai1(x,y),ssi1(x,y),shr1(x,y)
real p(z),t(z),td(z),fp(z)
real te(120),ta(120)
data p/1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100/
real CAPE,CIN,LI,H0,p0,shr,SSI,ki,ai,tti
!-----------------------------------------------
!-------------------------从源文件读入数据
open(1,file='D:\learn\work\diagnosis\myNeed\0127_00\0127_00.dat',status='unknown',form='unformatted',access='direct',recl=x*y)
do k=1,z
!read(1,rec=k)  (((tk(i,j,k),i=1,x),j=1,y),k=1,z)
	read(1,rec=k)  ((tk(i,j,k),i=1,x),j=1,y)
	read(1,rec=z+k)  ((h(i,j,k),i=1,x),j=1,y)
	read(1,rec=2*z+k)  ((u(i,j,k),i=1,x),j=1,y)
	read(1,rec=3*z+k)  ((v(i,j,k),i=1,x),j=1,y)
	read(1,rec=4*z+k)  ((f(i,j,k),i=1,x),j=1,y)
enddo
read(1,rec=5*z+2)  ((cape1(i,j),i=1,x),j=1,y)
!read(1,rec=5*z+6)  ((shr1(i,j),i=1,x),j=1,y)
close(1)
!-----------------------------------计算 thse & td
do i=1,x
	do j=1,y
		do k=1,z
			call tem_se(f(i,j,k),tk(i,j,k),p(k),tdk(i,j,k),q(i,j,k),citase(i,j,k))
		enddo
!		print*, tk(i,j,10),h(i,j,10),tdk(i,j,10),cape1(i,j),shr1(i,j)
!		pause
	enddo
enddo
!--------------------calculate easy indexes
do i=1,x
	do j=1,y
		do k=1,z
		t(k)=tk(i,j,k)-T0
		td(k)=tdk(i,j,k)-T0
		enddo
!	print*, 't=',t,'td=',td	
	H0=3600
	p0=1000
	cape=cape1(i,j)
!	shr=shr1(i,j)
	call CAL_SSI(H0,h(i,j,:),U(i,j,:),V(i,j,:),z,CAPE,shr,SSI)
!	SSI=100*(2+0.276*ALOG(abs(SHR))+2.011*CAPE*0.0001)
	call cal_tti(t,td,tti)
	call cal_ki(t,td,ki)
	call cal_ai(t,td,ai)
	ssi1(i,j)=ssi
	tti1(i,j)=tti
	ki1(i,j)=ki
	ai1(i,j)=ai
	print*, 'cape=',cape1(i,j),'ssi=',ssi1(i,j),'tti=',tti1(i,j),'ki=',ki1(i,j),'ai=',ai1(i,j)
!	pause
	enddo
enddo
!---------------把诊断量写入目标文件
open(3,file='D:\learn\work\diagnosis\myNeed\0127_00\index0127.dat',status='replace',form='unformatted',access='direct',recl=x*y)
write(3,rec=1)	((cape1(i,j),i=1,x),j=1,y)
write(3,rec=2)	((ssi1(i,j),i=1,x),j=1,y)
write(3,rec=3)	((tti1(i,j),i=1,x),j=1,y)
write(3,rec=4)	((ki1(i,j),i=1,x),j=1,y)
write(3,rec=5)	((ai1(i,j),i=1,x),j=1,y)
do k=1,z
	write(3,rec=5+k) ((citase(i,j,k),i=1,x),j=1,y)
	write(3,rec=5+z+k) ((q(i,j,k),i=1,x),j=1,y)
enddo
close(3)
stop 
end program
!---------------计算假相当位温和露点的子程序
SUBROUTINE tem_se(f,t,p,td,q,se)
      T0=273.16
      d=4.9283
      c=6764.9
      ET=6.11*((T0/t)**d)*exp(c/T0-c/t)
      qs=0.622*ET/(p-0.378*ET)
      q=qs*f/100
             call tem_td(q,t,p,td)
                c=0.28586*alog((1000.0/p))
               d=q/(338.52-0.24*t+1.24*td)
               se=t*exp((c+2500*d))
return
end
SUBROUTINE tem_td(q,t,p,td)
      real q,t,p,td,qs,ET
      td=t
        qs=q+1
        do 20, while(q.le.qs)
        td=td-0.02
        if (td<200) goto 10
      T0=273.16
        d=4.9283
        c=6764.9
      ET=6.11*((T0/td)**d)*exp(c/T0-c/td)
      qs=0.622*ET/(p-0.378*ET)
20     continue
10      return
 end
!-----------------------主程序结束
!------------------计算修正的全总指数、K指数、A指数
subroutine cal_tti(t,td,tti)
integer,parameter :: x=101,y=81,z=21
real t(z),td(z)
real tave,tdave,tti
tave=0.
tdave=0.
do i=1,6
	tave=tave+t(i)
	tdave=tdave+t(i)
enddo
tave=tave/6
tdave=tdave/6
tti=tave+tdave-2*t(13)
return
end
subroutine cal_ki(t,td,ki)
integer,parameter :: x=101,y=81,z=21
real t(z),td(z)
real ki
ki=(t(1)+t(6))/2+(td(1)+td(6))/2-t(13)-(t(9)-td(9))
return
end
subroutine cal_ai(t,td,ai)
integer,parameter :: x=101,y=81,z=21
real t(z),td(z)
real ai
ai=t(6)-t(13)-((t(6)-td(6))+(t(9)-td(9))+(t(13)-td(13)))
return
end