Subversion Repositories lagranto.wrf

Rev

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

Rev 11 Rev 21
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
c	Sebastian: Why 9999 -> conflicts with Pressure if it is in Pa
229
c	Sebastian: Why 9999 -> conflicts with Pressure if it is in Pa
230
               do j=5,ncol
230
               do j=5,ncol
231
c                  if ( abs(tra(n,i,j)).gt.9999.) then   
231
c                  if ( abs(tra(n,i,j)).gt.9999.) then   
232
		   if ( abs(tra(n,i,j)).gt.999999.) then                  
232
		   if ( abs(tra(n,i,j)).gt.999999.) then                  
233
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
233
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
234
                     tra(n,i,j) = -999.99
234
                     tra(n,i,j) = -999.99
235
                  endif
235
                  endif
236
               enddo
236
               enddo
237
 
237
 
238
               write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)') 
238
               write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)') 
239
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
239
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
240
     >               nint(tra(n,i,4)),               ! p
240
     >               nint(tra(n,i,4)),               ! p
241
     >               (tra(n,i,j),j=5,ncol)           ! fields
241
     >               (tra(n,i,j),j=5,ncol)           ! fields
242
            enddo
242
            enddo
243
         enddo
243
         enddo
244
 
244
 
245
c     Write ascii mode, sorted by time (mode=2)
245
c     Write ascii mode, sorted by time (mode=2)
246
      elseif (mode.eq.2) then
246
      elseif (mode.eq.2) then
247
         do i=1,ntim
247
         do i=1,ntim
248
            write(fid,*)
248
            write(fid,*)
249
            do n=1,ntra
249
            do n=1,ntra
250
 
250
 
251
c              Avoid ugly *s or missing space in output
251
c              Avoid ugly *s or missing space in output
252
               do j=5,ncol
252
               do j=5,ncol
253
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
253
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
254
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
254
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
255
                     tra(n,i,j) = -999.99
255
                     tra(n,i,j) = -999.99
256
                  endif
256
                  endif
257
               enddo
257
               enddo
258
 
258
 
259
               write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)') 
259
               write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)') 
260
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
260
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
261
     >               nint(tra(n,i,4)),               ! p
261
     >               nint(tra(n,i,4)),               ! p
262
     >               (tra(n,i,j),j=5,ncol)           ! fields
262
     >               (tra(n,i,j),j=5,ncol)           ! fields
263
            enddo
263
            enddo
264
         enddo
264
         enddo
265
 
265
 
266
c     Write fortran mode (mode=3)
266
c     Write fortran mode (mode=3)
267
      elseif (mode.eq.3) then
267
      elseif (mode.eq.3) then
268
         write(fid) tra                              
268
         write(fid) tra                              
269
 
269
 
270
c     Write netcdf mode (mode=4)
270
c     Write netcdf mode (mode=4)
271
      elseif (mode.eq.4) then
271
      elseif (mode.eq.4) then
272
         call getvars(fid,nvars,vars,ierr)
272
         call getvars(fid,nvars,vars,ierr)
273
         do i=1,ntim
273
         do i=1,ntim
274
            time=tra(1,i,1)
274
            time=tra(1,i,1)
275
            do j=2,ncol
275
            do j=2,ncol
276
               do n=1,ntra
276
               do n=1,ntra
277
                  arr(n)=tra(n,i,j)
277
                  arr(n)=tra(n,i,j)
278
               enddo
278
               enddo
279
               call putdat(fid,vars(j),time,0,arr,ierr)
279
               call putdat(fid,vars(j),time,0,arr,ierr)
280
            enddo
280
            enddo
281
         enddo
281
         enddo
282
      endif
282
      endif
283
 
283
 
284
      end
284
      end
285
 
285
 
286
 
286
 
287
c     ----------------------------------------------------------------
287
c     ----------------------------------------------------------------
288
c     Read header from trajectory file
288
c     Read header from trajectory file
289
c     ----------------------------------------------------------------
289
c     ----------------------------------------------------------------
290
 
