#if defined(ROW_LAND)
#define SEA_P .true.
#define SEA_U .true.
#define SEA_V .true.
#elif defined(ROW_ALLSEA)
#define SEA_P allip(j).or.ip(i,j).ne.0
#define SEA_U alliu(j).or.iu(i,j).ne.0
#define SEA_V alliv(j).or.iv(i,j).ne.0
#else
#define SEA_P ip(i,j).ne.0
#define SEA_U iu(i,j).ne.0
#define SEA_V iv(i,j).ne.0
#endif
      subroutine convch(m,n)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
c
c --- hycom version 1.0
      implicit none
c
      integer m,n
c
cdiag real    sigup,uup,vup,siglo,ulo,vlo,coluin(idm),colout(idm)
      real    q,tem,sal,thet,trc(mxtrcr)
      real    dthet,plev,alfadt,betads
      integer i,iter,j,k,k1,ks,kp,ktr,margin
      logical llayer
c
      integer, parameter :: itmax=5
c
      include 'stmt_fns.h'
c
c --- ---------------------
c --- convective adjustment
c --- ---------------------
c
 103  format (i9,2i5,a/(33x,i3,2f8.3,f8.3,f8.2,f8.1))
cdiag if (itest.gt.0 .and. jtest.gt.0) then
cdiag   write (lp,103) nstep,itest+i0,jtest+i0,
cdiag&  '  entering convec:  temp    saln    dens    thkns    dpth',
cdiag&  (k,temp(itest,jtest,k,n),saln(itest,jtest,k,n),
cdiag&  th3d(itest,jtest,k,n)+thbase,dp(itest,jtest,k,n)*qonem,
cdiag&  p(itest,jtest,k+1)*qonem,k=1,kk)
cdiag  endif
c
      call xctilr(th3d(1-nbdy,1-nbdy,1,n),1,kk, 1,1, halo_ps)
c
c --- set thstar, note there is no thermobaric correction
c
      margin = 1
c
!$OMP PARALLEL DO PRIVATE(j,k,i,q)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
      do k=2,kk
      k1=k-1
      do i=1-margin,ii+margin
      if (SEA_P) then
      p(i,j,k)=p(i,j,k-1)+dp(i,j,k,n)
      if (.not.mxlkta .and. p(i,j,k).lt.dpmixl(i,j,n)) then
        nmlb(i,j,n)=k
      endif
      if (locsig) then
        p(i,j,k)=p(i,j,k-1)+dp(i,j,k,n)
        if (k.eq.2) then
          thstar(i,j,k,1)=th3d(i,j,k1,n)
        else
          if (k1.gt.nmlb(i,j,n)) then
            alfadt=dsiglocdt(ahalf*(temp(i,j,k1,n)+
     &                              temp(i,j,k ,n) ),
     &                       ahalf*(saln(i,j,k1,n)+
     &                              saln(i,j,k ,n) ),p(i,j,k))*
     &                             (temp(i,j,k1,n)-
     &                              temp(i,j,k, n) )
            betads=dsiglocds(ahalf*(temp(i,j,k1,n)+
     &                              temp(i,j,k ,n) ),
     &                       ahalf*(saln(i,j,k1,n)+
     &                              saln(i,j,k ,n) ),p(i,j,k))*
     &                             (saln(i,j,k1,n)-
     &                              saln(i,j,k, n) )
            thstar(i,j,k,1)=thstar(i,j,k1,1)-alfadt-betads
          else
            thstar(i,j,k,1)=thstar(i,j,k1,1)
          endif
        endif
      else
        thstar(i,j,k,1)=th3d(i,j,k,n)
      endif
      endif !ip
      enddo !i
      enddo !k
      enddo !j
!$OMP END PARALLEL DO
c
c --- convective adjustment
c
!$OMP PARALLEL DO PRIVATE(j,k,i,q)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1,jj
c
c --- convection of u
c
      do k=2,kk
      do i=1,ii
      if (SEA_U) then
      pu(i,j,k+1)=pu(i,j,k)+dpu(i,j,k,n)
      if (pu(i,j,k+1).lt.depthu(i,j)-onemm .and.
     &    thstar(i,j,k  ,1)+thstar(i-1,j,k  ,1).lt.
     &    thstar(i,j,k-1,1)+thstar(i-1,j,k-1,1)    ) then
