Subversion Repositories lagranto.um

Rev

Rev 3 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3 Rev 10
1
c     ****************************************************************
1
c     ****************************************************************
2
c     * This package provides IO routines for trajectories. A file   *
2
c     * This package provides IO routines for trajectories. A file   *
3
c     * is characterised by the filename <filename> and the file     *
3
c     * is characterised by the filename <filename> and the file     *
4
c     * identifier <fid>. Different modes <mode> are supported:      *
4
c     * identifier <fid>. Different modes <mode> are supported:      *
5
c     *     mode=1: ascii, sorted by trajectory;                     *
5
c     *     mode=1: ascii, sorted by trajectory;                     *
6
c     *     mode=2: ascii, sorted by time;                           *
6
c     *     mode=2: ascii, sorted by time;                           *
7
c     *     mode=3: fortran (unformatted)                            * 
7
c     *     mode=3: fortran (unformatted)                            * 
8
c     *     mode=4: IVE netcdf (for compatibiltzy reasons)           *
8
c     *     mode=4: IVE netcdf (for compatibiltzy reasons)           *
9
c     * A trajectory set is given as 3d array <tra(ntra,ntim,ncol)>  *
9
c     * A trajectory set is given as 3d array <tra(ntra,ntim,ncol)>  *
10
c     * where <ntra> is the number of trajectories, <ntim> the       *
10
c     * where <ntra> is the number of trajectories, <ntim> the       *
11
c     * number of times of each trajectory and <ncol> the number of  *
11
c     * number of times of each trajectory and <ncol> the number of  *
12
c     * columns of the trajectory. The first 4 columns are: time,    *
12
c     * columns of the trajectory. The first 4 columns are: time,    *
13
c     * longitude, latitude, pressure. The other columns are traced  *
13
c     * longitude, latitude, pressure. The other columns are traced  *
14
c     * fields. The complete list of all columns is given in the     *
14
c     * fields. The complete list of all columns is given in the     *
15
c     * array <vars(ncol)>. Finally, the reference date is given in  *
15
c     * array <vars(ncol)>. Finally, the reference date is given in  *
16
c     * the array <time(6)=year,month,day,hour,time length of the    *
16
c     * the array <time(6)=year,month,day,hour,time length of the    *
17
c     * trajectory (hour,min)>.                                      *
17
c     * trajectory (hour,min)>.                                      *
18
c     *                                                              *
18
c     *                                                              *
19
c     * Author: Michael Sprenger, September 2008                     *
19
c     * Author: Michael Sprenger, September 2008                     *
20
c     ****************************************************************
20
c     ****************************************************************
21
 
21
 
22
c     ----------------------------------------------------------------
22
c     ----------------------------------------------------------------
23
c     Open a trajectory file for reading
23
c     Open a trajectory file for reading
24
c     ----------------------------------------------------------------
24
c     ----------------------------------------------------------------
25
      
25
      
26
      subroutine ropen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)
26
      subroutine ropen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)
27
 
27
 
28
      implicit none
28
      implicit none
29
      
29
      
30
c     Declaration of subroutine parameters
30
c     Declaration of subroutine parameters
31
      integer      fid
31
      integer      fid
32
      character*80 filename
32
      character*80 filename
33
      integer      mode
33
      integer      mode
34
      integer      ntra,ntim,ncol
34
      integer      ntra,ntim,ncol
35
      integer      time(6)
35
      integer      time(6)
36
      character*80 vars(ncol)
36
      character*80 vars(ncol)
37
 
37
 
38
c     Auxiliary variables
38
c     Auxiliary variables
39
      integer      vardim(4)
39
      integer      vardim(4)
40
      real         varmin(4),varmax(4),stag(4)
40
      real         varmin(4),varmax(4),stag(4)
41
      real         mdv
41
      real         mdv
42
      character*80 cfn
42
      character*80 cfn
43
      integer      ierr
43
      integer      ierr
44
      integer      i
44
      integer      i
45
      integer      nvars
45
      integer      nvars
46
 
46
 
47
c     Open file
47
c     Open file
48
      if (mode.eq.1) then
48
      if (mode.eq.1) then
49
         fid = 10
49
         fid = 10
50
         open(fid,file=filename)
50
         open(fid,file=filename)
51
      elseif (mode.eq.2) then
51
      elseif (mode.eq.2) then
52
         fid = 10
52
         fid = 10
53
         open(fid,file=filename)
53
         open(fid,file=filename)
54
      elseif (mode.eq.3) then
54
      elseif (mode.eq.3) then
55
         open(fid,file=filename,form='unformatted')
55
         open(fid,file=filename,form='unformatted')
56
      elseif (mode.eq.4) then
56
      elseif (mode.eq.4) then
57
         call cdfopn(filename,fid,ierr)
57
         call cdfopn(filename,fid,ierr)
58
      endif
58
      endif
59
 
59
 
60
c     Read header information
60
c     Read header information
61
      call read_hea(fid,time,vars,ntra,ntim,ncol,mode)
61
      call read_hea(fid,time,vars,ntra,ntim,ncol,mode)
62
 
62
 
63
      end
63
      end
64
      
64
      
65
 
65
 
