      program calccaldat
      implicit none

c -----
c data file IO
      integer ID
      parameter (ID=21)

c data is valid after TMIN
      real*8  TMIN
      parameter (TMIN=25.0d0)

c maximum number of calibration table data
      integer NTMAX
      parameter (NTMAX=100)

c -----
      real*8  t
      real*8  WTs, WTa, CTs, CTa

      integer ndat
      real*8  sWTs, sWTa, sCTs, sCTa

      integer nt
      real*8  Tct(NTMAX), Tst(NTMAX)
      real*8  WTc, CTc

      integer i

c -----
c open data file
      open(ID,file="datafile")

c skip header
      read(ID,*)
      read(ID,*)
      read(ID,*)
      read(ID,*)
      read(ID,*)
      read(ID,*)
      read(ID,*)
      read(ID,*)

c read data and caulculate
      ndat = 0
      sWTs = 0.0d0
      sWTa = 0.0d0
      sCTs = 0.0d0
      sCTa = 0.0d0
 200  continue
      read(ID,*,end=210) t, WTs, WTa, CTs, CTa
      if (t.ge.TMIN) then
          ndat = ndat + 1
          sWTs = sWTs + WTs
          sWTa = sWTa + WTa
          sCTs = sCTs + CTs
          sCTa = sCTa + CTa
      end if
      go to 200
 210  continue
      close(ID)
      if (ndat.ne.0) then
          WTs = sWTs / (ndat*1.0d0)
          WTa = sWTa / (ndat*1.0d0)
          CTs = sCTs / (ndat*1.0d0)
          CTa = sCTa / (ndat*1.0d0)
        else
          write(*,*) 'no data'
          stop
      end if

c output
c      write(*,*) CTs, CTa
c      write(*,*) WTs, WTa



c open calibration table
      open(ID,file="caltable")
c read table
      read(ID,*)
      do 310 i = 1, NTMAX
        read(ID,*,end=320) Tct(i), Tst(i)
 310  continue
      write(*,*) 'data length is too long'
      stop
 320  continue
      nt = i - 1

c calculate correct temperature
      call calTs(nt,Tct,Tst,CTs,CTc)
      call calTs(nt,Tct,Tst,WTs,WTc)

c output
      write(*,*) CTc, CTa
      write(*,*) WTc, WTa


      stop
      end



      subroutine calTs(nt,Tct,Tst,Ts,Tc)
      implicit none

c -----
c maximum number of calibration table data
      integer NTMAX
      parameter (NTMAX=100)

c -----
      integer nt
      real*8  Tct(NTMAX), Tst(NTMAX)

      real*8  Ts

c -----
      real*8  Tc

c -----
      real*8  dT

      integer i

c -----
      if (nt.eq.0) then
          write(*,*) 'calibration table is empty'
          stop
        else if (nt.eq.1) then
          dT = Tct(1) - Tst(1)
          Tc = Ts + dT
        else
          do 110 i = 2, nt
            if (Ts.le.Tst(i)) then
                Tc = ( (Tct(i)-Tct(i-1))/(Tst(i)-Tst(i-1)) )
     &                  *(Ts-Tst(i-1)) + Tct(i-1)
                go to 120
            end if
 110      continue
          Tc = ( (Tct(nt)-Tct(nt-1))/(Tst(nt)-Tst(nt-1)) )
     &            *(Ts-Tst(nt-1)) + Tct(nt-1)
 120      continue
      end if

      return
      end
