Subversion Repositories lagranto.icon

Rev

Rev 3 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
c     ****************************************************************
2
c     * This package provides IO routines for trajectories. A file   *
3
c     * is characterised by the filename <filename> and the file     *
4
c     * identifier <fid>. Different modes <mode> are supported:      *
5
c     *     mode=1: ascii, sorted by trajectory;                     *
6
c     *     mode=2: ascii, sorted by time;                           *
7
c     *     mode=3: fortran (unformatted)                            * 
8
c     *     mode=4: IVE netcdf (for compatibility reasons)           *
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       *
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,    *
13
c     * longitude, latitude, pressure. The other columns are traced  *
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  *
16
c     * the array <time(6)=year,month,day,hour,time length of the    *
17
c     * trajectory (hour,min)>.                                      *
18
c     *                                                              *
19
c     * Author: Michael Sprenger, September 2008                     *
20
c     ****************************************************************
21
 
22
c     ----------------------------------------------------------------
23
c     Open a trajectory file for reading
24
c     ----------------------------------------------------------------
25
 
26
      subroutine ropen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)
27
 
28
      use netcdf
29
 
30
      implicit none
31
 
32
c     Declaration of subroutine parameters
33
      integer      fid
34
      character*80 filename
35
      integer      mode
36
      integer      ntra,ntim,ncol
37
      integer      time(6)
38
      character*80 vars(ncol)
39
 
40
c     Auxiliary variables
41
      integer      vardim(4)
42
      real         varmin(4),varmax(4),stag(4)
43
      real         mdv
44
      character*80 cfn
45
      integer      ierr
46
      integer      i
47
      integer      nvars
48
 
49
c     Open file
50
      if (mode.eq.1) then
51
         fid = 10
52
         open(fid,file=filename)
53
 
54
      elseif (mode.eq.2) then
55
         fid = 10
56
         open(fid,file=filename)
57
 
58
      elseif (mode.eq.3) then
59
         open(fid,file=filename,form='unformatted')
60
 
61
      elseif (mode.eq.4) then
62
 
63
         ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
64
         IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
65
 
66
      elseif (mode.eq.5) then
67
         ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
68
         IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
69
 
70
      elseif (mode.eq.6) then
71
         print*,' ERROR: Reading KML not supported'
72
         stop
73
      endif
74
 
75
c     Read header information
76
      call read_hea(fid,time,vars,ntra,ntim,ncol,mode)
77
 
78
      end
79
 
80
 
81
c     ----------------------------------------------------------------
82
c     Open a trajectory file for wrinting
83
c     ----------------------------------------------------------------
84
 
85
      subroutine wopen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)
86
 
87
      use netcdf
88
 
89
      implicit none
90
 
91
c     Declaration of subroutine parameters
92
      integer      fid
93
      character*80 filename
94
      integer      mode
95
      integer      ntra,ntim,ncol
96
      integer      time(6) 
97
      character*80 vars(ncol)
98
 
99
c     Auxiliary variables
100
      integer      vardim(4)
101
      real         varmin(4),varmax(4),stag(4)
102
      real         mdv
103
      character*80 cfn
104
      integer      ierr
105
      integer      i
106
      character*80 varname
107
      real         rtime(6)
108
 
109
c     Open file
110
      if (mode.eq.1) then
111
         fid = 10
112
         open(fid,file=filename)
113
 
114
      elseif (mode.eq.2) then
115
         fid = 10
116
         open(fid,file=filename)
117
 
118
      elseif (mode.eq.3) then
119
         open(fid,file=filename,form='unformatted')
120
 
121
      elseif (mode.eq.4) then
122
         ierr = nf90_create(path  = filename, 
123
     >                      cmode = nf90_clobber + nf90_64bit_offset,
124
     >                      ncid  = fid)
125
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
126
 
127
      elseif (mode.eq.5) then
128
         print*,' ERROR: writing online format not supported'
129
         stop
130
 
131
      elseif (mode.eq.6) then
132
         fid = 10
133
         open(fid,file=filename)
134
 
135
      endif
136
 
137
c     Write header
138
      call write_hea(fid,time,vars,ntra,ntim,ncol,mode)
139
 
140
      end
141
 
142
 
