#include "options.inc"
    

c=======================================================================
c      account for static instabilities
c=======================================================================

      subroutine convect
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      K_b(:,js:je,:) = 0.0
      do j=js,je
       do k=1,km-1
        do i=2,imt
         if (b(i,j,k+1,taum1) > b(i,j,k,taum1)) 
     &         K_b(i,j,k)=1000.*maskW(i,j,k)
        enddo
       enddo
      enddo
      call implicit_tracer_mix(b,K_b)
c      call convect_expl
      end subroutine convect


      subroutine implicit_tracer_mix(var,K_vert)
      use cpflame_module
      implicit none
      real :: var(imt,jmt,km,0:2)
      real :: K_vert(imt,jmt,km)
      integer :: j,k,js,je
      real :: a(imt,km),bb(imt,km),c(imt,km),bet(imt)
      real :: pu(imt,km),gam(imt,km),fxa, r(imt,km)

      js=max(2,js_pe); je = min(je_pe,jmt-1)
      fxa = c2dt/dz**2

      do j=js,je

       bb(:,1) = 1+fxa*K_vert(:,j,1)
       c(:,1)  =  -fxa*K_vert(:,j,1)
       do k=2,km-1
        a(:,k)  =  -fxa *   K_vert(:,j,k-1)
        bb(:,k) = 1+fxa * ( K_vert(:,j,k)+K_vert(:,j,k-1) )
        c(:,k)  =  -fxa *   K_vert(:,j,k)
       enddo
       a(:,km)  =  -fxa * K_vert(:,j,km-1)
       bb(:,km) = 1+fxa * K_vert(:,j,km-1) 

       pu=0.0;gam=0.0
       r=var(:,j,:,taup1)
       r(:,km-1)=r(:,km-1)
c     &           +aidif_tracer*surf_flux(:,j)*c2dt/dz
       bet=bb(:,1)
       where (bet/=0.0) pu(:,1)=r(:,1)/bet
       do k=2,km
        where (bet/=0.0) gam(:,k)=c(:,k-1)/bet
        bet=bb(:,k)-a(:,k)*gam(:,k)
        where(bet/=0.0) pu(:,k)=(r(:,k)-a(:,k)*pu(:,k-1))/bet
       enddo
       do k=km-1,1,-1
        pu(:,k)=pu(:,k)-gam(:,k+1)*pu(:,k+1)
       enddo
       var(:,j,:,taup1)=pu
      enddo
      end subroutine implicit_tracer_mix





      subroutine convect_expl
c-----------------------------------------------------------------------
c      explicit convection
c-----------------------------------------------------------------------
      use cpflame_module
      implicit none
      integer :: i,j,k,js,je,n
      real :: dztsum,bmix,bsum
      integer :: l,l1,kcon,lcon,lcona,lconb,lmix
      integer :: lctot(imt,jmt),lcven(imt,jmt)
c
c     check each row column by column; note that 'goto 1310'
c     finishes a particular column and moves to the next one.
c
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js,je
        do i=2,imt
          kcon     = kmt(i,j)
c
c         1. initial search for uppermost unstable pair; if none is
c            found, move on to next column
c
          lcon       = 0
          do k=kcon+1,km
            if (b(i,j,k,taup1) > b(i,j,k-1,taup1)) lcon = k
          enddo
          if (lcon==0) cycle

1319      lcona = lcon
          lconb = lcon - 1
c
c         2. mix the first two unstable layers
c
          dztsum = 2*dz 
          bsum = b(i,j,lcona,taup1)*dz +b(i,j,lconb,taup1)*dz
          bmix = bsum / dztsum
          b(i,j,lconb:lcona,taup1) = bmix /dztsum
c
c         3. test layer below lconb
c
1306      continue
          if (lconb==kcon) goto 1308

          if (b(i,j,lconb,taup1) > b(i,j,lconb-1,taup1)) then
            lconb = lconb-1
            dztsum = dztsum + dz
            bsum = bsum + b(i,j,lconb,taup1)*dz
            bmix = bsum / dztsum
            b(i,j,lconb:lcona,taup1) = bmix
            goto 1306
          end if
c
c         4. test layer above lcona
c
1308      continue
          if (lcona < km-1) then
            if (b(i,j,lcona+1,taup1) > b(i,j,lcona,taup1)) then
              lcona = lcona+1
              dztsum = dztsum + dz
              bsum = bsum + b(i,j,lcona,taup1)*dz
              bmix = bsum / dztsum 
              b(i,j,lconb:lcona,taup1) = bmix
              goto 1306
            end if
          end if
c
c         6. resume search if step 3. and 4. have been passed and this
c            unstable part of the water column has thus been removed,
c            i.e. find further unstable areas further down the column
c
          if (lconb == kcon) cycle
          lcon = lconb
c
1302      continue
          lcon = lcon - 1
          if (lcon == kcon) cycle
          if (b(i,j,lcon,taup1) <= b(i,j,lcon-1,taup1)) goto 1302
          goto 1319

1310      continue
        enddo
      enddo
      end subroutine convect_expl