cdiag   if (i.eq.itest .and. j.eq.jtest) then
cdiag   uup=u(i,j,k-1,n)
cdiag   ulo=u(i,j,k,  n)
cdiag   endif
        q=1.0-max(dpu(i,j,k,n),0.)/
     &   (max(dpu(i,j,k,n),0.)+max(dpu(i,j,k-1,n),onemm))
        u(i,j,k,  n)=u(i,j,k,n)+q*(u(i,j,k-1,n)-u(i,j,k,n))
        u(i,j,k-1,n)=u(i,j,k,n)
cdiag   if (i.eq.itest .and. j.eq.jtest) then
cdiag     write (lp,100) nstep,i+i0,j+j0,k,
cdiag&      1,'  upr,lwr,final u:',uup,ulo,u(i,j,k,n),q
cdiag   endif
      end if
      endif !ip
      enddo !i
      enddo !k
c
c --- convection of v
c
      do k=2,kk
      do i=1,ii
      if (SEA_V) then
      pv(i,j,k+1)=pv(i,j,k)+dpv(i,j,k,n)
      if (pv(i,j,k+1).lt.depthv(i,j)-onemm .and.
     &    thstar(i,j,k  ,1)+thstar(i,j-1,k  ,1).lt.
     &    thstar(i,j,k-1,1)+thstar(i,j-1,k-1,1)    ) then
cdiag   vup=v(i,j,k-1,n)
cdiag   vlo=v(i,j,k,  n)
        q=1.0-max(dpv(i,j,k,n),0.)/
     &   (max(dpv(i,j,k,n),0.)+max(dpv(i,j,k-1,n),onemm))
        v(i,j,k,  n)=v(i,j,k,n)+q*(v(i,j,k-1,n)-v(i,j,k,n))
        v(i,j,k-1,n)=v(i,j,k,n)
cdiag   if (i.eq.itest .and. j.eq.jtest) then
cdiag      write (lp,100) nstep,i+i0,j+i0,k,
cdiag&       1,'  upr,lwr,final v:',vup,vlo,v(i,j,k,n),q
cdiag   endif
      end if
      endif !ip
      enddo !i
      enddo !k
      enddo !j
!$OMP END PARALLEL DO
c
c --- convection of thermodynamical variables and tracer
c
!$OMP PARALLEL DO PRIVATE(j,i,k,ktr,iter,ks,kp,
!$OMP&                    q,tem,sal,thet,trc)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1,jj
      do i=1,ii
      if (SEA_P) then
c
ccc      coluin(i)=0.
ccc      colout(i)=0.
c
      do 12 k=1,kk
ccc      coluin(i)=coluin(i)+temp(i,j,k,n)*dp(i,j,k,n)
      p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n)
 12   continue
c
      do 6 iter=1,itmax
c
      klist(i,j)=1
      util1(i,j)=dp(i,j,1,n)
c
      do 6 k=2,kk
      k1=k-1
