Subversion Repositories lagranto.icon

Rev

Go to most recent revision | Details | 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
 
219
         call input_getvars(fid,nvars,vars,ierr)
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
 
245
         call input_getvars(fid,nvars,vars,ierr)
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
363
 
364
         call input_getvars(fid,nvars,vars,ierr)
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
 
536
         call  input_getvars(fid,ncol,vars,ierr)
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
 
554
         call  input_getvars(fid,ncol,vars,ierr)
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
674
 
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)
696
 
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)
708
 
709
         call input_getvars(fid,nvars,vars_on_file,ierr)
710
         if ( nvars.ne.ncol ) then
711
            print*,' ERROR: wrong number of fields on netCDF....'
712
            stop
713
         endif
714
         do i=1,nvars
715
            if ( vars(i).ne.vars_on_file(i) ) then
716
               print*,' ERROR: Field order wrong on netCDF...'
717
               stop
718
            endif
719
         enddo
720
 
721
c     Write COSMO online netcdf format (mode=5)
722
      elseif (mode.eq.5) then
723
         print*," ERROR: writing online format not supported"
724
         stop
725
 
726
c     Write KML format (mode=5)
727
      elseif (mode.eq.6) then
728
 
729
      write(fid,"(A)") '<?xml version="1.0" encoding="UTF-8"?>'
730
      write(fid,"(A)") '<kml xmlns="http://www.opengis.net/kml/2.2">'
731
      write(fid,"(A)") '<Document>'
732
      write(fid,"(A)") '<name>Paths</name>'
733
      write(fid,"(A)") '<Style id="yellowLineGreenPoly">'
734
      write(fid,"(A)") '<LineStyle>'
735
c      write(fid,*) '<color>7f00ffff</color>'    ! Yellow
736
      write(fid,"(A)") '<color>500A0A0A</color>'     ! Black
737
      write(fid,"(A)") '<width>4</width>'
738
      write(fid,"(A)") '</LineStyle>'
739
      write(fid,"(A)") '<PolyStyle>'
740
      write(fid,"(A)") '<color>7f00ff00</color>'
741
      write(fid,"(A)") '</PolyStyle>'
742
      write(fid,"(A)") '</Style>'
743
 
744
      endif
745
 
746
      end
747
 
748
 
749
c     ----------------------------------------------------------------
750
c     Close a trajectory file
751
c     ----------------------------------------------------------------
752
 
753
      subroutine close_tra(fid,mode)
754
 
755
      use netcdf
756
      implicit none
757
 
758
c     Declaration of subroutine parameters
759
      integer      fid
760
      integer      mode
761
 
762
c     Auxiliary variables
763
      integer      ierr
764
 
765
c     Close file
766
      if (mode.eq.1) then
767
         close(abs(fid))
768
      elseif (mode.eq.2) then
769
         close(abs(fid))
770
      elseif (mode.eq.3) then
771
         close(fid)
772
      elseif (mode.eq.4) then
773
         ierr = NF90_CLOSE(fid)
774
         IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
775
      elseif (mode.eq.5) then
776
         ierr = NF90_CLOSE(fid)
777
         IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
778
      elseif (mode.eq.6) then
779
         write(fid,"(A)") '</Document>'
780
         write(fid,"(A)") '</kml>'
781
         close(abs(fid))
782
      endif
783
 
784
      end
785
 
786
 
787
c     ----------------------------------------------------------------
788
c     Determine the mode of a trajectory file
789
c     ----------------------------------------------------------------
790
 
791
      subroutine mode_tra(mode,filename)
792
 
793
      use netcdf
794
      implicit none
795
 
796
c     Declaration of subroutine parameters
797
      integer        mode
798
      character*80   filename
799
 
800
c     Auxiliary variables
801
      integer        len
802
      character      char0,char1,char2,char3,char4
803
 
804
c     Get mode
805
      mode=-1
806
 
807
      len  = len_trim(filename)
808
 
809
      char0 = filename((len-1):(len-1))
810
      char1 = filename(len:len)
811
 
812
      if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1
813
      if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2
814
      if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3
815
      if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4
816
      if ( (char0.eq.'.').and.(char1.eq.'5') ) mode=5
817
      if ( (char0.eq.'.').and.(char1.eq.'6') ) mode=6
818
 
819
      if ( mode.gt.0 ) return
820
 
821
c     Mode specified by appendix
822
      char0 = filename((len-3):(len-3))
823
      char1 = filename((len-2):(len-2))
824
      char2 = filename((len-1):(len-1))
825
      char3 = filename(len:len)
826
      if ( (char1.eq.'.').and.(char2.eq.'l').and.(char3.eq.'s') ) mode=1
827
      if ( (char1.eq.'.').and.(char2.eq.'t').and.(char3.eq.'i') ) mode=2
828
      if ( (char1.eq.'.').and.(char2.eq.'d').and.(char3.eq.'u') ) mode=3
829
      if ( (char1.eq.'.').and.(char2.eq.'n').and.(char3.eq.'c') ) mode=4
830
      if ( (char1.eq.'.').and.(char2.eq.'o').and.(char3.eq.'l') ) mode=5
831
 
832
      if ( (char0.eq.'.').and.(char1.eq.'k').and.
833
     >                        (char2.eq.'m').and.
834
     >                        (char3.eq.'l') ) mode = 6
