c 観測報デコードデータファイルのうち
c 3110または3160データから気圧・温度・露点・風のデータなどを
c 各々出力するプログラム
c 2000/12/27 by A.Tsubota
c 2001/07/17 特異点出力追加
c コンパイル時 要-byteswapio

      integer*2 idata(4000)
      integer*2 idt1(14),idt2(10),idt3(20),idt4(4000),idt5(4000)
      integer*2 no(4),miss,ti1,tdi1,wti1,wsi1
      integer l,iunit,tunit,tdunit,wiunit,no1,no2,flag
      integer miunit,maunit,tilunit,tirunit,irr,irrflag,numt,numtd
      integer tunit2,tdunit2
      real ti,tdi,wti,wsi,min,max,tmin,tmax,tdmin,tdmax,lat,lon
      real pi(100),pimin,pscan,pmin,wt(100),ws(100),wto,p(100)
      real t(100),td(100),p3,t3,td3

      miss=-32768
      flag=0
      irrflag=0
      pimin=1100.
      numt=0
      numtd=0

      iunit=10
      open(iunit,file='inputdat',access='direct',recl=2)
      tunit=11
      open(tunit,file='tdatav.dat',status='new')
      tdunit=12
      open(tdunit,file='tddatav.dat',status='new')
      wiunit=13
      open(wiunit,file='wind.dat',status='new')
      miunit=15
      open(miunit,file='min.dat',status='new')
      maunit=16
      open(maunit,file='max.dat',status='new')
      tilunit=17
      open(tilunit,file='tileft.dat',status='new')
      tirunit=18
      open(tirunit,file='tiright.dat',status='new')
      read(*,*) l,irr

      nc=1
 10   call get_record(iunit,idata,ndata)
      if (ndata.gt.0) then
*     write(*,*) nc,ndata,(idata(i),i=1,7)
         call get_data(idata,ndata,idt1,ndt1,idt2,ndt2,idt3,ndt3,
     &                 idt4,ndt4,idt5,ndt5)
*     write(*,*) nc,idt1

         else
            go to 9000
      end if

      if(nc.eq.l.or.nc.eq.l+1.and.irrflag.eq.1) then
      if(irrflag.eq.1.and.idt1(8).ne.3120.and.idt1(8).ne.3170) then
        go to 9000
      end if
      if(irrflag.ne.1) then
        do i=1,4
          no1=256
 20       no2=idt3(i)-no1
          if(no2.ge.256) then
            no1=no1+256
            go to 20
          end if
          no1=no1/256
          no2=no2*256
          no(i)=no1+no2
        end do
        write(tilunit,'(4a2)') no(1),no(2),no(3),no(4)
        write(tilunit,'(i4)') idt2(6)
	write(tilunit,*) int(idt2(5)/100)
	write(tilunit,*) idt2(5)-int(idt2(5)/100)*100

        if(idt1(9).ge.0) then
           lat=idt1(9)/100.
           write(tirunit,'(f6.2)') lat
           write(tirunit,*) 'N'
           else
             lat=abs(idt1(9))/100.
             write(tirunit,'(f6.2)') lat
             write(tirunit,*) 'S'
        end if
        if(idt1(10).ge.0) then
           lon=idt1(10)/100.
           write(tirunit,'(f6.2)') lon
           write(tirunit,*) 'E'
           else
              lon=abs(idt1(10))/100.
              write(tirunit,'(f6.2)') lon
              write(tirunit,*) 'W'
        end if
        end if

        do i=1,idt4(2)
          pscan=idt4(2+(i-1)*10+3)/10.
	  if(pscan.lt.pimin) then
	    pimin=pscan
	  end if
	end do

        do i=1,idt4(2)
          pi(i)=idt4(2+(i-1)*10+3)/10.
          if(flag.ne.1) then
            ti1=idt4(2+(i-1)*10+6)
            tdi1=idt4(2+(i-1)*10+7)
            wti1=idt4(2+(i-1)*10+9)
            wsi1=idt4(2+(i-1)*10+10)
            if(ti1.ne.miss) then
              ti=ti1/10.-273.15
              write(tunit,*) ti,pi(i)
              numt=numt+1
            end if
            if(tdi1.ne.miss) then
              tdi=tdi1/10.-273.15
              write(tdunit,*) tdi,pi(i)
              numtd=numtd+1
            end if
            wt(i)=wti1
            ws(i)=wsi1/10.
            call minmax(tmin,tmax,ti)
            call minmax(tdmin,tdmax,tdi)
          end if
          if(pi(i).eq.pimin) then
             flag=1
          end if
        end do

        if(tmax.ge.tdmax) then
          max=int(tmax/10)*10+10
          else
            max=int(tdmax/10)*10+10
        endif
        if(tmin.le.tdmin) then
          min=int(tmin/10)*10-20
          else
            min=int(tdmin/10)*10-20
        endif
        write(miunit,*) int(min)
        write(maunit,*) int(max)

        flag=0

        do i=1,idt4(2)
          if(flag.ne.1) then
            if(wt(i).ne.miss.and.ws(i).ne.miss) then
              wto=90.-wt(i)
              write(wiunit,*) max-8,pi(i),wto,ws(i)
            end if
          end if
          if(pi(i).eq.pimin) then
             flag=1
          end if
        end do
      if(irr.ne.1.or.irrflag.eq.1) then
      go to 9000
      else
      irrflag=1
      flag=0
      pimin=1100.
      end if
      end if
      nc=nc+1
      go to 10

 9000 close(iunit)
      close(tunit)
      close(tdunit)
      close(wiunit)
      close(munit)
      close(tilunit)
      close(tirunit)

      tunit=11
      open(tunit,file='tdatav.dat')
      tdunit=12
      open(tdunit,file='tddatav.dat')
      tunit2=19
      open(tunit2,file='tdata.dat',status='new')
      tdunit2=20
      open(tdunit2,file='tddata.dat',status='new')

