#include "options.inc"


      subroutine stab
      use spflame_module
      implicit none
      real tdig,cflcrt,epsln
      integer maxcfl,i,j,k,numcfl,io
      real cl,cosur,f1,f2,f3
      real cflu,cflv,cflwu,cflwt
      real umax,vmax,wmax
      real pcflu,pcflv,pcflwu,pcflwt
      integer is,ie,js,je,ks,ke,kk


      real :: ramb,rahb,rame,ramn,rahe,rahn
      real :: reyx,reyy,reyz
      real :: reynx=0.,reyny=0.,reynz=0.
      real :: reynu,reynv,reynw
      real :: reynmu,reynmv,reynmw
      integer :: ireynx,jreynx,kreynx
      integer :: ireyny,jreyny,kreyny
      integer :: ireynz,jreynz,kreynz
      real :: pecx,pecy,pecz
      real :: peclx=0.,pecly=0.,peclz=0.
      real :: peclu,peclv,peclw
      real :: peclmu,peclmv,peclmw
      integer :: ipeclx,jpeclx,kpeclx
      integer :: ipecly,jpecly,kpecly
      integer :: ipeclz,jpeclz,kpeclz

c-----------------------------------------------------------------------
c     test for various measures of stability  
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c     SPFLAME version : c.eden
c-----------------------------------------------------------------------
c

      tdig     = 1.e-4
      cflcrt   = 1.5
      maxcfl   = 3
      epsln=1e-12
      io=6

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c      Checking vertical diffusive criteria
c-----------------------------------------------------------------------
c
      do j=js,je
       do k=1,km
        do i=is,ie
         if (.not. enable_implicit_vert_fric) then
          if ((dt*visc_cbu(i,k,j))/dzt(k)**2 >= .25) 
     &     print*,' vertical diffusive criterium for momentum ',
     &      'exceeded at (i,j,k)=',i,j,k,' visc_cbu(i,k,j)*dt/dz^2=',
     &      visc_cbu(i,k,j)*dt/dzt(k)**2,' use implicit vertical mixing'
         endif
         if (.not. enable_implicit_vert_diff) then
          if ((dt*diff_cbt(i,k,j))/dzt(k)**2 >= .25) 
     &     print*,' vertical diffusive criterium for tracer ',
     &      'exceeded at (i,j,k)=',i,j,k,' diff_cbt(i,k,j)*dt/dz^2=',
     &      diff_cbt(i,k,j)*dt/dzt(k)**2,' use implicit vertical mixing'
         endif
        enddo
       enddo
      enddo
c
c-----------------------------------------------------------------------
c     scan for CFL violations. save locations of closest approach to
c     local CFL limit.
c-----------------------------------------------------------------------
c
      do j=js,je
       cl = cflcrt * 0.5
       cosur = max(csur(j),epsln)
       do k=1,km
        f1    = dt*dyur(j)
        f2    = dt*dzwr(k)
	f3    = dt*cosur
        do i=is,ie
          cflu  = abs(f3*dxur(i)*u(i,k,j,1,tau))
          cflv  = abs(f1*u(i,k,j,2,tau))
          cflwu = abs(f2*adv_vbu(i,k,j))*umask(i,k,j)
          cflwt = abs(f2*adv_vbt(i,k,j))*tmask(i,k,j)
          if (cflu .ge. cl .or. cflv .ge. cl .or. cflwu .ge. cl .or.
     &        cflwt .ge. cl) then
c          if (cflu .ge. cl .or. cflv .ge. cl) then
            umax  = 0.5*csu(j)*dxu(i)/dt
            pcflu = abs(100.0*u(i,k,j,1,tau)/umax)
            vmax  = 0.5*dyu(j)/dt
            pcflv = abs(100.0*u(i,k,j,2,tau)/vmax)
            wmax  = 0.5*dzw(k)/dt
            pcflwu= abs(100.0*adv_vbu(i,k,j)/wmax)
            pcflwt= abs(100.0*adv_vbt(i,k,j)/wmax)
            write (io,
     &  '(/,a,i4,a1,i3,a1,i3,a13,f6.3,/,4(/,a,f8.2,a,g15.8,a,/))')
     &       ' ==> CFL exceeded at coordinate (i,j,k) = (',i,
     &      ',',j,',',k,') by factor =',cflcrt,
     &       ' u reached   ', pcflu,' % of the CFL limit (',umax,')',
     &       ' v reached   ', pcflv,' % of the CFL limit (',vmax,')',
     &     ' adv_vbu reached', pcflwu,' % of the CFL limit (',wmax,')',
     &     ' adv_vbt reached', pcflwt,' % of the CFL limit (',wmax,')'
          endif
        enddo
      enddo

