** This program reads in the force files created by Sharon Vadas' ray trace program.
** Fx=force in the zonal direction=forcefxy(i,j,k,it,1).
** Fy=force in the meridional direction=forcefxy(i,j,k,it,2).
** Here, the indices are: i=zonal, j=meridional, k=altitude, and it=time.
** It also reads in the corresponding longitude (lon), latitude (lat), altitude (zzzkm),
** and timeUT (in min from the beginning of Oct 01).  The body forces here only last from
** 18 UT to 27 UT.  
      implicit none
      integer Nxx,Nyy,Nzz,Nt,Nf,iword
      integer i,j,k,it
      character*100 filexyzt,filefxy
!     parameter(Nxx=111,Nyy=55,Nzz=51,Nt=54,Nf=2,iword=1)
      parameter(Nxx=111,Nyy=55,Nzz=51,Nt=54,Nf=2,iword=4)
      integer iunit,Nxxre,Nyyre,Nzzre,Ntre
      real*4 lon(Nxx),lat(Nyy),zzzkm(Nzz),timeUTmn(Nt)
      real*4 forcefxy(Nxx,Nyy,Nzz,Nt,Nf)


!     filexyzt='data/LonLatZT.txt.com'
!     filefxy='data/forcefxy.com'

      filexyzt='LonLatZT.txt.com'
      filefxy='forcefxy.com.bin'

      iunit=10
      call checkparams(filexyzt,iunit,
     *     Nxx,Nyy,Nzz,Nt, Nxxre,Nyyre,Nzzre,Ntre)


** read in each of the files
      CALL readinfiles(lon,lat,zzzkm,timeUTmn,forcefxy,
     *     Nxx,Nyy,Nzz,Nt,Nf,filexyzt,filefxy,iunit,iword)


      i=47
      j=22
      k=23
      it=10
      print *,'At the longitude=',lon(i),', latitude of',lat(j),
     *', altitude of',zzzkm(k),', and time of',timeUTmn(it),'(min),'
      print *,'the thermospheric body force has the following compts:'
      print*,'Force in the zonal direction=',forcefxy(i,j,k,it,1)
      print*,'Force in the meridional direction=',forcefxy(i,j,k,it,2)
      print *,' '



c 77   format(20(1x,1pg12.5))
      end
*///////////////////////**********************//////////////////////////
      SUBROUTINE readinfiles(lon,lat,zzzkm,timeUTmn,forcefxy,
     *     Nxx,Nyy,Nzz,Nt,Nf,filexyzt,filefxy,iunit,iword)
********* global variables **********
      implicit none
      integer Nxx,Nyy,Nzz,Nt,Nf,iunit,iword
      character*100 filexyzt,filefxy
      real*4 lon(Nxx),lat(Nyy),zzzkm(Nzz),timeUTmn(Nt),
     *     forcefxy(Nxx,Nyy,Nzz,Nt,Nf)
********* local variables **********
      integer i,Nxxre,Nyyre,Nzzre,Ntre


      print *,'in readinfiles, reading in files'


      open(unit=iunit,file=filexyzt,status='unknown')
      print *,'opening the file for reading: ',filexyzt
      read(iunit,*) Nxxre,Nyyre,Nzzre,Ntre
      read(iunit,*) lon
      read(iunit,*) lat
      read(iunit,*) zzzkm
      read(iunit,*) timeUTmn
      close(iunit)

      print *,'lon',lon
      print *,'lat',lat
      print *,'zzzkm',zzzkm
      print *,'timeUTmn',timeUTmn


      open(unit=iunit,access='direct',form='unformatted',
     *     file=filefxy,status='unknown',recl=Nxx*iword)
      print *,'opening the file for reading: ',filefxy,' iword=',
     |  iword,' Nxx=',Nxx,' iunit=',iunit

      do i=1,Nf
         CALL read_bin(iunit,(i-1),forcefxy, Nxx,Nyy,Nzz,Nt,Nf, i)
      enddo
      close(iunit)



      RETURN
      END
*///////////////////////**********************//////////////////////////
      SUBROUTINE read_bin(iunit,offset,fxn,nx,ny,nz,nt,Nuv,iu)
** This routines reads in the function fxn in binary form into file
** with unit #=iunit.
********* global variables **********
      implicit none
      integer iunit,offset,nx,ny,nz,nt,Nuv,iu
      real*4 fxn(nx,ny,nz,nt,Nuv)
********* local variables **********
      integer i,j,k,it,irecnum


** write out fxn:
      do it=1,nt
         do k=1,nz
            do j=1,ny
               irecnum= j + (k-1)*ny + (it-1)*ny*nz + offset*ny*nz*nt
               read(iunit,rec=irecnum)( fxn(i,j,k,it,iu),
     *              i=1,nx)
            enddo
         enddo
      enddo
      
      RETURN
      END
*///////////////////////**********************//////////////////////////
      SUBROUTINE checkparams(fileuvwshort,iunit,
     *     Nxx,Nyy,Nzz,Nt, Nxxre,Nyyre,Nzzre,Ntre)
*********global variables **********
      implicit none
      character*100 fileuvwshort
      integer iunit,Nxx,Nyy,Nzz,Nt,Nxxre,Nyyre,Nzzre,Ntre
********* local variables **********
      integer iaj
         
      open(unit=iunit,file=fileuvwshort,status='unknown')
      read(iunit,*) Nxxre,Nyyre,Nzzre,Ntre
      IF(Nxxre.ne.Nxx.or.Nyyre.ne.Nyy.or.Nzzre.ne.Nzz.or.
     *     Ntre.ne.Nt) then
         print *,'params are not equal',Nxx,Nxxre,Nyy,Nyyre,
     *        Nzz,Nzzre,Nt,Ntre
         stop
      ENDIF
      close(iunit)

      print *,'params checked out ok'

      RETURN
      END
*///////////////////////**********************//////////////////////////