c=====================================================================
      subroutine crceng(day,data,edata)
c---------------------------------------------------------------------
      implicit NONE
c  input
      real day,data(19),edata(19)
c  var
      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 data2(19),edata2(19),egcrc1(16),egcrc2(16)
      real gain,engch0,factor,factor2
      integer i1,i11,i2,i21,i3,ich
c  common
      common /engrng/ engch1,engch2
c  function
      real calgin
      real 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 day
c  output
      real calgin
c  const
      real P1, P2, P3, P4
      parameter ( P1=7.849, P2=-3.1320e-4,
     &            P3=-2.1046e-3, P4=-2.1426e-4 )
c  var
      real 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 e1,e2
c  output
      real effic
c  function
      real 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