290
 
291
      subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)
291
      subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)
292
 
292
 
293
      implicit none
293
      implicit none
294
      
294
      
295
c     Declaration of subroutine parameters
295
c     Declaration of subroutine parameters
296
      integer       fid
296
      integer       fid
297
      integer       time(6)
297
      integer       time(6)
298
      integer       ntra,ntim,ncol
298
      integer       ntra,ntim,ncol
299
      character*80  vars(ncol)
299
      character*80  vars(ncol)
300
      integer       mode
300
      integer       mode
301
 
301
 
302
c     Auxiliary variables
302
c     Auxiliary variables
303
      integer       i
303
      integer       i
304
      character     ch(500)
304
      character     ch(500)
305
      character*500 str
305
      character*500 str
306
      integer       ich(500)
306
      integer       ich(500)
307
      integer       isstr,ileft,iright
307
      integer       isstr,ileft,iright
308
      character*80  varname
308
      character*80  varname
309
      real          rtime(6)
309
      real          rtime(6)
310
      integer       ierr
310
      integer       ierr
311
      integer       nvars
311
      integer       nvars
312
      character*15  str1
312
      character*15  str1
313
      character     str2
313
      character     str2
314
      character*13  str3
314
      character*13  str3
315
      character*4   str4
315
      character*4   str4
316
      character*80  linestr
316
      character*80  linestr
317
      integer       itmp1,itmp2
317
      integer       itmp1,itmp2
318
 
318
 
319
c     Read ascii format (mode=1,2)
319
c     Read ascii format (mode=1,2)
320
      if ( (mode.eq.1).or.(mode.eq.2) ) then
320
      if ( (mode.eq.1).or.(mode.eq.2) ) then
321
 
321
 
322
c        Read the time specification (old and new format)
322
c        Read the time specification (old and new format)
323
         read(fid,'(a80)') linestr
323
         read(fid,'(a80)') linestr
324
 
324
 
325
         if ( linestr(1:14).eq.'Reference date' ) then
325
         if ( linestr(1:14).eq.'Reference date' ) then
326
            read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)') 
326
            read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)') 
327
     >           str1,
327
     >           str1,
328
     >           time(1),time(2),time(3),str2,time(4),time(5),
328
     >           time(1),time(2),time(3),str2,time(4),time(5),
329
     >           str3,time(6),str4
329
     >           str3,time(6),str4
330
                       
330
                       
331
         elseif ( linestr(1:11).eq.'time period' ) then
331
         elseif ( linestr(1:11).eq.'time period' ) then
332
            read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)') 
332
            read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)') 
333
     >           str1,
333
     >           str1,
334
     >           time(1),time(2),time(3),str2,time(4),
334
     >           time(1),time(2),time(3),str2,time(4),
335
     >           str3,itmp1,str3,itmp2,str4
335
     >           str3,itmp1,str3,itmp2,str4
336
            time(5) = 0
336
            time(5) = 0
337
            time(6) = itmp1 * 60 + itmp2
337
            time(6) = itmp1 * 60 + itmp2
338
 
338
 
339
         endif
339
         endif
340
 
340
 
341
c        Skip the empty line and read field names
341
c        Skip the empty line and read field names
342
         read(fid,*)
342
         read(fid,*)
343
         read(fid,'(a500)',end=100) str
343
         read(fid,'(a500)',end=100) str
344
         do i=1,500
344
         do i=1,500
345
            ch(i)=str(i:i)
345
            ch(i)=str(i:i)
346
         enddo
346
         enddo
347
 
347
 
348
c        Split the input string
348
c        Split the input string
349
         isstr=0
349
         isstr=0
350
         nvars=0
350
         nvars=0
351
         do i=1,500
351
         do i=1,500
352
            if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then
352
            if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then
353
               isstr=1
353
               isstr=1
354
               ileft=i