c 温度データのソート

      do i=1,numt
         read(tunit,*) t(i),p(i)
      end do
      do j=1,numt-1
         do i=j,numt
          if(p(i).gt.p(j)) then
             p3=p(i)
             p(i)=p(j)
             p(j)=p3
             t3=t(i)
             t(i)=t(j)
             t(j)=t3
          end if
         end do
      end do
      do i=1,numt
         write(tunit2,*) t(i),p(i)
      end do

c 露点データのソート

      do i=1,numtd
         read(tdunit,*) td(i),p(i)
      end do
      do j=1,numtd-1
         do i=j,numtd
          if(p(i).gt.p(j)) then
             p3=p(i)
             p(i)=p(j)
             p(j)=p3
             td3=td(i)
             td(i)=td(j)
             td(j)=td3
          end if
         end do
      end do
      do i=1,numtd
         write(tdunit2,*) td(i),p(i)
      end do

      close(tunit)
      close(tdunit)
      close(tunit2)
      close(tdunit2)

      stop
      end
c
      subroutine get_data(idata,ndata,idt1,ndt1,idt2,ndt2,idt3,ndt3,
     &                    idt4,ndt4,idt5,ndt5)
      integer*2 idata(*),idt1(*),idt2(*),idt3(*),idt4(*),idt5(*)

      ndt1=idata(2)
      ndt2=idata(3)
      ndt3=idata(4)
      ndt4=idata(5)
      ndt5=idata(6)

      nn=0
      do 10 n=1,ndt1
         nn=nn+1
         idt1(n)=idata(nn)
 10   continue
      do 20 n=1,ndt2
         nn=nn+1
         idt2(n)=idata(nn)
 20   continue
      do 30 n=1,ndt3
         nn=nn+1
         idt3(n)=idata(nn)
 30   continue
      do 40 n=1,ndt4
         nn=nn+1
         idt4(n)=idata(nn)
 40   continue
      do 50 n=1,ndt5
         nn=nn+1
         idt5(n)=idata(nn)
 50   continue
      return
      end
c
      subroutine get_record(iunit,idata,ndata)
      integer*2 idata(*)
      data nrec/1/

      call get_1_data(iunit,nrec,idata(1),ist)
      if (ist.eq.-1) then
         ndata=ist
         return
         else
            nrec=nrec+1
            ndata=1
      end if

      do 10 n=2,idata(1)
         call get_1_data(iunit,nrec,idata(n),ist)
         if (ist.eq.-1) then
            ndata=ist
            return
            else
               nrec=nrec+1
               ndata=n
         end if
 10   continue
      return
      end
c
      subroutine get_1_data(i_unit,n_rec,i_data,i_st)
      integer*2 i_data

      read(i_unit,rec=n_rec,end=9000,err=8000) i_data
      i_st=0
      return
 8000 write(*,*) 'Error!'
      i_st=-1
      return
 9000 i_st=-1
      return
      end

      subroutine minmax(min,max,in)
      real min,max,in

      if(min.eq.0) then
         min=in
      end if
      if(max.eq.0) then
         max=in
      end if
      if(in.lt.min) then
         min=in
      end if
      if(in.gt.max) then
         max=in
      end if

      return
      end
