energy.for

计算气象场中扰动能量和rmse随时间演变

      program main
	parameter (m=81,n=71,z=22,tl=17,int=3)
	real u1(m,n,z,tl),v1(m,n,z,tl),t1(m,n,z,tl),u2(m,n,z,tl),
     &     v2(m,n,z,tl),t2(m,n,z,tl),de(m,n,z,tl),temp(m,n),tde(z,tl)
    	character time*2,day*8,dend1*8,dend2*8,dend3*8
	day='20030703'
	dend1='20030704'
	dend2='20030705'
	dend3='20030706'
      
	open (10,file='d:\energy\0703\de.dat',form='unformatted',
     &      access='direct',recl=z)
	
	do itl=1,tl
	   it=(itl-1)*int
	   if (it.ge.24 .and. it.lt.48) then
	      it=it-24
	      day=dend1
	    else if (it.ge.48 .and. it.lt.72) then
            it=it-48
	      day=dend2
	    else if (it.eq.72) then
	      it=it-72
	      day=dend3
	   end if
      
      call n2a(it,time)
      print *, time,itl 
      open (20,file='d:\outlinux\prornot\0703con48h\'
     &     //'H'//day//time//'.dat',form='unformatted',
     &     access='direct',recl=m*n)
      open (30,file='d:\outlinux\prornot\0703ran1h48\'
     &     //'H'//day//time//'.dat',form='unformatted',
     &     access='direct',recl=m*n)
      
	input=20
      iru=79
	irv=104
	irt=29
	call rdata(u1,temp,m,n,z,tl,itl,iru,input)
      call rdata(v1,temp,m,n,z,tl,itl,irv,input)
	call rdata(t1,temp,m,n,z,tl,itl,irt,input)
	input=30
	iru=79
	irv=104
	irt=29
	call rdata(u2,temp,m,n,z,tl,itl,iru,input)
	call rdata(v2,temp,m,n,z,tl,itl,irv,input)
	call rdata(t2,temp,m,n,z,tl,itl,irt,input)
	end do !t loop is over
      close (20)
	close (30)
      
	do itl=1,tl	
	 do k=1,z
	  do j=11,51    !here i,j are used to choose the calculated region
	   do i=31,61
	      de(i,j,k,itl)=(u2(i,j,k,itl)-u1(i,j,k,itl))**2+
     &		  (v2(i,j,k,itl)-v1(i,j,k,itl))**2+
     &	      (1004/287)*((t2(i,j,k,itl)-t1(i,j,k,itl))**2)
	      tde(k,itl)=tde(k,itl)+de(i,j,k,itl)
	   end do
	  end do
        end do
	end do


      ir=1
	do itl=1,tl
	   write(10,rec=ir) (tde(k,itl),k=1,z)
	   ir=ir+1
	end do
	
	stop
	end program main
	
	subroutine n2a(num,tim)
	integer num
	character a*1,b*1,tim*2
	in=num/10          !higher order num
	im=mod(num,10)     !lower order num
      a=char(in+48)
	b=char(im+48)
	tim=a//b
	return
	end subroutine n2a
	
	subroutine rdata(var,temp,m,n,z,tl,itl,irec,input)
	real var(m,n,z,tl),temp(m,n)
	integer it,irec,input
	 do k=1,z 
	    read(input,rec=irec) ((temp(i,j),i=1,m),j=1,n)
	   do j=1,n
	    do i=1,m
		 var(i,j,k,itl)=temp(i,j)
	   end do	   
	   end do
	   irec=irec+1
	end do 
	end subroutine rdata