Subversion Repositories lagranto.ecmwf

Rev

Rev 3 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
      PROGRAM trainfo
2
 
3
c     ***********************************************************************
4
c     * Get infos for a trajectory file                                     *
5
c     * Michael Sprenger / Spring, summer 2010                              *
6
c     ***********************************************************************
7
 
8
      implicit none
9
 
10
c     ----------------------------------------------------------------------
11
c     Declaration of variables
12
c     ----------------------------------------------------------------------
13
 
14
c     Input file
15
      character*80                           inpfile     ! Input filename
16
      character*80                           mode        ! Mode
17
 
18
c     Trajectories
19
      integer                                ntra        ! Number of trajectories
20
      integer                                ntim        ! Number of times
21
      integer                                ncol        ! Number of columns
22
      character*80                           vars(100)   ! Variable names
23
      integer                                refdate(6)  ! Reference date
24
      real,allocatable, dimension (:,:,:) :: tra         ! Trajectories (ntra,ntim,ncol)
25
 
26
c     Auxiliary variables
27
      integer                                inpmode
28
      integer                                stat
29
      integer                                fid
30
      integer                                i,j,n
31
      integer                                old(5)
32
      integer                                new(5)
33
      integer                                hour,min
34
      character*20                           datestr
35
      character*120                          str
36
      character*4                            str1
37
      character*2                            str2,str3,str4,str5,str6
5 michaesp 38
      integer                                isok
39
      integer                                nleaving
3 michaesp 40
 
41
c     ----------------------------------------------------------------------
42
c     Do the reformating
43
c     ----------------------------------------------------------------------
44
 
45
c     Read parameters
46
      open(10,file='trainfo.param')
47
       read(10,*) inpfile
48
       read(10,*) mode
49
      close(10)
50
 
51
c     Determine the formats
52
      call mode_tra(inpmode,inpfile)
53
      if (inpmode.eq.-1) inpmode=1
54
 
55
c     Get the dimension of the trajectory file
56
      call info_tra(inpfile,ntra,ntim,ncol,inpmode)
57
 
58
c     Get haeder information
59
      call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
60
      call close_tra(fid,inpmode)
61
 
62
c     Write dimensions
63
      if ( (mode.eq.'dim').or.(mode.eq.'all') ) then
64
         print*,ntra,ntim,ncol
65
      endif
66
 
67
c     Write single dimensions
68
      if ( (mode.eq.'ntra').or.(mode.eq.'all') ) then
69
         print*,ntra
70
      endif
71
      if ( (mode.eq.'ntim').or.(mode.eq.'all') ) then
72
         print*,ntim
73
      endif
74
      if ( (mode.eq.'ncol').or.(mode.eq.'all') ) then
75
         print*,ncol
76
      endif
77
 
78
c     Write variable names
79
      if ( (mode.eq.'vars').or.(mode.eq.'all') ) then
