	subroutine readd(fname,day,data,edata,ndata,cmnt,iret)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
ce
	real*4 data1(19),edata1(19)
	integer   nchk(nlim)
ce	character cmnt*20,fname*25
	character cmnt*20,fname*60
	character dummy*80
	logical exflag
*------- clear data ------------
 2001	do 2000 i21=1,nlim
 2101	  do 2100 i22=1,19
	    data(i21,i22)=0.0
	    edata(i21,i22)=0.0
 2100	  continue
 2000	continue
	iret=0
*------- inquire file ----------
	print *,"Input file name=",fname
	inquire(file=fname,exist=exflag)
	if(exflag)then
*------- open file -------------
	open(7,file=fname,status='old')
*------- read comment ----------
	read(7,'(a20)') cmnt
	print *,cmnt
	ndata=0
 1000	continue
*------- read data -------------
	read(7,'(a80)',end=999) dummy
	read(dummy,'(i4)') no
	if(no.eq.0) goto 1000

	read(dummy,'(4x,i4,4f8.4,2f8.3)') ich,d,e1,e2,ed,da,el

	if(ich.eq.108) then
	  mch=17
	elseif(ich.eq.916) then
	  mch=18
	elseif(ich.eq.116) then
	  mch=19
	else
	  mch=int(float(ich)/101)
	endif
	if(ndata.le.0) goto 1200
*------- search previous day ---
 1101	do 1100 i11=1,ndata
	  if(day(i11).eq.da)then
	    data(i11,mch)=d
	    edata(i11,mch)=ed
	    goto 1000
	  endif
 1100	continue
*------- new data --------------
 1200	continue
	ndata=ndata+1
	data(ndata,mch)=d
	edata(ndata,mch)=ed
	day(ndata)=da

*	day(ndata)=da+46795.5
*	    ---- to julian day

*------- check ndata -----------
	if(ndata.ge.nlim) then
	  write(*,'(1h,'' Too many data '')')
	  goto 999
	endif
*------- next data -------------
	goto 1000
  999 	continue
	close(7)
*------- for all channel -------

	if(ndata.le.0) goto 990
 1301   do 1300 i13=1,ndata
	  if(data(i13,19).eq.0.0) then
	    if(data(i13,17).ne.0.0.and.data(i13,18).ne.0.0) then
	      data(i13,19)=data(i13,17)+data(i13,18)
	      edata(i13,19)=sqrt(edata(i13,17)**2
     &	                              +edata(i13,18)**2)
	    else
	      data(i13,19)=0.0
	      edata(i13,19)=0.0
	    endif
	  endif
 1300	continue
ce------ Check bad data --------
 1401   do 1400 i14=1,ndata
           nchk(i14)=1
 1411      do 1410 i141=1,18
              if(data(i14,i141).le.-0.3) then
                 nchk(i14)=0
                 goto 1400
              elseif(edata(i14,i141).ge.0.5) then
                 nchk(i14)=0
                 goto 1400
              elseif(edata(i14,i141).le.2e-3) then
                 nchk(i14)=0
                 goto 1400
              elseif(data(i14,i141).le.0) then
                 zero = data(i14,i141) + edata(i14,i141)*3
                 if(zero.lt.0) then
                    nchk(i14)=0
                    goto 1400
                 endif
              endif
 1410      continue
 1400   continue


*------- Check consistency -----
 1501	do 1500 i15=1,ndata
	  csum1=0.0
	  esum1=0.0
	  csum2=0.0
	  esum2=0.0
	  rsum1=0.0
	  ersm1=0.0
	  rsum2=0.0
	  ersm2=0.0
ce	  nchk(i15)=1
	  if(edata(i15,17).gt.0.0) then
 1511	    do 1510 i151=1,8
	      csum1=csum1+data(i15,i151)
	      esum1=sqrt(esum1**2+edata(i15,i151)**2)
 1510	    continue
	    if(data(i15,17).eq.0.0) then
	      rsum1=csum1+1.0
	      esum1=0.0
	    else
	      rsum1=csum1/data(i15,17)
	      ersm1=sqrt((csum1*edata(i15,17))**2
     &             +(esum1*data(i15,17))**2)/(data(i15,17)**2)
	    endif
	    if(data(i15,17).lt.0.0) rsum1=abs(rsum1)
	    if(rsum1+ersm1.lt.0.45.or.rsum1-ersm1.gt.1.20)
     &        nchk(i15)=0
