var/convertor/kmabufr/read_bufr.F

   1       PROGRAM BUFR
   2 C
   3       IMPLICIT LOGICAL(L,O,G), CHARACTER*8(C,H,Y)
   4 C
   5       PARAMETER(JSUP =   9,JSEC0=   3,JSEC1= 40,JSEC2= 64 ,JSEC3=    4,
   6      1        JSEC4=   2,JELEM=80000,JSUBS=400,JCVAL=150 ,JBUFL=512000,
   7 #ifdef JBPW_64
   8      2        JBPW =  64,JTAB =1000,JCTAB=120,JCTST=1800,JCTEXT=1200,
   9 #else
  10      2        JBPW =  32,JTAB =1000,JCTAB=120,JCTST=1800,JCTEXT=1200,
  11 #endif
  12      3        JWORK=4096000,JKEY=46,JBYTE=2048000)
  13 C
  14       PARAMETER (KELEM=40000)
  15       PARAMETER (KVALS=4096000)
  16 C
  17       DIMENSION KBUFF(JBUFL)
  18       DIMENSION KBUFR(JBUFL)
  19       DIMENSION KSUP(JSUP)  ,KSEC0(JSEC0),KSEC1(JSEC1)
  20       DIMENSION KSEC2(JSEC2),KSEC3(JSEC3),KSEC4(JSEC4)
  21       DIMENSION KEY  (JKEY),KREQ(2)
  22       DIMENSION NREQUEST(2)
  23 C
  24       REAL*8 VALUES(KVALS),VALUE(KVALS)
  25       REAL*8 VALS(KVALS)
  26       REAL*8 RQV(KELEM)
  27       REAL*8 RVIND,EPS
  28 C
  29       DIMENSION KTDLST(JELEM),KTDEXP(JELEM),KRQ(KELEM)
  30       DIMENSION KDATA(200),KBOXR(JELEM*4)
  31
  32 C
  33       CHARACTER*256 CFIN,COUT,CARG(4)
  34       CHARACTER*64 CNAMES(KELEM),CBOXN(JELEM*4)
  35       CHARACTER*24 CUNITS(KELEM),CBOXU(JELEM*4)
  36       CHARACTER*80 CVALS(kelem)
  37       CHARACTER*80 CVAL(kelem)
  38 C
  39 C
  40 C     ------------------------------------------------------------------
  41 C*          1. INITIALIZE CONSTANTS AND VARIABLES.
  42 C              -----------------------------------
  43  100  CONTINUE
  44 C
  45 C     MISSING VALUE INDICATOR
  46 C
  47       RVIND=1.7E38
  48       NVIND=2147483647
  49 C
  50       NBYTPW=JBPW/8
  51       IOBS=0
  52       EPS=1.E-8
  53       NPACK=0
  54       N=0
  55       OO=.FALSE.
  56       KKK=0
  57 C
  58       ict = 0    ! report counter
  59       open(3,file='littler',status='unknown',form='formatted')
  60       open(12,file='flist',status='old',form='formatted')
  61    87 continue
  62       read(12,'(a256)',end=988,err=987) cfin
  63
  64 c determine fm number of the data based on file name
  65 c
  66       iln = index(cfin,'_') - 1
  67       if (cfin(1:iln) .eq. 'temp') then
  68         ifm = 35
  69       elseif (cfin(1:iln) .eq. 'pilot') then
  70         ifm = 32
  71       elseif (cfin(1:iln) .eq. 'airep') then
  72         ifm = 96
  73       elseif (cfin(1:iln) .eq. 'acars') then
  74         ifm = 96
  75       elseif (cfin(1:iln) .eq. 'ship') then
  76         ifm = 13
  77       elseif (cfin(1:iln) .eq. 'buoy') then
  78         ifm = 18
  79       elseif (cfin(1:iln) .eq. 'synop') then
  80         ifm = 12
  81       elseif (cfin(1:iln) .eq. 'aws') then
  82         ifm = 15
  83       elseif (cfin(1:iln) .eq. 'satob') then
  84         ifm = 88
  85       elseif (cfin(1:iln) .eq. 'satem') then
  86         ifm = 87
  87       else
  88         write(6,*) 'Observation type ',cfin(1:iln),' is not supported.
  89      &Skipping.'
  90         go to 87
  91       endif
  92
  93       ILN=INDEX(CFIN,' ')
  94       ILN=ILN-1
  95
  96 c     write(6,*) 'Processing FM-',ifm,' from file ',cfin(1:iln)
  97 C
  98 C     SET REQUEST FOR PARTIAL EXPANSION
  99 C
 100       KRQL=0
 101       NR=0
 102       KREQ(1)=0
 103       KREQ(2)=0
 104       DO 103 I=1,KELEM
 105       RQV(I)=RVIND
 106       KRQ(I)=NVIND
 107  103  CONTINUE
 108 C
 109 C*          1.2 OPEN FILE CONTAINING BUFR DATA.
 110 C               -------------------------------
 111  120  CONTINUE
 112 C
 113       IRET=0
 114       CALL PBOPEN(IUNIT,CFIN(1:ILN),'R',IRET)
 115       write(6,*) 'opening ',cfin(1:iln)
 116       IF(IRET.EQ.-1) STOP 'OPEN FAILED'
 117       IF(IRET.EQ.-2) STOP 'INVALID FILE NAME'
 118       IF(IRET.EQ.-3) STOP 'INVALID OPEN MODE SPECIFIED'
 119
 120
 121 C
 122 C
 123 C     -----------------------------------------------------------------
 124 C*          2. SET REQUEST FOR EXPANSION.
 125 C              --------------------------
 126  200  CONTINUE
 127 C
 128       OPRT=.FALSE.
 129       OENC=.FALSE.
 130 c     WRITE(*,'(A,$)') ' DO YOU WANT TO PRINT( Y/N ) : '
 131          OPRT=.TRUE.
 132       ICODE=0
 133 c     WRITE(*,'(A,$)') ' CODE TABLES TO BE PRINTED ( Y/N ) : '
 134          ICODE=1
 135 c     WRITE(*,'(A,$)') ' RECORD NUMBER TO START FROM : '
 136       NR = 1
 137 C
 138  201  CONTINUE
 139 C
 140       OSEC3=.TRUE.
 141 C
 142 C*          2.1 SET REQUEST FOR PARTIAL EXPANSION.
 143 C               ----------------------------------
 144  210  CONTINUE
 145 C
 146 c     KERR=0
 147       CALL BUSRQ(KREQ,KRQL,KRQ,RQV,KERR)
 148 C
 149 C     SET VARIABLE TO PACK BIG VALUES AS MISSING VALUE INDICATOR
 150 C
 151       KPMISS=1
 152       KPRUS=0
 153       KOKEY=0
 154       CALL BUPRQ(KPMISS,KPRUS,KOKEY)
 155 C
 156 C     -----------------------------------------------------------------
 157 C*          3.  READ BUFR MESSAGE.
 158 C               ------------------
 159  300  CONTINUE
 160 C
 161       IERR=0
 162       KBUFL=0
 163 C
 164       IRET=0
 165       CALL PBBUFR(IUNIT,KBUFF,JBYTE,KBUFL,IRET)
 166       IF(IRET.EQ.-1) THEN
 167          PRINT*,'NUMBER OF SUBSETS     ',IOBS
 168          PRINT*,'NUMBER OF MESSAGES    ',N
 169          CALL PBCLOSE(IUNIT,IRET)
 170          GO TO 101
 171       END IF
 172       IF(IRET.EQ.-2) STOP 'FILE HANDLING PROBLEM'
 173       IF(IRET.EQ.-3) STOP 'ARRAY TOO SMALL FOR PRODUCT'
 174 C
 175       N=N+1
 176       PRINT*,'----------------------------------',N,' ',KBUFL
 177       KBUFL=KBUFL/NBYTPW+1
 178       IF(N.LT.NR) GO TO 300
 179 C
 180 C     -----------------------------------------------------------------
 181 C*          4. EXPAND BUFR MESSAGE.
 182 C              --------------------
 183  400  CONTINUE
 184 C
 185       CALL BUS012(KBUFL,KBUFF,KSUP,KSEC0,KSEC1,KSEC2,KERR)
 186       IF(KERR.NE.0) THEN
 187          PRINT*,'ERROR IN BUS012: ',KERR
 188          PRINT*,' BUFR MESSAGE NUMBER ',N,' CORRUPTED.'
 189          KERR=0
 190          GO TO 300
 191       END IF
 192 C
 193       KEL=KVALS/KSUP(6)
 194       IF(KEL.GT.JELEM) KEL=JELEM
 195 C
 196          CALL BUFREX(KBUFL,KBUFF,KSUP,KSEC0 ,KSEC1,KSEC2 ,KSEC3 ,KSEC4,
 197      1            KEL,CNAMES,CUNITS,KVALS,VALUES,CVALS,IERR)
 198 C
 199       IF(IERR.NE.0) THEN
 200          IF(IERR.EQ.39) GO TO 300
 201          CALL EXIT(2)
 202       END IF
 203 C
 204 C
 205       IOBS=IOBS+KSEC3(3)
 206 C
 207       NPACK=NPACK+1
 208 C
 209       CALL BUSEL(KTDLEN,KTDLST,KTDEXL,KTDEXP,KERR)
 210       IF(KERR.NE.0) CALL EXIT(2)
 211
 212       iyr = ksec1(9)
 213       imo = ksec1(10)
 214       idy = ksec1(11)
 215       ihr = ksec1(12)
 216       imn = ksec1(13)
 217 c     write(6,*) 'date = ',iyr,imo,idy,ihr,imn
 218 C
 219 C
 220 C*          4.1 PRINT CONTENT OF EXPANDED DATA.
 221 C               -------------------------------
 222  410  CONTINUE
 223 C
 224       IF(.NOT.OPRT) GO TO 500
 225       IF(.NOT.OSEC3) GO TO 450
 226 C
 227 C*          4.2 PRINT SECTION ZERO OF BUFR MESSAGE.
 228 C               -----------------------------------
 229  420  CONTINUE
 230 C
 231
 232       CALL BUPRS0(KSEC0)
 233 C
 234 C*          4.3 PRINT SECTION ONE OF BUFR MESSAGE.
 235 C               -----------------------------------
 236  430  CONTINUE
 237 C
 238       CALL BUPRS1(KSEC1)
 239 C
 240 C
 241 C*          4.4 PRINT SECTION TWO OF BUFR MESSAGE.
 242 C               -----------------------------------
 243  440  CONTINUE
 244 C
 245 C*          4.5 PRINT SECTION 3 OF BUFR MESSAGE.
 246 C               -----------------------------------
 247  450  CONTINUE
 248 C
 249 C               FIRST GET DATA DESCRIPTORS
 250 C
 251       CALL BUSEL(KTDLEN,KTDLST,KTDEXL,KTDEXP,KERR)
 252 C     IF(KERR.NE.0) CALL EXIT(2)
 253 C
 254 C               PRINT  CONTENT
 255 C
 256       IF(OSEC3) THEN
 257          CALL BUPRS3(KSEC3,KTDLEN,KTDLST,KTDEXL,KTDEXP,KEL,CNAMES)
 258       END IF
 259
 260 c     write(6,*) 'ktdexl = ',ktdexl
 261 c     do jj = 1,9
 262 c       write(6,*) jj,' ksup = ',ksup(jj)
 263 c     enddo
 264 c     write(6,*) 'kel = ', KVALS/KSUP(6)
 265 c     write(6,*) 'jelem = ',jelem
 266 c     write(6,*) 'KVALS = ',kvals
 267
 268       do ij = 1, ksup(6)    ! loop over reports
 269         do kj = 1,ktdexl
 270           jj = kj + (ij-1)*kel
 271 c         write(6,*) 'sta = ',ij,' val = ',values(jj),' cnames = ',
 272 c    &        cnames(kj)
 273 c         write(6,*) 'sta = ',ij,' ktdexp = ',ktdexp(kj),' cunits = ',
 274 c    &        cunits(kj)
 275 c         write(6,789) ij,ktdexp(kj),values(jj),cunits(kj),cnames(kj)
 276         enddo
 277       enddo
 278   789 format (i6,i8,f20.6,1x,a20,a80)
 279       call wlittler (values,kvals,ktdexp,jelem,ktdexl,kel,ksup(6),
 280      &   ict,3,ifm,cvals)
 281 C
 282 C*         4.6 PRINT SECTION 4 (DATA).
 283 C              -----------------------
 284  460  CONTINUE
 285 C
 286       ist = 1
 287       iend = ktdexl
 288 C
 289 C              PRINT DATA
 290 C
 291       ICODE=0
 292
 293 c        CALL BUPRT(ICODE,IST,IEND,KEL,CNAMES,CUNITS,CVALS,
 294 c    1              KVALS,VALUES,KSUP,KSEC1,IERR)
 295 C
 296 C     -----------------------------------------------------------------
 297 C*          5. COLLECT DATA FOR REPACKING.
 298 C              ---------------------------
 299  500  CONTINUE
 300  600  CONTINUE
 301 C
 302       NPACK=0
 303 C
 304       GO TO 300
 305 C     -----------------------------------------------------------------
 306 C
 307  101  CONTINUE
 308 C
 309       go to 87
 310   987 write(6,*) 'error reading flist'
 311       stop
 312   988 write(6,*) 'end of file in flist'
 313       END
 314 c-------------------------------------------------------
 315       subroutine wlittler (v,kvals,kt,jelem,ktdexl,kel,nsta,
 316      &    ict,iunit,ifm,cvals)
 317       parameter (xmis = -888888., mkx=300)
 318       real*8 v(kvals)
 319       real lat, lon, elev
 320       real p(mkx), z(mkx), t(mkx), td(mkx), spd(mkx), dir(mkx)
 321       integer kt(jelem)
 322       character ctime*12, ch2*2, ch3*3, ch4*4
 323       character*80 cvals(kel)
 324       character*40 id, name, platform, source
 325       logical is_sound
 326       id     = '                                        '
 327       name   = '                                        '
 328       source = '                                        '
 329       kx = 1
 330       if (ifm .eq. 12) then
 331         platform = 'FM-12 SYNOP                             '
 332         is_sound = .false.
 333       else if (ifm .eq. 13) then
 334         platform = 'FM-13 SHIP                              '
 335         is_sound = .false.
 336       else if (ifm .eq. 15) then
 337         platform = 'FM-15 METAR                             '
 338         is_sound = .false.
 339       else if (ifm .eq. 32) then
 340         platform = 'FM-32 PILOT                             '
 341         is_sound = .true.
 342       else if (ifm .eq. 35) then
 343         platform = 'FM-35 TEMP                              '
 344         is_sound = .true.
 345       else if (ifm .eq. 88) then
 346         platform = 'FM-88 SATOB                             '
 347         is_sound = .true.
 348       else if (ifm .eq. 96) then
 349         platform = 'FM-96 AIREP                             '
 350         is_sound = .true.
 351       endif
 352       do 100 i = 1, nsta
 353         ict = ict + 1
 354         iw = 0
 355         call initzero(p, z, t, td, spd, dir,
 356      &      slp, psfc, elev ,lat, lon, xmis, kx)
 357         kx = 0
 358         do k = 1, ktdexl
 359           j = k + (i-1)*kel
 360 c         write(6,*) i,kt(k),v(j)
 361           if (kt(k) .eq. 1001) then
 362             write(ch2,'(i2.2)') int(v(j))
 363             id(1:2) = ch2
 364           endif
 365           if (kt(k) .eq. 1002) then
 366             write(ch3,'(i3.3)') int(v(j))
 367             id(3:5) = ch3
 368           endif
 369           if (kt(k) .eq. 1006) then    ! aircraft
 370             m = int(v(j)) / 1000
 371             ml = int(v(j)) - m*1000
 372             id(1:ml) = cvals(m)(1:ml)
 373             elev = xmis
 374           endif
 375           if (kt(k) .eq. 1006) then    ! satellite  not implemented by kma
 376             write(ch2,'(i2.2)') int(v(j))
 377             id = 'SAT'//ch2
 378           endif
 379           if (kt(k) .eq. 1011) then
 380             m = int(v(j)) / 1000
 381             ml = int(v(j)) - m*1000
 382             id(1:ml) = cvals(m)(1:ml)
 383             elev = 0.    ! watch out for Great Lakes
 384           endif
 385           if (kt(k) .eq. 4001) then
 386             write(ch4,'(i4.4)') int(v(j))
 387             ctime(1:4) = ch4
 388           endif
 389           if (kt(k) .eq. 4002) then
 390             write(ch2,'(i2.2)') int(v(j))
 391             ctime(5:6) = ch2
 392           endif
 393           if (kt(k) .eq. 4003) then
 394             write(ch2,'(i2.2)') int(v(j))
 395             ctime(7:8) = ch2
 396           endif
 397           if (kt(k) .eq. 4004) then
 398             write(ch2,'(i2.2)') int(v(j))
 399             ctime(9:10) = ch2
 400           endif
 401           if (kt(k) .eq. 4005) then
 402             write(ch2,'(i2.2)') int(v(j))
 403             ctime(11:12) = ch2
 404           endif
 405           if (kt(k) .eq. 5001) lat = verifi(v(j),-90.,90.,xmis)
 406           if (kt(k) .eq. 6001) lon = verifi(v(j),-180.,180.,xmis)
 407           if (kt(k) .eq. 5002) lat = verifi(v(j),-90.,90.,xmis)
 408           if (kt(k) .eq. 6002) lon = verifi(v(j),-180.,180.,xmis)
 409           if (kt(k) .eq. 7001) elev = verifi(v(j),-200.,9100.,xmis)
 410           if (kt(k) .eq. 7030) elev = verifi(v(j),-200.,9100.,xmis)
 411           if (kt(k) .eq. 7004) then
 412             kx = kx + 1
 413             p(kx) = verifi(v(j),100.,110000.,xmis)
 414           endif
 415           if (kt(k) .eq. 10004) then
 416             psfc = verifi(v(j),100.,1100.,xmis)
 417             if (psfc .gt. 0.) psfc = psfc * 100.
 418             p(1) = psfc
 419             z(1) = elev   ! for sfc stations,
 420             kx = kx + 1
 421           endif
 422           if (kt(k) .eq. 10051) then
 423             slp = verifi(v(j),100.,1100.,xmis)
 424             if (slp .gt. 0.) slp = slp * 100.
 425           endif
 426           if (kt(k) .eq.  7002) then
 427             kx = kx + 1
 428             z(kx) = verifi(v(j),-200.,55000.,xmis)  ! m
 429           endif
 430           if (kt(k) .eq. 10009) z(kx) = verifi(v(j),-200.,55000.,xmis)  ! gpm
 431           if (kt(k) .eq. 11011) dir(1) = verifi(v(j),0.,360.,xmis)  ! 10m
 432           if (kt(k) .eq. 11012) spd(1) = verifi(v(j),0.,200.,xmis)  ! 10m
 433           if (kt(k) .eq. 11001) dir(kx) = verifi(v(j),0.,360.,xmis)
 434           if (kt(k) .eq. 11002) spd(kx) = verifi(v(j),0.,200.,xmis)
 435           if (kt(k) .eq. 12004) t(1) = verifi(v(j),0.,360.,xmis)    ! 2m
 436           if (kt(k) .eq. 12006) td(1) = verifi (v(j),0.,360.,xmis)  ! 2m
 437           if (kt(k) .eq. 12001) t(kx) = verifi(v(j),0.,360.,xmis)   ! aircraft temp
 438           if (kt(k) .eq. 12101) t(kx) = verifi(v(j),0.,360.,xmis)
 439           if (kt(k) .eq. 12003) td(kx) = verifi (v(j),0.,360.,xmis) ! aircraft td
 440           if (kt(k) .eq. 12103) td(kx) = verifi (v(j),0.,360.,xmis)
 441           enddo
 442           if (lat .eq. xmis .or. lon .eq. xmis) then
 443             ict = ict - 1
 444             go to 100
 445           endif
 446           if (is_sound .and. z(1) .eq. elev) then
 447             psfc = p(1)
 448           endif
 449           if (is_sound .and. (ifm .ne. 96 .and. ifm .ne. 88 )) then
 450             call cmprsnd (p, z, t, td, spd, dir, kx, xmis)
 451           endif
 452             call write_obs (p, z,
 453      &      t, td, spd, dir,
 454      &      slp, psfc, elev ,lat, lon,
 455      &      ctime , kx ,
 456      &      id, name, platform, source,
 457      &      is_sound, bogus, ict, iunit)
 458   100 continue
 459       end
 460 c-------------------------------------
 461       SUBROUTINE write_obs ( p, z, t, td, spd, dir,
 462      &                      slp, psfc, ter, xlat, xlon, cdate, kx,
 463      & string1, string2, string3, string4, is_sound, bogus,
 464      & iseq_num, iunit)
 465
 466       dimension p(kx), z(kx),t(kx),td(kx),spd(kx),dir(kx)
 467
 468       character *20 date_char
 469       character *40 string1, string2 , string3 , string4
 470       CHARACTER *84  rpt_format
 471       CHARACTER *22  meas_format
 472       CHARACTER *14  end_format
 473       character cdate*12, cmin*4
 474       logical is_sound,bogus
 475
 476       rpt_format =  ' ( 2f20.5 , 2a40 , '
 477      *             // ' 2a40 , 1f20.5 , 5i10 , 3L10 , '
 478      *             // ' 2i10 , a20 ,  13( f13.5 , i7 ) ) '
 479       meas_format =  ' ( 10( f13.5 , i7 ) ) '
 480       end_format = ' ( 3 ( i7 ) ) '
 481
 482       date_char(17:20)='0000'
 483       date_char(1:6)='      '
 484       date_char(7:18) = cdate
 485
 486       WRITE ( UNIT = iunit, ERR = 19, FMT = rpt_format )
 487      *        xlat,xlon, string1 , string2 ,
 488      *        string3, string4, ter, kx*6, 0, 0, iseq_num, 0,
 489      *        is_sound,bogus,.false.,
 490      *         -888888, -888888, date_char,
 491      *         slp,0,-888888.,0, -888888.,0, -888888.,0, psfc,0,
 492      *               -888888.,0,
 493      *               -888888.,0, -888888.,0, -888888.,0, -888888.,0,
 494      *               -888888.,0,
 495      *               -888888.,0, -888888.,0
 496
 497       do 100 k = 1 , kx
 498          WRITE ( UNIT = iunit , ERR = 19 , FMT = meas_format )
 499      *          p(k), 0, z(k),0, t(k),0, td(k),0,
 500      *          spd(k),0, dir(k),0,
 501      *          -888888.,0, -888888.,0,-888888.,0, -888888.,0
 502 100   continue
 503       WRITE ( UNIT = iunit , ERR = 19 , FMT = meas_format )
 504      * -777777.,0, -777777.,0,float(kx),0,
 505      * -888888.,0, -888888.,0, -888888.,0,
 506      * -888888.,0, -888888.,0, -888888.,0,
 507      * -888888.,0
 508       WRITE ( UNIT = iunit , ERR = 19 , FMT = end_format )  kx, 0, 0
 509
 510       return
 511 19    continue
 512       print *,'troubles writing little_r observation'
 513       stop 19
 514       END
 515 c------------------------------------------
 516       subroutine subdat (ccyymmddhh, dh, idate)
 517       INTEGER ccyymmddhh,ccyy,mmddhh,mm,dd,hh,dh
 518       character*10 cmin
 519
 520       ccyy   = ccyymmddhh / 1000000
 521       mmddhh = MOD ( ccyymmddhh , 1000000 )
 522       mm     = mmddhh / 10000
 523       dd     = MOD ( mmddhh , 10000 ) / 100
 524       hh     = MOD ( mmddhh , 100 )
 525
 526       hh = hh + dh
 527    10 IF ( hh .LT. 0 ) THEN
 528          hh = hh + 24
 529          CALL change_date ( ccyy, mm, dd, -1 )
 530       ELSEIF ( hh .GT. 23 ) THEN
 531          hh = hh - 24
 532          CALL change_date ( ccyy, mm, dd, 1 )
 533       ELSE
 534          WRITE (cmin,'(I4.4,3I2.2)') ccyy,mm,dd,hh
 535          read (cmin,'(i10)') idate
 536          return
 537       ENDIF
 538       GOTO 10
 539       END
 540 C
 541       SUBROUTINE change_date ( ccyy, mm, dd, delta )
 542       INTEGER ccyy, mm, dd, delta
 543       INTEGER mmday(12)
 544       DATA    mmday/31,28,31,30,31,30,31,31,30,31,30,31/
 545
 546       mmday(2) = 28
 547       IF ( MOD(ccyy,4) .EQ. 0 ) THEN
 548          IF     ( MOD(ccyy,400) .EQ. 0 ) THEN
 549             mmday(2) = 29
 550          ELSEIF ( MOD(ccyy,100) .NE. 0 ) THEN
 551             mmday(2) = 29
 552          ENDIF
 553       ENDIF
 554
 555       dd = dd + delta
 556       IF ( dd .EQ. 0 ) THEN
 557          mm = mm - 1
 558          IF ( mm .EQ. 0 ) THEN
 559             mm = 12
 560             ccyy = ccyy - 1
 561          ENDIF
 562          dd = mmday(mm)
 563       ELSEIF ( dd .GT. mmday(mm) ) THEN
 564          dd = 1
 565          mm = mm + 1
 566          IF ( mm .GT. 12 ) THEN
 567             mm = 1
 568             ccyy = ccyy + 1
 569          ENDIF
 570       ENDIF
 571       RETURN
 572       END
 573 c-----------------------------------------------------------------
 574       subroutine jdate (iyr, jday, ihr, jtime)
 575       character tmp*8
 576       integer mon(12)
 577       data mon/31,28,31,30,31,30,31,31,30,31,30,31/
 578       write(tmp,'(i4)') iyr
 579       read(tmp,'(2x,i2)') jyr
 580       if (mod(jyr,4) .eq. 0) mon(2) = 29
 581       m = 0
 582       do i = 1, 12
 583         m = m + mon(i)
 584         if (jday .le. m ) go to 10
 585       enddo
 586    10 continue
 587       if ( i .gt. 1 ) then
 588         idy = jday - ( m - mon(i))
 589       else
 590         idy = jday
 591       endif
 592       write(tmp,'(4(i2.2))') jyr, i, idy, ihr
 593       read(tmp,'(i8)') jtime
 594       end
 595 c-----------------------------------------------------------------
 596       real function verifi (x,xmin,xmax,xmis)
 597       real*8 x
 598       if (x .lt. xmin .or. x .gt. xmax) then
 599         verifi = xmis
 600       else
 601         verifi = x
 602       endif
 603       end
 604 c-----------------------------------------------------------------
 605       subroutine initzero(p, z, t, td, spd, dir,
 606      &      slp, psfc, elev ,lat, lon, xmis, kx)
 607       real p(kx), z(kx), t(kx), td(kx), spd(kx), dir(kx)
 608       real slp, psfc, elev, lat, lon, xmis
 609       do i = 1, kx
 610         p(i) = xmis
 611         z(i) = xmis
 612         t(i) = xmis
 613         td(i) = xmis
 614         spd(i) = xmis
 615         dir(i) = xmis
 616       enddo
 617       slp = xmis
 618       psfc = xmis
 619       elev = xmis
 620       lat = xmis
 621       lon = xmis
 622       return
 623       end
 624 c-----------------------------------------------------------------
 625       subroutine cmprsnd (p, z, t, td, spd, dir, kx, xmis)
 626 c kma files have levels with all missing data, so delete them
 627       real p(kx), z(kx), t(kx), td(kx), spd(kx), dir(kx)
 628       il = 1
 629       ih = kx
 630    10 continue
 631 c     write(6,*) 'begin loop, il = ',il,' ih = ',ih
 632       do k = il, ih
 633         if (p(k) .eq. xmis .and. z(k) .eq. xmis .and. t(k) .eq. xmis
 634      &      .and. td(k) .eq. xmis .and. spd(k) .eq. xmis .and.
 635      &      dir(k) .eq. xmis) go to 20
 636       enddo
 637       kx = ih
 638       return
 639    20 continue
 640       do j = k, ih-1
 641         p(j) = p(j+1)
 642         z(j) = z(j+1)
 643         t(j) = t(j+1)
 644         td(j) = td(j+1)
 645         spd(j) = spd(j+1)
 646         dir(j) = dir(j+1)
 647       enddo
 648       il = min0(k,kx)
 649       ih = max0(ih - 1,1)
 650       if ( il .eq. 1 .and. ih .eq. 1 ) then
 651         kx = 1    ! if all levels are missing, return just one.
 652         return
 653       else
 654         goto 10
 655       endif
 656       end