Subversion Repositories lagranto.um

      PROGRAM geo2rot
      
c     ***********************************************************************
c     * Rotate lat/lon coordinates                                          *
c     * Michael Sprenger / Spring, summer 2010                              *
c     ***********************************************************************

      implicit none
      
c     ----------------------------------------------------------------------
c     Declaration of variables
c     ----------------------------------------------------------------------

c     Mode for geo2rot
      character*80                           mode

c     Input and output format for trajectories (see iotra.f)
      character*80                           inpfile     ! Input filename
      character*80                           outfile     ! Output filename

c     Trajectories
      integer                                ntra        ! Number of trajectories
      integer                                ntim        ! Number of times
      integer                                ncol        ! Number of columns
      real,allocatable, dimension (:,:,:) :: tra         ! Trajectories (ntra,ntim,ncol)
      character*80                           vars(100)   ! Variable names
      integer                                refdate(6)  ! Reference date

c     Pole position
      real                                   pollon      ! Rotated pole longitude
      real                                   pollat      ! Rotated pole latitude

c     Auxiliary variables
      integer                                inpmode
      integer                                outmode
      integer                                stat
      integer                                fid
      integer                                i,j,k
      real                                   lat,lon,rlat,rlon

c     Numerical epsilon
      real                                   eps
      parameter                              (eps=0.0001)

c     Externals
      real                                   lmtolms        
      real                                   phtophs    
      real                                   lmstolm
      real                                   phstoph        
      external                               lmtolms,phtophs
      external                               lmstolm,phstoph

c     ----------------------------------------------------------------------
c     Read parameters
c     ----------------------------------------------------------------------

c     Read parameters
      open(10,file='geo2rot.param')
       read(10,*) mode

       if ( mode.eq.'trajectory' ) then
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
         
       elseif ( mode.eq.'position' ) then
          read(10,*) lon
          read(10,*) lat
          
       endif

       read(10,*) pollon,pollat
       
      close(10)

c     ----------------------------------------------------------------------
c     Rotate a single position
c     ----------------------------------------------------------------------
      
      if ( mode.ne.'position') goto 100
      
      if ( abs(pollat-90.).gt.eps) then
         rlon = lmtolms(lat,lon,pollat,pollon)
         rlat = phtophs(lat,lon,pollat,pollon)
      else
         rlon = lon
         rlat = lat
      endif

      print*,rlon,rlat

 100  continue

c     ----------------------------------------------------------------------
c     Rotate a trajectory file
c     ----------------------------------------------------------------------
      
      if ( mode.ne.'trajectory') goto 200

c     Determine the formats
      call mode_tra(inpmode,inpfile)
      if (inpmode.eq.-1) inpmode=1
      call mode_tra(outmode,outfile)
      if (outmode.eq.-1) outmode=1

c     Allocate memory
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra      ***' 

c     Read inpufile
      call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
      call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
      call close_tra(fid,inpmode)

c     Do the coordinate rotation
      do i=1,ntra
         do j=1,ntim

            lon = tra(i,j,2)
            lat = tra(i,j,3)
            
            if ( abs(pollat-90.).gt.eps) then
               rlon = lmtolms(lat,lon,pollat,pollon)
               rlat = phtophs(lat,lon,pollat,pollon)
            else
               rlon = lon
               rlat = lat
            endif
            
            tra(i,j,2) = rlon
            tra(i,j,3) = rlat
            
         enddo
      enddo

c     Write output file
      call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
      call write_tra(fid,tra,ntra,ntim,ncol,outmode)
      call close_tra(fid,outmode)
      
 200  continue

      end