*	    if(rsum1.lt.0.1.or.rsum1.gt.10.0) nchk(i15)=0
	  endif
	  if(edata(i15,18).gt.0.0) then
 1521       do 1520 i152=9,16
              csum2=csum2+data(i15,i152)
              esum2=sqrt(esum2**2+edata(i15,i152)**2)
 1520       continue
	    if(data(i15,18).eq.0.0) then
	      rsum2=csum2+1.0
	      esum2=0.0
	    else
              rsum2=csum2/data(i15,18)
              ersm2=sqrt((csum2*edata(i15,18))**2
     &             +(esum2*data(i15,18))**2)/(data(i15,18)**2)
	    endif
	    if(data(i15,18).lt.0.0) rsum2=abs(rsum2)
            if(rsum2+ersm2.lt.0.45.or.rsum2-ersm2.gt.1.20)
     &         nchk(i15)=0
*	    if(rsum2.lt.0.1.or.rsum2.gt.10.0) nchk(i15)=0
	  endif
 1500	continue



*------- Sort data -------------

	mdata=0
 1531	do 1530 i153=1,ndata
	  if(nchk(i153).eq.0) goto 1530
	  mdata=mdata+1
	  day(mdata)=day(i153)
 1541	  do 1540 i154=1,19
	    data(mdata,i154) =data(i153,i154)
	    edata(mdata,i154)=edata(i153,i154)
 1540	  continue
 1530	continue
	ndata=mdata

ce------ gain crrection --------

	call setengch

 3001   do 3000 i3=1,ndata
 3101      do 3100 i31=1,19
              data1(i31)=data(i3,i31)
              edata1(i31)=edata(i3,i31)
 3100      continue
           call crceng(day(i3),data1,edata1)
 3201      do 3200 i32=1,19
              data(i3,i32)=data1(i32)
              edata(i3,i32)=edata1(i32)
 3200      continue
 3000   continue



*------- return ----------------
	else
	  write(*,*) 'File cannot be found! Input once again!'
	  write(*,*)
	  iret=9
	endif
  990 	continue
	return
	end
c=====================================================================
      subroutine crceng(day,data,edata)
c---------------------------------------------------------------------
      implicit NONE
c  input
      real*4 day,data(19),edata(19)
c  var
ce--- real*4 engch1(16)/1.000,1.709,2.255,2.861,
*    &                  3.466,4.084,4.737,5.329,
*    &                  5.967,7.124,8.320,9.505,
*    &                  10.698,12.979,15.393,17.898/
*     real*4 engch2(16)/1.709,2.255,2.861,3.466,
*    &                  4.084,4.737,5.329,5.967,
*    &                  7.124,8.320,9.505,10.698,
ce---&                  12.979,15.393,17.898,20.270/

*      real*4 engch1(16)/1.000,1.640,2.152,2.717,
*     &                  3.281,3.857,4.469,5.024,
*     &                  5.634,6.743,7.899,9.055,
*     &                  10.223,12.468,14.860,17.369/
*      real*4 engch2(16)/1.640,2.152,2.717,3.281,
*     &                  3.857,4.469,5.024,5.634,
*     &                  6.743,7.899,9.055,10.223,
*     &                  12.468,14.860,17.369,19.827/

	real*4 engch1(16), engch2(16)
      real*4 data2(19),edata2(19),egcrc1(16),egcrc2(16)
      real*4 gain,engch0,factor,factor2
      integer i1,i11,i2,i21,i3,ich
c  common
      common /engrng/ engch1,engch2
c  function
      real*4 calgin
      real*4 effic
c  begin

 1001 do 1000 i1=1,19
         data2(i1)=0
         edata2(i1)=0
 1000 continue

c  calculate gain
      gain = calgin(day)

 1101 do 1100 i11=1,16
         egcrc1(i11) = engch1(i11) / gain
         egcrc2(i11) = engch2(i11) / gain
 1100 continue

c  correction of each enegy ch
 2001 do 2000 i2=1,16

         if ( i2 .eq. 1 ) then

            factor = effic(engch1(1),engch1(2))
     &             / effic(egcrc1(1),egcrc1(2))
            data2(i2) = factor * data(1)
            edata2(i2) = factor * edata(1)
            goto 2000
         elseif ( i2 .eq. 2 ) then
            factor = effic(engch1(2),egcrc2(1))
     &             / effic(egcrc1(1),egcrc2(1))
            factor2 = effic(egcrc1(2),engch2(2))
     &              / effic(egcrc1(2),egcrc2(2))
            data2(i2) = factor * data(1) + factor2 * data(2)
            edata2(i2) = sqrt( ( factor * edata(1) )**2
     &                       + ( factor2 * edata(2) )**2 )
            goto 2000
         endif

         engch0 = engch1(i2)
         ich = 0
 2101    do 2100 i21=1,16
            if ( engch1(i2) .ge. egcrc1(i21)
     &           .and. engch1(i2) .lt. egcrc2(i21) ) then
               ich = i21
               goto 2200
            endif
 2100    continue

 2200    continue
         if ( engch2(i2) .le. egcrc2(ich) ) then
            factor = ( engch2(i2) - engch0 ) 
     &               / ( egcrc2(ich) - egcrc1(ich) )
            data2(i2) = data2(i2) + factor * data(ich)
            edata2(i2) = sqrt( edata2(i2)**2
     &                      + ( factor*edata(ich) )**2 )
         else
            factor = ( egcrc2(ich) - engch0 )
     &               / ( egcrc2(ich) - egcrc1(ich) )
            data2(i2) = data2(i2) + factor * data(ich)
            edata2(i2) = sqrt( edata2(i2)**2
     &                      + ( factor*edata(ich) )**2 )
            engch0 = egcrc2(ich)
            ich = ich + 1
            goto 2200
         endif

 2000 continue

