c    @(#)bndheat.f	1.1
      subroutine bndheat
c
c
c
c leonard walstad
c
c-----------------------------------------------------------------------
c
c
      include'mladv.inc'
c
c----------------------------------------------------------------------
c
      common/parm/kz,m,mm1,mm2,mp1,msq,tim,tmax,tstart,hy,rea,dt,
     +   xbasin,yy(mx0),icalc,kzp1,kzm1,iftop,ifbot,iftvv,ifrel,theta,
     +   ifsbl,ifeva
                

c
c----------------------------------------------------------------------
c
      data icall/0/
c
c
c----------------------------------------------------------------------
c
      if(icall.ne.0)go to 99
        icall=1
        do 10 k=1,l
          z1=zml(k)+dzml(k)/2.0
          z2=z1-dzml(k)
          absrp(k)=.01*(exp(decay*z1) - exp(decay*z2))/dzml(k)
          if(absrp(k).lt. 0.0001*decay)go to 10
          kmax=k
   10   continue
c
c
   99 continue
c
c
      do 110 k=1,kmax
        abs=absrp(k)
        f1=abs*qgts
        f2=ag*f1
        do 105 ibnd=1,mbnd
          ip=ipofb(ibnd)
          bndtemp(ibnd,k)=bndtemp(ibnd,k)+f1*qsw(ip)
          bndbuoy(ibnd,k)=bndbuoy(ibnd,k)+f2*qsw(ip)
  105   continue
  110 continue
c
c
c
      f1=qgts/100.
      f3t=ag*f1
      dz=dzml(1)
      dz2=dz+.5*dzml(2)
      do 120 ibnd=1,mbnd
        ip=ipofb(ibnd)
         h=bnddpth(ibnd)
         dzt=dz
         if(h.lt.dz2)dzt=h
         f2=f1/dzt
         f3=f3t/dzt
         bndtemp(ibnd,1)=bndtemp(ibnd,1)+f2*qsurf(ip)
         bndbuoy(ibnd,1)=bndbuoy(ibnd,1)+f3*qsurf(ip)+
     .                    f1*wflux(ip)*bbar(ip)
  120 continue
      return
      end
