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