subroutine gfidi_sig(lon_dim,lons_lat,lat,
1 dg,tg,zg,ug,vg,rqg,dphi,dlam,qg,
1 rcl,del,rdel2,ci,tov,spdmax,deltim,sl,nvcn,xvcn,
2 dtdf,dtdl,drdf,drdl,dudl,dvdl,dudf,dvdf,
2 dqdt,dtdt,drdt,dudt,dvdt)
cc
USE MACHINE , ONLY : kind_phys
cc
use resol_def
use physcons, rerth => con_rerth, rd => con_rd, cp => con_cp
&, omega => con_omega
implicit none
cc
integer j,k,n,nvcn
real(kind=kind_evod) fnor,rcl,rcl2,rk,sinra,deltim,xvcn
c
c input variables
c
integer lon_dim,lat
integer lons_lat
c
real(kind=kind_evod)
1 dg(lon_dim,levs),tg(lon_dim,levs), zg(lon_dim,levs),
2 ug(lon_dim,levs),vg(lon_dim,levs),rqg(lon_dim,levs,ntrac),
3 dphi(lon_dim),dlam(lon_dim),qg(lon_dim)
c
real(kind=kind_evod)
1 dtdf(lon_dim,levs),dtdl(lon_dim,levs),
1 drdf(lon_dim,levs,ntrac),drdl(lon_dim,levs,ntrac),
1 dudl(lon_dim,levs),dvdl(lon_dim,levs),
1 dudf(lon_dim,levs),dvdf(lon_dim,levs)
c
c output variables
c
real(kind=kind_evod) spdmax(levs),
1 dudt(lon_dim,levs),dvdt(lon_dim,levs),
1 dtdt(lon_dim,levs),drdt(lon_dim,levs,ntrac),
1 dqdt(lon_dim)
c
c constant arrays
c
real(kind=kind_evod)
1 sl(levs),
1 del(levs),rdel2(levs),
2 ci(levp1),tov(levs)
c
c local variables
c
real(kind=kind_evod)
1 cg (lon_dim,levs), db(lon_dim,levs),cb(lon_dim,levs),
2 dot(lon_dim,levp1),dup(lon_dim,levs),dvp(lon_dim,levs),
3 dum(lon_dim,levs ),dvm(lon_dim,levs), ek(lon_dim,levs),
4 rmu(levs ),rnu(levs),rho(levs),si(levp1),
5 x1(levs ), x2(levs), x3(levs),x4(levs)
c
real(kind=kind_evod) cons0,cons0p5,cons1,cons2 !constant
c
cons0 = 0.d0 !constant
cons0p5 = 0.5d0 !constant
cons1 = 1.d0 !constant
cons2 = 2.d0 !constant
c
rk= rd /cp
sinra=sqrt(cons1-cons1/rcl) !constant
fnor=cons2*omega*sinra !constant
sinra=sinra/rerth
if(lat.gt.latg2) then
fnor=-fnor
sinra=-sinra
endif
c
si(1)=cons1 !constant
do 4 k=1,levs
si(k+1)=si(k)-del(k)
4 continue
c
do 1 k=1,levm1
rho(k)= log(si(k)/si(k+1))
1 continue
rho(levs)=cons0 !constant
c
do 2 k=1,levs
rmu(k)=cons1-si(k+1)*rho(k)/del(k) !constant
2 continue
c
do 3 k=1,levm1
rnu(k+1)=-cons1+si(k)*rho(k)/del(k) !constant
3 continue
rnu(1)=cons0 !constant
c
do 20 k=1,levs
x1(k)=rmu(k)*(cons1-rk*rnu(k))/(rmu(k)+rnu(k)) !constant
x2(k)=cons1-x1(k) !constant
x3(k)=(cons1+rk*rmu(k))/(cons1-rk*rnu(k)) !constant
x4(k)=cons1/x3(k) !constant
20 continue
c
do 1234 k=1,levs
spdmax(k)=cons0 !constant
1234 continue
c
rcl2=cons0p5*rcl !constant
c
do 140 k=1,levs
do 140 j=1,lons_lat
ek(j,k)=(ug(j,k)*ug(j,k)+vg(j,k)*vg(j,k))*rcl
140 continue
c
do 10 k=1,levs
do 10 j=1,lons_lat
if (ek(j,k) .gt. spdmax(k)) spdmax(k)=ek(j,k)
10 continue
c
c compute c=v(true)*del(ln(ps)).divide by cos for del, cos for v
c
do 150 j=1,lons_lat
dphi(j)=dphi(j)*rcl
dlam(j)=dlam(j)*rcl
150 continue
c
do 180 k=1,levs
do 180 j=1,lons_lat
cg(j,k)=ug(j,k)*dlam(j)+vg(j,k)*dphi(j)
180 continue
c
do 190 j=1,lons_lat
db(j,1)=del(1)*dg(j,1)
cb(j,1)=del(1)*cg(j,1)
190 continue
c
do 210 k=1,levm1
do 210 j=1,lons_lat
db(j,k+1)=db(j,k)+del(k+1)*dg(j,k+1)
cb(j,k+1)=cb(j,k)+del(k+1)*cg(j,k+1)
210 continue
c
c store integral of cg in dlax
c
do 220 j=1,lons_lat
dqdt(j)= -cb(j,levs)
220 continue
c
c sigma dot computed only at interior interfaces.
c
do 230 j=1,lons_lat
dot(j,1)=cons0 !constant
dvm(j,1)=cons0 !constant
dum(j,1)=cons0 !constant
dot(j,levp1)=cons0 !constant
dvp(j,levs )=cons0 !constant
dup(j,levs )=cons0 !constant
c
230 continue
c
do 240 k=1,levm1
do 240 j=1,lons_lat
dot(j,k+1)=dot(j,k)+
1 del(k)*(db(j,levs)+cb(j,levs)-
2 dg(j,k)-cg(j,k))
240 continue
c
c
c implicitly filter input profiles
c if vertical advection would be numerically unstable
c
call vcnfil(lons_lat,lon_dim,levs,ntrac,
x deltim,del,sl,si,tov,dot(1,1),
x ug(1,1),vg(1,1),tg(1,1),rqg(1,1,1),nvcn,xvcn)
c
do 260 k=1,levm1
do 260 j=1,lons_lat
dvp(j,k )=vg(j,k+1)-vg(j,k)
dup(j,k )=ug(j,k+1)-ug(j,k)
dvm(j,k+1)=vg(j,k+1)-vg(j,k)
dum(j,k+1)=ug(j,k+1)-ug(j,k)
260 continue
c
do j=1,lons_lat
dphi(j)=dphi(j)/rcl
dlam(j)=dlam(j)/rcl
enddo
c
do k=1,levs
do j=1,lons_lat
dudt(j,k)=-ug(j,k)*dudl(j,k)-vg(j,k)*dudf(j,k)
1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k))
2 -rd*tg(j,k)*dlam(j)
c
dvdt(j,k)=-ug(j,k)*dvdl(j,k)-vg(j,k)*dvdf(j,k)
1 -rdel2(k)*(dot(j,k+1)*dvp(j,k)+dot(j,k)*dvm(j,k))
2 -rd*tg(j,k)*dphi(j)
enddo
enddo
c
do k=1,levs
do j=1,lons_lat
dudt(j,k)=dudt(j,k)+vg(j,k)*fnor
c
dvdt(j,k)=dvdt(j,k)-ug(j,k)*fnor-sinra*ek(j,k)
enddo
enddo
c
do k=1,levs
do j=1,lons_lat
dudt(j,k)=dudt(j,k)*rcl
dvdt(j,k)=dvdt(j,k)*rcl
enddo
enddo
c
c
do 280 k=1,levm1
do 280 j=1,lons_lat
cecmwf:
dup(j,k )=tg(j,k+1)+tov(k+1)-tg(j,k)-tov(k)
x +cons2*rk*rnu(k+1)*(tg(j,k)+tov(k)) !constant
cecmwf:
dum(j,k+1)=tg(j,k+1)+tov(k+1)-tg(j,k)-tov(k)
x +cons2*rk*rmu(k+1)*(tg(j,k+1)+tov(k+1)) !constant
cecmwf:
cecmwf:
280 continue
c
c
do k=1,levs
do j=1,lons_lat
c dtdt(j,k)=
dtdt(j,k)=-ug(j,k)*dtdl(j,k)-vg(j,k)*dtdf(j,k)
1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k))
enddo
enddo
c
do k=1,levs
do j=1,lons_lat
dtdt(j,k)=dtdt(j,k)
1 +rk*(tov(k)+tg(j,k))*(cg(j,k)-cb(j,levs)-db(j,levs))
enddo
enddo
c
do 330 n=1,ntrac
do 300 k=1,levm1
do 300 j=1,lons_lat
dup(j,k )=rqg(j,k+1,n)-rqg(j,k,n)
dum(j,k+1)=rqg(j,k+1,n)-rqg(j,k,n)
300 continue
c
do 310 j=1,lons_lat
dup(j,levs)=cons0 !constant
310 continue
c
do 320 k=1,levs
do 320 j=1,lons_lat
drdt(j,k,n)=-ug(j,k)*drdl(j,k,n)-vg(j,k)*drdf(j,k,n)
1 -rdel2(k)*(dot(j,k+1)*dup(j,k)+dot(j,k)*dum(j,k))
320 continue
330 continue
c
return
end
c-----------------------------------------------------------------------
subroutine vcnfil(im,ix,km,nt,deltim,del,sl,si,tov,dot,
& u,v,t,q,nvcn,xvcn)
c$$$ subprogram documentation block
c . . . .
c subprogram: vcnfil vertical advection instability filter
c prgmmr: iredell org: w/nmc23 date: 91-05-07
c
c abstract: prefilters fields appearing in the nonlinear terms
c in the dynamics tendency equation in order to ensure stability
c when the vertical velocity exceeds the cfl criterion.
c the vertical velocity in this case is sigmadot.
c for simple second-order centered eulerian advection,
c filtering is needed when vcn=sigmadot*deltim/dsigma>1.
c the maximum eigenvalue of the linear advection equation
c with second-order implicit filtering on the tendencies
c is less than one for all resolvable wavenumbers (i.e. stable)
c if the nondimensional filter parameter is nu=(vcn**2-1)/4.
c
c program history log:
c 97-07-30 iredell
c 98-02-20 iredell increase margin of error in filter
c by starting filter when vcn>0.9.
c
c usage: call vcnfil(im,ix,km,nt,deltim,del,sl,si,tov,dot,
c & u,v,t,q,nvcn,xvcn)
c
c input argument list:
c im - integer number of gridpoints to filter
c ix - integer first dimension of dot,u,v,t,q
c km - integer number of vertical levels
c nt - integer number of tracers in q
c deltim - real timestep in seconds
c del - real (km) sigma thicknesses
c sl - real (km) full sigma values
c si - real (km+1) interface sigma values
c tov - real (km) temperature base
c dot - real (ix,km+1) sigmadot in 1/seconds
c u - real (ix,km) zonal wind
c v - real (ix,km) meridional wind
c t - real (ix,km) virtual temperature deviation
c q - real (ix,km,nt) tracers
c
c output argument list:
c u - real (ix,km) zonal wind
c v - real (ix,km) meridional wind
c t - real (ix,km) virtual temperature deviation
c q - real (ix,km,nt) tracers
c nvcn - integer number of points requiring filtering
c xvcn - real maximum vertical courant number
c
c subprograms called:
c tridim - tridiagonal matrix solver
c
c$$$
cfpp$ noconcur r
cc
USE MACHINE , ONLY : kind_phys
use resol_def
use physcons, rerth => con_rerth, rocp => con_rocp
implicit none
cc
integer i,im,ix,j,k,km,n,nt,nvcn
real(kind=kind_evod) cvcn,deli,deltim,rnu,t0term,t1term
! real(kind=kind_evod) cvcn,deli,deltim,rnu,rocp,t0term,t1term
real(kind=kind_evod) tterm,xvcn
cc
! parameter(rocp=rd/cp)
parameter(cvcn=0.9d0) !constant
real(kind=kind_evod) del(km),sl(km),si(km+1),tov(km)
real(kind=kind_evod) dot(ix,km+1)
real(kind=kind_evod) u(ix,km),v(ix,km),t(ix,km),q(ix,km,nt)
logical lvcn(im)
integer ivcn(im)
real(kind=kind_evod) vcn(ix,km-1)
real(kind=kind_evod) cm(ix,km),cu(ix,km-1),cl(ix,km-1)
real(kind=kind_evod) rr(ix,km,3+nt)
cc
real(kind=kind_evod) cons0,cons0p5,cons1,cons4 !constant
cc
cons0 = 0.d0 !constant
cons0p5 = 0.5d0 !constant
cons1 = 1.d0 !constant
cons4 = 4.d0 !constant
cc
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c compute vertical courant number
nvcn=0
xvcn=cons0 !constant
do i=1,im
lvcn(i)=.false.
enddo
do k=1,km-1
do i=1,im
deli=sl(k)-sl(k+1)
vcn(i,k)=abs(dot(i,k+1))*deltim/deli
lvcn(i)=lvcn(i).or.vcn(i,k).gt.cvcn
xvcn=max(xvcn,vcn(i,k))
enddo
enddo
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c determine points requiring filtering
if(xvcn.gt.cvcn) then
do i=1,im
if(lvcn(i)) then
ivcn(nvcn+1)=i
nvcn=nvcn+1
endif
enddo
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c compute tridiagonal matrix
do j=1,nvcn
cm(j,1)=cons1 !constant
enddo
do k=1,km-1
deli=sl(k)-sl(k+1)
do j=1,nvcn
i=ivcn(j)
if(vcn(i,k).gt.cvcn) then
rnu=(vcn(i,k)**2-cvcn**2)/cons4 !constant
cu(j,k)=-rnu*deli/del(k)
cl(j,k)=-rnu*deli/del(k+1)
cm(j,k)=cm(j,k)-cu(j,k)
cm(j,k+1)=cons1-cl(j,k) !constant
else
cu(j,k)=cons0 !constant
cl(j,k)=cons0 !constant
cm(j,k+1)=cons1 !constant
endif
enddo
enddo
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c fill fields to be filtered
do k=1,km
do j=1,nvcn
i=ivcn(j)
rr(j,k,1)=u(i,k)
rr(j,k,2)=v(i,k)
rr(j,k,3)=t(i,k)+tov(k)
enddo
do n=1,nt
do j=1,nvcn
i=ivcn(j)
rr(j,k,3+n)=q(i,k,n)
enddo
enddo
enddo
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c adjust temperature for adiabatic change
do k=1,km-1
deli=sl(k)-sl(k+1)
t1term=cons0p5*rocp*deli/si(k+1) !constant
t0term=t1term*(tov(k)+tov(k+1))
do j=1,nvcn
i=ivcn(j)
tterm=t0term+t1term*(t(i,k)+t(i,k+1))
rr(j,k,3)=rr(j,k,3)-cu(j,k)*tterm
rr(j,k+1,3)=rr(j,k+1,3)+cl(j,k)*tterm
enddo
enddo
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c solve tridiagonal system
call tridim(nvcn,ix,km,km,3+nt,cl,cm,cu,rr,cu,rr)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c replace filtered fields
do k=1,km
do j=1,nvcn
i=ivcn(j)
u(i,k)=rr(j,k,1)
v(i,k)=rr(j,k,2)
t(i,k)=rr(j,k,3)-tov(k)
enddo
do n=1,nt
do j=1,nvcn
i=ivcn(j)
q(i,k,n)=rr(j,k,3+n)
enddo
enddo
enddo
endif
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end
c-----------------------------------------------------------------------
subroutine tridim(l,lx,n,nx,m,cl,cm,cu,r,au,a)
c$$$ subprogram documentation block
c . . . .
c subprogram: tridim solves tridiagonal matrix problems.
c prgmmr: iredell org: w/nmc23 date: 91-05-07
c
c abstract: this routine solves multiple tridiagonal matrix problems
c with multiple right-hand-side and solution vectors for every matrix.
c the solutions are found by eliminating off-diagonal coefficients,
c marching first foreward then backward along the matrix diagonal.
c the computations are vectorized around the number of matrices.
c no checks are made for zeroes on the diagonal or singularity.
c
c program history log:
c 97-07-30 iredell
c
c usage: call tridim(l,lx,n,nx,m,cl,cm,cu,r,au,a)
c
c input argument list:
c l - integer number of tridiagonal matrices
c lx - integer first dimension (lx>=l)
c n - integer order of the matrices
c nx - integer second dimension (nx>=n)
c m - integer number of vectors for every matrix
c cl - real (lx,2:n) lower diagonal matrix elements
c cm - real (lx,n) main diagonal matrix elements
c cu - real (lx,n-1) upper diagonal matrix elements
c (may be equivalent to au if no longer needed)
c r - real (lx,nx,m) right-hand-side vector elements
c (may be equivalent to a if no longer needed)
c
c output argument list:
c au - real (lx,n-1) work array
c a - real (lx,nx,m) solution vector elements
c
c attributes:
c language: fortran 77.
c machine: cray.
c
c$$$
cfpp$ noconcur r
cc
use machine
implicit none
c
integer i,j,k,l,lx,m,n,nx
real(kind=kind_evod) fk
cc
real(kind=kind_evod) cl(lx,2:n),cm(lx,n),cu(lx,n-1),r(lx,nx,m),
& au(lx,n-1),a(lx,nx,m)
cc
real(kind=kind_evod) cons1 !constant
cc
cons1 = 1.d0 !constant
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c march up
do i=1,l
fk=cons1/cm(i,1) !constant
au(i,1)=fk*cu(i,1)
enddo
do j=1,m
do i=1,l
fk=cons1/cm(i,1) !constant
a(i,1,j)=fk*r(i,1,j)
enddo
enddo
do k=2,n-1
do i=1,l
fk=cons1/(cm(i,k)-cl(i,k)*au(i,k-1)) !constant
au(i,k)=fk*cu(i,k)
enddo
do j=1,m
do i=1,l
fk=cons1/(cm(i,k)-cl(i,k)*au(i,k-1)) !constant
a(i,k,j)=fk*(r(i,k,j)-cl(i,k)*a(i,k-1,j))
enddo
enddo
enddo
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c march down
do j=1,m
do i=1,l
fk=cons1/(cm(i,n)-cl(i,n)*au(i,n-1)) !constant
a(i,n,j)=fk*(r(i,n,j)-cl(i,n)*a(i,n-1,j))
enddo
enddo
do k=n-1,1,-1
do j=1,m
do i=1,l
a(i,k,j)=a(i,k,j)-au(i,k)*a(i,k+1,j)
enddo
enddo
enddo
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end