#ifdef notdef
      do k=1,km
        vmax  = 0.5*dyu(j)/dt
        wmax  = 0.5*dzw(k)/dt
        do i=is,ie
          umax  = 0.5*csu(j)*dxu(i)/dt
          if (abs(100.0*u(i,k,j,1,tau)/umax) .gt. cflup) then
            cflup = abs(100.0*u(i,k,j,1,tau)/umax)
            cflum = u(i,k,j,1,tau)
            icflu = i
            jcflu = j
            kcflu = k
          endif
          if (abs(100.0*u(i,k,j,2,tau)/vmax) .gt. cflvp) then
            cflvp = abs(100.0*u(i,k,j,2,tau)/vmax)
            cflvm = u(i,k,j,2,tau)
            icflv = i
            jcflv = j
            kcflv = k
          endif
          if (abs(100.0*umask(i,k,j)*adv_vbu(i,k,j)/wmax) .gt. cflwup)
     &      then
            cflwup = abs(100.0*adv_vbu(i,k,j)/wmax)
            cflwum = adv_vbu(i,k,j)
            icflwu = i
            jcflwu = j
            kcflwu = k
          endif
          if (abs(100.0*tmask(i,k,j)*adv_vbt(i,k,j)/wmax) .gt. cflwtp)
     &      then
            cflwtp = abs(100.0*adv_vbt(i,k,j)/wmax)
            cflwtm = adv_vbt(i,k,j)
            icflwt = i
            jcflwt = j
            kcflwt = k
          endif
       enddo
      enddo
#endif
      enddo

c
c-----------------------------------------------------------------------
c     look for max peclet numbers using velocities at "tau".
c     look for max reynolds numbers using velocities at "tau".
c-----------------------------------------------------------------------
c

      do j=js,je
       do k=1,km
        do i=is,ie
	  ramb = 1./(visc_cbu(i,k,j) + epsln)
	  rahb = 1./(diff_cbt(i,k,j) + epsln)
          rame = 1./(am*am_scale(i,j)+ epsln)
          ramn = 1./(am*am_scale(i,j)+ epsln)
          rahe = 1./(ah*ah_scale(i,j)+ epsln)
          rahn = 1./(ah*ah_scale(i,j)+ epsln)

          reyx = abs(u(i,k,j,1,tau)*dxu(i))*rame
          if (reyx .gt. reynx) then
            ireynx = i ; jreynx = j ; kreynx = k
            reynx  = reyx
            reynu  = u(i,k,j,1,tau)
            reynmu = 1./rame
          endif
          reyy = abs(u(i,k,j,2,tau)*dyu(j))*ramn
          if (reyy .gt. reyny) then
            ireyny = i ; jreyny = j ; kreyny = k
            reyny  = reyy
            reynv  = u(i,k,j,2,tau)
            reynmv = 1./ramn
          endif
	  kk = min(k+1,km)
	  if (k .ge. kmu(i,j)) then
	    reyz = 0.
	  else 
            reyz =umask(i,kk,j)*abs(adv_vbu(i,k,j)*dzw(k))*ramb
          endif
	  if (reyz .gt. reynz) then
            ireynz = i ; jreynz = j ; kreynz = k
            reynz  = reyz
            reynw  = adv_vbu(i,k,j)
            reynmw = 1./ramb
          endif


          pecx = abs(u(i,k,j,1,tau)*dxu(i))*rahe
          if (pecx .gt. peclx) then
            ipeclx = i ; jpeclx = j ; kpeclx = k
            peclx  = pecx
            peclu  = u(i,k,j,1,tau)
            peclmu = 1./rahe
          endif
          pecy = abs(u(i,k,j,2,tau)*dyu(j))*rahn
          if (pecy .gt. pecly) then
            ipecly = i ; jpecly = j ; kpecly = k
            pecly  = pecy
            peclv  = u(i,k,j,2,tau)
            peclmv = 1./rahn
          endif
	  kk = min(k+1,km)
	  if (k .ge. kmt(i,j)) then
	    pecz = 0.0
	  else 
            pecz =tmask(i,kk,j)*abs(adv_vbt(i,k,j)*dzw(k))*rahb
          endif
	  if (pecz .gt. peclz) then
            ipeclz = i ; jpeclz = j ; kpeclz = k
            peclz  = pecz
            peclw  = adv_vbt(i,k,j)
            peclmw = 1./rahb
          endif
        enddo
       enddo
      enddo