66
c     ----------------------------------------------------------------
66
c     ----------------------------------------------------------------
67
c     Open a trajectory file for wrinting
67
c     Open a trajectory file for wrinting
68
c     ----------------------------------------------------------------
68
c     ----------------------------------------------------------------
69
      
69
      
70
      subroutine wopen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)
70
      subroutine wopen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)
71
      
71
      
72
      implicit none
72
      implicit none
73
      
73
      
74
c     Declaration of subroutine parameters
74
c     Declaration of subroutine parameters
75
      integer      fid
75
      integer      fid
76
      character*80 filename
76
      character*80 filename
77
      integer      mode
77
      integer      mode
78
      integer      ntra,ntim,ncol
78
      integer      ntra,ntim,ncol
79
      integer      time(6) 
79
      integer      time(6) 
80
      character*80 vars(ncol)
80
      character*80 vars(ncol)
81
 
81
 
82
c     Auxiliary variables
82
c     Auxiliary variables
83
      integer      vardim(4)
83
      integer      vardim(4)
84
      real         varmin(4),varmax(4),stag(4)
84
      real         varmin(4),varmax(4),stag(4)
85
      real         mdv
85
      real         mdv
86
      character*80 cfn
86
      character*80 cfn
87
      integer      ierr
87
      integer      ierr
88
      integer      i
88
      integer      i
89
      character*80 varname
89
      character*80 varname
90
      real         rtime(6)
90
      real         rtime(6)
91
 
91
 
92
c     Open file
92
c     Open file
93
      if (mode.eq.1) then
93
      if (mode.eq.1) then
94
         fid = 10
94
         fid = 10
95
         open(fid,file=filename)
95
         open(fid,file=filename)
96
      elseif (mode.eq.2) then
96
      elseif (mode.eq.2) then
97
         fid = 10
97
         fid = 10
98
         open(fid,file=filename)
98
         open(fid,file=filename)
99
      elseif (mode.eq.3) then
99
      elseif (mode.eq.3) then
100
         open(fid,file=filename,form='unformatted')
100
         open(fid,file=filename,form='unformatted')
101
      elseif (mode.eq.4) then
101
      elseif (mode.eq.4) then
102
         vardim(1)=ntra
102
         vardim(1)=ntra
103
         vardim(2)=1
103
         vardim(2)=1
104
         vardim(3)=1
104
         vardim(3)=1
105
         vardim(4)=1
105
         vardim(4)=1
106
         cfn      =trim(filename)//'_cst'
106
         cfn      =trim(filename)//'_cst'
107
         mdv      =-999.98999
107
         mdv      =-999.98999
108
         call crecdf(filename,fid,varmin,varmax,3,cfn,ierr)
108
         call crecdf(filename,fid,varmin,varmax,3,cfn,ierr)
109
      endif
109
      endif
110
 
110
 
111
c     Write header information
111
c     Write header information
112
      call write_hea(fid,time,vars,ntra,ntim,ncol,mode)
112
      call write_hea(fid,time,vars,ntra,ntim,ncol,mode)
113
 
113
 
114
      end
114
      end
115
 
115
 
116
 
116
 
117
c     ----------------------------------------------------------------
117
c     ----------------------------------------------------------------
118
c     Read a trajectory
118
c     Read a trajectory
119
c     ----------------------------------------------------------------
119
c     ----------------------------------------------------------------
120
 
120
 
121
      subroutine read_tra(fid,tra,ntra,ntim,ncol,mode)
121
      subroutine read_tra(fid,tra,ntra,ntim,ncol,mode)
122
 
122
 
123
      implicit none
123
      implicit none
124
 
124
 
125
c     Declaration of subroutine parameters
125
c     Declaration of subroutine parameters
126
      integer   fid
126
      integer   fid
127
      integer   ntim
127
      integer   ntim
128
      integer   ncol
128
      integer   ncol
129
      integer   ntra
129
      integer   ntra
130
      real      tra(ntra,ntim,ncol)
130
      real      tra(ntra,ntim,ncol)
131
      integer   mode
131
      integer   mode
132
 
132
 
133
c     Auxiliary variables
133
c     Auxiliary variables
134
      integer      i,j,n
134
      integer      i,j,n
135
      real         arr(ntra)
135
      real         arr(ntra)
136
      integer      ntimes
136
      integer      ntimes
137
      real         times(1000)
137
      real         times(1000)
138
      integer      ierr
138
      integer      ierr
139
      character*80 vars(ncol+2)
139
      character*80 vars(ncol+2)
140
      integer      nvars
140
      integer      nvars
141
 
141
 
142
c     Read ascii mode, sorted by trajectory (mode=1)
142
c     Read ascii mode, sorted by trajectory (mode=1)
143
      if (mode.eq.1) then
143
      if (mode.eq.1) then
144
         read(fid,*,end=100)
144
         read(fid,*,end=100)
145
         do n=1,ntra
145
         do n=1,ntra
146
            do i=1,ntim
146
            do i=1,ntim
147
               read(fid,*,end=110) (tra(n,i,j),j=1,ncol)
147
               read(fid,*,end=110) (tra(n,i,j),j=1,ncol)
148
            enddo
148
            enddo
149
         enddo
149
         enddo
150
 
150
 
