#include "options.inc"


      subroutine checks()
      use spflame_module
      implicit none
      integer :: is,ie,js,je,i,j,k
      real :: dxsq,dysq,h1,h2,clix,hx,dfdy
c
c-----------------------------------------------------------------------
c     some checks
c     SPFLAME version:   c.eden
c-----------------------------------------------------------------------
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      call barrier
      call sub_flush(6)

c
c     check consistency of the vertical mixing setup
c
      if (.not. enable_const_vert_mixing .and.
     &    .not. enable_cgh_vert_mixing   .and.
     &    .not. enable_tkemix ) then
        if (my_pe==0) 
     &   print*,' you have to specify one of the vertical mixing',
     &   ' schemes in the namelist input'
        call sub_flush(6)
        call halt_stop(' in checks')
      endif

      if (enable_const_vert_mixing .and.
     &    enable_cgh_vert_mixing) then
       if (my_pe==0) then
        print*,' specify either constant or CGH vertical mixing'
        print*,' not both '
       endif
       call sub_flush(6)
       call halt_stop(' in checks')
      endif

      if (enable_impl_convection .and.
     &    enable_expl_convection) then
       if (my_pe==0) then
        print*,' specify either implicit of exlicit convection'
        print*,' not both '
       endif
       call sub_flush(6)
       call halt_stop(' in checks')
      endif

      if (enable_impl_convection .and.
     &    .not. enable_const_vert_mixing) then
       if (my_pe==0) then
        print*,' implicit mixing depends on constant parameters'
        print*,' switch on enable_const_vert_mixing  '
       endif
       call sub_flush(6)
       call halt_stop(' in checks')
      endif

      if (enable_const_vert_mixing .and.
     &    enable_tkemix) then
       if (my_pe==0) then
        print*,' specify either constant or TKE vertical mixing'
        print*,' not both '
       endif
       call sub_flush(6)
       call halt_stop(' in checks')
      endif

      if (enable_cgh_vert_mixing .and.
     &    enable_tkemix) then
       if (my_pe==0) then
        print*,' specify either CGH or TKE vertical mixing'
        print*,' not both '
       endif
       call sub_flush(6)
       call halt_stop(' in checks')
      endif

      if (enable_ktmix .and. .not. enable_expl_convection) then
       if (my_pe==0) then
        print*,'Kraus Turner wind induced stirring depends on a',
     & ' stable stratification and thus on explicit convection'
        print*,' enable explicit convection scheme with KTmix '
       endif
       call sub_flush(6)
       call halt_stop(' in checks')
      endif
c
c     check consistency of the horizontal mixing setup
c
      if (enable_diffusion_isoneutral .and.
     &    enable_diffusion_isopycnic)  then
       if (my_pe==0) then
        print*,' specify either isoneutral mixing or isopynic mixing'
        print*,' not both '
       endif
       call sub_flush(6)
       call halt_stop(' in checks')
      endif

      if ((enable_diffusion_isoneutral .or.
     &     enable_diffusion_isopycnic)
     &     .and. .not. enable_implicit_vert_diff)  then
       if (my_pe==0) then
        print*,' use implicit vertical mixing with isoneutral mixing '
        print*,' or isopynic mixing (option enable_implicit_vert_diff) '
       endif
       call sub_flush(6)
       call halt_stop(' in checksp')
      endif

#ifdef partial_cell

      if (enable_diffusion_isopycnic) then
       if (my_pe==0) then
        print*,' partial cells with isopycnic mixing '
        print*,' not yet implemented, use isoneutral mixing'
       endif
       call sub_flush(6)
       call halt_stop(' in checks')
      endif

c      if (enable_bbl) then
c       if (my_pe==0) then
c        print*,' partial cells with Beckmann/Doescher BBL '
c        print*,' not yet implemented'
c       endif
c       call sub_flush(6)
c       call halt_stop(' in checks')
c      endif

#endif
c
c     Check some numerical criteria
c
      if (my_pe==0) print*,' Checking Killworths topographic critrium'

      do j=js,je
        do i=is,ie
         dysq=dyt(j)**2

         if ((kmt_big(i+1,j-1)/=0).and.(kmt_big(i+1,j)/=0) ) then 

          dxsq = (dxt(i)*cst(j))**2
          clix = am*dt/dxsq
          h1   = zt(kmt_big(i+1,j-1))
          h2   = zt(kmt_big(i+1,j) )
          hx   = (8.0*h1*h2/(h1+h2)**2+dxsq/dysq)/(4.0 + dxsq/dysq)

          if (clix .ge. hx .and. enable_friction_harmonic ) then
           print*,' criterium exceeded at i=',i,' j=',j,
     &             'kmt(i,j)=',kmt_big(i,j),' h=',zt(kmt_big(i,j)),' cm'
           print*,' clix=',clix,' hx=',hx
           call sub_flush(6)
           call halt_stop(' in checks')
          endif
         endif
        enddo
      enddo

      if (my_pe==0) print*,' ok'

      call barrier; call sub_flush(6)