c
      ks=klist(i,j)
      if (locsig) then
        plev=0.5*(p(i,j,ks+1)+p(i,j,k))
        alfadt=0.5*
     &        (dsiglocdt(temp(i,j,ks,n),saln(i,j,ks,n),plev)+
     &         dsiglocdt(temp(i,j,k ,n),saln(i,j,k ,n),plev))*
     &        (temp(i,j,ks,n)-temp(i,j,k,n))
        betads=0.5*
     &        (dsiglocds(temp(i,j,ks,n),saln(i,j,ks,n),plev)+
     &         dsiglocds(temp(i,j,k ,n),saln(i,j,k ,n),plev))*
     &        (saln(i,j,ks,n)-saln(i,j,k,n))
        dthet=alfadt+betads
      else
        dthet=th3d(i,j,ks,n)-th3d(i,j,k,n)
      endif
      if (dthet.gt.0.0) then
        if (dp(i,j,k,n).lt.onemm) then
          if (locsig) then
            alfadt=dsiglocdt(ahalf*(temp(i,j,k1,n)+
     &                              temp(i,j,k ,n) ),
     &                       ahalf*(saln(i,j,k1,n)+
     &                              saln(i,j,k ,n) ),p(i,j,k))*
     &                             (temp(i,j,k1,n)-
     &                              temp(i,j,k, n) )
            betads=dsiglocds(ahalf*(temp(i,j,k1,n)+
     &                              temp(i,j,k ,n) ),
     &                       ahalf*(saln(i,j,k1,n)+
     &                              saln(i,j,k ,n) ),p(i,j,k))*
     &                             (saln(i,j,k1,n)-
     &                              saln(i,j,k, n) )
            dthet=alfadt+betads
          else
            dthet=th3d(i,j,k1,n)-th3d(i,j,k,n)
          endif
          if(dthet.ge.0.0) then
            saln(i,j,k,n)=saln(i,j,k1,n)
            temp(i,j,k,n)=temp(i,j,k1,n)
            th3d(i,j,k,n)=sig(temp(i,j,k,n),saln(i,j,k,n))-thbase
            do ktr= 1,ntracr
              tracer(i,j,k,n,ktr)=tracer(i,j,k1,n,ktr)
            enddo
          end if
        else				!  dp > onemm
          if (iter.eq.itmax) then
            if (locsig) then
              plev=0.5*(p(i,j,ks+1)+p(i,j,k))
              alfadt=0.5*
     &              (dsiglocdt(temp(i,j,ks,n),saln(i,j,ks,n),plev)+
     &               dsiglocdt(temp(i,j,k ,n),saln(i,j,k ,n),plev))*
     &              (temp(i,j,ks,n)-temp(i,j,k,n))
              betads=0.5*
     &              (dsiglocds(temp(i,j,ks,n),saln(i,j,ks,n),plev)+
     &               dsiglocds(temp(i,j,k ,n),saln(i,j,k ,n),plev))*
     &              (saln(i,j,ks,n)-saln(i,j,k,n))
              dthet=alfadt+betads
            else
              dthet=th3d(i,j,ks,n)-th3d(i,j,k,n)
            endif
            if (dthet.gt.sigjmp .and. dp(i,j,k,n).gt.onem) then
              ! only write out the lowest unstable layer
              llayer = k.eq.kk
              if     (.not.llayer) then
                llayer =   dp(i,j,k+1,n).lt.onem .or.
     &                   th3d(i,j,k+1,n).ge.th3d(i,j,ks,n)-sigjmp
              endif
              if     (llayer) then
!$OMP           CRITICAL
                write (lp,'(i9,2i5,i3,a,i3,a,i3,a,2f10.4)')
     &           nstep,i+i0,j+j0,k,
     &           ' colmn unstbl (wrt',ks,') after',
     &           iter-1,' its',
     &           th3d(i,j,ks,n)+thbase,th3d(i,j,k,n)+thbase
!$OMP           END CRITICAL
              endif
            endif
          else				!  it < itmax
cdiag       sigup=th3d(i,j,ks,n)
cdiag       siglo=th3d(i,j,k, n)
            util1(i,j)=util1(i,j)+dp(i,j,k,n)
            q=1.0-max(dp(i,j,k,n),0.5*onemm)/max(util1(i,j),onemm)
            tem=temp(i,j,k,n)+q*(temp(i,j,ks,n)-temp(i,j,k,n))
            sal=saln(i,j,k,n)+q*(saln(i,j,ks,n)-saln(i,j,k,n))
            do ktr= 1,ntracr
              trc(ktr)=tracer(i,j,k,n,ktr)+q*(tracer(i,j,ks,n,ktr)-
     &                                        tracer(i,j,k,n,ktr))
            enddo
            thet=sig(tem,sal)-thbase
            do 10 kp=ks,k
            temp(i,j,kp,n)=tem
            saln(i,j,kp,n)=sal
            th3d(i,j,kp,n)=thet
            do ktr= 1,ntracr
              tracer(i,j,kp,n,ktr)=trc(ktr)
            enddo
 10         continue
c
cdiag if (i.eq.itest .and. j.eq.jtest) then
cduag   write (lp,100) nstep,i+i0,j+j0,k,iter,
cdiag&    '  upr,lwr,final dens:',(sigup+thbase),
cdiag&    (siglo+thbase),(th3d(i,j,k,n)+thbase),q
cdiag endif
 100    format (i9,2i5,i3,'  it',i2,a,3f8.3,f5.2)
c
          end if
        end if
      else				!  th3d(kn) > th3d(ksn)
        klist(i,j)=k
        util1(i,j)=dp(i,j,k,n)
      end if
 6    continue
c
ccc      do k=1,kk
ccc      colout(i)=colout(i)+temp(i,j,k,n)*dp(i,j,k,n)
ccc      enddo !k
ccc      if (abs((colout(i)-coluin(i))/coluin(i)).gt.1.e-6)
ccc     .  write (lp,'(i9,2i5,a/1p,3e14.6)') nstep,i,j,
ccc     .  '  column integral not conserved in convec:',
ccc     .  coluin(i),colout(i),(colout(i)-coluin(i))/coluin(i)
      endif !ip
      enddo !i
      enddo !j
