c
c    @(#)mix.f	1.4
c
      subroutine mix
c
c
c
c leonard walstad
c
c
c
c----------------------------------------------------------------------
c
c
      include'mladv.inc'
c
c
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
      common /extpar/ rlat0,rlng0,t0,v0,dhor,ht,time0,r0,f0star
c
      data icall/0/
c
c----------------------------------------------------------------------
c
c
c
      h1=-zml(1)
      hb=-zml(l)-.10
      entmiss=0.0
c
c
        if(h.lt.h1)h=h1
        if(h.ge.hb)h=hb
        if(hmx.lt.h1)hmx=h1
    9   continue
        kml=0
   10   continue
          kml=kml+1
          dkml=-zml(kml)
          dkmlp1=-zml(kml+1)
          if(h.ge.dkmlp1)go to 10
        if(h.gt.hmx)go to 110
c
        btot=entr
        if(btot.lt.0.0)btot=0.0
c
        dz=.5*dzml(kml)+(h-dkml)
        mltint=dz*t1d(kml)
        mlbint=dz*b1d(kml)
        do 11 i=1,ntrace
           trcint(i)=dz*trc1d(kml,i)
  11     continue
c
c
c
        if(kml.eq.1)go to 20
c
        do 15 k=1,kml-1
          dz=dzml(k)
          mltint=mltint+dz*t1d(k)
          mlbint=mlbint+dz*b1d(k)
          do 14 i=1,ntrace
            trcint(i)=trcint(i)+dz*trc1d(k,i)
   14      continue
   15   continue
c
   20   continue
        mltint=mltint/h
        mlbint=mlbint/h
        do 21 i=1,ntrace
          trcint(i)=trcint(i)/h
   21    continue
c
   25   continue
          kmlp1=kml+1
          dkmlp1=-zml(kmlp1)
          dnext=dkmlp1+.5*dzml(kmlp1)
          dmax=min(dnext,hb,hmx)
          dzmax=dmax-h
c
          delt=mltint-t1d(kmlp1)
          delb=mlbint-b1d(kmlp1)
          do 26 i=1,ntrace
            deltrc(i)=trcint(i)-trc1d(kmlp1,i)
   26      continue
c
          dz=dzmax
          if(delb.gt.0.0)dz=btot/delb
          dz=min(dz,dzmax)
c
          hold=h
          h=h+dz
          mltint=mltint-dz*delt/h
          mlbint=mlbint-dz*delb/h
          do 27 i=1,ntrace
             trcint(i)=trcint(i)-dz*deltrc(i)/h
   27      continue
          btot=btot-dz*delb
c
c
          if(h.ge.dkmlp1)go to 90
            if(btot.le.0.0)btot=0.0
            if((h+.001).gt.hmx)btot=0.0
            entmiss=btot
            go to 130
c
   90     continue
          kml=kml+1
          if(dz.lt.dzmax)go to 95
          if(hmx.le.dnext)go to 94
          if(hb.le.dnext)go to 94
c
          goto 100
c
   94     continue
          if(hb.le.hmx)go to 95
          btot=0.0
   95     continue
            kmlp1=kml+1
            dzp1=dzml(kmlp1)
            dzl=dnext-h
            dtot=dzl+dzp1
            t1d(kmlp1)=(dzl*t1d(kml) + dzp1*t1d(kmlp1))/dtot 
            b1d(kmlp1)=(dzl*b1d(kml) + dzp1*b1d(kmlp1))/dtot 
            do 96 i=1,ntrace
               trc1d(kmlp1,i)=(dzl*trc1d(kml,i) + 
     .                         dzp1*trc1d(kmlp1,i))/dtot
   96        continue
            if(btot.le.0.0)btot=0.0
            entmiss=btot
            go to 130
c
  100    continue
         go to 25
c
c detrainment
c
  110   continue
        entr=0.0
        hold=h
        h=hmx
        kmlp1=kml+1
        dkml=-zml(kml)
        hkml=dkml-.5*dzml(kml)
        dkmlp1=-zml(kmlp1)
        hkmlp1=dkmlp1-.5*dzml(kmlp1)
        hkmlp2=dkmlp1+.5*dzml(kmlp1)
c
        if(h.lt.dkml)go to 120
          dh=hold-h
          dzl=hkmlp2-hold
          dtot=dzl+dh
          t1d(kmlp1)=(dh*t1d(kml) + dzl*t1d(kmlp1))/dtot
          b1d(kmlp1)=(dh*b1d(kml) + dzl*b1d(kmlp1))/dtot
          do 111 i=1,ntrace
             trc1d(kmlp1,i)=(dh*trc1d(kml,i) + dzl*trc1d(kmlp1,i))/dtot
  111      continue
          go to 125
c
  120   continue
        if(hold.lt.hkmlp1)go to 122
          dzl=hold-hkmlp1
          dh=hkmlp2-hold
          dtot=dh+dzl
          t1d(kmlp1) = (dzl*t1d(kml)+ dh*t1d(kmlp1))/dtot
          b1d(kmlp1) = (dzl*b1d(kml)+ dh*b1d(kmlp1))/dtot
          do 121 i=1,ntrace
             trc1d(kmlp1,i)=(dzl*trc1d(kml,i) + dh*trc1d(kmlp1,i))/dtot
  121      continue
          go to 125
c
  122   continue
        dzl=hkmlp1-hold
        dh=min(hold-h,hold-hkml)
        dtot=dzl+dh
        t1d(kml)=(dzl*t1d(kmlp1) + dh*t1d(kml))/dtot 
        b1d(kml)=(dzl*b1d(kmlp1) + dh*b1d(kml))/dtot 
        do 123 i=1,ntrace
           trc1d(kml,i)=(dzl*trc1d(kmlp1,i) + dh*trc1d(kml,i))/dtot
  123    continue
        go to 125
c
c
c
  125   continue
        goto 9
c
c
  130   continue
        do 140 k=1,kml
          b1d(k)=mlbint
          t1d(k)=mltint
          do 131 i=1,ntrace
            trc1d(k,i)=trcint(i)
  131      continue
  140   continue
        dmix=h
        kmix=kml
        kmixp1=kmix+1
        if(abs(mltint-t1d(kmixp1)).ge.tthresh)go to 160
  150   continue
          kmix=kmix+1
          dmix=-zml(kmix)
          if(kmix.ge.l)go to 160
          kmixp1=kmix+1
          delt=t1d(kmix)-t1d(kmixp1)
          dzmax=-zml(kmixp1)-dmix
          dz=dzmax
          sgn=1.0
          if(delt.lt.0.)sgn=-1.0
          if(delt.ne.0.0)dz=dz*(t1d(kmix)-mltint+sgn*tthresh)/delt
c
          dz=min(dz,dzmax)
          if(dz.eq.dzmax)go to 150
        dmix=dmix+dz
  160   continue
  200 continue

      return
      end
