WebSVN – lagranto.um – Diff – /tags/2.1/goodies/rot2geo.f


      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
 
 
 
 

Rev 3	Rev 12
1	`PROGRAM rot2geo`	1	`PROGRAM rot2geo`
2		2
3	`c ***********************************************************************`	3	`c ***********************************************************************`
4	`c * Rotate lat/lon coordinates *`	4	`c * Rotate lat/lon coordinates *`
5	`c * Michael Sprenger / Spring, summer 2010 *`	5	`c * Michael Sprenger / Spring, summer 2010 *`
6	`c ***********************************************************************`	6	`c ***********************************************************************`
7		7
8	`implicit none`	8	`implicit none`
9		9
10	`c ----------------------------------------------------------------------`	10	`c ----------------------------------------------------------------------`
11	`c Declaration of variables`	11	`c Declaration of variables`
12	`c ----------------------------------------------------------------------`	12	`c ----------------------------------------------------------------------`
13		13
14	`c Mode for geo2rot`	14	`c Mode for geo2rot`
15	`character*80 mode`	15	`character*80 mode`
16		16
17	`c Input and output format for trajectories (see iotra.f)`	17	`c Input and output format for trajectories (see iotra.f)`
18	`character*80 inpfile ! Input filename`	18	`character*80 inpfile ! Input filename`
19	`character*80 outfile ! Output filename`	19	`character*80 outfile ! Output filename`
20		20
21	`c Trajectories`	21	`c Trajectories`
22	`integer ntra ! Number of trajectories`	22	`integer ntra ! Number of trajectories`
23	`integer ntim ! Number of times`	23	`integer ntim ! Number of times`
24	`integer ncol ! Number of columns`	24	`integer ncol ! Number of columns`
25	`real,allocatable, dimension (:,:,:) :: tra ! Trajectories (ntra,ntim,ncol)`	25	`real,allocatable, dimension (:,:,:) :: tra ! Trajectories (ntra,ntim,ncol)`
26	`character*80 vars(100) ! Variable names`	26	`character*80 vars(100) ! Variable names`
27	`integer refdate(6) ! Reference date`	27	`integer refdate(6) ! Reference date`
28		28
29	`c Pole position`	29	`c Pole position`
30	`real pollon ! Rotated pole longitude`	30	`real pollon ! Rotated pole longitude`
31	`real pollat ! Rotated pole latitude`	31	`real pollat ! Rotated pole latitude`
32		32
33	`c Auxiliary variables`	33	`c Auxiliary variables`
34	`integer inpmode`	34	`integer inpmode`
35	`integer outmode`	35	`integer outmode`
36	`integer stat`	36	`integer stat`
37	`integer fid`	37	`integer fid`
38	`integer i,j,k`	38	`integer i,j,k`
39	`real lat,lon,rlat,rlon`	39	`real lat,lon,rlat,rlon`
40		40
41	`c Numerical epsilon`	41	`c Numerical epsilon`
42	`real eps`	42	`real eps`
43	`parameter (eps=0.0001)`	43	`parameter (eps=0.0001)`
44		44
45	`c Externals`	45	`c Externals`
46	`real lmtolms`	46	`real lmtolms`
47	`real phtophs`	47	`real phtophs`
48	`real lmstolm`	48	`real lmstolm`
49	`real phstoph`	49	`real phstoph`
50	`external lmtolms,phtophs`	50	`external lmtolms,phtophs`
51	`external lmstolm,phstoph`	51	`external lmstolm,phstoph`
52		52
53	`c ----------------------------------------------------------------------`	53	`c ----------------------------------------------------------------------`
54	`c Read parameters`	54	`c Read parameters`
55	`c ----------------------------------------------------------------------`	55	`c ----------------------------------------------------------------------`
56		56
57	`c Read parameters`	57	`c Read parameters`
58	`open(10,file='rot2geo.param')`	58	`open(10,file='rot2geo.param')`
59	`read(10,*) mode`	59	`read(10,*) mode`
60		60
61	`if ( mode.eq.'trajectory' ) then`	61	`if ( mode.eq.'trajectory' ) then`
62	`read(10,*) inpfile`	62	`read(10,*) inpfile`
63	`read(10,*) outfile`	63	`read(10,*) outfile`
64	`read(10,*) ntra,ntim,ncol`	64	`read(10,*) ntra,ntim,ncol`
65		65
66	`elseif ( mode.eq.'position' ) then`	66	`elseif ( mode.eq.'position' ) then`
67	`read(10,*) rlon`	67	`read(10,*) rlon`
68	`read(10,*) rlat`	68	`read(10,*) rlat`
69		69
70	`endif`	70	`endif`
71		71
72	`read(10,*) pollon,pollat`	72	`read(10,*) pollon,pollat`
73		73
74	`close(10)`	74	`close(10)`
75		75
76	`c ----------------------------------------------------------------------`	76	`c ----------------------------------------------------------------------`
77	`c Rotate a single position`	77	`c Rotate a single position`
78	`c ----------------------------------------------------------------------`	78	`c ----------------------------------------------------------------------`
79		79
80	`if ( mode.ne.'position') goto 100`	80	`if ( mode.ne.'position') goto 100`
81		81
82	`if ( abs(pollat-90.).gt.eps) then`	82	`if ( abs(pollat-90.).gt.eps) then`
83	`lon = lmstolm(rlat,rlon,pollat,pollon)`	83	`lon = lmstolm(rlat,rlon,pollat,pollon)`
84	`lat = phstoph(rlat,rlon,pollat,pollon)`	84	`lat = phstoph(rlat,rlon,pollat,pollon)`
85	`else`	85	`else`
86	`lon = rlon`	86	`lon = rlon`
87	`lat = rlat`	87	`lat = rlat`
88	`endif`	88	`endif`
89		89
90	`print*,lon,lat`	90	`print*,lon,lat`
91		91
92	`100 continue`	92	`100 continue`
93		93
94	`c ----------------------------------------------------------------------`	94	`c ----------------------------------------------------------------------`
95	`c Rotate a trajectory file`	95	`c Rotate a trajectory file`
96	`c ----------------------------------------------------------------------`	96	`c ----------------------------------------------------------------------`
97		97
98	`if ( mode.ne.'trajectory') goto 200`	98	`if ( mode.ne.'trajectory') goto 200`
99		99
100	`c Determine the formats`	100	`c Determine the formats`
101	`call mode_tra(inpmode,inpfile)`	101	`call mode_tra(inpmode,inpfile)`
102	`if (inpmode.eq.-1) inpmode=1`	102	`if (inpmode.eq.-1) inpmode=1`
103	`call mode_tra(outmode,outfile)`	103	`call mode_tra(outmode,outfile)`
104	`if (outmode.eq.-1) outmode=1`	104	`if (outmode.eq.-1) outmode=1`
105		105
106	`c Allocate memory`	106	`c Allocate memory`
107	`allocate(tra(ntra,ntim,ncol),stat=stat)`	107	`allocate(tra(ntra,ntim,ncol),stat=stat)`
108	`if (stat.ne.0) print,' error allocating array tra *'`	108	`if (stat.ne.0) print,' error allocating array tra *'`
109		109
110	`c Read inpufile`	110	`c Read inpufile`
111	`call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)`	111	`call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)`
112	`call read_tra (fid,tra,ntra,ntim,ncol,inpmode)`	112	`call read_tra (fid,tra,ntra,ntim,ncol,inpmode)`
113	`call close_tra(fid,inpmode)`	113	`call close_tra(fid,inpmode)`
114		114
115	`c Do the coordinate rotation`	115	`c Do the coordinate rotation`
116	`do i=1,ntra`	116	`do i=1,ntra`
117	`do j=1,ntim`	117	`do j=1,ntim`
118		118
119	`rlon = tra(i,j,2)`	119	`rlon = tra(i,j,2)`
120	`rlat = tra(i,j,3)`	120	`rlat = tra(i,j,3)`
121		121
122	`if ( abs(pollat-90.).gt.eps) then`	122	`if ( abs(pollat-90.).gt.eps) then`
123	`lon = lmstolm(rlat,rlon,pollat,pollon)`	123	`lon = lmstolm(rlat,rlon,pollat,pollon)`
124	`lat = phstoph(rlat,rlon,pollat,pollon)`	124	`lat = phstoph(rlat,rlon,pollat,pollon)`
125	`else`	125	`else`
126	`lon = rlon`	126	`lon = rlon`
127	`lat = rlat`	127	`lat = rlat`
128	`endif`	128	`endif`
129		129
130	`tra(i,j,2) = lon`	130	`tra(i,j,2) = lon`
131	`tra(i,j,3) = lat`	131	`tra(i,j,3) = lat`
132		132
133	`enddo`	133	`enddo`
134	`enddo`	134	`enddo`
135		135
136	`c Write output file`	136	`c Write output file`
137	`call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)`	137	`call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)`
138	`call write_tra(fid,tra,ntra,ntim,ncol,outmode)`	138	`call write_tra(fid,tra,ntra,ntim,ncol,outmode)`
139	`call close_tra(fid,outmode)`	139	`call close_tra(fid,outmode)`
140		140
141	`200 continue`	141	`200 continue`
142		142
143	`end`	143	`end`
144		144
145		145
146		146
147		147

Subversion Repositories lagranto.um

(root)/tags/2.1/goodies/rot2geo.f – Rev 3 → 12