151
c     Read ascii mode, sorted by time (mode=2)
151
c     Read ascii mode, sorted by time (mode=2)
152
      elseif (mode.eq.2) then
152
      elseif (mode.eq.2) then
153
         read(fid,*,end=100)
153
         read(fid,*,end=100)
154
         do i=1,ntim
154
         do i=1,ntim
155
            do n=1,ntra
155
            do n=1,ntra
156
               read(fid,*,end=100) (tra(n,i,j),j=1,ncol)
156
               read(fid,*,end=100) (tra(n,i,j),j=1,ncol)
157
            enddo
157
            enddo
158
         enddo
158
         enddo
159
 
159
 
160
c     Read fortran mode (mode=3)
160
c     Read fortran mode (mode=3)
161
      elseif (mode.eq.3) then
161
      elseif (mode.eq.3) then
162
         read(fid) tra
162
         read(fid) tra
163
 
163
 
164
c     Read IVE netcdf mode (mode=4)
164
c     Read IVE netcdf mode (mode=4)
165
      elseif (mode.eq.4) then
165
      elseif (mode.eq.4) then
166
          call gettimes(fid,times,ntimes,ierr)
166
          call gettimes(fid,times,ntimes,ierr)
167
          call getvars(fid,nvars,vars,ierr)
167
          call getvars(fid,nvars,vars,ierr)
168
          do i=1,ntim
168
          do i=1,ntim
169
             do j=1,ncol
169
             do j=1,ncol
170
                if (j.eq.1) then
170
                if (j.eq.1) then
171
                   do n=1,ntra
171
                   do n=1,ntra
172
                      tra(n,i,1)=times(i)
172
                      tra(n,i,1)=times(i)
173
                   enddo
173
                   enddo
174
                else
174
                else
175
                   call getdat(fid,vars(j),times(i),0,arr,ierr)
175
                   call getdat(fid,vars(j),times(i),0,arr,ierr)
176
                   do n=1,ntra
176
                   do n=1,ntra
177
                      tra(n,i,j)=arr(n)
177
                      tra(n,i,j)=arr(n)
178
                   enddo
178
                   enddo
179
                endif
179
                endif
180
            enddo
180
            enddo
181
         enddo
181
         enddo
182
 
182
 
183
      endif
183
      endif
184
 
184
 
185
      return
185
      return
186
 
186
 
187
c     End of file has been reached: set negative <fid>
187
c     End of file has been reached: set negative <fid>
188
 100  fid=-fid
188
 100  fid=-fid
189
      return
189
      return
190
 
190
 
191
c     Error: incomplete trajectory
191
c     Error: incomplete trajectory
192
 110  print*,'<read_tra>: Incomplete trajectory... Stop'
192
 110  print*,'<read_tra>: Incomplete trajectory... Stop'
193
      stop
193
      stop
194
      
194
      
195
      end
195
      end
196
 
196
 
197
 
197
 
198
c     ----------------------------------------------------------------
198
c     ----------------------------------------------------------------
199
c     Write a trajectory
199
c     Write a trajectory
200
c     ----------------------------------------------------------------
200
c     ----------------------------------------------------------------
201
 
201
 
202
      subroutine write_tra(fid,tra,ntra,ntim,ncol,mode)
202
      subroutine write_tra(fid,tra,ntra,ntim,ncol,mode)
203
 
203
 
204
      implicit none
204
      implicit none
205
 
205
 
206
c     Declaration of subroutine parameters
206
c     Declaration of subroutine parameters
207
      integer   fid
207
      integer   fid
208
      integer   ntim
208
      integer   ntim
209
      integer   ncol
209
      integer   ncol
210
      integer   ntra
210
      integer   ntra
211
      real      tra(ntra,ntim,ncol)
211
      real      tra(ntra,ntim,ncol)
212
      integer   mode
212
      integer   mode
213
 
213
 
214
c     Auxiliary variables
214
c     Auxiliary variables
215
      integer      i,j,n
215
      integer      i,j,n
216
      real         arr(ntra)
216
      real         arr(ntra)
217
      integer      ierr
217
      integer      ierr
218
      real         time
218
      real         time
219
      character*80 vars(ncol+2)
219
      character*80 vars(ncol+2)
220
      integer      nvars
220
      integer      nvars
221
 
221
 
222
c     Write ascii mode, sorted by trajectory (mode=1)
222
c     Write ascii mode, sorted by trajectory (mode=1)
223
      if (mode.eq.1) then
223
      if (mode.eq.1) then
224
         do n=1,ntra
224
         do n=1,ntra
225
            write(fid,*)
225
            write(fid,*)
226
            do i=1,ntim
226
            do i=1,ntim
227
 
227
 
228
c              Avoid ugly *s or missing space in output
228
c              Avoid ugly *s or missing space in output
229
               do j=5,ncol
229
               do j=5,ncol
230
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
230
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
231
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
231
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
232
                     tra(n,i,j) = -999.99
232
                     tra(n,i,j) = -999.99
233
                  endif
233
                  endif
234
               enddo
234
               enddo
235
 
235
 
236
               write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)') 
236
               write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)') 
237
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
237
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
238
     >               nint(tra(n,i,4)),               ! p
238
     >               nint(tra(n,i,4)),               ! p
239
     >               (tra(n,i,j),j=5,ncol)           ! fields