!$OMP END PARALLEL DO
c
cdiag if (itest.gt.0 .and. jtest.gt.0) then
cdiag   write (lp,103) nstep,itest+i0,jtest+j0,
cdiag&  '  exiting  convec:  temp    saln    dens    thkns    dpth',
cdiag&  (k,temp(itest,jtest,k,n),saln(itest,jtest,k,n),
cdiag&  th3d(itest,jtest,k,n)+thbase,dp(itest,jtest,k,n)*qonem,
cdiag&  p(itest,jtest,k+1)*qonem,k=1,kk)
cdiag endif
c
      return
      end
c
      subroutine convcm(m,n)
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
c
c --- hycom version 1.0 (adapted from micom version 2.8)
      implicit none
c
      integer m,n
c
      integer i,j,k,k1
      real    dthet,delp,alfadt,betads
c
      include 'stmt_fns.h'
c
 103  format (i9,2i5,a/(33x,i3,2f8.3,3p,f8.3,0p,f8.2,f8.1))
cdiag if (itest.gt.0 .and. jtest.gt.0) then
cdiag   write (lp,103) nstep,itest+i0,jtest+j0,
cdiag&  '  entering convec:  temp    saln    dens    thkns    dpth',
cdiag&  (k,temp(itest,jtest,k,n),saln(itest,jtest,k,n),
cdiag&  th3d(itest,jtest,k,n)+thbase,dp(itest,jtest,k,n)*qonem,
cdiag&  p(itest,jtest,k+1)*qonem,k=1,kk)
cdiag endif
c
c --- -------------------------------------------------------------
c --- entrain water lighter than mixed-layer water into mixed layer
c --- -------------------------------------------------------------
c
!$OMP PARALLEL DO PRIVATE(j,k,i,dthet,delp)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1,jj
      do i=1,ii
      if (SEA_P) then
      klist(i,j)=0
      p(i,j,2)=dp(i,j,1,n)
      dpo(i,j,1,n)=dp(i,j,1,n)
c
      do k=2,kk
      p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n)
      dpo(i,j,k,n)=0.
      enddo !k
c
      do k=2,kk
      k1=k-1
c
      p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n)
      if (dp(i,j,k,n).le.0. .or.
     &   (klist(i,j).gt.0 .and. k.gt.klist(i,j)+1)) exit !do k
      if (locsig) then
        alfadt=dsiglocdt(ahalf*(temp(i,j,k1,n)+
     &                          temp(i,j,k ,n) ),
     &                   ahalf*(saln(i,j,k1,n)+
     &                          saln(i,j,k ,n) ),p(i,j,k))*
     &                         (temp(i,j,k1,n)-
     &                          temp(i,j,k, n) )
        betads=dsiglocds(ahalf*(temp(i,j,k1,n)+
     &                          temp(i,j,k ,n) ),
     &                   ahalf*(saln(i,j,k1,n)+
     &                          saln(i,j,k ,n) ),p(i,j,k))*
     &                         (saln(i,j,k1,n)-
     &                          saln(i,j,k, n) )
        dthet=-alfadt-betads
      else
        dthet=th3d(i,j,k,n)-th3d(i,j,1,n)
      endif
      if (dthet.lt.0. .and. p(i,j,k+1).le.p(i,j,kk+1)) then
c
cdiag   if (i.eq.itest .and. j.eq.jtest) then
cdiag   write (lp,100) nstep,i+i0,j+j0,' convec',1,th3d(i,j,1,n)
cdiag&  +thbase,dp(i,j,1,n)*qonem,temp(i,j,1,n),saln(i,j,1,n),k,
cdiag&  th3d(i,j,k,n)+thbase,dp(i,j,k,n)*qonem,temp(i,j,k,n),saln(i,j,k,n)
cdiag endif
 100    format (i9,2i5,a,i3,'  th3d,dp,t,s =',3pf7.3,0pf7.1,2f8.3
     &    /26x,i3,15x,3pf7.3,0pf7.1,2f8.3)
c
c --- layer -k- contains mass less dense than mixed layer. entrain it.
        delp=dp(i,j,1,n)+dp(i,j,k,n)
        saln(i,j,1,n)=(saln(i,j,1,n)*dp(i,j,1,n)
     &                +saln(i,j,k, n)*dp(i,j,k, n))/delp
        temp(i,j,1,n)=(temp(i,j,1,n)*dp(i,j,1,n)
     &                +temp(i,j,k, n)*dp(i,j,k, n))/delp
        th3d(i,j,1,n)=sig(temp(i,j,1,n),saln(i,j,1,n))-thbase
