program gravity
      PARAMETER (kgit=56,igit=64,nb=36)
      include 'com_main.fc'
! Input parameters for egwd
      REAL*8 :: vy1d(kgit)     !v-wind
      REAL*8 :: vx1d(kgit)     !u-wind
      REAL*8 :: temp1d(kgit)   ! 
      REAL*8 :: cp1d(kgit)     ! Cp
      REAL*8 :: hm(kgit)        !geopotential height (in m)
      REAL*8 :: pres(kgit)     !pressure (in Pa)
      REAL*8 :: rho(kgit)      !density
! Output variables for egwd
      REAL*8 :: ut_gwd(kgit)
      REAL*8 :: vt_gwd(kgit)
     &         ,gwh_ir(kgit), gwh_dif(kgit)
     &         ,var_tot(kgit), flux_tot(kgit)
      integer nvar
      integer :: im2, im1, ip1, ip2
      parameter (nvar=3)
      real uvt(igit,nb,kgit,nvar),acc(igit,nb,kgit,6)
      real accsm(igit,nb,kgit,6)
c-------------------- 
      open (72,file='uvt300_Jan_tst.dx',
     +          status='old',form='unformatted',
     +          access='direct', recl=nb*kgit*igit*3)
      open (73,file='grav_DE.dx',
     +       status='unknown',form='unformatted',
     +         access='direct', recl=igit*nb*kgit*6)
      a = 6366197.
        zoben =135000.
      pi=2.*asin(1.)
C------------------------------- (g/sm^3) 
      hh=7000.
	dich0  =1.273e-3
      rhos   =1.273  
      RgSGS=8.31441e7
      RgSI =8314.41
      G    =9.81
      z_48 =135000.  !!!recommended!!!
      dz    = z_48/(48.-0.5)
c--------------------     zunten  is the lowest level of dynamic fields 
c-----     (!!!vertical wind is calculated at ground, not at zunten!!!)
      zunten=dz*0.5
      do k=1,kgit
         rm(k)=28.
         z(k)=zunten/1000.+float(k-1)*dz/1000.
         rou(k) = exp(-z(k)*1000./hh)
      end do      
      do l=101,105
         print*,'l=',l
         read(72,rec=l) uvt
        do k=1,kgit
           do j=1,nb
              do i=1,igit
                 do m=1,3
                   an1(j,k,i,m)=uvt(i,j,k,m)
                 enddo
              enddo
           enddo
        enddo
       do i=1,igit
         do j=1,nb
c	     rhos=1.e5/(RGAS*an1(j,1,i,3))    !rhos depends on local T
           do k=1,kgit
             vx1d(k)=an1(j,k,i,1)
             vy1d(k)=an1(j,k,i,2)
             temp1d(k)=an1(j,k,i,3)

      !AM  Below must be the GEOPOTENTIAL HEIGHT!! We put the log-pressure
	!    height for trial only. Changde this!!! 
	       hm(k)=z(k)*1000.
             cp1d(k)=1005.0
             pres(k)=1.e5*exp(-z(k)/7.)
             rho(k)=rhos*rou(k)
           enddo
	
           call EYGwave(rho, vx1d, vy1d, temp1d, cp1d, hm, 
     &       ut_gwd, vt_gwd, pres, gwh_ir, gwh_dif,
     &       var_tot, flux_tot)
	     
           do k=1,kgit
	           acc(i,j,k,1)=real(ut_gwd(k))*86400.
		   acc(i,j,k,2)=real(vt_gwd(k))*86400.
		   acc(i,j,k,3)=real(gwh_ir(k)+gwh_dif(k))*86400.
		   acc(i,j,k,4)=real(vx1d(k))
		   acc(i,j,k,5)=real(vy1d(k))
		   acc(i,j,k,6)=real(temp1d(k)) 
	   enddo
	    
         enddo
       enddo
c-------------------------------------------------------------------
       do j=1,nb
        do k=1,kgit
         do i=1,igit
           ip1=i+1
           ip2=i+2
           im1=i-1
           im2=i-2
           if(i.eq.1) then
            im1=igit
            im2=igit-1
           end if
           if(i.eq.2)im2=igit
           if(i.eq.igit) then
            ip1=1
            ip2=2
           end if
           if(i.eq.igit-1)ip2=1
           accsm(i,j,k,1)= (acc(im2,j,k,1)+acc(im1,j,k,1)+
     &      acc(i,j,k,1)+acc(ip2,j,k,1)+acc(ip1,j,k,1))/5.  
           accsm(i,j,k,2)= (acc(im2,j,k,2)+acc(im1,j,k,2)+
     &      acc(i,j,k,2)+acc(ip2,j,k,2)+acc(ip1,j,k,2))/5. 
           accsm(i,j,k,3)= (acc(im2,j,k,3)+acc(im1,j,k,3)+
     &      acc(i,j,k,3)+acc(ip2,j,k,3)+acc(ip1,j,k,3))/5. 
           accsm(i,j,k,4)= acc(i,j,k,4) 
           accsm(i,j,k,5)= acc(i,j,k,5)
           accsm(i,j,k,6)= acc(i,j,k,6)  
         enddo
        enddo
       enddo
c-------------------------------------------------------------------
      write(73,rec=l-100) accsm
      enddo !rec
      close (72)
      close (73)
      end