239
     >               (tra(n,i,j),j=5,ncol)           ! fields
240
            enddo
240
            enddo
241
         enddo
241
         enddo
242
 
242
 
243
c     Write ascii mode, sorted by time (mode=2)
243
c     Write ascii mode, sorted by time (mode=2)
244
      elseif (mode.eq.2) then
244
      elseif (mode.eq.2) then
245
         do i=1,ntim
245
         do i=1,ntim
246
            write(fid,*)
246
            write(fid,*)
247
            do n=1,ntra
247
            do n=1,ntra
248
 
248
 
249
c              Avoid ugly *s or missing space in output
249
c              Avoid ugly *s or missing space in output
250
               do j=5,ncol
250
               do j=5,ncol
251
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
251
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
252
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
252
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
253
                     tra(n,i,j) = -999.99
253
                     tra(n,i,j) = -999.99
254
                  endif
254
                  endif
255
               enddo
255
               enddo
256
 
256
 
257
               write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)') 
257
               write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)') 
258
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
258
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
259
     >               nint(tra(n,i,4)),               ! p
259
     >               nint(tra(n,i,4)),               ! p
260
     >               (tra(n,i,j),j=5,ncol)           ! fields
260
     >               (tra(n,i,j),j=5,ncol)           ! fields
261
            enddo
261
            enddo
262
         enddo
262
         enddo
263
 
263
 
264
c     Write fortran mode (mode=3)
264
c     Write fortran mode (mode=3)
265
      elseif (mode.eq.3) then
265
      elseif (mode.eq.3) then
266
         write(fid) tra                              
266
         write(fid) tra                              
267
 
267
 
268
c     Write netcdf mode (mode=4)
268
c     Write netcdf mode (mode=4)
269
      elseif (mode.eq.4) then
269
      elseif (mode.eq.4) then
270
         call getvars(fid,nvars,vars,ierr)
270
         call getvars(fid,nvars,vars,ierr)
271
         do i=1,ntim
271
         do i=1,ntim
272
            time=tra(1,i,1)
272
            time=tra(1,i,1)
273
            do j=2,ncol
273
            do j=2,ncol
274
               do n=1,ntra
274
               do n=1,ntra
275
                  arr(n)=tra(n,i,j)
275
                  arr(n)=tra(n,i,j)
276
               enddo
276
               enddo
277
               call putdat(fid,vars(j),time,0,arr,ierr)
277
               call putdat(fid,vars(j),time,0,arr,ierr)
278
            enddo
278
            enddo
279
         enddo
279
         enddo
280
      endif
280
      endif
281
 
281
 
282
      end
282
      end
283
 
283
 
284
 
284
 
285
c     ----------------------------------------------------------------
285
c     ----------------------------------------------------------------
286
c     Read header from trajectory file
286
c     Read header from trajectory file
287
c     ----------------------------------------------------------------
287
c     ----------------------------------------------------------------
288
 
288
 
289
      subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)
289
      subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)
290
 
290
 
291
      implicit none
291
      implicit none
292
      
292
      
293
c     Declaration of subroutine parameters
293
c     Declaration of subroutine parameters
294
      integer       fid
294
      integer       fid
295
      integer       time(6)
295
      integer       time(6)
296
      integer       ntra,ntim,ncol
296
      integer       ntra,ntim,ncol
297
      character*80  vars(ncol)
297
      character*80  vars(ncol)
298
      integer       mode
298
      integer       mode
299
 
299
 
300
c     Auxiliary variables
300
c     Auxiliary variables
301
      integer       i
301
      integer       i
302
      character     ch(120)
302
      character     ch(120)
303
      character*120 str
303
      character*120 str
304
      integer       ich(120)
304
      integer       ich(120)
305
      integer       isstr,ileft,iright
305
      integer       isstr,ileft,iright
306
      character*80  varname
306
      character*80  varname
307
      real          rtime(6)
307
      real          rtime(6)
308
      integer       ierr
308
      integer       ierr
309
      integer       nvars
309
      integer       nvars
310
      character*15  str1
310
      character*15  str1
311
      character     str2
311
      character     str2
312
      character*13  str3
312
      character*13  str3
313
      character*4   str4
313
      character*4   str4
314
      character*80  linestr
314
      character*80  linestr
315
      integer       itmp1,itmp2
315
      integer       itmp1,itmp2
316
 
316
 
317
c     Read ascii format (mode=1,2)
317
c     Read ascii format (mode=1,2)
318
      if ( (mode.eq.1).or.(mode.eq.2) ) then
318
      if ( (mode.eq.1).or.(mode.eq.2) ) then
319
 
319
 
320
c        Read the time specification (old and new format)
320
c        Read the time specification (old and new format)
321
         read(fid,'(a80)') linestr
321
         read(fid,'(a80)') linestr
322
         
322
         
323
         if ( linestr(1:14).eq.'Reference date' ) then
323
         if ( linestr(1:14).eq.'Reference date' ) then
324
            read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)') 
324
            read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)') 
325
     >           str1,
325
     >           str1,
326
     >           time(1),time(2),time(3),str2,time(4),time(5),
326
     >           time(1),time(2),time(3),str2,time(4),time(5),