354
               ileft=i
355
            elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then
355
            elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then
356
               isstr=0
356
               isstr=0
357
               iright=i-1
357
               iright=i-1
358
               nvars=nvars+1
358
               nvars=nvars+1
359
               vars(nvars)=str(ileft:iright)
359
               vars(nvars)=str(ileft:iright)
360
            endif
360
            endif
361
         enddo
361
         enddo
362
 
362
 
363
c        Skip the empty line
363
c        Skip the empty line
364
         read(fid,*,end=100)
364
         read(fid,*,end=100)
365
 
365
 
366
c     Read fortran mode (mode=3)
366
c     Read fortran mode (mode=3)
367
      elseif (mode.eq.3) then
367
      elseif (mode.eq.3) then
368
         read(fid) ntra,ntim,ncol
368
         read(fid) ntra,ntim,ncol
369
         read(fid) time
369
         read(fid) time
370
         read(fid) vars
370
         read(fid) vars
371
 
371
 
372
c     Read IVE netcdf mode (mode=4)
372
c     Read IVE netcdf mode (mode=4)
373
      elseif (mode.eq.4) then
373
      elseif (mode.eq.4) then
374
         call getvars(fid,nvars,vars,ierr)  
374
         call getvars(fid,nvars,vars,ierr)  
375
         varname='BASEDATE'
375
         varname='BASEDATE'
376
         call getdat(fid,varname,0.,0,rtime,ierr)
376
         call getdat(fid,varname,0.,0,rtime,ierr)
377
         do i=1,6
377
         do i=1,6
378
            time(i)=nint(rtime(i))
378
            time(i)=nint(rtime(i))
379
         enddo
379
         enddo
380
      endif
380
      endif
381
 
381
 
382
      return
382
      return
383
 
383
 
384
c     End of file has been reached
384
c     End of file has been reached
385
 100  fid=-fid
385
 100  fid=-fid
386
      return
386
      return
387
 
387
 
388
c     Excetion handling
388
c     Excetion handling
389
 110  print*,'<read_hea>: Unexspected time format.... Stop'
389
 110  print*,'<read_hea>: Unexspected time format.... Stop'
390
      stop
390
      stop
391
   
391
   
392
      end
392
      end
393
      
393
      
394
 
394
 
395
c     ----------------------------------------------------------------
395
c     ----------------------------------------------------------------
396
c     Write header to trajectory file (in ascii mode)
396
c     Write header to trajectory file (in ascii mode)
397
c     ----------------------------------------------------------------
397
c     ----------------------------------------------------------------
398
 
398
 
399
      subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)
399
      subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)
400
 
400
 
401
      implicit none
401
      implicit none
402
      
402
      
403
c     Declaration of subroutine parameters
403
c     Declaration of subroutine parameters
404
      integer       fid
404
      integer       fid
405
      integer       time(6)
405
      integer       time(6)
406
      integer       ntra,ntim,ncol
406
      integer       ntra,ntim,ncol
407
      character*80  vars(ncol)
407
      character*80  vars(ncol)
408
      integer       mode
408
      integer       mode
409
 
409
 
410
c     Auxiliary variables
410
c     Auxiliary variables
411
      integer       i
411
      integer       i
412
      character*500 str
412
      character*500 str
413
      character*4   str1
413
      character*4   str1
414
      character*2   str2,str3,str4,str5,str6
414
      character*2   str2,str3,str4,str5,str6
415
      integer       vardim(4)
415
      integer       vardim(4)
416
      real          varmin(4),varmax(4),stag(4)
416
      real          varmin(4),varmax(4),stag(4)
417
      real          mdv
417
      real          mdv
418
      integer       ierr
418
      integer       ierr
419
      character*80  varname
419
      character*80  varname
420
      real          rtime(6)
420
      real          rtime(6)
421
      integer       nvars
421
      integer       nvars
422
 
422
 
