cape1.for

此程序用于计算单站对流有效位能CAPE值。

**********************************************************************
****************单站探空资料计算cape(对流有效位能)*******************************
*********************************************************************
	parameter(irec=522)
      real pc(irec),tc(irec),td(irec),pz(irec),tz(irec),tvi(irec),
     $	 tve(irec),cape,e1(irec),e2(irec)
		  
        !      Rd=287J/kg*k
         rd=287
         cape=0
 	   aa1=7.5
	   bb1=237.7
	   aa2=9.5
	   bb2=265.5


	open(1,file='d:\exercise\ztjifenstep070808.txt')
	open(2,file='d:\exercise\cjjifenstep070808.txt')
	
      read(1,*)
	do ii=1, irec
	read(1,*) pz(ii),tz(ii)
	enddo

	! 读入对流有效位能状态曲线积分步长
	 
	read(2,*)
	do ii=1, irec
	read(2,*) pc(ii),tc(ii),td(ii)
	enddo

	! 读入对流有效位能层结曲线积分步长
 

      do ii=1,irec-1
	if (td(ii).gt.0)  then
	   e1(ii)=6.1078*exp(aa1*td(ii)/(td(ii)+bb1))    
      	else
	   e1(ii)=6.1078*exp(aa2*td(ii)/(td(ii)+bb2))
	endif
	
	if (tz(ii).gt.0) then
	   e2(ii)=6.1078*exp(aa1*tz(ii)/(tz(ii)+bb1))
	else
	   e2(ii)=6.1078*exp(aa2*tz(ii)/(tz(ii)+bb2))
	endif
	 !  tc、tz分别为环境、气块的温度
	 !	e1、e2分别为的环境、气块水汽压

         tve(ii)=(tc(ii)+273.16)*(1+0.378*e1(ii)/pc(ii))
         tvi(ii)=(tz(ii)+273.16)*(1+0.378*e2(ii)/pz(ii))
	 !  tve、tvi分别为环境和气块的虚温
	 
         cape=cape+((-1)*rd*(tvi(ii)-tve(ii))
     $            *log((pz(ii+1)/pz(ii))))

	 !  积分计算对流有效位能 
	
      enddo
	write(*,*) 'cape',cape
      end