327
     >           str3,time(6),str4
327
     >           str3,time(6),str4
328
                       
328
                       
329
         elseif ( linestr(1:11).eq.'time period' ) then
329
         elseif ( linestr(1:11).eq.'time period' ) then
330
            read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)') 
330
            read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)') 
331
     >           str1,
331
     >           str1,
332
     >           time(1),time(2),time(3),str2,time(4),
332
     >           time(1),time(2),time(3),str2,time(4),
333
     >           str3,itmp1,str3,itmp2,str4
333
     >           str3,itmp1,str3,itmp2,str4
334
            time(5) = 0
334
            time(5) = 0
335
            time(6) = itmp1 * 60 + itmp2
335
            time(6) = itmp1 * 60 + itmp2
336
 
336
 
337
         endif
337
         endif
338
 
338
 
339
c        Skip the empty line and read field names
339
c        Skip the empty line and read field names
340
         read(fid,*)
340
         read(fid,*)
341
         read(fid,'(a120)',end=100) str
341
         read(fid,'(a120)',end=100) str
342
         do i=1,120
342
         do i=1,120
343
            ch(i)=str(i:i)
343
            ch(i)=str(i:i)
344
         enddo
344
         enddo
345
 
345
 
346
c        Split the input string
346
c        Split the input string
347
         isstr=0
347
         isstr=0
348
         nvars=0
348
         nvars=0
349
         do i=1,120
349
         do i=1,120
350
            if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then
350
            if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then
351
               isstr=1
351
               isstr=1
352
               ileft=i
352
               ileft=i
353
            elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then
353
            elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then
354
               isstr=0
354
               isstr=0
355
               iright=i-1
355
               iright=i-1
356
               nvars=nvars+1
356
               nvars=nvars+1
357
               vars(nvars)=str(ileft:iright)
357
               vars(nvars)=str(ileft:iright)
358
            endif
358
            endif
359
         enddo
359
         enddo
360
 
360
 
361
c        Skip the empty line
361
c        Skip the empty line
362
         read(fid,*,end=100)
362
         read(fid,*,end=100)
363
 
363
 
364
c     Read fortran mode (mode=3)
364
c     Read fortran mode (mode=3)
365
      elseif (mode.eq.3) then
365
      elseif (mode.eq.3) then
366
         read(fid) ntra,ntim,ncol
366
         read(fid) ntra,ntim,ncol
367
         read(fid) time
367
         read(fid) time
368
         read(fid) vars
368
         read(fid) vars
369
 
369
 
370
c     Read IVE netcdf mode (mode=4)
370
c     Read IVE netcdf mode (mode=4)
371
      elseif (mode.eq.4) then
371
      elseif (mode.eq.4) then
372
         call getvars(fid,nvars,vars,ierr)  
372
         call getvars(fid,nvars,vars,ierr)  
373
         varname='BASEDATE'
373
         varname='BASEDATE'
374
         call getdat(fid,varname,0.,0,rtime,ierr)
374
         call getdat(fid,varname,0.,0,rtime,ierr)
375
         do i=1,6
375
         do i=1,6
376
            time(i)=nint(rtime(i))
376
            time(i)=nint(rtime(i))
377
         enddo
377
         enddo
378
      endif
378
      endif
379
 
379
 
380
      return
380
      return
381
 
381
 
382
c     End of file has been reached
382
c     End of file has been reached
383
 100  fid=-fid
383
 100  fid=-fid
384
      return
384
      return
385
 
385
 
386
c     Excetion handling
386
c     Excetion handling
387
 110  print*,'<read_hea>: Unexspected time format.... Stop'
387
 110  print*,'<read_hea>: Unexspected time format.... Stop'
388
      stop
388
      stop
389
   
389
   
390
      end
390
      end
391
      
391
      
392
 
392
 
393
c     ----------------------------------------------------------------
393
c     ----------------------------------------------------------------
394
c     Write header to trajectory file (in ascii mode)
394
c     Write header to trajectory file (in ascii mode)
395
c     ----------------------------------------------------------------
395
c     ----------------------------------------------------------------
396
 
396
 
397
      subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)
397
      subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)
398
 
398
 
399
      implicit none
399
      implicit none
400
      
400
      
401
c     Declaration of subroutine parameters
401
c     Declaration of subroutine parameters
402
      integer       fid
402
      integer       fid
403
      integer       time(6)
403
      integer       time(6)
404
      integer       ntra,ntim,ncol
404
      integer       ntra,ntim,ncol
405
      character*80  vars(ncol)
405
      character*80  vars(ncol)
406
      integer       mode
406
      integer       mode
407
 
407
 
408
c     Auxiliary variables
408
c     Auxiliary variables
409
      integer       i
409
      integer       i
410
      character*120 str
410
      character*120 str
411
      character*4   str1
411
      character*4   str1
412
      character*2   str2,str3,str4,str5,str6
412
      character*2   str2,str3,str4,str5,str6
413
      integer       vardim(4)
413
      integer       vardim(4)
414
      real          varmin(4),varmax(4),stag(4)
414
      real          varmin(4),varmax(4),stag(4)
415
      real          mdv
415
      real          mdv
416
      integer       ierr
416
      integer       ierr
417
      character*80  varname
417
      character*80  varname