143
c     ----------------------------------------------------------------
144
c     Read a trajectory
145
c     ----------------------------------------------------------------
146
 
147
      subroutine read_tra(fid,tra,ntra,ntim,ncol,mode)
148
 
149
      use netcdf
150
 
151
      implicit none
152
 
153
c     Declaration of subroutine parameters
154
      integer   fid
155
      integer   ntim
156
      integer   ncol
157
      integer   ntra
158
      real      tra(ntra,ntim,ncol)
159
      integer   mode
160
 
161
c     Remember format of trajectory file
162
      integer oldtrajform
163
      common oldtrajform
164
 
165
c     Auxiliary variables
166
      integer       i,j,n,d
167
      real          arr(ntra)
168
      integer       ntimes
169
      real          times(ntim)
170
      integer       ierr
171
      character*80  vars(ncol+3)
172
      integer       nvars
173
      real          pollon
174
      real          pollat
175
      integer       tsid
176
      integer       polid
177
      integer       varid
178
      integer       hh,mm
179
      real          frac
180
 
181
c     Read ascii mode, sorted by trajectory (mode=1)
182
      if (mode.eq.1) then
183
         read(fid,*,end=100)
184
         do n=1,ntra
185
            do i=1,ntim
186
               read(fid,*,end=110) (tra(n,i,j),j=1,ncol)
187
            enddo
188
         enddo
189
 
190
c        Adjust time from fractional time to hh.mm time; this is a little
191
c        ugly because the old and new Lagranto assume different time formats
192
         if ( oldtrajform.eq.1 ) then
193
            do n=1,ntra
194
              do i=1,ntim
195
                 frac       = tra(n,i,1)
196
                 hh         = int(frac)
197
                 mm         = nint(60. * (frac-real(int(frac))))
198
                 tra(n,i,1) = real(hh) + 0.01 * real(mm)
199
              enddo
200
            enddo
201
         endif
202
 
203
c     Read ascii mode, sorted by time (mode=2)
204
      elseif (mode.eq.2) then
205
         read(fid,*,end=100)
206
         do i=1,ntim
207
            do n=1,ntra
208
               read(fid,*,end=100) (tra(n,i,j),j=1,ncol)
209
            enddo
210
         enddo
211
 
212
c     Read fortran mode (mode=3)
213
      elseif (mode.eq.3) then
214
         read(fid) tra
215
 
216
c     Read IVE netcdf mode (mode=4)
217
      elseif (mode.eq.4) then
218
 
5 michaesp 219
         call input_getvars(fid,vars,nvars,ierr)
3 michaesp 220
 
221
         ierr = NF90_INQ_VARID(fid,'time',varid)
