Blame | Last modification | View Log | Download | RSS feed
PROGRAM rot2geo
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='rot2geo.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,*) rlon
read(10,*) rlat
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
lon = lmstolm(rlat,rlon,pollat,pollon)
lat = phstoph(rlat,rlon,pollat,pollon)
else
lon = rlon
lat = rlat
endif
print*,lon,lat
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
rlon = tra(i,j,2)
rlat = tra(i,j,3)
if ( abs(pollat-90.).gt.eps) then
lon = lmstolm(rlat,rlon,pollat,pollon)
lat = phstoph(rlat,rlon,pollat,pollon)
else
lon = rlon
lat = rlat
endif
tra(i,j,2) = lon
tra(i,j,3) = lat
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