C NCLFORTSTART
      subroutine ES (t,h,w,r,buio,ffmo,
     1                       xisio,fwio,dmco,dco)
C        integer tlen
      dimension t(183),h(183)
      dimension w(183),r(183)
      dimension ffmo(183),xisio(183)
      dimension fwio(183),dmco(183)
      dimension buio(183),dco(183)
C        integer indate,yr,mon,day,hr
C        character*10 timez
C NCLEND
C        write(timez,'(i10)') indate
C        read(timez,'(i4,3i2)') yr,mon,day,hr


      call canada(t,h,w,r,buio,ffmo,xisio,
     1           fwio,dmco,dco)

      return
      end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine canada(t,h,w,r,buio,ffmo,xisio,
     1           fwio,dmco,dco)

      dimension lmon(12), el(12), fl(12)
      dimension t(183),h(183),w(183),r(183)
      dimension tin(183),hin(183),win(183),rin(183)
      dimension ffmc(183),dmc(183),dc(183)
     1          ,xisi(183),bui(183),fwi(183),dsr(183)
      dimension ffmo(183),dmco(183),dco(183),buio(183)
      dimension fwio(183),xisio(183)
c data statement
      data lmon / 31,28,21,30,31,30,31,31,30,31,30,31/
      data el /6.5,7.5,9.0,12.8,13.9,13.9,12.4,10.9,9.4,8.0,7.0,
     1         6.0/
      data fl /-1.6,-1.6,-1.6,.9,3.8,5.8,6.4,5.0,2.4,.4,-1.6,
     1         -1.6/
c
c&& enter in the inital values here
c 3 fuel moisture codes, initialize
       ffmo(1) =85.5
       dmco(1)=6.0
       dco(1)=15.0

c keep track of julian day
            iday=90
c start on april 1 and go to september 30 
        do im=4,9
        do jday=1,lmon(im)
            iday = iday +1
c move the data from the array to the variable for the calculation on a given day
            t = tin(iday)
            h = hin(iday)
            w = win(iday)
            r = rin(iday)
c
c  convert m/s to km/hr
c      w(i,j) = w(i,j)*3.6
c
c     FINE FUEL MOISTURE CODE
c
      wmo = 147.2*(101-ffmo(iday))/(59.5+ffmo(iday))
      if(r(iday)>0.5) then
        ra=r(iday)-0.5
        if (ra.lt.0.1) ra = 0.1
        if(wmo.gt.150) then
          wmo= wmo+ 42.5*ra*exp(-100.0/(251-wmo))*(1.0-exp(-6.93/ra))
     1      +0.0015*(wmo-150)*(wmo-150)*sqrt(ra)
        else
          wmo=wmo+42.5*ra*exp(-100.0/(251-wmo))*(1.0-exp(-6.93/ra))
        endif
      endif
      if (wmo.gt.250.) wmo = 250.
      ed=0.942 * h(iday)**0.679 +
     1          (11.0*exp((h(iday)-100.0)/10.0)) +
     1    0.18 * (21.1-t(iday)) * (1.0-1.0/exp(h(iday)*0.115))
      ew=0.618 * h(iday)**0.753 +
     1       (10.0*exp((h(iday)-100.0)/10.0)) +
     1    0.18 * (21.1-t(iday)) * (1.0-1.0/exp(h(iday)*0.115))
      if(wmo.lt.ed.and.wmo.lt.ew) then
        z = 0.424 * ( 1.0 - ((100.0-h(iday))/100.0)**1.7 ) +
     1     0.0694 * sqrt(w(iday)) * (1.0 -
     1  ((100.0-h(iday))/100.0)**8.0 )
        x = z * 0.581 * exp(0.0365*t(iday))
        wm = ew-(ew-wmo) / (10.0**x)
      elseif (wmo.gt.ed) then
        z = 0.424 * ( 1.0 - (h(iday)/100.)**1.7 ) + 0.0694 *
     1    sqrt(w(iday)) * ( 1 - (h(iday)/100)**8.0 )
        x = z * 0.581 * exp(0.0365*t(iday))
        wm = ed + ( wmo - ed) / (10.0**x)
      else
        wm = wmo
      endif
      if (wm.lt.0.) wm = 0.0
      ffmc(iday) = 59.5 * (250.0 - wm) / (147.2 + wm)
      if (ffmc(iday).gt.101.0) ffmc(iday) = 101.0
      if (ffmc(iday).lt.0.) ffmc(iday) = 0.