222
         IF (ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
223
 
224
         ierr = nf90_get_var(fid, varid, tra(1,:,1)  )
225
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
226
         do i=1,ntra
227
            do j=1,ntim
228
               tra(i,:,1) = tra(1,:,1)
229
            enddo
230
         enddo
231
 
232
         do i=2,ncol
233
 
234
            ierr = NF90_INQ_VARID(fid,vars(i),varid)
235
            IF (ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
236
 
237
            ierr = nf90_get_var(fid, varid, tra(:,:,i)  )
238
            IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
239
 
240
         enddo
241
 
242
c     Read COSMO online netcdf mode (mode=5)
243
      elseif (mode.eq.5) then
244
 
5 michaesp 245
         call input_getvars(fid,vars,nvars,ierr)
3 michaesp 246
 
247
         ierr = NF90_INQ_VARID(fid,'time',varid)
248
         IF (ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
249
 
250
         ierr = nf90_get_var(fid, varid, tra(1,:,1)  )
251
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
252
 
253
         do i=2,ncol
254
 
255
            ierr = NF90_INQ_VARID(fid,vars(i),varid)
256
            IF (ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
257
 
258
            ierr = nf90_get_var(fid, varid, tra(:,:,i)  )
259
            IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
260
 
261
         enddo
262
 
263
         do i=1,ntim
264
            hh = int( tra(1,i,1) / 3600. )
265
            mm = int( tra(1,i,1) /   60. ) - 60 * hh
266
            print*,tra(1,i,1), hh,mm
267
            do j=1,ntra
268
               tra(j,i,1) = real(hh) + 0.01 * real(mm)
269
            enddo
270
         enddo
271
 
272
      endif
273
 
274
      return
275
 
276
c     End of file has been reached: set negative <fid>
277
 100  fid=-fid
278
      return
279
 
280
c     Error: incomplete trajectory
281
 110  print*,'<read_tra>: Incomplete trajectory... Stop'
282
      stop
283
 
284
      end
285
 
286
 
287
c     ----------------------------------------------------------------
288
c     Write a trajectory
289
c     ----------------------------------------------------------------
290
 
291
      subroutine write_tra(fid,tra,ntra,ntim,ncol,mode)
292
 
293
      use netcdf
294
 
295
      implicit none
296
 
297
c     Declaration of subroutine parameters
298
      integer   fid
299
      integer   ntim
300
      integer   ncol
301
      integer   ntra
302
      real      tra(ntra,ntim,ncol)
303
      integer   mode
304
 
305
c     Auxiliary variables
306
      integer      i,j,n
307
      real         arr(ntra)
308
      integer      ierr
309
      real         time
310
      character*80 vars(ncol+4)
311
      integer      nvars
312
      integer      varid
313
      character*80 outstr,lonstr,levstr,latstr
314
 
315
c     Write ascii mode, sorted by trajectory (mode=1)
316
      if (mode.eq.1) then
317
         do n=1,ntra
318
            write(fid,*)
319
            do i=1,ntim
320
 
321
c              Avoid ugly *s or missing space in output
322
               do j=5,ncol
323
                  if ( abs(tra(n,i,j)).gt.99999.) then
324
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
325
                     tra(n,i,j) = -999.99
326
                  endif
327
               enddo
328
 
329
               write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)') 
330
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
331
     >               nint(tra(n,i,4)),               ! z
332
     >               (tra(n,i,j),j=5,ncol)           ! fields
333
            enddo
334
         enddo
335
 
336
c     Write ascii mode, sorted by time (mode=2)
337
      elseif (mode.eq.2) then
338
         do i=1,ntim
339
            write(fid,*)
340
            do n=1,ntra
341
 
342
c              Avoid ugly *s or missing space in output
343
               do j=5,ncol
344
                  if ( abs(tra(n,i,j)).gt.99999.) then
345
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
346
                     tra(n,i,j) = -999.99
347
                  endif
348
               enddo
349
 
350
               write(fid,'(1f7.2,f10.3,f9.3,i6,100f10.3)') 
351
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
352
     >               nint(tra(n,i,4)),               ! z
353
     >               (tra(n,i,j),j=5,ncol)           ! fields
354
            enddo
355
         enddo
356
 
357
c     Write fortran mode (mode=3)
358
      elseif (mode.eq.3) then
359
         write(fid) tra                              
360
 
361
c     Write netcdf mode (mode=4)
362
      elseif (mode.eq.4) then
5 michaesp 363
 
364
         call input_getvars(fid,vars,nvars,ierr)
3 michaesp 365
 
366
         ierr = NF90_INQ_VARID(fid,'time',varid)
367
         IF (ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
368
 
369
         ierr = nf90_put_var(fid, varid, tra(1,:,1)  )
370
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
371
 
372
         do i=2,ncol
373
 
374
            ierr = NF90_INQ_VARID(fid,vars(i),varid)
375
            IF (ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
376
 
377
            ierr = nf90_put_var(fid, varid, tra(:,:,i)  )
378
            IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
379
 
380
         enddo
381
 
382
c     Write COSMO online netcdf mode (mode=5)
383
      elseif (mode.eq.5) then
384
 
385
         print*,' ERROR: writing online format not supported...'
386
         stop
387
 
388
c     Write KML mode (mode=5)
389
      elseif (mode.eq.6) then
390
 
391
         do n=1,ntra
392
           write(fid,"(A)") '<Placemark>'
393
           write(fid,"(A)") '<name>Absolute Extruded</name>'
394
           write(fid,"(A)") '<styleUrl>#yellowkLineGreenPoly</styleUrl>'
395
           write(fid,"(A)") '<LineString>'
396
           write(fid,"(A)") '<extrude>1</extrude>'
397
           write(fid,"(A)") '<tessellate>1</tessellate>'
398
           write(fid,"(A)") '<altitudeMode>absolute</altitudeMode>'
399
           write(fid,"(A)") '<coordinates>'
400
 
401
           do i=1,ntim
402
             write(lonstr,*) tra(n,i,2)
403
             write(latstr,*) tra(n,i,3)
404
             write(levstr,*) tra(n,i,4)
405
 
406
             outstr = trim(adjustl(lonstr))//','//
407
     >                trim(adjustl(latstr))//','//
408
     >                trim(adjustl(levstr))
409
 
410
             write(fid,"(A)") outstr
411
 
412
           enddo
413
 
414
           write(fid,*) '</coordinates>'
415
           write(fid,*) '</LineString>'
416
           write(fid,*) '</Placemark>'
417
         enddo
418
 
419
 
420
 
421
      endif
422
 
423
      end
424
 
425
 
426
c     ----------------------------------------------------------------
427
c     Read header from trajectory file
428
c     ----------------------------------------------------------------
429
 
430
      subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)
431
 
432
      use netcdf
433
 
434
      implicit none
435
 
436
c     Declaration of subroutine parameters
437
      integer       fid
438
      integer       time(6)
439
      integer       ntra,ntim,ncol
440
      character*80  vars(ncol+3)
441
      character*80  tmp(ncol)
442
      integer       mode
443
 
444
c     Remember format of trajectory file
445
      integer oldtrajform
446
      common oldtrajform
447
 
448
c     Auxiliary variables
449
      integer       i
450
      character     ch(200)
451
      character*200 str
452
      integer       ich(200)
453
      integer       isstr,ileft,iright
454
      character*80  varname
455
      real          rtime(6)
456
      integer       ierr
457
      integer       nvars
458
      character*15  str1
459
      character     str2
460
      character*13  str3
461
      character*4   str4
462
      character*80  linestr
463
      integer       itmp1,itmp2
464
      character*80  vars_on_file(100)
465
 
466
c     Read ascii format (mode=1,2)
467
      if ( (mode.eq.1).or.(mode.eq.2) ) then
468
 
469
c        Read the time specification (old and new format)
470
         read(fid,'(a80)') linestr
471
 
472
         if ( linestr(1:14).eq.'Reference date' ) then
473
            read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)') 
474
     >           str1,
475
     >           time(1),time(2),time(3),str2,time(4),time(5),
476
     >           str3,time(6),str4
477
            oldtrajform = 0
478
 
479
         elseif ( linestr(1:11).eq.'time period' ) then
480
            read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)') 
481
     >           str1,
482
     >           time(1),time(2),time(3),str2,time(4),
483
     >           str3,itmp1,str3,itmp2,str4
484
            time(5) = 0
485
            time(6) = itmp1 * 60 + itmp2
486
            oldtrajform = 1
487
         endif
488
 
489
c        Skip the empty line and read field names
490
         read(fid,*)
491
         read(fid,'(a200)',end=100) str
492
         do i=1,200
493
            ch(i)=str(i:i)
494
         enddo
495
 
496
c        Split the input string
497
         isstr=0
498
         nvars=0
499
         do i=1,200
500
            if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then
501
               isstr=1
502
               ileft=i
503
            elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then
504
               isstr=0
505
               iright=i-1
506
               nvars=nvars+1
507
               vars(nvars)=str(ileft:iright)
508
            endif
509
         enddo
510
 
511
c        Skip the empty line
512
         read(fid,*,end=100)
513
 
514
c     Read fortran mode (mode=3)
515
      elseif (mode.eq.3) then
516
         read(fid) ntra,ntim,ncol
517
         read(fid) time
518
         read(fid) vars
519
 
520
c     Read IVE netcdf mode (mode=4)
521
      elseif (mode.eq.4) then
522
 
523
         ierr  = nf90_get_att(fid,nf90_global,"ref_year",  time(1) )
524
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
525
         ierr  = nf90_get_att(fid,nf90_global,"ref_month", time(2) )
526
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
527
         ierr  = nf90_get_att(fid,nf90_global,"ref_day",   time(3) )
528
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
529
         ierr  = nf90_get_att(fid,nf90_global,"ref_hour",  time(4) )
530
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
531
         ierr  = nf90_get_att(fid,nf90_global,"ref_min",   time(5) )
532
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
533
         ierr  = nf90_get_att(fid,nf90_global,"duration",  time(6) )
534
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
535
 
5 michaesp 536
         call  input_getvars(fid,vars,nvars,ierr)
3 michaesp 537
 
538
c     Read COSMO online netcdf mode (mode=5)
539
      elseif (mode.eq.5) then
540
 
541
         ierr  = nf90_get_att(fid,nf90_global,"ref_year",  time(1) )
542
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
543
         ierr  = nf90_get_att(fid,nf90_global,"ref_month", time(2) )
544
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
545
         ierr  = nf90_get_att(fid,nf90_global,"ref_day",   time(3) )
546
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
547
         ierr  = nf90_get_att(fid,nf90_global,"ref_hour",  time(4) )
548
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
549
         ierr  = nf90_get_att(fid,nf90_global,"ref_min",   time(5) )
550
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
551
         ierr  = nf90_get_att(fid,nf90_global,"duration",  time(6) )
552
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
553
 
5 michaesp 554
         call  input_getvars(fid,vars,nvars,ierr)
3 michaesp 555
 
556
         if ( vars(2).eq.'longitude' ) then
557
            vars(2) = 'lon'
558
         else
559
            print*,' ERROR: lon missing on netCDF file'
560
            stop
561
         endif
562
 
563
         if ( vars(3).eq.'latitude' ) then
564
            vars(3) = 'lat'
565
         else
566
            print*,' ERROR: lat missing on netCDF file'
567
            stop
568
         endif
569
 
570
         if ( vars(4).ne.'z' ) then
571
            print*,' ERROR: lat missing on netCDF file'
572
            stop
573
         endif
574
 
575
      endif
576
 
577
      return
578
 
579
c     End of file has been reached
580
 100  fid=-fid
581
      return
582
 
583
c     Excetion handling
584
 110  print*,'<read_hea>: Unexspected time format.... Stop'
585
      stop
586
 
587
      end
588
 
589
 
590
c     ----------------------------------------------------------------
591
c     Write header to trajectory file (in ascii mode)
592
c     ----------------------------------------------------------------
593
 
594
      subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)
595
 
596
      use netcdf
597
 
598
      implicit none
599
 
600
c     Declaration of subroutine parameters
601
      integer       fid
602
      integer       time(6)
603
      integer       ntra,ntim,ncol
604
      character*80  vars(ncol)
605
      integer       mode
606
 
607
c     Auxiliary variables
608
      integer       i
609
      character*200 str
610
      character*4   str1
611
      character*2   str2,str3,str4,str5,str6
612
      integer       vardim(4)
613
      real          varmin(4),varmax(4),stag(4)
614
      real          mdv
615
      integer       ierr
616
      character*80  varname
617
      real          rtime(6)
618
      integer       nvars
619
      character*80  str80
620
      integer       ntraDimID,ntimDimID,ncolDimID
621
      integer       dimids(2)
622
      integer       varid
623
      character*80  vars_on_file(100)
624
 
625
c     Write ascii format (mode=1,2)
626
      if ( (mode.eq.1).or.(mode.eq.2) ) then
627
 
628
c        Get the strings for output
629
         write(str1,'(i4)') time(1)
630
         write(str2,'(i2)') time(2)
631
         write(str3,'(i2)') time(3)
632
         write(str4,'(i2)') time(4)
633
         write(str5,'(i2)') time(5)
634
         if (time(2).eq. 0) str2(1:1)='0'
635
         if (time(3).eq. 0) str3(1:1)='0'
636
         if (time(4).eq. 0) str4(1:1)='0'
637
         if (time(5).eq. 0) str5(1:1)='0'
638
         if (time(2).lt.10) str2(1:1)='0'
639
         if (time(3).lt.10) str3(1:1)='0'
640
         if (time(4).lt.10) str4(1:1)='0'
641
         if (time(5).lt.10) str5(1:1)='0'
642
 
643
c        Write the time specification
644
         write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)') 