423
c     Write ascii format (mode=1,2)
423
c     Write ascii format (mode=1,2)
424
      if ( (mode.eq.1).or.(mode.eq.2) ) then
424
      if ( (mode.eq.1).or.(mode.eq.2) ) then
425
 
425
 
426
c        Get the strings for output
426
c        Get the strings for output
427
         write(str1,'(i4)') time(1)
427
         write(str1,'(i4)') time(1)
428
         write(str2,'(i2)') time(2)
428
         write(str2,'(i2)') time(2)
429
         write(str3,'(i2)') time(3)
429
         write(str3,'(i2)') time(3)
430
         write(str4,'(i2)') time(4)
430
         write(str4,'(i2)') time(4)
431
         write(str5,'(i2)') time(5)
431
         write(str5,'(i2)') time(5)
432
         if (time(2).eq. 0) str2(1:1)='0'
432
         if (time(2).eq. 0) str2(1:1)='0'
433
         if (time(3).eq. 0) str3(1:1)='0'
433
         if (time(3).eq. 0) str3(1:1)='0'
434
         if (time(4).eq. 0) str4(1:1)='0'
434
         if (time(4).eq. 0) str4(1:1)='0'
435
         if (time(5).eq. 0) str5(1:1)='0'
435
         if (time(5).eq. 0) str5(1:1)='0'
436
         if (time(2).lt.10) str2(1:1)='0'
436
         if (time(2).lt.10) str2(1:1)='0'
437
         if (time(3).lt.10) str3(1:1)='0'
437
         if (time(3).lt.10) str3(1:1)='0'
438
         if (time(4).lt.10) str4(1:1)='0'
438
         if (time(4).lt.10) str4(1:1)='0'
439
         if (time(5).lt.10) str5(1:1)='0'
439
         if (time(5).lt.10) str5(1:1)='0'
440
 
440
 
441
c        Write the time specification
441
c        Write the time specification
442
         write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)') 
442
         write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)') 
443
     >          'Reference date ',
443
     >          'Reference date ',
444
     >           str1,str2,str3,'_',str4,str5,
444
     >           str1,str2,str3,'_',str4,str5,
445
     >          ' / Time range',time(6), ' min'
445
     >          ' / Time range',time(6), ' min'
446
         write(fid,*)
446
         write(fid,*)
447
 
447
 
448
c        Write variable names
448
c        Write variable names
449
         str=''
449
         str=''
450
         do i=1,ncol
450
         do i=1,ncol
451
            str=trim(str)//trim(vars(i))
451
            str=trim(str)//trim(vars(i))
452
         enddo
452
         enddo
453
         write(fid,'(a6,a9,a8,a6,100a10)') (trim(vars(i)),i=1,ncol)
453
         write(fid,'(a6,a10,a9,a6,100a10)') (trim(vars(i)),i=1,ncol)
454
         write(fid,'(a6,a9,a8,a6,100a10)') 
454
         write(fid,'(a6,a10,a9,a6,100a10)') 
455
     >              '------','---------','--------','------',
455
     >              '------','----------','---------','------',
456
     >              ('----------',i=5,ncol)
456
     >              ('----------',i=5,ncol)
457
 
457
 
458
c     Write fortran mode (mode=3)
458
c     Write fortran mode (mode=3)
459
      elseif (mode.eq.3) then
459
      elseif (mode.eq.3) then
460
         write(fid) ntra,ntim,ncol
460
         write(fid) ntra,ntim,ncol
461
         write(fid) time
461
         write(fid) time
462
         write(fid) vars
462
         write(fid) vars
463
 
463
 
464
c     Write IVE netcdf format (mode=4)
464
c     Write IVE netcdf format (mode=4)
465
      elseif (mode.eq.4) then
465
      elseif (mode.eq.4) then
466
         vardim(1)=ntra
466
         vardim(1)=ntra
467
         vardim(2)=1
467
         vardim(2)=1
468
         vardim(3)=1
468
         vardim(3)=1