835
 
836
      end
837
 
838
 
839
c     ----------------------------------------------------------------
840
c     Get dimension of a trajectory file
841
c     ----------------------------------------------------------------
842
 
843
      subroutine info_tra(filename,ntra,ntim,ncol,mode)
844
 
845
      use netcdf
846
      implicit none
847
 
848
c     Declaration of subroutine parameters
849
      integer      fid
850
      character*80 filename
851
      integer      mode
852
      integer      ntra,ntim,ncol
853
 
854
c     Auxiliary variables
855
      integer       vardim(4)
856
      real          varmin(4),varmax(4),stag(4)
857
      real          mdv
858
      character*80  cfn
859
      integer       ierr
860
      integer       i,ndim
861
      character*80  vars(100)
862
      integer       nvars
863
      integer       ntimes
864
      real          times(100)
865
      character*500 str
866
      integer       nline0,nline1,nline2
867
      integer       isstr,isok
868
      character     ch
869
      integer       ntraID,ntimID
870
 
871
c     Open file
872
      if (mode.eq.1) then
873
         fid=10
874
         open(fid,file=filename)
875
      elseif (mode.eq.2) then
876
         fid=10
877
         open(fid,file=filename)
878
      elseif (mode.eq.3) then
879
         fid=10
880
         open(fid,file=filename,form='unformatted')
881
      elseif (mode.eq.4) then
882
         ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
883
         IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
884
      elseif (mode.eq.5) then
885
         ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
886
         IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
887
      endif
888
 
889
c     Get dimension information
890
      if ( (mode.eq.1).or.(mode.eq.2) ) then
891
         read(fid,*)
892
         read(fid,*)
893
         read(fid,'(a500)') str
894
         read(fid,*)
895
 
896
c        Get the number of columns
897
         isstr=0
898
         ncol =0
899
         do i=1,500
900
            ch = str(i:i)
901
            if ( (isstr.eq.0).and.(ch.ne.' ') ) then
902
               isstr=1
903
            elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
904
               isstr=0
905
               ncol=ncol+1
906
            endif
907
         enddo
908
 
909
c        Get the first data block
910
         nline0  = 5
911
         nline1  = 5
912
         read(fid,*)
913
 100     read(fid,'(a500)',end=110) str
914
         if (str.ne.'') then
915
            nline1 = nline1 + 1
916
            goto 100
917
         endif
918
 110     continue
919
 
920
c        Get the total numbers of lines in the data block
921
         nline2 = nline1
922
 120     read(fid,*,end=130)
923
         nline2 = nline2 + 1
924
         goto 120
925
 130     nline2 = nline2 + 1
926
 
927
c        Set the dimensions
928
         if (mode.eq.1) then
929
            ntim = nline1 - nline0
930
            ntra = (nline2-nline0+1)/(ntim+1)
931
         else
932
            ntra = nline1 - nline0
933
            ntim = (nline2-nline0+1)/(ntra+1)
934
         endif
935
 
936
      elseif (mode.eq.3) then
937
         read(fid) ntra,ntim,ncol
938
 
939
      elseif (mode.eq.4) then
940
 
941
         call input_getvars(fid,ncol,vars,ierr)
942
 
943
         ierr = nf90_inq_dimid(fid,'ntra', ntraID)
944
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
945
         ierr = nf90_inquire_dimension(fid, ntraID, len = ntra)
946
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
947
 
948
         ierr = nf90_inq_dimid(fid,'ntim', ntimID)
949
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
950
         ierr = nf90_inquire_dimension(fid, ntimID, len = ntim)
951
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
952
 
953
      elseif (mode.eq.5) then
954
 
955
         call input_getvars(fid,ncol,vars,ierr)
956
 
957
         ierr = nf90_inq_dimid(fid,'id', ntraID)
958
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
959
         ierr = nf90_inquire_dimension(fid, ntraID, len = ntra)
960
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
961
 
962
         ierr = nf90_inq_dimid(fid,'time', ntimID)
963
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
964
         ierr = nf90_inquire_dimension(fid, ntimID, len = ntim)
965
         IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
966
 
967
      endif
968
 
969
c     Close file
970
      if (mode.eq.1) then
971
         close(fid)
972
      elseif (mode.eq.2) then
973
         close(fid)
974
      elseif (mode.eq.3) then
975
         close(fid)
976
      elseif (mode.eq.4) then
977
          ierr = NF90_CLOSE(fid)
978
          IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
979
      elseif (mode.eq.5) then
980
         ierr = NF90_CLOSE(fid)
981
         IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
982
      endif
983
 
984
      end
985
 
986
c     ------------------------------------------------------------
987
c     Get a list of variables on netCDF file
988
c     ------------------------------------------------------------
989
 
990
      subroutine input_getvars(fid,vnam,nvars)
991
 
992
c     List of variables on netCDF file
993
 
994
      use netcdf
995
      implicit none
996
 
997
c     Declaration of subroutine parameters
998
      integer      fid
999
      integer      nvars
1000
      character*80 vnam(200)
1001
 
1002
c     Auxiliary variables
1003
      integer ierr
1004
      integer i
1005
      integer nDims,nGlobalAtts,unlimdimid
1006
 
1007
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
1008
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
1009
 
1010
      do i=1,nVars
1011
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
1012
      enddo
1013
 
1014
      end
1015