418
      real          rtime(6)
418
      real          rtime(6)
419
      integer       nvars
419
      integer       nvars
420
 
420
 
421
c     Write ascii format (mode=1,2)
421
c     Write ascii format (mode=1,2)
422
      if ( (mode.eq.1).or.(mode.eq.2) ) then
422
      if ( (mode.eq.1).or.(mode.eq.2) ) then
423
 
423
 
424
c        Get the strings for output
424
c        Get the strings for output
425
         write(str1,'(i4)') time(1)
425
         write(str1,'(i4)') time(1)
426
         write(str2,'(i2)') time(2)
426
         write(str2,'(i2)') time(2)
427
         write(str3,'(i2)') time(3)
427
         write(str3,'(i2)') time(3)
428
         write(str4,'(i2)') time(4)
428
         write(str4,'(i2)') time(4)
429
         write(str5,'(i2)') time(5)
429
         write(str5,'(i2)') time(5)
430
         if (time(2).eq. 0) str2(1:1)='0'
430
         if (time(2).eq. 0) str2(1:1)='0'
431
         if (time(3).eq. 0) str3(1:1)='0'
431
         if (time(3).eq. 0) str3(1:1)='0'
432
         if (time(4).eq. 0) str4(1:1)='0'
432
         if (time(4).eq. 0) str4(1:1)='0'
433
         if (time(5).eq. 0) str5(1:1)='0'
433
         if (time(5).eq. 0) str5(1:1)='0'
434
         if (time(2).lt.10) str2(1:1)='0'
434
         if (time(2).lt.10) str2(1:1)='0'
435
         if (time(3).lt.10) str3(1:1)='0'
435
         if (time(3).lt.10) str3(1:1)='0'
436
         if (time(4).lt.10) str4(1:1)='0'
436
         if (time(4).lt.10) str4(1:1)='0'
437
         if (time(5).lt.10) str5(1:1)='0'
437
         if (time(5).lt.10) str5(1:1)='0'
438
 
438
 
439
c        Write the time specification
439
c        Write the time specification
440
         write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)') 
440
         write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)') 
441
     >          'Reference date ',
441
     >          'Reference date ',
442
     >           str1,str2,str3,'_',str4,str5,
442
     >           str1,str2,str3,'_',str4,str5,
443
     >          ' / Time range',time(6), ' min'
443
     >          ' / Time range',time(6), ' min'
444
         write(fid,*)
444
         write(fid,*)
445
 
445
 
446
c        Write variable names
446
c        Write variable names
447
         str=''
447
         str=''
448
         do i=1,ncol
448
         do i=1,ncol
449
            str=trim(str)//trim(vars(i))
449
            str=trim(str)//trim(vars(i))
450
         enddo
450
         enddo
451
         write(fid,'(a6,a9,a8,a6,100a10)') (trim(vars(i)),i=1,ncol)
451
         write(fid,'(a6,a9,a8,a6,100a10)') (trim(vars(i)),i=1,ncol)
452
         write(fid,'(a6,a9,a8,a6,100a10)') 
452
         write(fid,'(a6,a9,a8,a6,100a10)') 
453
     >              '------','---------','--------','------',
453
     >              '------','---------','--------','------',
454
     >              ('----------',i=5,ncol)
454
     >              ('----------',i=5,ncol)
455
 
455
 
456
c     Write fortran mode (mode=3)
456
c     Write fortran mode (mode=3)
457
      elseif (mode.eq.3) then
457
      elseif (mode.eq.3) then
458
         write(fid) ntra,ntim,ncol
458
         write(fid) ntra,ntim,ncol
459
         write(fid) time
459
         write(fid) time
460
         write(fid) vars
460
         write(fid) vars
461
 
461
 
462
c     Write IVE netcdf format (mode=4)
462
c     Write IVE netcdf format (mode=4)
463
      elseif (mode.eq.4) then
463
      elseif (mode.eq.4) then
464
         vardim(1)=ntra
464
         vardim(1)=ntra
465
         vardim(2)=1
465
         vardim(2)=1
466
         vardim(3)=1
466
         vardim(3)=1
467
         vardim(4)=1
467
         vardim(4)=1
468
         mdv      =-999.98999
468
         mdv      =-999.98999
469
         do i=2,ncol
469
         do i=2,ncol
470
            call putdef(fid,vars(i),4,mdv,vardim,
470
            call putdef(fid,vars(i),4,mdv,vardim,
471
     >                  varmin,varmax,stag,ierr)
471
     >                  varmin,varmax,stag,ierr)
472
         enddo
472
         enddo
473
         varname='BASEDATE'
473
         varname='BASEDATE'
474
         vardim(1)=6
474
         vardim(1)=6
475
         call putdef(fid,varname,4,mdv,vardim,
475
         call putdef(fid,varname,4,mdv,vardim,
476
     >               varmin,varmax,stag,ierr)
476
     >               varmin,varmax,stag,ierr)
477
         do i=1,6
477
         do i=1,6
478
            rtime(i)=real(time(i))
478
            rtime(i)=real(time(i))
479
         enddo
479
         enddo
480
         call putdat(fid,varname,0.,0,rtime,ierr)
480
         call putdat(fid,varname,0.,0,rtime,ierr)