469
         vardim(4)=1
469
         vardim(4)=1
470
         mdv      =-999.98999
470
         mdv      =-999.98999
471
         do i=2,ncol
471
         do i=2,ncol
472
            call putdef(fid,vars(i),4,mdv,vardim,
472
            call putdef(fid,vars(i),4,mdv,vardim,
473
     >                  varmin,varmax,stag,ierr)
473
     >                  varmin,varmax,stag,ierr)
474
         enddo
474
         enddo
475
         varname='BASEDATE'
475
         varname='BASEDATE'
476
         vardim(1)=6
476
         vardim(1)=6
477
         call putdef(fid,varname,4,mdv,vardim,
477
         call putdef(fid,varname,4,mdv,vardim,
478
     >               varmin,varmax,stag,ierr)
478
     >               varmin,varmax,stag,ierr)
479
         do i=1,6
479
         do i=1,6
480
            rtime(i)=real(time(i))
480
            rtime(i)=real(time(i))
481
         enddo
481
         enddo
482
         call putdat(fid,varname,0.,0,rtime,ierr)
482
         call putdat(fid,varname,0.,0,rtime,ierr)
483
 
483
 
484
      endif
484
      endif
485
 
485
 
486
      end
486
      end
487
      
487
      
488
      
488
      
489
c     ----------------------------------------------------------------
489
c     ----------------------------------------------------------------
490
c     Close a trajectory file
490
c     Close a trajectory file
491
c     ----------------------------------------------------------------
491
c     ----------------------------------------------------------------
492
      
492
      
493
      subroutine close_tra(fid,mode)
493
      subroutine close_tra(fid,mode)
494
 
494
 
495
      implicit none
495
      implicit none
496
      
496
      
497
c     Declaration of subroutine parameters
497
c     Declaration of subroutine parameters
498
      integer      fid
498
      integer      fid
499
      integer      mode
499
      integer      mode
500
      
500
      
501
c     Auxiliary variables
501
c     Auxiliary variables
502
      integer      ierr
502
      integer      ierr
503
 
503
 
504
c     Close file
504
c     Close file
505
      if (mode.eq.1) then
505
      if (mode.eq.1) then
506
         close(abs(fid))
506
         close(abs(fid))
507
      elseif (mode.eq.2) then
507
      elseif (mode.eq.2) then
508
         close(abs(fid))
508
         close(abs(fid))
509
      elseif (mode.eq.3) then
509
      elseif (mode.eq.3) then
510
         close(fid)
510
         close(fid)
511
      elseif (mode.eq.4) then
511
      elseif (mode.eq.4) then
512
         call clscdf(fid,ierr)
512
         call clscdf(fid,ierr)
513
      endif
513
      endif
514
 
514
 
515
      end
515
      end
516
      
516
      
517
 
517
 
518
c     ----------------------------------------------------------------
518
c     ----------------------------------------------------------------
519
c     Determine the mode of a trajectory file
519
c     Determine the mode of a trajectory file
520
c     ----------------------------------------------------------------
520
c     ----------------------------------------------------------------
521
 
521
 
522
      subroutine mode_tra(mode,filename)
522
      subroutine mode_tra(mode,filename)
523
      
523
      
524
      implicit none
524
      implicit none
525
 
525
 
526
c     Declaration of subroutine parameters
526
c     Declaration of subroutine parameters
527
      integer        mode
527
      integer        mode
528
      character*80   filename
528
      character*80   filename
529
 
529
 
530
c     Auxiliary variables
530
c     Auxiliary variables
531
      integer        len
531
      integer        len
532
      character      char0,char1
532
      character      char0,char1
533
 
533
 
534
c     Get mode
534
c     Get mode
535
      mode=-1
535
      mode=-1
536
      
536
      
537
      len  = len_trim(filename)
537
      len  = len_trim(filename)
538
 
538
 
539
      char0 = filename((len-1):(len-1))