645
     >          'Reference date ',
646
     >           str1,str2,str3,'_',str4,str5,
647
     >          ' / Time range',time(6), ' min'
648
         write(fid,*)
649
 
650
c        Write variable names
651
         str=''
652
         do i=1,ncol
653
            if ( len_trim(vars(i)).ge.10 ) then
654
               print*,' WARNING: Field name too long... taking 9 chars'
655
               str80   = vars(i)
656
               vars(i) = trim( str80(1:9) )
657
               print*,'             ',trim(str80),' -> ',trim(vars(i)) 
658
            endif
659
            str=trim(str)//trim(vars(i))
660
         enddo
661
         write(fid,'(a7,a10,a9,a6,100a10)') (trim(vars(i)),i=1,ncol)
662
         write(fid,'(a7,a10,a9,a6,100a10)')
663
     >              '-------','----------','---------','------',
664
     >              ('----------',i=5,ncol)
665
 
666
c     Write fortran mode (mode=3)
667
      elseif (mode.eq.3) then
668
         write(fid) ntra,ntim,ncol
669
         write(fid) time
670
         write(fid) vars
671
 
672
c     Write IVE netcdf format (mode=4)
673
      elseif (mode.eq.4) then
5 michaesp 674
 
3 michaesp 675
         ierr  = nf90_put_att(fid,nf90_global,"ref_year",  time(1) )