c
c     DUFF MOISTURE CODE
c
      t1=t(iday)
      if (t(iday).lt.-1.1) t1 = -1.1
      rk = 1.894 * (t1 + 1.1) * (100.0 - h(iday)) * el(im) * 0.0001
      if(r(iday).le.1.5) then
        pr = dmco(iday)
      else
        ra = r(iday)
        rw = 0.92 * ra - 1.27
        wmi = 20.0 + 280.0 / exp(0.023 * dmco(iday))
        if(dmco(iday).le.33) then
          b = 100.0 / (0.5 + 0.3 * dmco(iday))
        elseif (dmco(iday).le.65) then
          b = 14.0 - 1.3 * log(dmco(iday))
        else
          b = 6.2 * log(dmco(iday)) - 17.2
        endif
        wmr = wmi + 1000.0 * rw / (48.77 + b * rw)
        pr = 43.43 * (5.6348 - log(wmr-20))
      endif
      if (pr.lt.0.) pr = 0.
      dmc(iday) = pr + rk
      if (dmc(iday).lt.0.) dmc(iday) = 0.
c
c     DROUGHT CODE

      t2 = t(iday)
      if (t(iday).lt.-2.8) t2 = -2.81
      pe = (.36 * (t2 + 2.8) + fl(mon)) / 2.0
      if (pe.lt.0.) pe = 0.0
c the fix for winter negative DC change
      if(r(iday).le.2.801) then
        dr = dco(iday)
      else
       ra = r(iday)
       rw = 0.83 * ra - 1.27
       smi = 800 * exp(-dco(iday) / 400)
       qr = smi + 3.937 * rw
       dr = 400 * log(800 / qr)
       if (dr.lt.0.) dr = 0.
      endif
      v = 0.36 * (t2 + 2.8) + fl(im)
      if (v.lt.0.) v = 0.
      dc(iday) = dr + v * 0.5
      if (dc(iday).lt.0.) dc(iday)=0.
c
c     INITIAL SPREAD INDEX
c
      fm = 147.2 * (101.0 - ffmc(iday)) / (59.5 + ffmc(iday))
      sf = 19.115 * exp(-0.1386 * fm) * (1.0 + fm**5.31 / 4.93e07)
      xisi(iday) = sf * exp(0.05039 * w(iday))
c
c     BUILD UP INDEX
c
      if(dmc(iday).eq.0.and.dc(iday).eq.0.) then
        bui(iday)=0
      elseif (dmc(iday).le.0.4*dc(iday)) then
        bui(iday) = 0.8 * dc(iday) * dmc(iday)
     1   / (dmc(iday) + 0.4 * dc(iday))
      else
        bui(iday) = dmc(iday) - (1 - 0.8 * dc(iday)
     1   / (dmc(iday) + 0.4 *
     1    dc(iday))) * (0.92 + (0.0114 * dmc(iday))** 1.7)
      endif
      if (bui(iday).lt.0.) bui(iday) = 0.0
c
c     FIRE WEATHER INDEX AND DSR
c
      if(bui(iday).gt.80.) then
        bb = 0.1 * xisi(iday) *
     1    (1000.0 / (25.0 + 108.64 / exp(0.023*bui(iday))))
      else
        bb = 0.1 * xisi(iday) * (0.626 * bui(iday)**0.809 + 2.0)
      endif
      if (bb.le.1.) then
        fwi(iday) = bb
      else
        fwi(iday) = exp(2.72 * (0.434 * log(bb))**0.647 )
      endif
      dsr(iday) = 0.0272 * fwi(iday)**1.77
c
      ffmo(iday)=ffmc(iday)
      dmco(iday)=dmc(iday)
      dco(iday)=dc(iday)
      buio(iday)=bui(iday)
      fwio(iday)=fwi(iday)
      xisio(iday)=xisi(iday)
c
       enddo
       enddo
       return
      end