#ifdef notdef
      print*,' max x-Pecl.# for PE=',my_pe,' at (i,j,k)=(',
     &       ipeclx,jpeclx,kpeclx,') = ',peclx
      print*,' max y-Pecl.# for PE=',my_pe,' at (i,j,k)=(',
     &       ipecly,jpecly,kpecly,') = ',pecly
      print*,' max z-Pecl.# for PE=',my_pe,' at (i,j,k)=(',
     &       ipeclz,jpeclz,kpeclz,') = ',peclz

      print*,' max x-Reyn.# for PE=',my_pe,' at (i,j,k)=(',
     &       ireynx,jreynx,kreynx,') = ',reynx
      print*,' max y-Reyn.# for PE=',my_pe,' at (i,j,k)=(',
     &       ireyny,jreyny,kreyny,') = ',reyny
      print*,' max z-Reyn.# for PE=',my_pe,' at (i,j,k)=(',
     &       ireynz,jreynz,kreynz,') = ',reynz
#endif

#ifdef notdef
c
c-----------------------------------------------------------------------
c     look for ficticious creation of local extremum for tracers
c     by finding local min and max tracer at "tau" and comparing to
c     tracer at "tau+1"
c-----------------------------------------------------------------------
c
      call getunit (iostab, 'iostab', 'fsa')
      ks = max(2,kscfl)
      ke = min(km-1,kecfl)
      is = max(2,iscfl)
      ie = min(imt-1,iecfl)
      do n=1,nt
        do k=ks,ke
          do i=is,ie
	    if (tmask(i,k,j) .ne. c0) then
	      tbig = max(t(i,k,j,n,tau),t(i,k,j,n,taum1))
	      tsml = min(t(i,k,j,n,tau),t(i,k,j,n,taum1))
	      do jj=j-1,j+1,2
	        if (tmask(i,k,jj) .ne. c0) then
		  tbig = max(tbig,t(i,k,jj,n,tau),t(i,k,jj,n,taum1))
		  tsml = min(tsml,t(i,k,jj,n,tau),t(i,k,jj,n,taum1))
		endif
	      enddo
c
	      do ii=i-1,i+1,2
	        if (tmask(ii,k,j) .ne. c0) then
		  tbig = max(tbig,t(ii,k,j,n,tau),t(ii,k,j,n,taum1))
		  tsml = min(tsml,t(ii,k,j,n,tau),t(ii,k,j,n,taum1))
		endif
	      enddo
c
	      do kk=k-1,k+1,2
	        if (tmask(i,kk,j) .ne. c0) then
		  tbig = max(tbig,t(i,kk,j,n,tau),t(i,kk,j,n,taum1))
		  tsml = min(tsml,t(i,kk,j,n,tau),t(i,kk,j,n,taum1))
		endif
	      enddo
c
              tcrit = tdig*abs(t(i,k,j,n,taup1))
              if (tmask(i,k,j) .ne. c0 .and. 
     &          ((t(i,k,j,n,taup1) .gt. tbig + tcrit)
     &	        .or. (t(i,k,j,n,taup1) .lt. tsml - tcrit))) then
		write (iostab,'(i4, i4, i4, i2, 3g14.7)')
     &           i, k, j, n, t(i,k,j,n,taup1), tsml, tbig
              endif
	    endif
	  enddo
	enddo
      enddo
      close (iostab)
#endif
c      
9100  format(1x,a12,1x,'ts=',i10,1x,',j=',i3,', lat=',f6.2
     &,', lon:',f6.2,' ==> ',f6.2,', depth(m):',f6.1,' ==> ',f6.1
     &,', scaling=',1pg10.3)
      end subroutine stab