676
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
677
         ierr  = nf90_put_att(fid,nf90_global,"ref_month", time(2) )
678
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
679
         ierr  = nf90_put_att(fid,nf90_global,"ref_day",   time(3) )
680
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
681
         ierr  = nf90_put_att(fid,nf90_global,"ref_hour",  time(4) )
682
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
683
         ierr  = nf90_put_att(fid,nf90_global,"ref_min",   time(5) )
684
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
685
         ierr  = nf90_put_att(fid,nf90_global,"duration",  time(6) )
686
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
687
         ierr  = nf90_put_att(fid,nf90_global,"pollon",  0.        )
688
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr) 
689
         ierr  = nf90_put_att(fid,nf90_global,"pollat",  90.       )
690
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
691
 
692
         ierr=nf90_def_dim(fid,'ntra',ntra, dimids(1) )
693
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
694
         ierr=nf90_def_dim(fid,'ntim',ntim, dimids(2) )
695
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
5 michaesp 696
 
3 michaesp 697
         ierr = nf90_def_var(fid,'time',
698
     >                       NF90_FLOAT,dimids(2),varid)
699
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
700
         do i=2,ncol
701
            ierr = nf90_def_var(fid,vars(i),
702
     >                           NF90_FLOAT,dimids,varid)
703
            IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