539
      char0 = filename((len-1):(len-1))
540
      char1 = filename(len:len)
540
      char1 = filename(len:len)
541
 
541
 
542
      if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1
542
      if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1
543
      if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2
543
      if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2
544
      if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3
544
      if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3
545
      if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4
545
      if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4
546
 
546
 
547
      end
547
      end
548
 
548
 
549
 
549
 
550
c     ----------------------------------------------------------------
550
c     ----------------------------------------------------------------
551
c     Get dimension of a trajectory file
551
c     Get dimension of a trajectory file
552
c     ----------------------------------------------------------------
552
c     ----------------------------------------------------------------
553
    
553
    
554
      subroutine info_tra(filename,ntra,ntim,ncol,mode)
554
      subroutine info_tra(filename,ntra,ntim,ncol,mode)
555
 
555
 
556
      implicit none
556
      implicit none
557
      
557
      
558
c     Declaration of subroutine parameters
558
c     Declaration of subroutine parameters
559
      integer      fid
559
      integer      fid
560
      character*80 filename
560
      character*80 filename
561
      integer      mode
561
      integer      mode
562
      integer      ntra,ntim,ncol
562
      integer      ntra,ntim,ncol
563
 
563
 
564
c     Auxiliary variables
564
c     Auxiliary variables
565
      integer       vardim(4)
565
      integer       vardim(4)
566
      real          varmin(4),varmax(4),stag(4)
566
      real          varmin(4),varmax(4),stag(4)
567
      real          mdv
567
      real          mdv
568
      character*80  cfn
568
      character*80  cfn
569
      integer       ierr
569
      integer       ierr
570
      integer       i,ndim
570
      integer       i,ndim
571
      character*80  vars(100)
571
      character*80  vars(100)
572
      integer       nvars
572
      integer       nvars
573
      integer       ntimes
573
      integer       ntimes
574
      real          times(100)
574
      real          times(100)
575
      character*500 str
575
      character*500 str
576
      integer       nline0,nline1,nline2
576
      integer       nline0,nline1,nline2
577
      integer       isstr,isok
577
      integer       isstr,isok
578
      character     ch
578
      character     ch
579
 
579
 
580
c     Open file
580
c     Open file
581
      if (mode.eq.1) then
581
      if (mode.eq.1) then
582
         fid=10
582
         fid=10
583
         open(fid,file=filename)
583
         open(fid,file=filename)
584
      elseif (mode.eq.2) then
584
      elseif (mode.eq.2) then
585
         fid=10
585
         fid=10
586
         open(fid,file=filename)
586
         open(fid,file=filename)
587
      elseif (mode.eq.3) then
587
      elseif (mode.eq.3) then
588
         fid=10
588
         fid=10
589
         open(fid,file=filename,form='unformatted')
589
         open(fid,file=filename,form='unformatted')
590
      elseif (mode.eq.4) then
590
      elseif (mode.eq.4) then
591
         call cdfopn(filename,fid,ierr)
591
         call cdfopn(filename,fid,ierr)
592
      endif
592
      endif
593
 
593
 
594
c     Get dimension information
594
c     Get dimension information
595
      if ( (mode.eq.1).or.(mode.eq.2) ) then
595
      if ( (mode.eq.1).or.(mode.eq.2) ) then
596
         read(fid,*)
596
         read(fid,*)
597
         read(fid,*)
597
         read(fid,*)
598
         read(fid,'(a500)') str
598
         read(fid,'(a500)') str
599
         read(fid,*)
599
         read(fid,*)
600
 
600
 
601
c        Get the number of columns
601
c        Get the number of columns
602
         isstr=0
602
         isstr=0
603
         ncol =0
603
         ncol =0
604
         do i=1,500
604
         do i=1,500
605
            ch = str(i:i)
605
            ch = str(i:i)
606
 
606
 
607
            if ( (isstr.eq.0).and.(ch.ne.' ') ) then