c  correction of enegy ch 1-8, 9-16, 1-16
      factor = ( engch2(8)-engch1(1) ) / ( egcrc2(8)-egcrc1(1) )
      data2(17) = factor * data(17)
      edata2(17) = factor * edata(17)

      factor = ( egcrc2(8)-engch1(9) ) / ( egcrc2(8)-egcrc1(1) )
      factor2 = ( engch2(16)-egcrc1(9) ) / ( egcrc2(16)-egcrc1(9) )
      data2(18) = factor * data(17) + factor2 * data(18)
      edata2(18) = sqrt( ( factor*edata(17) )**2
     &                + ( factor2*edata(18) )**2 )

      factor = ( engch2(16)-engch1(1) ) / ( egcrc2(16)-egcrc1(1) )
      data2(19) = factor * data(19)
      edata2(19) = factor * edata(19)

 3001 do 3000 i3=1,19
         data(i3) = data2(i3)
         edata(i3) = edata2(i3)
 3000 continue


      return
      end
c=====================================================================
      function calgin(day)
c---------------------------------------------------------------------
      implicit NONE
c  input
      real*4 day
c  output
      real*4 calgin
c  const
      real*4 P1, P2, P3, P4
      parameter ( P1=7.849, P2=-3.1320e-4,
     &            P3=-2.1046e-3, P4=-2.1426e-4 )
c  var
      real*4 peak, offset
c  begin
      if ( day .lt. 1300.0 ) then
         peak = P1 + P2*day
      elseif ( day .lt. 1400.0 ) then
         peak = P1 + P2*1300.0 + P3*(day-1300.0)
      else
         peak = P1 + P2*1300.0 + P3*100.0 + P4*(day-1400.0)
      endif
      offset = P1 / ( 5.896 - 1.000 )
      calgin = ( peak + offset ) / ( P1 + offset )
      return
      end
c=====================================================================
      function effic(e1,e2)
c---------------------------------------------------------------------
      implicit NONE
c  input
      real*4 e1,e2
c  output
      real*4 effic
c  function
      real*4 eff,e,eff1,eff5,p1,p2
      eff(e) = 1.0 - exp( -p1 * e**p2 )
      eff1 = 1.0 - 0.001
      eff5 = 1.0 - 0.9
      p1 = -log(eff1)
      p2 = log( -log(eff5) / p1 ) / log(5.0)
c  begin

      if ( e2 .lt. e1 ) then
         write(*,*) 'ERROR e2 < e1 in function effic.'
         stop
      endif
      effic = ( eff(e1) + eff(e2) ) * ( e2 - e1 ) / 2.0

      return
      end

***********************************************************
	subroutine setengch

	real*4 engch1(16)
	real*4 engch2(16)
      common /engrng/ engch1,engch2

        engch1(1)=1.000
        engch1(2)=1.640
        engch1(3)=2.152   
        engch1(4)=2.717
        engch1(5)=3.281
        engch1(6)=3.857
        engch1(7)=4.469
        engch1(8)=5.024
        engch1(9)=5.634
        engch1(10)=6.743
        engch1(11)=7.899
        engch1(12)=9.055
        engch1(13)=10.223
        engch1(14)=12.468
        engch1(15)=14.860
        engch1(16)=17.369

        engch2(1)=1.640
        engch2(2)=2.152
        engch2(3)=2.717
        engch2(4)=3.281
        engch2(5)=3.857
        engch2(6)=4.469
        engch2(7)=5.024
        engch2(8)=5.634
        engch2(9)=6.743
        engch2(10)=7.899
        engch2(11)=9.055
        engch2(12)=10.223
        engch2(13)=12.468
        engch2(14)=14.860
        engch2(15)=17.369
        engch2(16)=19.827

	return
	end