704
         enddo
705
 
706
         ierr = nf90_enddef(fid)
707
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
5 michaesp 708
         print*,1
709
         call input_getvars(fid,vars_on_file,nvars,ierr)
710
         print*,2,nvars
3 michaesp 711
         if ( nvars.ne.ncol ) then
712
            print*,' ERROR: wrong number of fields on netCDF....'
713
            stop
714
         endif
715
         do i=1,nvars
716
            if ( vars(i).ne.vars_on_file(i) ) then
717
               print*,' ERROR: Field order wrong on netCDF...'
718
               stop
719
            endif
720
         enddo
5 michaesp 721
 
3 michaesp 722
c     Write COSMO online netcdf format (mode=5)
723
      elseif (mode.eq.5) then
724
         print*," ERROR: writing online format not supported"
725
         stop
726
 
727
c     Write KML format (mode=5)
728
      elseif (mode.eq.6) then
729
 
730
      write(fid,"(A)") '<?xml version="1.0" encoding="UTF-8"?>'
731
      write(fid,"(A)") '<kml xmlns="http://www.opengis.net/kml/2.2">'
732
      write(fid,"(A)") '<Document>'
733
      write(fid,"(A)") '<name>Paths</name>'
734
      write(fid,"(A)") '<Style id="yellowLineGreenPoly">'
735
      write(fid,"(A)") '<LineStyle>'
736
c      write(fid,*) '<color>7f00ffff</color>'    ! Yellow
737
      write(fid,"(A)") '<color>500A0A0A</color>'     ! Black
738
      write(fid,"(A)") '<width>4</width>'
739
      write(fid,"(A)") '</LineStyle>'
740
      write(fid,"(A)") '<PolyStyle>'
741
      write(fid,"(A)") '<color>7f00ff00</color>'
742
      write(fid,"(A)") '</PolyStyle>'
743
      write(fid,"(A)") '</Style>'
744
 
745
      endif
746
 
747
      end