481
 
481
 
482
      endif
482
      endif
483
 
483
 
484
      end
484
      end
485
      
485
      
486
      
486
      
487
c     ----------------------------------------------------------------
487
c     ----------------------------------------------------------------
488
c     Close a trajectory file
488
c     Close a trajectory file
489
c     ----------------------------------------------------------------
489
c     ----------------------------------------------------------------
490
      
490
      
491
      subroutine close_tra(fid,mode)
491
      subroutine close_tra(fid,mode)
492
 
492
 
493
      implicit none
493
      implicit none
494
      
494
      
495
c     Declaration of subroutine parameters
495
c     Declaration of subroutine parameters
496
      integer      fid
496
      integer      fid
497
      integer      mode
497
      integer      mode
498
      
498
      
499
c     Auxiliary variables
499
c     Auxiliary variables
500
      integer      ierr
500
      integer      ierr
501
 
501
 
502
c     Close file
502
c     Close file
503
      if (mode.eq.1) then
503
      if (mode.eq.1) then
504
         close(abs(fid))
504
         close(abs(fid))
505
      elseif (mode.eq.2) then
505
      elseif (mode.eq.2) then
506
         close(abs(fid))
506
         close(abs(fid))
507
      elseif (mode.eq.3) then
507
      elseif (mode.eq.3) then
508
         close(fid)
508
         close(fid)
509
      elseif (mode.eq.4) then
509
      elseif (mode.eq.4) then
510
         call clscdf(fid,ierr)
510
         call clscdf(fid,ierr)
511
      endif
511
      endif
512
 
512
 
513
      end
513
      end
514
      
514
      
515
 
515
 
516
c     ----------------------------------------------------------------
516
c     ----------------------------------------------------------------
517
c     Determine the mode of a trajectory file
517
c     Determine the mode of a trajectory file
518
c     ----------------------------------------------------------------
518
c     ----------------------------------------------------------------
519
 
519
 
520
      subroutine mode_tra(mode,filename)
520
      subroutine mode_tra(mode,filename)
521
      
521
      
522
      implicit none
522
      implicit none
523
 
523
 
524
c     Declaration of subroutine parameters
524
c     Declaration of subroutine parameters
525
      integer        mode
525
      integer        mode
526
      character*80   filename
526
      character*80   filename
527
 
527
 
528
c     Auxiliary variables
528
c     Auxiliary variables
529
      integer        len
529
      integer        len
530
      character      char0,char1
530
      character      char0,char1
531
 
531
 
532
c     Get mode
532
c     Get mode
533
      mode=-1
533
      mode=-1
534
      
534
      
535
      len  = len_trim(filename)
535
      len  = len_trim(filename)
536
 
536
 
537
      char0 = filename((len-1):(len-1))
537
      char0 = filename((len-1):(len-1))
538
      char1 = filename(len:len)
538
      char1 = filename(len:len)
539
 
539
 
540
      if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1
540
      if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1
541
      if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2
541
      if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2
542
      if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3
542
      if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3
543
      if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4
543
      if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4
544
 
544
 
545
      end
545
      end
546
 
546
 
547
 
547
 
548
c     ----------------------------------------------------------------
548
c     ----------------------------------------------------------------
549
c     Get dimension of a trajectory file
549
c     Get dimension of a trajectory file
550
c     ----------------------------------------------------------------
550
c     ----------------------------------------------------------------
551
    
551
    
552
      subroutine info_tra(filename,ntra,ntim,ncol,mode)
552
      subroutine info_tra(filename,ntra,ntim,ncol,mode)
553
 
553
 
554
      implicit none
554
      implicit none
555
      
555
      
556
c     Declaration of subroutine parameters
556
c     Declaration of subroutine parameters
557
      integer      fid
557
      integer      fid
558
      character*80 filename
558
      character*80 filename
559
      integer      mode
559
      integer      mode
560
      integer      ntra,ntim,ncol
560
      integer      ntra,ntim,ncol
561
 
561
 
562
c     Auxiliary variables
562
c     Auxiliary variables
563
      integer       vardim(4)
563
      integer       vardim(4)
564
      real          varmin(4),varmax(4),stag(4)
564
      real          varmin(4),varmax(4),stag(4)
565
      real          mdv
565
      real          mdv
566
      character*80  cfn
566
      character*80  cfn
567
      integer       ierr
567
      integer       ierr
568
      integer       i,ndim
568
      integer       i,ndim
569
      character*80  vars(100)
569
      character*80  vars(100)
570
      integer       nvars
570
      integer       nvars
571
      integer       ntimes
571
      integer       ntimes
572
      real          times(100)
572
      real          times(100)
573
      character*120 str
573
      character*120 str
574
      integer       nline0,nline1,nline2
574
      integer       nline0,nline1,nline2
575
      integer       isstr,isok
575
      integer       isstr,isok
576
      character     ch
576
      character     ch
577
 
577
 
578
c     Open file
578
c     Open file
579
      if (mode.eq.1) then
579
      if (mode.eq.1) then
580
         fid=10
580
         fid=10
581
         open(fid,file=filename)
581
         open(fid,file=filename)
582
      elseif (mode.eq.2) then
582
      elseif (mode.eq.2) then
583
         fid=10
