1      PROGRAM BUFR
2C
3C**** *BUFR*
4C
5C
6C     PURPOSE.
7C     --------
8C         Repacks Bufr multisubset data into single Bufr
9C         messages creating appropriate RDB key.
10C
11C
12C**   INTERFACE.
13C     ----------
14C
15C          NONE.
16C
17C     METHOD.
18C     -------
19C
20C          NONE.
21C
22C
23C     EXTERNALS.
24C     ----------
25C
26C         CALL BUSEL
27C         CALL BUFREX
28C         CALL BUFREN
29C         CALL BUPRS0
30C         CALL BUPRS1
31C         CALL BUPRS2
32C         CALL BUPRS3
33C         CALL BUPRT
34C         CALL BUUKEY
35C
36C     REFERENCE.
37C     ----------
38C
39C          NONE.
40C
41C     AUTHOR.
42C     -------
43C
44C          M. DRAGOSAVAC    *ECMWF*       15/02/95.
45C
46C
47C     MODIFICATIONS.
48C     --------------
49C
50C          NONE.
51C
52C
53      IMPLICIT LOGICAL(L,O,G), CHARACTER*8(C,H,Y)
54C
55      PARAMETER(JSUP =   9,JSEC0=   3,JSEC1= 40,JSEC2=4096 ,JSEC3= 4,
56     1       JSEC4=   2,JELEM=160000,JSUBS=400,JCVAL=150 ,JBUFL=512000,
57#ifdef JBPW_64
58     2          JBPW =  64,JTAB =3000,JCTAB=120,JCTST=1800,JCTEXT=1200,
59#else
60     2          JBPW =  32,JTAB =3000,JCTAB=120,JCTST=1800,JCTEXT=1200,
61#endif
62     3          JWORK=4096000,JKEY=46,JBYTE=512000)
63C
64      PARAMETER (KELEM=20000)
65      PARAMETER (KVALS=360000)
66C
67      DIMENSION KBUFF(JBUFL)
68      DIMENSION KBUFR(JBUFL)
69      DIMENSION KSUP(JSUP)  ,KSEC0(JSEC0),KSEC1(JSEC1)
70      DIMENSION KSEC2(JSEC2),KSEC3(JSEC3),KSEC4(JSEC4)
71      DIMENSION KEY  (JKEY),KREQ(2)
72C
73      REAL*8 VALUES(KVALS),VALUE(KVALS)
74      REAL*8 RQV(KELEM)
75      REAL*8 RVIND, EPS
76      REAL*8 RLA,RLO,ZEN_ANG
77      DIMENSION KTDLST(KELEM),KTDEXP(KELEM),KRQ(KELEM)
78      DIMENSION KTDLST1(KELEM)
79      DIMENSION KDATA(200)
80      DIMENSION IOUT(12800)
81C
82      CHARACTER*256 CF(100),COUT,CFIN
83      CHARACTER*64 CNAMES(KELEM)
84      CHARACTER*24 CUNITS(KELEM)
85      CHARACTER*80 CVALS(KVALS)
86      CHARACTER*80 YENC
87      CHARACTER*256 CARG(10)
88c
89cs      EXTERNAL GETARG,sat_zenith_angle
90      EXTERNAL sat_zenith_angle
91C
92C     ------------------------------------------------------------------
93C*          1. INITIALIZE CONSTANTS AND VARIABLES.
94C              -----------------------------------
95 100  CONTINUE
96C
97C     Missing value indicator
98C
99      itlen=6400
100      itl=0
101      jz=0
102      nbytes=jbpw/8
103      RVIND=1.7D38
104      EPS=1.0D08
105      nvind=2147483647
106      iobs=0
107      KRQL=0
108      NR=0
109      KREQ(1)=0
110      KREQ(2)=0
111      do 102 i=1,kelem
112      rqv(i)=rvind
113      krq(i)=nvind
114 102  continue
115c
116C     Input file names
117C
118      narg=iargc()
119      if(narg.lt.4) then
120         print*,'Usage -- bufr_repack -i infiles -o outfile'
121         stop
122      end if
123      nfile=narg
124c
125      do 104 j=1,narg
126      call getarg(j,carg(j))
127 104  continue
128
129      ii=0
130      io=0
131      do 105 j=1,narg
132      if(carg(j).eq.'-i') then
133         in=j
134      elseif(carg(j).eq.'-o') then
135         io=j
136      end if
137 105  continue
138      if(io.eq.0.or.in.eq.0) then
139          print*,'Usage -- bufr_repack -i infiles -o outfile'
140          stop
141      end if
142c
143      cout=carg(io+1)
144c
145      if(io.lt.in) then
146         ist=in+1
147         iend=narg
148      else
149         ist=in+1
150         iend=io-1
151      end if
152c
153      jj=index(cout,' ')
154c
155      CALL PBOPEN(IUNIT1,cout(1:jj),'w',IRET)
156      IF(IRET.EQ.-1) STOP 'open failed on bufr.dat'
157      IF(IRET.EQ.-2) STOP 'Invalid file name'
158      IF(IRET.EQ.-3) STOP 'Invalid open mode specified'
159c
160      do 101 ii=ist,iend
161      cfin=carg(ii)
162      iln=index(cfin,' ')
163C
164C*          1.2 OPEN FILE CONTAINING BUFR DATA.
165C               -------------------------------
166 120  CONTINUE
167C
168      IRET=0
169      CALL PBOPEN(IUNIT,CFIN(1:iln),'r',IRET)
170      IF(IRET.EQ.-1) STOP 'open failed'
171      IF(IRET.EQ.-2) STOP 'Invalid file name'
172      IF(IRET.EQ.-3) STOP 'Invalid open mode specified'
173c
174C     -----------------------------------------------------------------
175C*          2. SET REQUEST FOR EXPANSION.
176C              --------------------------
177 200  CONTINUE
178C
179      OPRT=.FALSE.
180      OENC=.true.
181      NCOM=1
182      OCOMP=.FALSE.
183      NR=1
184      OSEC3=.FALSE.
185C
186C*          2.1 SET REQUEST FOR PARTIAL EXPANSION.
187C               ----------------------------------
188 210  CONTINUE
189C
190C     set variable to pack big values as missing value indicator
191C
192C      KPMISS=1
193C      KPRUS=0
194C      CALL BUPRQ(KPMISS,KPRUS)
195C
196C     -----------------------------------------------------------------
197C*          3.  READ BUFR MESSAGE.
198C               ------------------
199 300  CONTINUE
200C
201      IERR=0
202      KBUFL=0
203C
204      IRET=0
205      CALL PBBUFR(IUNIT,KBUFF,JBYTE,KBUFL,IRET)
206      IF(IRET.EQ.-1) THEN
207          go to 900
208      END IF
209      IF(IRET.EQ.-2) STOP 'File handling problem'
210      IF(IRET.EQ.-3) STOP 'Array too small for product'
211C
212      N=N+1
213      ikbufl=kbufl
214      KBUFL=KBUFL/nbytes+1
215      IF(N.LT.NR) GO TO 300
216C
217C     -----------------------------------------------------------------
218C*          4. EXPAND BUFR MESSAGE.
219C              --------------------
220 400  CONTINUE
221C
222      CALL BUS012(KBUFL,KBUFF,KSUP,KSEC0,KSEC1,KSEC2,KERR)
223      IF(KERR.NE.0) THEN
224         PRINT*,'Error in BUS012: ',KERR
225         PRINT*,' BUFR MESSAGE NUMBER ',N,' CORRUPTED.'
226         KERR=0
227         GO TO 300
228      END IF
229C
230      IF(KSUP(6).GT.1) THEN
231         KEL=KVALS/KSUP(6)
232         IF(KEL.GT.KELEM) KEL=KELEM
233      ELSE
234         KEL=KELEM
235      END IF
236C
237      CALL BUFREX(KBUFL,KBUFF,KSUP,KSEC0 ,KSEC1,KSEC2 ,KSEC3 ,KSEC4,
238     1            KEL,CNAMES,CUNITS,KVALS,VALUES,CVALS,IERR)
239C
240      IF(IERR.NE.0) then
241         if(ierr.eq.45) go to 300
242         call exit(2)
243      end if
244      iobs=iobs+ksec3(3)
245C
246      CALL BUSEL(KTDLEN,KTDLST,KTDEXL,KTDEXP,KERR)
247      IF(KERR.NE.0) CALL EXIT(2)
248C
249C*          4.1 PRINT CONTENT OF EXPANDED DATA.
250C               -------------------------------
251 410  CONTINUE
252C
253      IF(.NOT.OPRT) GO TO 500
254      IF(.NOT.OSEC3) GO TO 450
255C
256C*          4.2 PRINT SECTION ZERO OF BUFR MESSAGE.
257C               -----------------------------------
258 420  CONTINUE
259C
260
261      CALL BUPRS0(KSEC0)
262C
263C*          4.3 PRINT SECTION ONE OF BUFR MESSAGE.
264C               -----------------------------------
265 430  CONTINUE
266C
267      CALL BUPRS1(KSEC1)
268C
269C
270C*          4.4 PRINT SECTION TWO OF BUFR MESSAGE.
271C               -----------------------------------
272 440  CONTINUE
273c
274C              AT ECMWF SECTION 2 CONTAINS RDB KEY.
275C              SO UNPACK KEY
276C
277      CALL BUUKEY(KSEC1,KSEC2,KEY,KSUP,KERR)
278C
279C              PRINT KEY
280C
281      CALL BUPRS2(KSUP ,KEY)
282C
283C*          4.5 PRINT SECTION 3 OF BUFR MESSAGE.
284C               -----------------------------------
285 450  CONTINUE
286C
287C               FIRST GET DATA DESCRIPTORS
288C
289      CALL BUSEL(KTDLEN,KTDLST,KTDEXL,KTDEXP,KERR)
290      IF(KERR.NE.0) CALL EXIT(2)
291C
292C               PRINT  CONTENT
293C
294      IF(OSEC3) THEN
295         CALL BUPRS3(KSEC3,KTDLEN,KTDLST,KTDEXL,KTDEXP,KEL,CNAMES)
296      END IF
297C
298C*         4.6 PRINT SECTION 4 (DATA).
299C              -----------------------
300 460  CONTINUE
301C
302C          IN THE CASE OF MANY SUBSETS DEFINE RANGE OF SUBSETS
303C
304      IF(.NOT.OO) THEN
305      WRITE(*,'(a,$)') ' STARTING SUBSET TO BE PRINTED : '
306      READ(*,'(BN,I4)')   IST
307      WRITE(*,'(a,$)') ' ENDING SUBSET TO BE PRINTED : '
308      READ(*,'(BN,I4)')   IEND
309      OO=.false.
310      END IF
311C
312C              PRINT DATA
313C
314      ICODE=0
315      CALL BUPRT(ICODE,IST,IEND,KEL,CNAMES,CUNITS,CVALS,
316     1           KVALS,VALUES,KSUP,KSEC1,IERR)
317c
318c              Resolve bit maps
319C
320      if(iend.gt.ksec3(3)) iend=ksec3(3)
321c
322      DO 461 IK=IST,IEND
323C
324      CALL BUBOX(IK,KSUP,KEL,KTDEXP,CNAMES,CUNITS,KVALS,VALUES,
325     1           KBOX,KAPP,KLEN,KBOXR,VALS,CBOXN,CBOXU,IERR)
326C
327      CALL BUPRTBOX(KBOX,KAPP,KLEN,KBOXR,VALS,CBOXN,CBOXU)
328C
329 461  CONTINUE
330C
331C     -----------------------------------------------------------------
332C*          5. COLLECT DATA FOR REPACKING.
333C              ---------------------------
334 500  CONTINUE
335C
336      IF(.NOT.OENC) GO TO 300
337C
338C               FIRST GET DATA DESCRIPTORS
339C
340      CALL BUSEL(KTDLEN,KTDLST,KTDEXL,KTDEXP,KERR)
341      IF(KERR.NE.0) CALL EXIT(2)
342C
343C     -----------------------------------------------------------------
344C*          6. PACK BUFR MESSAGE BACK INTO BUFR.
345C              ---------------------------------
346 600  CONTINUE
347C
348
349      KKK=0
350      KBUFL=JBUFL
351C
352C     GET REPLICATION FACTORS
353C
354      KK=0
355      DO 601 K=1,KSUP(5)
356      IF(KTDEXP(K).EQ.31001.OR.KTDEXP(K).EQ.31002.OR.
357     1   KTDEXP(K).EQ.31000) THEN
358         KK=KK+1
359         KDATA(KK)=NINT(VALUES(K))
360      END IF
361 601  CONTINUE
362C
363      KDLEN=2
364      IF(KK.NE.0) KDLEN=KK
365C
366C           6.1 Convert subtype 88 into 89
367C
368 610  CONTINUE
369C
370      do i=1,KSUP(6)
371      do j=1,222
372         iv=j+(i-1)*222
373         value(iv)=rvind
374      end do
375      end do
376c
377      do i=1,KSUP(6)
378      do j=1,KSUP(5)
379        ij=1+(i-1)*kel
380        iv=1+(i-1)*222
381        id=nint(values(ij))
382        value(iv)=values(ij)
383c
384        ij=3+(i-1)*kel
385        iv=6+(i-1)*222
386        iyear=nint(values(ij))
387        value(iv)=values(ij)
388c
389        ij=4+(i-1)*kel
390        iv=7+(i-1)*222
391        value(iv)=values(ij)
392c
393        ij=5+(i-1)*kel
394        iv=8+(i-1)*222
395        value(iv)=values(ij)
396c
397        ij=6+(i-1)*kel
398        iv=9+(i-1)*222
399        value(iv)=values(ij)
400c
401        ij=7+(i-1)*kel
402        iv=10+(i-1)*222
403        value(iv)=values(ij)
404c
405        ij=8+(i-1)*kel
406        iv=11+(i-1)*222
407        value(iv)=values(ij)
408c
409        ij=9+(i-1)*kel
410        iv=12+(i-1)*222
411        rla=values(ij)
412        value(iv)=values(ij)
413c
414        ij=10+(i-1)*kel
415        iv=13+(i-1)*222
416        rlo=values(ij)
417        value(iv)=values(ij)
418c
419        ij=12+(i-1)*kel
420        iv=66+(i-1)*222
421        value(iv)=values(ij)
422c
423        iv=14+(i-1)*222
424        call sat_zenith_angle(id,iyear,rla,rlo,zen_ang)
425        value(iv)=zen_ang
426
427      end do
428      end do
429c
430c     set bit map
431c
432      do i=1,ksup(6)
433      do j=109,214
434        iv=j+(i-1)*222
435        value(iv)=1.0
436        if(j.eq.174) value(iv)=0.0
437      end do
438      end do
439c
440      do j=1,ksup(6)
441        iv=215+(j-1)*222
442        value(iv)=254.
443        iv=216+(j-1)*222
444        value(iv)=1.
445        iv=217+(j-1)*222
446        value(iv)=100.
447        iv=220+(j-1)*222
448        value(iv)=254.
449        iv=221+(j-1)*222
450        value(iv)=1.
451        iv=222+(j-1)*222
452        value(iv)=0.0
453      end do
454c
455      ktdlst1( 1)=001007
456      ktdlst1( 2)=001031
457      ktdlst1( 3)=002196
458      ktdlst1( 4)=002221
459      ktdlst1( 5)=002222
460      ktdlst1( 6)=004001
461      ktdlst1( 7)=004002
462      ktdlst1( 8)=004003
463      ktdlst1( 9)=004004
464      ktdlst1(10)=004005
465      ktdlst1(11)=004006
466      ktdlst1(12)=005001
467      ktdlst1(13)=006001
468      ktdlst1(14)=007024
469      ktdlst1(15)=010002
470      ktdlst1(16)=002252
471      ktdlst1(17)=002023
472      ktdlst1(18)=007004
473      ktdlst1(19)=011001
474      ktdlst1(20)=011002
475      ktdlst1(21)=002197
476      ktdlst1(22)=002198
477      ktdlst1(23)=012193
478      ktdlst1(24)=002197
479      ktdlst1(25)=002198
480      ktdlst1(26)=020193
481      ktdlst1(27)=020194
482      ktdlst1(28)=020012
483      ktdlst1(29)=002197
484      ktdlst1(30)=002198
485      ktdlst1(31)=020193
486      ktdlst1(32)=020194
487      ktdlst1(33)=020012
488      ktdlst1(34)=002197
489      ktdlst1(35)=002198
490      ktdlst1(36)=020193
491      ktdlst1(37)=020194
492      ktdlst1(38)=020012
493      ktdlst1(39)=002197
494      ktdlst1(40)=002198
495      ktdlst1(41)=020193
496      ktdlst1(42)=020194
497      ktdlst1(43)=020012
498      ktdlst1(44)=002197
499      ktdlst1(45)=002198
500      ktdlst1(46)=020193
501      ktdlst1(47)=020194
502      ktdlst1(48)=020012
503      ktdlst1(49)=002197
504      ktdlst1(50)=002198
505      ktdlst1(51)=020193
506      ktdlst1(52)=020194
507      ktdlst1(53)=020012
508      ktdlst1(54)=002252
509      ktdlst1(55)=002199
510      ktdlst1(56)=007004
511      ktdlst1(57)=007004
512      ktdlst1(58)=013003
513      ktdlst1(59)=002252
514      ktdlst1(60)=002254
515      ktdlst1(61)=002251
516      ktdlst1(62)=002197
517      ktdlst1(63)=002198
518      ktdlst1(64)=012195
519      ktdlst1(65)=012196
520      ktdlst1(66)=012063
521      ktdlst1(67)=002252
522      ktdlst1(68)=002254
523      ktdlst1(69)=002251
524      ktdlst1(70)=002197
525      ktdlst1(71)=002198
526      ktdlst1(72)=012195
527      ktdlst1(73)=012196
528      ktdlst1(74)=012063
529      ktdlst1(75)=002252
530      ktdlst1(76)=002254
531      ktdlst1(77)=002251
532      ktdlst1(78)=002197
533      ktdlst1(79)=002198
534      ktdlst1(80)=012195
535      ktdlst1(81)=012196
536      ktdlst1(82)=012063
537      ktdlst1(83)=002252
538      ktdlst1(84)=002254
539      ktdlst1(85)=002251
540      ktdlst1(86)=002197
541      ktdlst1(87)=002198
542      ktdlst1(88)=012195
543      ktdlst1(89)=012196
544      ktdlst1(90)=012063
545      ktdlst1(91)=002252
546      ktdlst1(92)=002254
547      ktdlst1(93)=002251
548      ktdlst1(94)=002197
549      ktdlst1(95)=002198
550      ktdlst1(96)=012195
551      ktdlst1(97)=012196
552      ktdlst1(98)=012063
553      ktdlst1(99)=002252
554      ktdlst1(100)=002254
555      ktdlst1(101)=002251
556      ktdlst1(102)=002197
557      ktdlst1(103)=002198
558      ktdlst1(104)=012195
559      ktdlst1(105)=012196
560      ktdlst1(106)=012063
561      ktdlst1(107)=222000
562      ktdlst1(108)=236000
563      ktdlst1(109)=101106
564      ktdlst1(110)=031031
565      ktdlst1(111)=001031
566      ktdlst1(112)=001032
567      ktdlst1(113)=101001
568      ktdlst1(114)=033007
569      ktdlst1(115)=222000
570      ktdlst1(116)=237000
571      ktdlst1(117)=001031
572      ktdlst1(118)=001032
573      ktdlst1(119)=101001
574      ktdlst1(120)=033252
575c
576      ktdlen=120
577C
578C*          6.2 ENCODE DATA INTO BUFR MESSAGE.
579C               ------------------------------
580 620  CONTINUE
581C
582C     repacking multisubset data one by one
583C
584      CALL BUUKEY(KSEC1,KSEC2,KEY,KSUP,KERR)
585c
586      IF(KSUP(6).EQ.0) THEN
587         print*,'Zero subsets'
588         call exit(2)
589      END IF
590c
591      key(3)=89
592      ksec1(7)=89
593      ksec1(15)=13
594c
595c     Get information for RDB key
596c
597      CALL BUCREKEY(KTDEXP,KSUP,KSEC1,KEY,VALUES,CVALS,KERR)
598      IF(KERR.NE.0) THEN
599         print*,'Error in BUCREKEY.'
600         call exit(2)
601      end if
602c
603c     Pack new RDB key
604c
605      CALL BUPKEY(KEY,KSEC1,KSEC2,KERR)
606      IF(KERR.NE.0) CALL EXIT(2)
607c
608
609      kel1=222
610      if(ksec3(4).eq.0) ksec3(4)=192
611      CALL BUFREN( KSEC0,KSEC1,KSEC2,KSEC3,KSEC4,
612     1             KTDLEN,KTDLST1,KDLEN,KDATA,KEL1,        !ksup(5),
613     2             KVALS,VALUE,CVALS,KBUFL,KBUFR,KERR)
614C
615      IF(KERR.GT.0) THEN
616         PRINT*,'ERROR DURING ENCODING.'
617         CALL EXIT(2)
618      END IF
619C
620C           6.3 WRITE PACKED BUFR MESSAGE INTO FILE.
621C               ------------------------------------
622 630  CONTINUE
623C
624      ikbufl=KBUFL*4
625      CALL PBWRITE(IUNIT1,kbufr,ikbufl,IERR)
626      IF(IERR.LT.0) THEN
627         print*,'Error writing into target file.'
628         call exit(2)
629      END IF
630C
631      GO TO 300
632C     -----------------------------------------------------------------
633C
634 810  CONTINUE
635C
636      WRITE(*,'(1H ,A)') 'OPEN ERROR ON INPUT FILE'
637      GO TO 900
638C
639 800  CONTINUE
640C
641      IF(IRET.EQ.-1) THEN
642         print*,'Number of records processed ',n
643      ELSE
644         print*,' BUFR : error= ',ierr
645      END IF
646C
647 900  CONTINUE
648C
649      CALL PBCLOSE(IUNIT,IRET)
650 101  CONTINUE
651c
652      CALL PBCLOSE(IUNIT1,IRET)
653C
654      END
655      SUBROUTINE SAT_ZENITH_ANGLE(SAT_ID,YEAR,OBS_LAT,OBS_LON,ZEN_ANG)
656
657
658       real*8 obs_lat,obs_lon
659       integer sat_id,year
660       real*8 zen_ang
661
662       parameter ( dtr    = 0.017453292)
663       parameter ( radius = 6371.0)
664       parameter ( pi     = 3.14159)
665       parameter ( height = 35800.0)
666
667       real  sat_lat_rad,sat_lon_rad, sat_lat, sat_lon
668       real  obs_lat_rad,obs_lon_rad
669       real  xx,yy,zz
670       real  arcl, distance
671       real  ang, s
672
673c      begin main routine
674
675       sat_lat = 0.0
676       sat_lon = 0.0
677
678c     Meteosat 5 for INDOEX
679       if (sat_id .eq. 52 .and. year .gt. 1997) sat_lon = 63.0
680
681       sat_lat_rad = sat_lat * dtr
682       sat_lon_rad = sat_lon * dtr
683
684       obs_lat_rad = obs_lat * dtr
685       obs_lon_rad = obs_lon * dtr
686
687       xx = cos(obs_lat_rad)*cos(obs_lon_rad)
688     1       - cos(sat_lat_rad)*cos(sat_lon_rad)
689       yy = cos(obs_lat_rad)*sin(obs_lon_rad)
690     1       - cos(sat_lat_rad)*sin(sat_lon_rad)
691       zz = sin(obs_lat_rad) - sin(sat_lat_rad)
692
693       arcl = 2.0 * asin(sqrt(xx*xx + yy*yy + zz*zz)/2.0)
694       distance = arcl*radius
695
696       ang  = (360.0 * distance/(2.0 * pi * radius)) * dtr
697       s    = sqrt(radius**2.0 + (radius+height)**2.0 -
698     1        2.0 * radius * (radius+height) * cos(ang))
699
700       zen_ang = abs(asin((radius + height) * sin(ang)/s) / dtr )
701
702       return
703      end
704