c       what happens for biha ??
c
      if (my_pe==0) print*,' Checking horizontal diffusive criteria'

      do j=js,je
        do i=is,ie
         hx=(ah+1.e-5)  ! if ah is zero
     &                         *ah_scale(i,j)
         if (enable_diffusion_harmonic .and. 
     &         ( dt*tmask(i,1,j) > (dxt(i)*cst(j))**2/(2.*hx) ) ) then
           print*,' criterium exceeded for T grid at i=',i,' j=',j
           print*,' dt = ',dt
           print*,' dx^2/2/ah=', (dxt(i)*cst(j))**2/(2.*hx) 
           call sub_flush(6)
           call halt_stop(' in checks')
         endif
         hx=(am+1.e-5)*am_scale(i,j)
         if (enable_friction_harmonic .and. 
     &         ( dt*umask(i,1,j) > (dxu(i)*csu(j))**2/(2.*hx) ) ) then
           print*,' criterium exceeded for U grid at i=',i,' j=',j
           print*,' dt = ',dt
           print*,' dx^2/2/am=', (dxu(i)*csu(j))**2/(2.*hx) ,hx
           call sub_flush(6)
           call halt_stop(' in chekcs')
         endif
        enddo
      enddo
      if (my_pe==0) print*,' ok'

      call barrier; call sub_flush(6)


#ifdef notdef
      if (enable_friction_harmonic .and. 
     &   .not. enable_friction_biharmonic .and. 
     &   .not. enable_rotated_grid) then
      if (my_pe==0) print*,' Checking Munck criteria'
c
c   what happens if am=0?? and for biha?
c   and for a rotated grid ?
c
      do j=js,je
        do i=is,ie
         hx=(am+1e-5)*am_scale(i,j)
         dfdy   = (2.*omega*csu(j)/radius+1e-15)
         if  ( (hx/dfdy)**(1./3.) <= 0.5*dxu(i)*csu(j) ) then
          print*,' criterium exceeded for dxu at i=',i,' j=',j
          print*,' viscous boundary layer: ',  (hx/dfdy)**(1./3.) 
          print*,' dxu = ', 0.5*dxu(i)*csu(j) 
           call sub_flush(6)
          call halt_stop(' in checks'   )
         endif
         if  ( (hx/dfdy)**(1./3.) <= 0.5*dyu(j) ) then
          print*,' criterium exceeded for dyu at i=',i,' j=',j
          print*,' viscous boundary layer: ',  (hx/dfdy)**(1./3.) 
          print*,' dyu = ', 0.5*dyu(j) 
           call sub_flush(6)
          call halt_stop(' in checks'   )
         endif
c         print*, i,j, (hx/dfdy)**(1./3.) , 0.5*dyu(j) ,hx,dfdy
        enddo
      enddo
      if (my_pe==0) print*,' ok'
      elseif (my_pe==0) then
         print*,' do not checking Munck criteria '
      endif
#endif

      call barrier; call sub_flush(6)

      if (my_pe==0) then 
       print*,' Checking vertical diffusive criteria'

       do k=1,km

        if (.not. enable_implicit_vert_fric) then
         hx=kappa_m
         if (enable_cgh_vert_mixing .and. 
     &        enable_cgh_vert_momentum_mixing ) hx=max(hx,visc_cbu_cut)
         if (enable_tkemix)                      hx=max(hx,visc_cbu_cut)
         if ((dt*hx)/dzt(k)**2 >= .25) then
          print*,' momentum criterium exceeded at k=',k
          print*,' use option enable_implicit_vert_fric'
          call halt_stop(' in checks'   )
         endif
        endif

        if (.not. enable_implicit_vert_diff) then
         hx=kappa_h
         if (enable_cgh_vert_mixing ) hx=max(hx,diff_cbt_cut)
         if (enable_impl_convection) hx=max(hx,diff_cbt_cut)
         if (enable_cgh_vert_mixing .and. 
     &       enable_cgh_impl_convection ) hx=max(hx,1.0e6)
         if (enable_tkemix)       hx = max(hx,diff_cbt_cut)
         if ((dt*hx)/dzt(k)**2 >= .25) then
          print*,' tracer criterium exceeded at k=',k
          print*,' use option enable_implicit_vert_diff'
          call halt_stop(' in checks'   )
         endif
        endif

       enddo
       print*,' ok'
      endif

      call barrier; call sub_flush(6)

      if (enable_obc_south.or.enable_obc_north.or.
     &    enable_obc_west.or. enable_obc_east) then
       if (my_pe==0) then 
        print*,' Checking vertical diffusive criteria for obc'

        do k=1,km
         hx=kappa_h
         if ((dt*hx)/dzt(k)**2 >= .25) then
          print*,' tracer criterium exceeded at k=',k
          print*,' change kappa_h or time step'
          call halt_stop(' in checks'   )
         endif
        enddo

        if (kappa_h == 0.) then
          print*,' kappa_h is zero for open boundaries'
          print*,' this will not work'
          call halt_stop(' in checks'   )
        endif

        print*,' ok'
       endif
      endif
c
      print*,' pe # ',my_pe,' checks done '
      call sub_flush(6)

      end subroutine checks
