1 program tc_tracks
2!
3!
4!**** tc_tracks*
5!
6!
7!     Purpose.
8!     --------
9!
10!         Create tropical cyclone tracks
11!
12!
13!**   Interface.
14!     ----------
15!
16!
17!     Method.
18!     -------
19!
20!
21!
22!
23!
24!     Externals.
25!     ----------
26!
27!
28!
29!     Reference.
30!     ----------
31!
32!
33!
34!     Author.
35!     -------
36!
37!          M. Dragosavac    *ECMWF*      22/12/2001
38!
39!
40!     Modifications.
41!     --------------
42!
43!          NONE.
44!---------------------------------------------------------------------------
45 USE bufr_module
46
47 implicit none
48
49! Local variables
50
51  integer              :: iii,narg, iargc, n, nn
52  integer              :: j,ios,ij,ii
53  integer              :: iunit, iret,ierr, ilen
54  integer              :: kdlen
55  integer              :: second, subset
56
57  character(len=256)   :: input_file
58  character(len=256)   :: output_file
59  character(len=256)   :: det_file
60  character(len=256), dimension(10) :: carg
61  character(len=1024)   :: line
62  character(len=3)     :: tc_identifier
63  character(len=10)    :: tc_name
64
65
66  integer             :: increment
67  integer             :: member
68  integer             :: originating_centre
69  integer             :: number_of_steps
70  integer             :: is_missing_step(100)
71
72  real                :: lat_observed
73  real                :: lon_observed
74  real                :: lat_pressure(100)
75  real                :: lon_pressure(100)
76  real                :: lat_wind(100)
77  real                :: lon_wind(100)
78  integer             :: year
79  integer             :: month
80  integer             :: day
81  integer             :: hour
82  integer             :: minute
83  integer             :: key2(46)
84  integer             :: key3(46)
85  integer             :: kerr
86  integer             :: jj                             ! loop variable
87  real                :: wind(100)
88  real                :: pressure(100)
89  real                :: missing
90  integer             :: edition
91
92
93  integer, parameter :: kel=1000
94
95! Subroutine interface
96
97
98! ----------------------------------
99
100!  initialization of variables
101
102  missing=rvind
103  increment=0
104  number_of_steps=0
105  is_missing_step=0
106  lat_pressure=missing
107  lon_pressure=missing
108  lat_wind=missing
109  lon_wind=missing
110  wind=missing
111  pressure=missing
112
113!   Get program arguments
114
115
116    input_file=' '
117    output_file=' '
118
119    narg=iargc()
120
121    if(narg < 3) then
122       print*,'Usage -- tc_tracks_det -4 -i input_file -o det_file '
123       print*,'      -4 option to produce BUFR edition 4'
124       print*,'      input_file   -- input file name'
125       print*,'      det_file     -- determistic file name'
126       stop
127    end if
128
129    do  j=1,narg
130       call getarg(j,carg(j))
131    end do
132
133
134    edition=3
135    do j=1,narg
136     if(carg(j) == '-i') then
137       input_file=carg(j+1)
138     elseif(carg(j) == '-o') then
139       output_file=carg(j+1)
140     elseif(carg(j) == '-4') then
141       edition=4
142     else
143       cycle
144     end if
145    end do
146
147    write(0,*) 'Input  file name='//trim(input_file)
148    write(0,*) 'Deterministic file name='//trim(output_file)
149
150
151!-- Open input and output files
152!   ---------------------------
153
154
155!   Open input file
156
157    open(21,file=trim(input_file),iostat=ios)
158    if(ios /= 0) then
159      print*,'Open error on input file '
160      stop
161    end if
162
163!   Open output file
164
165    iret=0
166    call pbopen(iunit,trim(output_file),'w',iret)
167    if(iret == -1) then
168       print*,'Open failed'
169       stop
170    elseif(iret == -2) then
171       print*,'Invalid input file name'
172       stop
173    elseif(iret == -3) then
174       print*,'Invalid open mode specified'
175       stop
176    end if
177
178
179
180!   Read input file
181
182    subset=0
183    number_of_steps=0
184    values=rvind
185
186    do while( ios == 0 )
187
188    line=' '
189    read(21,'(a)', iostat=ios) line
190!   write(0,*) trim(line)
191
192    if( ios == 0 ) then
193      if(index(line,'originating centre') /= 0 ) then
194         read(line,'(47x,i3)',iostat=ios) originating_centre
195         write(*,'(47x,i3)',iostat=ios) originating_centre
196         if(ios /= 0 ) then
197            print*,'Internal read error for originating_centre'
198            call exit(2)
199         end if
200      elseif(index(line,'TC identifier') /= 0 ) then
201         read(line,'(16x,a3)',iostat=ios) tc_identifier
202         write(*,'(16x,a3)',iostat=ios) tc_identifier
203         if(ios /= 0 ) then
204            print*,'Internal read error for tc_identifier'
205            call exit(2)
206         end if
207      elseif(index(line,'TC name') /= 0 ) then
208         read(line,'(10x,a10)',iostat=ios) tc_name
209         write(*,'(10x,a10)',iostat=ios) tc_name
210         if(ios /= 0 ) then
211            print*,'Internal read error for tc_name'
212            call exit(2)
213         end if
214      elseif(index(line,'YYYYMMDDHHMM') /= 0 ) then
215         read(line,'(36x,i4,i2,i2,i2,i2)',iostat=ios) year, month, day, hour, minute
216         write(*,'(36x,i4,i2,i2,i2,i2)',iostat=ios) year, month, day, hour, minute
217         if(ios /= 0 ) then
218            print*,'Internal read error for date/time'
219            call exit(2)
220         end if
221      elseif(index(line,'obs lat') /= 0 ) then
222         read(line,'(14x,f6.2)',iostat=ios) lat_observed
223         write(*,'(14x,f6.2)',iostat=ios) lat_observed
224         if(ios /= 0 ) then
225            print*,'Internal read error for obs lat'
226            call exit(2)
227         end if
228      elseif(index(line,'obs lon') /= 0 ) then
229         read(line,'(13x,f7.2)',iostat=ios) lon_observed
230         write(*,'(13x,f7.2)',iostat=ios) lon_observed
231         if(ios /= 0 ) then
232            print*,'Internal read error for obs lon'
233            call exit(2)
234         end if
235      elseif(index(line,'member') /= 0 ) then
236         if (iii /= 0 ) then
237            print *,'Internal read error: too many members, only deterministic expected'
238            call exit(2)
239         endif
240         read(line,'(46x,i2)',iostat=ios) member
241         write(*,'(46x,i2)',iostat=ios) member
242         if(ios /= 0 ) then
243            print*,'Internal read error for member'
244            call exit(2)
245         end if
246      elseif(index(line,'number of steps') /= 0 ) then
247         read(line,'(46x,i2)',iostat=ios) number_of_steps
248         write(*,'(46x,i2)',iostat=ios) number_of_steps
249         if(ios /= 0 ) then
250            print*,'Internal read error for number of steps'
251            call exit(2)
252         end if
253
254
255      elseif(index(line,'increment') /= 0 ) then
256         read(line,'(46x,i2)',iostat=ios) increment
257         write(*,'(46x,i2)',iostat=ios) increment
258         if(ios /= 0 ) then
259            print*,'Internal read error for increment'
260            call exit(2)
261         end if
262      elseif(index(line,'missing') /= 0 ) then
263           iii=iii+1
264           is_missing_step(iii)=1
265      else
266         iii=iii+1
267         is_missing_step(iii)=0
268         read(line,'(4x,f6.2,3x,f7.2,3x,f7.2,4x,f6.2,4x,f6.2,3x,f7.2)',iostat=ios) lat_pressure(iii), &
269&             lon_pressure(iii),pressure(iii),wind(iii),lat_wind(iii), &
270&             lon_wind(iii)
271         write(*,'(4x,f6.2,3x,f7.2,3x,f7.2,4x,f6.2,4x,f6.2,3x,f7.2)',iostat=ios) lat_pressure(iii), &
272&             lon_pressure(iii),pressure(iii),wind(iii),lat_wind(iii), &
273&             lon_wind(iii)
274         pressure(iii)=pressure(iii)*100.
275!          print*,' ',iii,'----',lat(iii),lon(iii),pressure(iii),wind(iii)
276           if(ios /= 0 ) then
277              print*,'Internal read error for increment'
278              call exit(2)
279           end if
280      end if
281    else
282       print*,'END of data....................'
283       exit
284    end if
285    end do
286
287    if (number_of_steps .eq. 0) then
288      print *,'number of steps = 0, no output bufr written'
289      stop
290    end if
291
292    cvals(1)=tc_identifier
293    cvals(2)=tc_name
294
295
296    subset=subset+1
297
298    ij=(subset-1)*kel
299
300        n=1+ij
301        values(  n)=98.      ! Originating centre
302        n=2+ij
303        values(  n)=rvind    ! Originating subcentre
304        n=3+ij
305        values(  n)=1.       ! Generation application
306        n=4+ij
307        values(  n)=1003.    ! Storm identifier
308        n=5+ij
309        values(  n)=2010.    ! Storm name
310        n=6+ij
311        values(  n)=2.       ! singular vector
312        n=7+ij
313        values(  n)=member    ! Only deterministic member
314        n=8+ij
315        values(  n)=0        ! unperturbed high resol. control forecast
316!+++++++++++++++++++++
317        n=9+ij
318        values(  n)=year
319        n=10+ij
320        values(  n)=month
321        n=11+ij
322        values(  n)=day
323        n=12+ij
324        values(  n)=hour
325        n=13+ij
326        values(  n)=minute
327!+++++++++++++++++++++++++
328        n=14+ij
329        values(  n)=1.             ! Storm centre
330        n=15+ij
331        values(  n)=lat_observed
332        n=16+ij
333        values(  n)=lon_observed
334        n=17+ij
335        values(  n)=5. ! Location of the storm in the analysis
336        n=18+ij
337        if (is_missing_step(1) == 0) values(  n)=lat_pressure(1)
338        n=19+ij
339        if (is_missing_step(1) == 0) values(  n)=lon_pressure(1)
340        n=20+ij
341        if (is_missing_step(1) == 0) values(  n)=pressure(1)
342        n=21+ij
343        values(  n)=3.
344        n=22+ij
345        values(  n)=lat_wind(1)
346        n=23+ij
347        values(  n)=lon_wind(1)
348        n=24+ij
349        values(  n)=wind(1)
350!+++++++++++++++++++++++++
351        n=25+ij
352        values(  n)=number_of_steps-1       ! delayed replication up to any forecast steps
353
354        nn=25
355        do ii=2,number_of_steps
356          if (is_missing_step(ii) == 1) continue
357          nn=nn+1
358          n=nn+ij
359          values(  n)=4.         ! time significance  008021 forecast
360          nn=nn+1
361          n=nn+ij
362          values(  n)=(ii-1)*increment  ! time displacement
363          nn=nn+1
364          n=nn+ij
365          values(  n)=1.           ! 008005 storm centre
366          nn=nn+1
367          n=nn+ij
368          values(  n)=lat_pressure(ii)        !
369          nn=nn+1
370          n=nn+ij
371          values(  n)=lon_pressure(ii)
372          nn=nn+1
373          n=nn+ij
374          values(  n)=pressure(ii)
375          nn=nn+1
376          n=nn+ij
377          values(  n)=3.           ! 008005  location of max wind
378          nn=nn+1
379          n=nn+ij
380          values(  n)=lat_wind(ii)        !
381          nn=nn+1
382          n=nn+ij
383          values(  n)=lon_wind(ii)
384          nn=nn+1
385          n=nn+ij
386          values(  n)=wind(ii)
387        end do
388
389        ktdlst( 1)=001033
390        ktdlst( 2)=001034
391        ktdlst( 3)=001032
392        ktdlst( 4)=001025
393        ktdlst( 5)=001027
394        ktdlst( 6)=001090
395        ktdlst( 7)=001091
396        ktdlst( 8)=001092
397        ktdlst(09)=301011
398        ktdlst(10)=301012
399        ktdlst(11)=008005
400        ktdlst(12)=005002
401        ktdlst(13)=006002
402        ktdlst(14)=008005
403        ktdlst(15)=005002
404        ktdlst(16)=006002
405        ktdlst(17)=010051
406        ktdlst(18)=008005
407        ktdlst(19)=301023
408        ktdlst(20)=011012
409        ktdlst(21)=108000
410        ktdlst(22)=031001
411        ktdlst(23)=008021
412        ktdlst(24)=004024
413        ktdlst(25)=008005
414        ktdlst(26)=301023
415        ktdlst(27)=010051
416        ktdlst(28)=008005
417        ktdlst(29)=301023
418        ktdlst(30)=011012
419
420        ktdlen=30
421
422!          Create bufr message
423
424        if (edition == 3) then
425            ksec0( 3)=3
426            ksec1( 2)=3
427            ksec1( 1)=18
428        else if (edition == 4 ) then
429            ksec0( 3)=4
430            ksec1( 2)=4
431            ksec1( 1)=22
432        end if
433
434        ksec1( 3)=98
435        ksec1( 4)=1
436        ksec1( 5)=128        ! 128 section 2 exists
437        ksec1( 6)=7
438        ksec1( 7)=32
439        ksec1( 8)=0
440        if(values(9) > 2000.) then
441           ksec1( 9)=nint(values(9))-2000
442        else
443           ksec1( 9)=nint(values(9))-1900
444        end if
445        ksec1(10)=nint(values(10))
446        ksec1(11)=nint(values(11))
447        ksec1(12)=nint(values(12))
448        ksec1(13)=nint(values(13))
449        ksec1(14)=0
450        ksec1(15)=16
451        ksec1(16)=0
452        ksec1(17)=0
453        ksec1(18)=0
454
455        ksec2(1)=52 ! this is required to encode section 2
456
457! define section 2
458! initialise key2 variable
459        do jj=1,46
460          key2(jj)=0
461        enddo
462
463        key2(1)=52                 ! Length of section 2
464        key2(2)=8                  ! RDB type
465        key2(3)=32                 ! RDB subtype
466        key2(4)=nint(values(09))    ! Year
467        key2(5)=nint(values(10))   ! Month
468        key2(6)=nint(values(11))   ! Day
469        key2(7)=nint(values(12))   ! Hour
470        key2(8)=nint(values(13))   ! Minute
471        key2(9)=0                  ! Second
472
473        key2(10)=NINT(lon_observed*100000.+18000000)    ! Longitude 1 location of first storm centre
474        key2(11)=NINT(lat_observed*100000.+9000000)     ! Latitude 1 location of first storm centre
475        key2(12)=NINT(lon_observed*100000.+18000000)    ! Longitude 2 location of first storm centre
476        key2(13)=NINT(lat_observed*100000.+9000000)     ! Latitude 2 location of first storm centre
477        key2(14)=subset    ! no of observations/subsets
478
479! identifier
480        key2(16)=ichar(tc_identifier(1:1))           ! first character of short identifier
481        key2(17)=ichar(tc_identifier(2:2))           ! second character of short identifier
482        key2(18)=ichar(tc_identifier(3:3))           ! third character of short identifier
483        key2(19)=32
484        key2(20)=32
485        key2(21)=32
486        key2(22)=32
487        key2(23)=32
488        key2(24)=32
489
490        ksec3(1)=0
491        ksec3(2)=0
492        ksec3(3)=subset
493        ksec3(4)=0
494        if(ksec3(3).gt.1) ksec3(4)=64
495
496        kdlen=1
497        kdata(1)=number_of_steps-1   ! first step is pert. anal.
498        kbufl=jbufl
499
500! pack rdb key into ksec2 array
501
502        call bupkey(key2,ksec1,ksec2,kerr)
503
504        if(kerr.ne.0) then
505          print*,'BUPKEY: error',kerr
506          call exit(2)
507        end if
508
509! encode
510
511        ierr=0
512        call bufren(ksec0,ksec1,ksec2,ksec3,ksec4, &
513                    ktdlen,ktdlst,kdlen,kdata,kel, &
514                    kvals,values,cvals,kbufl,kbufr,ierr)
515        if(ierr > 0) then
516           print*,'bufren error ', ierr
517           call exit(2)
518        elseif(ierr < 0) then
519           print*,'Encoding return_code=',ierr
520        end if
521
522
523!          Write bufr message into output file
524
525        ilen=kbufl*4
526
527        call pbwrite(iunit,kbufr,ilen,ierr)
528        print*,'length=',ilen,' ierr=',ierr
529        print*,'Bufr message written into output file'
530
531 end program tc_tracks
532