748
 
749
 
750
c     ----------------------------------------------------------------
751
c     Close a trajectory file
752
c     ----------------------------------------------------------------
753
 
754
      subroutine close_tra(fid,mode)
755
 
756
      use netcdf
757
      implicit none
758
 
759
c     Declaration of subroutine parameters
760
      integer      fid
761
      integer      mode
762
 
763
c     Auxiliary variables
764
      integer      ierr
765
 
766
c     Close file
767
      if (mode.eq.1) then
768
         close(abs(fid))
769
      elseif (mode.eq.2) then
770
         close(abs(fid))
771
      elseif (mode.eq.3) then
772
         close(fid)
773
      elseif (mode.eq.4) then
774
         ierr = NF90_CLOSE(fid)
775
         IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
776
      elseif (mode.eq.5) then
777
         ierr = NF90_CLOSE(fid)
778
         IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
779
      elseif (mode.eq.6) then
780
         write(fid,"(A)") '</Document>'
781
         write(fid,"(A)") '</kml>'
782
         close(abs(fid))
783
      endif
784
 
785
      end
786
 
787
 
788
c     ----------------------------------------------------------------
789
c     Determine the mode of a trajectory file
790
c     ----------------------------------------------------------------
791
 
792
      subroutine mode_tra(mode,filename)
793
 
794
      use netcdf
795
      implicit none
796
 
797
c     Declaration of subroutine parameters
798
      integer        mode
799
      character*80   filename
800
 
801
c     Auxiliary variables
802
      integer        len
803
      character      char0,char1,char2,char3,char4
804
 
805
c     Get mode
806
      mode=-1
807
 
808
      len  = len_trim(filename)
809
 
810
      char0 = filename((len-1):(len-1))
811
      char1 = filename(len:len)
812
 
813
      if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1
814
      if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2
815
      if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3
816
      if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4
817
      if ( (char0.eq.'.').and.(char1.eq.'5') ) mode=5
818
      if ( (char0.eq.'.').and.(char1.eq.'6') ) mode=6
819
 
820
      if ( mode.gt.0 ) return
821
 
822
c     Mode specified by appendix
823
      char0 = filename((len-3):(len-3))
824
      char1 = filename((len-2):(len-2))
825
      char2 = filename((len-1):(len-1))
826
      char3 = filename(len:len)
827
      if ( (char1.eq.'.').and.(char2.eq.'l').and.(char3.eq.'s') ) mode=1
828
      if ( (char1.eq.'.').and.(char2.eq.'t').and.(char3.eq.'i') ) mode=2
829
      if ( (char1.eq.'.').and.(char2.eq.'d').and.(char3.eq.'u') ) mode=3
830
      if ( (char1.eq.'.').and.(char2.eq.'n').and.(char3.eq.'c') ) mode=4
831
      if ( (char1.eq.'.').and.(char2.eq.'o').and.(char3.eq.'l') ) mode=5
832
 
833
      if ( (char0.eq.'.').and.(char1.eq.'k').and.
834
     >                        (char2.eq.'m').and.
835
     >                        (char3.eq.'l') ) mode = 6
836
 
837
      end
838
 
839
 
840
c     ----------------------------------------------------------------
841
c     Get dimension of a trajectory file
842
c     ----------------------------------------------------------------
843
 
844
      subroutine info_tra(filename,ntra,ntim,ncol,mode)
845
 
846
      use netcdf
847
      implicit none
848
 
849
c     Declaration of subroutine parameters
850
      integer      fid
851
      character*80 filename
852
      integer      mode
853
      integer      ntra,ntim,ncol
854
 
855
c     Auxiliary variables
856
      integer       vardim(4)
857
      real          varmin(4),varmax(4),stag(4)
858
      real          mdv
859
      character*80  cfn
860
      integer       ierr
861
      integer       i,ndim
862
      character*80  vars(100)
863
      integer       nvars
864
      integer       ntimes
865
      real          times(100)
866
      character*500 str
867
      integer       nline0,nline1,nline2
868
      integer       isstr,isok
869
      character     ch
870
      integer       ntraID,ntimID
871
 
872
c     Open file
873
      if (mode.eq.1) then
874
         fid=10
875
         open(fid,file=filename)
876
      elseif (mode.eq.2) then
