c
c     @(#)mladv.f	1.14
c
      subroutine mladv
c
c leonard walstad
c
c
c
c main routine for the advected mixed layer model.  performs 
c time stepping and adjust pointers etc.
c
c----------------------------------------------------------------------
c
c
      include'mladv.inc'
c
      external bzero
      external btrace
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
                
      common /extpar/ rlat0, rlng0, t0, v0, dhor, ht, time0,
     +                r0,f0star

      common /scrat/r(mn0),r1(mn0),r2(mn0),r3(mn0)

      common/sbldef/qsw0,qsurf0,wflux0,u10m,angu10m,hdpth0
c
      data mlcur/1/, mlpast/2/, icall/0/, istep/12/
c
c
c
c
c
c
      if(icall.ne.0)go to 59
c
        icall=1
        mltime=t0*tim/86400. + time0
        call setup
        if (ifrst.eq.1) then
           icall = 2
           mlstep = qgts
           return
         endif
        call wget
        call getuv
        return
c
c
   59 continue
      if(icall.ne.1)go to 99
        icall=2
        mltime=t0*tim/86400. + time0
        mlstep=qgts+qgts
        call initial
        call boundary
        call mixall
        call bndmix
        call mlfilt
        do 70 itrace=1,ntrace
            if (ipdt(itrace).eq.ipdty) call posdef(itrace,mlcur)
  70     continue
        call wget
        call getuv
        call surflux
        call curl
        if ((ntrace .gt. 0) .and. (ibudf .eq. 1)) call budget
        if ((ntrace .gt. 0) .and. (ibudi .eq. 1)) call ibudget
c
        imlcnt=0
        return
c
c
c
   99 continue

      mltime=t0*tim/86400. + time0 - qgts/86400.

      do 71 itrace=1,ntrace
         if (ipdt(itrace).eq.ipdty) call posdef(itrace,mlcur)
  71   continue

      call surflux
      call mlout
      call wget
c
  100 continue
        if(imlcnt.eq.istep)imlcnt=0
        mladd=mlpast
        mlstep=qgts+qgts
        if(imlcnt.ne.0)go to 101
          mladd=mlcur
          mlstep=-qgts
c
  101 continue
      imlcnt=imlcnt+1
c
c
c

      if (ifvadv.eq.3) call cross

      call advset

      do 110 ip=1,msq
        tmp=bqgofk(ip,1)+buoy(ip,1,mlcur)
        r3(ip)=tmp
        r1(ip)=tmp*tauy(ip)/hdpth(ip,mlcur)/f0star
        r2(ip)=-tmp*taux(ip)/hdpth(ip,mlcur)/f0star
  110 continue
c
      call advectqg
c
      if (kmlmin .ne. kmlmax) then
         call advtrans(buoy(1,1,mlcur),1)
       endif
c
      call advect3d(buoy(1,1,mlcur),buoy(1,1,mladd),
     .                              buoy(1,1,mlpast),bzero)
c
c
      do 120 ip=1,msq
        tmp=tqgofk(ip,1)+temp(ip,1,mlcur)
        r3(ip)=tmp
        r1(ip)=tmp*tauy(ip)/hdpth(ip,mlcur)/f0star
        r2(ip)=-tmp*taux(ip)/hdpth(ip,mlcur)/f0star
  120 continue
c
      call advectqg
c
      if (kmlmin .ne. kmlmax) then
         call advtrans(temp(1,1,mlcur),1)
       endif
c
      call advect3d(temp(1,1,mlcur),temp(1,1,mladd),
     .                              temp(1,1,mlpast),bzero)
  130 continue
c
c

      do 150 itrace=1,ntrace

         icurtrace=itrace

         do 135 ip=1,msq
            tmp=trcqgofk(ip,1)+trace(ip,1,mlcur,itrace)
            r3(ip)=tmp
            r1(ip)=tmp*tauy(ip)/hdpth(ip,mlcur)/f0star
            r2(ip)=-tmp*taux(ip)/hdpth(ip,mlcur)/f0star
  135     continue
c
          call advectqg
c
          if (kmlmin .ne. kmlmax) then
             call advtrans(trace(1,1,mlcur,itrace),1)
           endif

          call advect3d(trace(1,1,mlcur,itrace),
     .                  trace(1,1,mladd,itrace),
     .                  trace(1,1,mlpast,itrace),btrace)

          if (wsink(itrace) .ne. 0.0) then
             call sink(trace(1,1,mlcur,itrace),
     .                 trace(1,1,mladd,itrace),
     .                 trace(1,1,mlpast,itrace)) 
            endif

  150   continue

c
c
      if (ifperh .eq. 0) then
         call advectml(hdpth(1,mlcur),hdpth(1,mladd),hdpth(1,mlpast))
        endif
c
c
      if (imlcnt.eq.1) then
         if (ifperh .eq. 1) then
            do 160 ip=1,msq
               hdpth(ip,mlcur)  = hdpth0
               hdpth(ip,mlpast) = hdpth0
  160        continue
           endif
         go to 100
       endif
c
c
      if (ifperh .eq. 1) goto 201

      do 200 ip=1,msq
        h=hdpth(ip,mlcur)
        call mlbflux(h,ip,usbs)
        ustar=sqrt(sqrt(taux(ip)**2 + tauy(ip)**2))
        call erate(ustar,usbs,h,hmax(ip),be(ip),ier)
        hmax(ip)=min(hmax(ip),maxh)
  200 continue

  201 continue

      call heat

      if (ifwbg .ne. 0) then
         call wgetbg
         call crossbg
         call advectbg
        endif

      call mixall

c     if (itphys .eq. 2) then
c        icurtrace = 1
c        do 205 ip=1,msq
c           out(ip) = trace(ip,1,mlcur,1) + trcqgofk(ip,1)
c 205     continue
c      endif

      if ((ntrace .gt. 0) .and. (itphys .gt. 0)) call trcphys(1)

c
      if (ifpers .ne. 0) goto 220
c
      if (ifperh .eq. 1) goto 211

        do 210 ibnd=1,mbnd
          h=bnddpth(ibnd)
          ip=ipofb(ibnd)
          call mlbflux(h,ip,usbs)
          ustar=sqrt(sqrt(taux(ip)**2 + tauy(ip)**2))
          call erate(ustar,usbs,h,hmax(ip),bndbe(ibnd),ier)
          hmax(ip)=min(hmax(ip),maxh)
          if(ier.ne.0)be(ip)=0.0
  210   continue

  211   continue

        call bndheat
        call bndmix
        if (ntrace .gt. 0) call bndtphys
c
c
  220 continue

      itmp=mlcur
      mlcur=mlpast
      mlpast=itmp
c
c
      call heat

      if (ifwbg .ne. 0) then
         call crossbg
         call advectbg
        endif

      call mixall

      if ((ntrace .gt. 0) .and. (itphys .gt. 0)) call trcphys(2)
c
      mltime=mltime+qgts/86400.
c
      if (ifpers .eq. 1) then

         call boundary

         do 230 ibnd=1,mbnd
            bndbe(ibnd)=0.0
            hmax(ipofb(ibnd))=maxh
  230     continue

         call bndmix

       endif

 
      if (ifpers .ge. 0) then
         call cfvn
        endif

      call mlfilt
      call stinst

      if ((ntrace .gt. 0) .and. (ibudf .eq. 1)) call budget
      if ((ntrace .gt. 0) .and. (ibudi .eq. 1)) call ibudget
c
      if(tim.gt.tmax)call mlout
c
c --- now finished with this time step, time to get read for the next
c
c

      call getuv

      end