80
         print*,(trim(vars(i))//' ',i=1,ncol)         
81
      endif
82
 
83
c     Write reference date 
84
      if ( (mode.eq.'refdate').or.(mode.eq.'all') ) then
85
 
86
c        Concatenate date string
87
         min = refdate(5)         
88
         call datestring(datestr,
89
     >             refdate(1),refdate(2),refdate(3),refdate(4) )
90
         if ( min.eq.0 ) then
91
            datestr = trim(datestr)//'00'
92
         elseif (min.lt.10) then
93
            datestr = trim(datestr)//'0'//char(ichar('0')+min)
94
         else
95
            datestr = trim(datestr)//
96
     >                char(ichar('0')+int(min/10))//
97
     >                char(ichar('0')+mod(min,10))
98
         endif
99
 
100
c        Write date string
101
         print*,trim(datestr)
102
 
103
      endif
104
 
105
c     Load all trajectory times if necessary
106
      if ( (mode.eq.'all'      ).or.
107
     >     (mode.eq.'times'    ).or.
108
     >     (mode.eq.'startdate').or.
109
     >     (mode.eq.'enddate'  ) )
110
     >then
111
         allocate(tra(ntra,ntim,ncol),stat=stat)
112
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
113
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
114
         call close_tra(fid,inpmode)
115
      endif
116
 
117
c     Write list of all times
118
      if ( (mode.eq.'times').or.(mode.eq.'all') ) then
119
         write(*,'(100f8.2)') (tra(1,i,1),i=1,ntim)
120
      endif
121
 
122
c     Write time range
123
      if ( (mode.eq.'timerange').or.(mode.eq.'all') ) then
124
         write(*,'(i6)') refdate(6)
125
      endif
126
 
127
c     Write firstdate (reference date + first time)
128
      if ( (mode.eq.'startdate').or.(mode.eq.'all') ) then
129
 
130
c        Set the time shift of first time relative to reference date
131
         hour = int(tra(1,1,1))
132
         min  = mod(nint(100.*tra(1,1,1)),100) + refdate(5)
133
         if (min.gt.60) then
134
            min  = min - 60
135
            hour = hour + 1
136
         endif
137
         if (min.lt.0) then
138
            min  = min + 60
139
            hour = hour - 1
140
         endif
141
         if (min.lt.0) then
142
            min  = min + 60
143
            hour = hour - 1
144
         endif
145
 
146
c        Get new date (hours and minutes)
147
         old(1) = refdate(1)
148
         old(2) = refdate(2)
149
         old(3) = refdate(3)
150
         old(4) = refdate(4)
151
         old(5) = 0
152
         call newdate(old,real(hour),new)
153
 
154
c        Concatenate the date string
155
         call datestring(datestr,
156
     >             new(1),new(2),new(3),new(4) )
157
         if ( min.eq.0 ) then
158
            datestr = trim(datestr)//'00'
159
         elseif (min.lt.10) then
160
            datestr = trim(datestr)//'0'//char(ichar('0')+min)
161
         else
162
            datestr = trim(datestr)//
163
     >                char(ichar('0')+int(min/10))//
164
     >                char(ichar('0')+mod(min,10))
165
         endif
166
 
167
c        Write date string
168
         print*,trim(datestr)
169
 
170
      endif
171
 
172
c     Write enddate (reference date + last time)
173
      if ( (mode.eq.'enddate').or.(mode.eq.'all') ) then
174
 
175
c        Set the time shift of first time relative to reference date
176
         hour = int(tra(1,ntim,1))
177
         min  = mod(nint(100.*tra(1,ntim,1)),100) + refdate(5)
178
         if (min.gt.60) then
179
            min  = min - 60
180
            hour = hour + 1
181
         endif
182
         if (min.lt.0) then
183
            min  = min + 60
184
            hour = hour - 1
185
         endif
186
         if (min.lt.0) then
187
            min  = min + 60
188
            hour = hour - 1
189
         endif
190
 
191
c        Get new date (hours and minutes)
192
         old(1) = refdate(1)
193
         old(2) = refdate(2)
194
         old(3) = refdate(3)
195
         old(4) = refdate(4)
196
         old(5) = 0
197
         call newdate(old,real(hour),new)
198
 
199
c        Concatenate the date string
200
         call datestring(datestr,
201
     >             new(1),new(2),new(3),new(4) )
202
         if ( min.eq.0 ) then
203
            datestr = trim(datestr)//'00'
204
         elseif (min.lt.10) then
205
            datestr = trim(datestr)//'0'//char(ichar('0')+min)
206
         else
207
            datestr = trim(datestr)//
208
     >                char(ichar('0')+int(min/10))//
209
     >                char(ichar('0')+mod(min,10))
210
         endif
211
 
212
c        Write date string
213
         print*,trim(datestr)
214
 
215
      endif
216
 
217
c     Write trajectories to screen
218
      if ( (mode.eq.'list') ) then
219
 
220
c        Read the complete trajectory file
221
         allocate(tra(ntra,ntim,ncol),stat=stat)
222
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
223
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
224
         call close_tra(fid,inpmode)
225
 
226
c        Get the strings for output
227
         write(str1,'(i4)') refdate(1)
228
         write(str2,'(i2)') refdate(2)
229
         write(str3,'(i2)') refdate(3)
230
         write(str4,'(i2)') refdate(4)
231
         write(str5,'(i2)') refdate(5)
232
         if (refdate(2).eq. 0) str2(1:1)='0'
233
         if (refdate(3).eq. 0) str3(1:1)='0'
234
         if (refdate(4).eq. 0) str4(1:1)='0'
235
         if (refdate(5).eq. 0) str5(1:1)='0'
236
         if (refdate(2).lt.10) str2(1:1)='0'
237
         if (refdate(3).lt.10) str3(1:1)='0'
238
         if (refdate(4).lt.10) str4(1:1)='0'
239
         if (refdate(5).lt.10) str5(1:1)='0'
240
 
241
c        Write the time specification
242
         write(*,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)') 
243
     >          'Reference date ',
244
     >           str1,str2,str3,'_',str4,str5,
245
     >          ' / Time range',refdate(6), ' min'
246
         write(*,*)
247
 
248
c        Write variable names
249
         str=''
250
         do i=1,ncol
251
            str=trim(str)//trim(vars(i))
252
         enddo
253
         write(*,'(a6,a9,a8,a6,100a10)') (trim(vars(i)),i=1,ncol)
254
         write(*,'(a6,a9,a8,a6,100a10)') 
255
     >              '------','---------','--------','------',
256
     >              ('----------',i=5,ncol)
257
 
258
         do n=1,ntra
259
            write(*,*)
260
            do i=1,ntim
261
               write(*,'(1f7.2,f9.2,f8.2,i6,100f10.3)') 
262
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
263
     >               nint(tra(n,i,4)),               ! p
264
     >               (tra(n,i,j),j=5,ncol)           ! fields
265
            enddo
266
         enddo
267
 
268
      endif
269
 
5 michaesp 270
c     Get number of trajectories leaving domain
271
      if ( (mode.eq.'leaving') ) then
3 michaesp 272
 
5 michaesp 273
c        Read the complete trajectory file
274
         allocate(tra(ntra,ntim,ncol),stat=stat)
275
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
276
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
277
         call close_tra(fid,inpmode)
278
 
279
c        Determine the number of trajectories leaving domain
280
         do i=1,ntra
281
            isok = 1
282
            do j=1,ntim
283
               if ( tra(i,j,4).lt.0. ) isok = 0
284
            enddo         
285
            if ( isok.eq.0 ) then
286
               nleaving = nleaving + 1
287
            endif
288
         enddo
289
 
290
c        Write output
291
	 print*,nleaving     
292
 
293
      endif
294
 
3 michaesp 295
      end
296
 
297
 
298
 
299
 
300