877
         fid=10
878
         open(fid,file=filename)
879
      elseif (mode.eq.3) then
880
         fid=10
881
         open(fid,file=filename,form='unformatted')
882
      elseif (mode.eq.4) then
883
         ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
884
         IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
885
      elseif (mode.eq.5) then
886
         ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
887
         IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
888
      endif
889
 
890
c     Get dimension information
891
      if ( (mode.eq.1).or.(mode.eq.2) ) then
892
         read(fid,*)
893
         read(fid,*)
894
         read(fid,'(a500)') str
895
         read(fid,*)
896
 
897
c        Get the number of columns
898
         isstr=0
899
         ncol =0
900
         do i=1,500
901
            ch = str(i:i)
902
            if ( (isstr.eq.0).and.(ch.ne.' ') ) then
903
               isstr=1
904
            elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
905
               isstr=0
906
               ncol=ncol+1
907
            endif
908
         enddo
909
 
910
c        Get the first data block
911
         nline0  = 5
912
         nline1  = 5
913
         read(fid,*)
914
 100     read(fid,'(a500)',end=110) str
915
         if (str.ne.'') then
916
            nline1 = nline1 + 1
917
            goto 100
918
         endif
919
 110     continue
920
 
921
c        Get the total numbers of lines in the data block
922
         nline2 = nline1
923
 120     read(fid,*,end=130)
924
         nline2 = nline2 + 1
925
         goto 120
926
 130     nline2 = nline2 + 1
927
 
928
c        Set the dimensions
929
         if (mode.eq.1) then
930
            ntim = nline1 - nline0
931
            ntra = (nline2-nline0+1)/(ntim+1)
932
         else
933
            ntra = nline1 - nline0
934
            ntim = (nline2-nline0+1)/(ntra+1)
935
         endif
936
 
937
      elseif (mode.eq.3) then
938
         read(fid) ntra,ntim,ncol
939
 
940
      elseif (mode.eq.4) then
941
 
5 michaesp 942
         call input_getvars(fid,vars,ncol,ierr)
3 michaesp 943
 
944
         ierr = nf90_inq_dimid(fid,'ntra', ntraID)
945
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
946
         ierr = nf90_inquire_dimension(fid, ntraID, len = ntra)
947
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
948
 
949
         ierr = nf90_inq_dimid(fid,'ntim', ntimID)
950
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
951
         ierr = nf90_inquire_dimension(fid, ntimID, len = ntim)
952
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
953
 
954
      elseif (mode.eq.5) then
955
 
5 michaesp 956
         call input_getvars(fid,vars,ncol,ierr)
3 michaesp 957
 
958
         ierr = nf90_inq_dimid(fid,'id', ntraID)
959
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
960
         ierr = nf90_inquire_dimension(fid, ntraID, len = ntra)
961
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
962
 
963
         ierr = nf90_inq_dimid(fid,'time', ntimID)
964
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
965
         ierr = nf90_inquire_dimension(fid, ntimID, len = ntim)
966
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
967
 
968
      endif
969
 
970
c     Close file
971
      if (mode.eq.1) then
972
         close(fid)
973
      elseif (mode.eq.2) then
974
         close(fid)
975
      elseif (mode.eq.3) then
976
         close(fid)
977
      elseif (mode.eq.4) then
978
          ierr = NF90_CLOSE(fid)
979
          IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
980
      elseif (mode.eq.5) then
981
         ierr = NF90_CLOSE(fid)
982
         IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
983
      endif
984
 
985
      end
986
 
987
c     ------------------------------------------------------------
988
c     Get a list of variables on netCDF file
989
c     ------------------------------------------------------------
990
 
991
      subroutine input_getvars(fid,vnam,nvars)
992
 
993
c     List of variables on netCDF file
994
 
995
      use netcdf
996
      implicit none
997
 
998
c     Declaration of subroutine parameters
999
      integer      fid
1000
      integer      nvars
1001
      character*80 vnam(200)
1002
 
1003
c     Auxiliary variables
1004
      integer ierr
1005
      integer i
1006
      integer nDims,nGlobalAtts,unlimdimid
1007
 
1008
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
1009
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
1010
 
1011
      do i=1,nVars
1012
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
1013
      enddo
1014
 
1015
      end
1016