c
        diaflx(i,j,1)=diaflx(i,j,1)+dp(i,j,k,n)			! diapyc.flux
        diaflx(i,j,k)=diaflx(i,j,k)-dp(i,j,k,n)			! diapyc.flux
c
c --- mass in layer -k- transferred to mixed layer is stored in -dpo-.
        dp(i,j,1,n)=delp
        dpo(i,j,k,n)=dp(i,j,k,n)
        dp(i,j,k,n)=0.
        klist(i,j)=k
      end if                 !  dthet < 0
c
      enddo !k
      endif !ip
      enddo !i
      enddo !j
!$OMP END PARALLEL DO
c
      do j=1,jj
        do i=1,ii
          if (ip(i,j).eq.1) then
            util1(i,j) = klist(i,j)
          else
            util1(i,j) = 0.0
          endif
        enddo !j
      enddo !i
      call xctilr(util1,                 1, 1, 1,1, halo_ps)
      call xctilr(dpo(1-nbdy,1-nbdy,1,n),1,kk, 1,1, halo_ps)
      do j=0,jj+1
        do i=0,ii+1
          klist(i,j) = util1(i,j)
        enddo !j
      enddo !i
c
!$OMP PARALLEL DO PRIVATE(j,k,i)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1,jj
c
      do i=1,ii
      if (SEA_P) then
      dpmixl(i,j,n)=dp(i,j,1,n)
      endif !ip
c
c --- entrain -u- momentum
c
      if (SEA_U) then
      util2(i,j)=min(0.5*(dpo(i,j,1,n)+dpo(i-1,j,1,n)),depthu(i,j))
      do k=2,max(klist(i,j),klist(i-1,j))
      util1(i,j)=util2(i,j)
      util2(i,j)=min(util2(i,j)+0.5*(dpo(i,j,k,n)+dpo(i-1,j,k,n)),
     &               depthu(i,j))
      u(i,j,1,n)=(u(i,j,1,n)*util1(i,j)
     &           +u(i,j,k,n)*(util2(i,j)-util1(i,j)))/util2(i,j)
      enddo !k
      endif !iu
c
c --- entrain -v- momentum
c
      if (SEA_V) then
      util2(i,j)=min(0.5*(dpo(i,j,1,n)+dpo(i,j-1,1,n)),
     &               depthv(i,j))
      do k=2,max(klist(i,j),klist(i,j-1))
      util1(i,j)=util2(i,j)
      util2(i,j)=min(util2(i,j)+0.5*(dpo(i,j,k,n)+dpo(i,j-1,k,n)),
     &               depthv(i,j))
      util2(i,j)=min(util2(i,j)+0.5*(dpo(i,j,k,n)+dpo(i,j-1,k,n)),
     &               depthv(i,j))
      v(i,j,1,n)=(v(i,j,1,n)*util1(i,j)
     &           +v(i,j,k,n)*(util2(i,j)-util1(i,j)))/util2(i,j)
      enddo !k
      endif !ip
      enddo !i
c
      enddo !j
!$OMP END PARALLEL DO
c
cdiag if (itest.gt.0 .and. jtest.gt.0) then
cdiag   write (lp,103) nstep,itest+i0,jtest+j0,
cdiag&  '  exiting  convec:  temp    saln    dens    thkns    dpth',
cdiag&  (k,temp(itest,jtest,k,n),saln(itest,jtest,k,n),
cdiag&  th3d(itest,jtest,k,n)+thbase,dp(itest,jtest,k,n)*qonem,
cdiag&  p(itest,jtest,k+1)*qonem,k=1,kk)
cdiag endif
c
      return
      end
c
c
c> Revision history:
c>
c> July 1997 - deleted final pressure calculation loop (moved to to thermf.f)
c> Apr. 1999 - added calculation of -th3d-
c> Aug. 2000 - convcm: adapted from micom 2.8 to run within hycom 1.0
c> Oct. 1999 - convch: convection of u and v added
c> Oct. 2009 - convcm: MPI bug fix.
c> Oct  2010 - replaced two calls to dsiglocdX with one call at mid-point
c> Aug  2011 - replaced dpold,dpoldm with dpo
c> May  2014 - use land/sea masks (e.g. ip) to skip land