583
         fid=10
584
         open(fid,file=filename)
584
         open(fid,file=filename)
585
      elseif (mode.eq.3) then
585
      elseif (mode.eq.3) then
586
         fid=10
586
         fid=10
587
         open(fid,file=filename,form='unformatted')
587
         open(fid,file=filename,form='unformatted')
588
      elseif (mode.eq.4) then
588
      elseif (mode.eq.4) then
589
         call cdfopn(filename,fid,ierr)
589
         call cdfopn(filename,fid,ierr)
590
      endif
590
      endif
591
 
591
 
592
c     Get dimension information
592
c     Get dimension information
593
      if ( (mode.eq.1).or.(mode.eq.2) ) then
593
      if ( (mode.eq.1).or.(mode.eq.2) ) then
594
         read(fid,*)
594
         read(fid,*)
595
         read(fid,*)
595
         read(fid,*)
596
         read(fid,'(a120)') str
596
         read(fid,'(a120)') str
597
         read(fid,*)
597
         read(fid,*)
598
 
598
 
599
c        Get the number of columns
599
c        Get the number of columns
600
         isstr=0
600
         isstr=0
601
         ncol =0
601
         ncol =0
602
         do i=1,120
602
         do i=1,120
603
            ch = str(i:i)
603
            ch = str(i:i)
604
            if ( (isstr.eq.0).and.(ch.ne.' ') ) then
604
            if ( (isstr.eq.0).and.(ch.ne.' ') ) then
605
               isstr=1
605
               isstr=1
606
            elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
606
            elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
607
               isstr=0
607
               isstr=0
608
               ncol=ncol+1
608
               ncol=ncol+1
609
            endif
609
            endif
610
         enddo
610
         enddo
611
 
611
 
612
c        Get the first data block
612
c        Get the first data block
613
         nline0  = 5
613
         nline0  = 5
614
         nline1  = 5
614
         nline1  = 5
615
         read(fid,*)
615
         read(fid,*)
616
 100     read(fid,'(a120)',end=110) str
616
 100     read(fid,'(a120)',end=110) str
617
         if (str.ne.'') then
617
         if (str.ne.'') then
618
            nline1 = nline1 + 1
618
            nline1 = nline1 + 1
619
            goto 100
619
            goto 100
620
         endif
620
         endif
621
 110     continue
621
 110     continue
622
         
622
         
623
c        Get the total numbers of lines in the data block
623
c        Get the total numbers of lines in the data block
624
         nline2 = nline1
624
         nline2 = nline1
625
 120     read(fid,*,end=130)
625
 120     read(fid,*,end=130)
626
         nline2 = nline2 + 1
626
         nline2 = nline2 + 1
627
         goto 120
627
         goto 120
628
 130     nline2 = nline2 + 1
628
 130     nline2 = nline2 + 1
629
 
629
 
630
c        Set the dimensions
630
c        Set the dimensions
631
         if (mode.eq.1) then
631
         if (mode.eq.1) then
632
            ntim = nline1 - nline0
632
            ntim = nline1 - nline0
633
            ntra = (nline2-nline0+1)/(ntim+1)
633
            ntra = (nline2-nline0+1)/(ntim+1)
634
         else
634
         else
635
            ntra = nline1 - nline0
635
            ntra = nline1 - nline0
636
            ntim = (nline2-nline0+1)/(ntra+1)
636
            ntim = (nline2-nline0+1)/(ntra+1)
637
         endif
637
         endif
638
 
638
 
639
      elseif (mode.eq.3) then
639
      elseif (mode.eq.3) then
640
         read(fid) ntra,ntim,ncol
640
         read(fid) ntra,ntim,ncol
641
 
641
 
642
      elseif (mode.eq.4) then
642
      elseif (mode.eq.4) then
643
         call gettimes(fid,times,ntimes,ierr)
643
         call gettimes(fid,times,ntimes,ierr)
644
         call getvars(fid,nvars,vars,ierr)
644
         call getvars(fid,nvars,vars,ierr)
645
         call getdef(fid,trim(vars(2)),ndim,mdv,vardim,
645
         call getdef(fid,trim(vars(2)),ndim,mdv,vardim,
646
     >               varmin,varmax,stag,ierr)
646
     >               varmin,varmax,stag,ierr)
647
         ntra = vardim(1)
647
         ntra = vardim(1)
648
         ntim = ntimes
648
         ntim = ntimes
649
         ncol = nvars-1
649
         ncol = nvars-1
650
 
650
 
651
      endif
651
      endif
652
 
652
 
653
c     Close file
653
c     Close file
654
      if (mode.eq.1) then
654
      if (mode.eq.1) then
655
         close(fid)
655
         close(fid)
656
      elseif (mode.eq.2) then
656
      elseif (mode.eq.2) then
657
         close(fid)
657
         close(fid)
658
      elseif (mode.eq.3) then
658
      elseif (mode.eq.3) then
659
         close(fid)
659
         close(fid)
660
      elseif (mode.eq.4) then
660
      elseif (mode.eq.4) then
661
         call clscdf(fid,ierr)
661
         call clscdf(fid,ierr)
662
      endif
662
      endif
663
 
663
 
664
 
664
 
665
      end
665
      end
666
      
666