subroutine cfl (j)
c
c=======================================================================
c                                                                    ===
c  This routine monitors the CFL (Courant-Friedrichs-Lewy) condition ===
c  necessary for stability: the Courant number be less than one.     ===
c                                                                    ===
c                         V dt/ds < 1                                ===
c                                                                    ===
c  (Adapted from MOM)                                                ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <scalar.h>
#include <fullwd.h>
#include <onedim.h>
#include <workspa.h>
#include <vertslabs.h>
#include <iounits.h>
#include <ccfl.h>
#if defined ext_tide & defined advtide
#  include <tidesp.h>
#endif
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,j,k
      FLOAT
     *      cflu,cflv,cflwt,cflwu,dtmax,f1,f2,f3t,f3u,umax,uwk,vmax,
     *      vwk,wtmax,wumax
#ifndef nocfl
     *      ,pcflu,pcflv,pcflwt,pcflwu
#endif
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
      dtmax=max(dtuv,dtts)
      f1=dtmax*csr(j)
      f2=dtmax*dyur(j)
      do 130 k=1,km
      do 130 i=2,imtm1
        f3t=dtmax/dzzqz(i,k,0)
        f3u=dtmax/dzzvqz(i,k,0)
#if !defined ext_tide | !defined advtide
        uwk =u(i,k)
        vwk =v(i,k)
# else
# ifndef advtide0
        uwk =u(i,k)+utide(i,k)
        vwk =v(i,k)+vtide(i,k)
# else
        uwk =u(i,k)+sadv*utide(i,k)
        vwk =v(i,k)+sadv*vtide(i,k)
# endif
#endif
        cflu=abs(f1*dxur(i)*uwk*gm(i,k))
        cflv=abs(f2*vwk*gm(i,k))
        cflwt=abs(f3t*w (i,k)*gm(i,k))
        cflwu=abs(f3u*wu(i,k)*gm(i,k))
#ifndef nocfl
# ifdef nowcfl
        if((cflu+cflv).ge.c1) then
# else
        if(((cflu+cflv+cflwt).ge.c1).or.((cflu+cflv+cflwu).ge.c1)) then
# endif
          write(stdout,100) i,j,k
 100      format(' CFL condition exceeded at coordinate (i,j,k) = ',
     *           i3,',',i3,',',i3)
          umax=dxu(i)*cs(j)/dtmax
          pcflu=abs(c100*uwk*gm(i,k)/umax)
          vmax=dyu(j)/dtmax
          pcflv=abs(c100*vwk*gm(i,k)/vmax)
          wtmax=dzzqz(i,k,0)/dtmax
          pcflwt=abs(c100*w (i,k)*gm(i,k)/wtmax)
          wumax=dzzvqz(i,k,0)/dtmax
          pcflwu=abs(c100*wu(i,k)*gm(i,k)/wumax)
          write(stdout,110) pcflu,umax,pcflv,vmax,pcflwt,wtmax,pcflwu,
     &                      wumax
 110      format(' U  reached ',f8.2,'% of the CFL limit (',g15.8,')'/
     *           ' V  reached ',f8.2,'% of the CFL limit (',g15.8,')'/
     *           ' W  reached ',f8.2,'% of the CFL limit (',g15.8,')'/
     *           ' WU reached ',f8.2,'% of the CFL limit (',g15.8,')')
          pcflwt=c100*(cflu+cflv+cflwt)
          pcflwu=c100*(cflu+cflv+cflwu)
          write(stdout,120) pcflwt,pcflwu
 120      format(' cflT reached ',f8.2,'% of the CFL limit'/
     *           ' cflV reached ',f8.2,'% of the CFL limit')
          call exitus('CFL')
        endif
#endif
 130  continue
c
c  Diagnose velocity which comes closest to its CFL criteria.
c
c
      if(diagts.and.eots) then
        dtmax=max(dtuv,dtts)
        vmax=dyu(j)/dtmax
        do 140 k=1,km
        do 140 i=2,imtm1
#if !defined ext_tide | !defined advtide
          uwk =u(i,k)
          vwk =v(i,k)
# else
          uwk =u(i,k)+utide(i,k)
          vwk =v(i,k)+vtide(i,k)
#endif
          wtmax=dzzqz(i,k,0)/dtmax
          wumax=dzzvqz(i,k,0)/dtmax
          umax=dxu(i)*cs(j)/dtmax
          if(abs(c100*uwk/umax).gt.cflup) then
            cflup=abs(c100*uwk/umax)
            cflum=uwk
            icflu=i
            jcflu=j
            kcflu=k
          endif
          if(abs(c100*vwk/vmax).gt.cflvp) then
            cflvp=abs(c100*vwk/vmax)
            cflvm=vwk
            icflv=i
            jcflv=j
            kcflv=k
          endif
          if(abs(c100*w(i,k)/wtmax).gt.cflwtp) then
            cflwtp=abs(c100*w(i,k)/wtmax)
            cflwtm=w(i,k)
            icflwt=i
            jcflwt=j
            kcflwt=k
          endif
          if(abs(c100*wu(i,k)/wumax).gt.cflwup) then
            cflwup=abs(c100*wu(i,k)/wumax)
            cflwum=wu(i,k)
            icflwu=i
            jcflwu=j
            kcflwu=k
          endif
 140    continue
      endif
      return
      end