#include <misc.h>
#include <params.h>
subroutine hordif(k,ztdt) 1,7
!-----------------------------------------------------------------------
!
! Purpose:
!
! Method:
! Horizontal diffusion of z,d,t,q
! 1. implicit del**2 form above level kmnhd4
! 2. implicit del**4 form at level kmnhd4 and below
! 3. courant number based truncation at level kmxhdc and above
! 4. increased del**2 coefficient at level kmxhd2 and above
!
! Computational note: this routine is multitasked by level, hence it
! is called once for each k
!
! Author:
! Original version: CCM1
! Standardized: J. Rosinski, June 1992
! Reviewed: B. Boville, J. Hack, August 1992
! Reviewed: B. Boville, April 1996
!
!-----------------------------------------------------------------------
!
! $Id: hordif.F90,v 1.3.2.7 2004/04/28 23:43:55 eaton Exp $
! $Author: eaton $
!
!-----------------------------------------------------------------------
use shr_kind_mod
, only: r8 => shr_kind_r8
use pmgrid
use pspect
use comspe
use time_manager
, only: get_step_size, is_first_step, get_nstep
use comhd
!-----------------------------------------------------------------------
implicit none
!------------------------------Commons----------------------------------
#include <comctl.h>
!------------------------------Arguments--------------------------------
!
! Input arguments
!
integer, intent(in) :: k ! level index
real(r8), intent(in) :: ztdt ! 2 times time step unless nstep=0
!
!---------------------------Local workspace-----------------------------
!
integer ir,ii ! spectral indices
integer lmr,lmc ! spectral indices
real(r8) dfac ! large coefficient on del^n multipliers to
! strongly damp waves req'd by Courant limiter
integer lm,m,n ! spectral indices
real(r8) ztodt ! 2 delta t
real(r8) zdt ! model time step
real(r8) dmpini ! used to compute divergence damp rate
real(r8) dmptim ! used to compute divergence damp rate
real(r8) dmprat ! divergence damping rate
real(r8) coef ! coeff. used to apply damping rate to divergence
real(r8) two
!
!-----------------------------------------------------------------------
!DIR$ NOSTREAM
two=2.
!
! Set the horizontal diffusion factors for each wavenumer at this level
! depending on: whether del^2 or del^4 diffusion is to be applied; and
! whether the courant number limit is to be applied.
!
if (k .ge. kmnhd4) then ! Del^4 diffusion factors
do n=1,pnmax
hdiftq(n,k) = hdfst4(n)
hdifzd(n,k) = hdfsd4(n)
end do
!
! Spectrally truncate selected levels (if courant number too large)
!
if (k.le. kmxhdc .and. nindex(k).le.pnmax) then
dfac = 1000.
do n=nindex(k),pnmax
hdiftq(n,k) = dfac*hdfst4(n)
hdifzd(n,k) = dfac*hdfsd4(n)
end do
end if
else ! Del^2 diffusion factors
if (k.le.kmxhd2) then
!
! Buggy sun compiler gives wrong answer for following line when
! using -Qoption f90comp -r8const flags
! dfac = 2.**(real(kmxhd2-k+1))
dfac = two**(real(kmxhd2-k+1))
else
dfac = 1.0
end if
do n=1,pnmax
hdiftq(n,k) = dfac*hdfst2(n)
hdifzd(n,k) = dfac*hdfsd2(n)
end do
!
! Spectrally truncate selected levels (if courant number too large)
!
if ((k.le.kmxhdc).and.(nindex(k).le.pnmax)) then
dfac = 1000.
do n=nindex(k),pnmax
hdiftq(n,k) = dfac*hdfst2(n)
hdifzd(n,k) = dfac*hdfsd2(n)
end do
end if
end if
!
! Define damping rate for divergence damper
!
zdt = get_step_size
()
ztodt = 2.*zdt
if (is_first_step()) ztodt = .5*ztodt
!
! Initial damping rate (e-folding time = zdt) and then linearly decrease
! to 0. over number of days specified by "divdampn".
!
coef = 1.
if (divdampn .gt. 0.0) then
dmpini = 1./(zdt)
dmptim = divdampn*86400.
dmprat = dmpini * (dmptim - float(get_nstep())*zdt) / dmptim
if (dmprat .gt. 0.0) coef = 1.0 / (1.0+ztodt*dmprat)
endif
!
! Compute time-split implicit factors for this level
!
do lm=1,numm(iam)
m=locm(lm,iam)
lmr = lnstart(lm)
lmc = 2*lmr
do n=1,nlen(m)
ir = lmc + 2*n - 1
ii = ir + 1
!
! time-split implicit factors
!
t(ir,k) = t(ir,k)/(1. + ztdt*hdiftq(n+m-1,k))
t(ii,k) = t(ii,k)/(1. + ztdt*hdiftq(n+m-1,k))
!
d(ir,k) = d(ir,k)*coef/(1. + ztdt*hdifzd(n+m-1,k))
d(ii,k) = d(ii,k)*coef/(1. + ztdt*hdifzd(n+m-1,k))
!
vz(ir,k) = vz(ir,k)/(1. + ztdt*hdifzd(n+m-1,k))
vz(ii,k) = vz(ii,k)/(1. + ztdt*hdifzd(n+m-1,k))
end do
end do
!
return
end subroutine hordif