607
            if ( (isstr.eq.0).and.(ch.ne.' ') ) then
608
               isstr=1
608
               isstr=1
609
            elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
609
            elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
610
               isstr=0
610
               isstr=0
611
               ncol=ncol+1
611
               ncol=ncol+1
612
            endif
612
            endif
613
         enddo
613
         enddo
614
c           backspace(fid)
614
c           backspace(fid)
615
c        Get the first data block
615
c        Get the first data block
616
         nline0  = 5
616
         nline0  = 5
617
         nline1  = 5
617
         nline1  = 5
618
         read(fid,*)
618
         read(fid,*)
619
 100     read(fid,'(a500)',end=110) str
619
 100     read(fid,'(a500)',end=110) str
620
         if (str.ne.'') then
620
         if (str.ne.'') then
621
            nline1 = nline1 + 1
621
            nline1 = nline1 + 1
622
            goto 100
622
            goto 100
623
         endif
623
         endif
624
 110     continue
624
 110     continue
625
           backspace(fid)
625
           backspace(fid)
626
         
626
         
627
c        Get the total numbers of lines in the data block
627
c        Get the total numbers of lines in the data block
628
         nline2 = nline1
628
         nline2 = nline1
629
 120     read(fid,*,end=130)
629
 120     read(fid,*,end=130)
630
         nline2 = nline2 + 1
630
         nline2 = nline2 + 1
631
         goto 120
631
         goto 120
632
 130     nline2 = nline2 + 1
632
 130     nline2 = nline2 + 1
633
 
633
 
634
c        Set the dimensions
634
c        Set the dimensions
635
         if (mode.eq.1) then
635
         if (mode.eq.1) then
636
            ntim = nline1 - nline0
636
            ntim = nline1 - nline0
637
            ntra = (nline2-nline0+1)/(ntim+1)
637
            ntra = (nline2-nline0+1)/(ntim+1)
638
         else
638
         else
639
            ntra = nline1 - nline0
639
            ntra = nline1 - nline0
640
            ntim = (nline2-nline0+1)/(ntra+1)
640
            ntim = (nline2-nline0+1)/(ntra+1)
641
         endif
641
         endif
642
 
642
 
643
      elseif (mode.eq.3) then
643
      elseif (mode.eq.3) then
644
         read(fid) ntra,ntim,ncol
644
         read(fid) ntra,ntim,ncol
645
 
645
 
646
      elseif (mode.eq.4) then
646
      elseif (mode.eq.4) then
647
         call gettimes(fid,times,ntimes,ierr)
647
         call gettimes(fid,times,ntimes,ierr)
648
         call getvars(fid,nvars,vars,ierr)
648
         call getvars(fid,nvars,vars,ierr)
649
         call getdef(fid,trim(vars(2)),ndim,mdv,vardim,
649
         call getdef(fid,trim(vars(2)),ndim,mdv,vardim,
650
     >               varmin,varmax,stag,ierr)
650
     >               varmin,varmax,stag,ierr)
651
         ntra = vardim(1)
651
         ntra = vardim(1)
652
         ntim = ntimes
652
         ntim = ntimes
653
         ncol = nvars-1
653
         ncol = nvars-1
654
 
654
 
655
      endif
655
      endif
656
 
656
 
657
c     Close file
657
c     Close file
658
      if (mode.eq.1) then
658
      if (mode.eq.1) then
659
         close(fid)
659
         close(fid)
660
      elseif (mode.eq.2) then
660
      elseif (mode.eq.2) then
661
         close(fid)
661
         close(fid)
662
      elseif (mode.eq.3) then
662
      elseif (mode.eq.3) then
663
         close(fid)
663
         close(fid)
664
      elseif (mode.eq.4) then
664
      elseif (mode.eq.4) then
665
         call clscdf(fid,ierr)
665
         call clscdf(fid,ierr)
666
      endif
666
      endif
667
 
667
 
668
 
668
 
669
      end
669
      end
670
      
670