#include <misc.h>
#include <params.h>


subroutine courlim (vmax2d, vmax2dt, vcour) 1,9

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Find out whether Courant limiter needs to be applied
! 
! Method: 
! 
! Author: 
! Original version:  CCM2
! 
!-----------------------------------------------------------------------
!
! $Id: courlim.F90,v 1.5.2.3 2003/12/15 18:52:50 hender Exp $
! $Author: hender $
!
!-----------------------------------------------------------------------

  use shr_kind_mod, only: r8 => shr_kind_r8
  use pmgrid
  use pspect
  use physconst, only: rga
  use time_manager, only: get_nstep, is_first_step
  use comhd
#ifdef TIMING_BARRIERS
  use mpishorthand
#endif

  implicit none

#include <comsta.h>

!
! Arguments
!
  real(r8), intent(inout) :: vmax2d(plev,plat)   ! Max. wind at each level, latitude
  real(r8), intent(inout) :: vmax2dt(plev,plat)  ! Max. truncated wind at each lvl,lat
  real(r8), intent(inout) :: vcour(plev,plat)    ! Maximum Courant number in slice
!
!--------------------------Local Variables------------------------------
!
  integer k,lat                ! Indices
  integer latarr(1)            ! Output from maxloc (needs to be array for conformability)
  integer :: nstep             ! Current timestep number

  real(r8) vcourmax            ! Max courant number in the vertical wind field
  real(r8) vmax1d(plev)        ! Sqrt of max wind speed 
  real(r8) vmax1dt(plev)       ! Sqrt of max wind speed
  real(r8) cn                  ! Estimate of truncated Courant number 
  real(r8) cnmax               ! Max. courant no. horiz. wind field
  real(r8) psurfsum            ! Summing variable - global mass
  real(r8) stqsum              ! Summing variable - global moisture
  real(r8) rmszsum             ! Summing variable - global vorticity
  real(r8) rmsdsum             ! Summing variable - global divergence
  real(r8) rmstsum             ! Summing variable - global temperature
  real(r8) stps                ! Global Mass integral
  real(r8) stqf                ! Global Moisture integral
  real(r8) rmszf               ! Global RMS Vorticity
  real(r8) rmsdf               ! Global RMS Divergence
  real(r8) rmstf               ! Global RMS Temperature
!
!-----------------------------------------------------------------------
!
#if ( defined SPMD )
#if ( defined TIMING_BARRIERS )
  call t_startf ('sync_realloc7')
  call mpibarrier (mpicom)
  call t_stopf ('sync_realloc7')
#endif
  call t_startf ('realloc7')
  call realloc7 (vmax2d, vmax2dt, vcour)
  call t_stopf ('realloc7')
#endif

  nstep = get_nstep()
!
! Compute maximum wind speed for each level
!
  do k=1,plev
     vmax1d(k)  = sqrt (maxval (vmax2d(k,:)))
     vmax1dt(k) = sqrt (maxval (vmax2dt(k,:)))
  end do
!
! Compute max. vertical Courant number (k is index to Model interfaces)
!
  vcourmax = maxval (vcour(2:,:))
!
! Determine whether the CFL limit has been exceeded for each level
! within the specified range (k<=kmxhdc). Set the truncation wave number
! (for each level independently) so that the CFL limit will not be
! violated and print a message (information only). The trunc wavenumber
! is used later in "hordif" to adjust the diffusion coefficients for
! waves beyond the limit. Store the maximum Courant number for printing
! on the stats line. Note that the max Courant number is not computed
! for the entire vertical domain, just the portion for which the limiter
! is actually applied.
!
  cnmax = 0.
  do k=1,kmxhdc
     cn = vmax1dt(k)*cnfac  ! estimate of truncated Courant number 
     cnmax = max(cnmax,cn)
     if (cn .gt. cnlim) then
        nindex(k) = int(nmaxhd*cnlim/cn + 1.)
        latarr = maxloc (vmax2dt(k,:))
        if (masterproc) write(6,800)k,latarr,cn,nindex(k)-1
     else
        nindex(k) = 2*nmaxhd
     endif
  end do
!
! Write out estimate of original Courant number if limit is exceeded
!
  do k=1,kmxhdc
     cn = vmax1d(k)*cnfac   ! estimate of original Courant number 
     if (cn .gt. cnlim) then
        latarr = maxloc (vmax2d(k,:))
        if (masterproc) write(6,805) k,latarr,cn
     end if
  end do
!
! Compute Max Courant # for whole atmosphere for diagnostic output
!
  cnmax = 0.
  do k=1,plev-1
     cn = vmax1dt(k)*cnfac  ! estimate of Courant number
     cnmax = max(cnmax,cn)
  end do
!
! Write out statisitics to standard output
!
  psurfsum = 0.
  stqsum = 0.
  rmszsum = 0.
  rmsdsum = 0.
  rmstsum = 0.

  do lat=1,plat
     psurfsum = psurfsum + psurf(lat)
     stqsum = stqsum + stq(lat)
     rmszsum = rmszsum + rmsz(lat)
     rmsdsum = rmsdsum + rmsd(lat)
     rmstsum = rmstsum + rmst(lat)
  end do

  stps = 0.5*psurfsum
  stqf = 0.5*rga*stqsum
  rmszf = sqrt(0.5*rmszsum)
  rmsdf = sqrt(0.5*rmsdsum)
  rmstf = sqrt(0.5*rmstsum)
  if (masterproc) then
     if (is_first_step()) write(6,810)
     write (6,820) nstep, rmszf, rmsdf, rmstf, stps, stqf, cnmax, vcourmax
  end if
!
  return
!
! Formats
!
800 format('COURLIM: *** Courant limit exceeded at k,lat=',2i3,  &
       ' (estimate = ',f6.3, '), solution has been truncated to ',  &
       'wavenumber ',i3,' ***')
805 format(' *** Original  Courant limit exceeded at k,lat=',2i3,  &
       ' (estimate = ',f6.3,')',' ***')
810 format(/109x,'COURANT'/10x,'NSTEP',4x,'RMSZ',19x,'RMSD',19x,  &
       'RMST',4x,'STPS',9x,'STQ',19x,'HOR   VERT')
820 format(' NSTEP =',i8,1x,1p2e23.15,0p1f8.3,1p1e13.5,e23.15,  &
       0p1f5.2,f6.2)
end subroutine courlim