subroutine advfct(iord,fld,u,v,scal,scali,dt,fco,fc) 8,3
c
c --- version 2.8 -- cyclic and noncyclic b.c. combined
c --- this is advem.f with MPDATA replaced by 2nd/4th order FCT
c --- if iord=1, scheme reduces to simple donor cell scheme
c
c fld - transported mixing ratio, e.g., salinity or temperature
c u,v - mass fluxes satisfying continuity equation
c scal - spatial increments (squared)
c scali - inverse of scal
c dt - temporal increment
c fco,fc - depth of the layer at previous and new time step
c
implicit none
include 'dimensions.h'
include 'dimension2.h'
c
real fld(idm,jdm),u(idm,jdm),v(idm,jdm),scal(idm,jdm),
. scali(idm,jdm),fco(idm,jdm),fc(idm,jdm)
real fmx(idm,jdm),fmn(idm,jdm),flp(idm,jdm),fln(idm,jdm),
. flx(idm,jdm),fly(idm,jdm),uan(idm,jdm),van(idm,jdm),
. flxdiv(idm,jdm),clipj(jdm),vlumj(jdm)
real dt,onemu,q,clip,vlume,amount,bfore,after,epsil
integer iord,itest,jtest,jaa
logical wrap,recovr
data recovr/.false./
c
parameter (epsil=1.e-11,onemu=.0098)
common/testpt/itest,jtest
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c --- optional code for checking conservation properties
ccc bfore=0.
ccc do 14 j=1,jj
ccc do 14 l=1,isp(j)
ccc do 14 i=ifp(j,l),ilp(j,l)
ccc 14 bfore=bfore+fld(i,j)*fco(i,j)*scal(i,j)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c --- compute low-order and antidiffusive (high-minus-low order) fluxes
c
call cpy_p
(fld)
c
c$OMP PARALLEL DO PRIVATE(ja,jaa,jb,q)
do 3 j=1,jj
ja=mod(j-2+jj,jj)+1
jb=mod(j ,jj)+1
jaa=mod(j-3+jj,jj)+1
c
do 2 l=1,isu(j)
do 2 i=ifu(j,l),ilu(j,l)
if (u(i,j).ge.0.) then
q=fld(i-1,j)
else
q=fld(i ,j)
end if
flx(i,j)=u(i,j)*q
q=fld(i,j)+fld(i-1,j) ! 2nd order
if (ip(i+1,j)+iu(i-1,j).eq.2)
. q=1.125*q-.125*(fld(i+1,j)+fld(i-2,j)) ! 4th order
2 uan(i,j)=.5*q*u(i,j)-flx(i,j)
c
do 3 l=1,isv(j)
do 3 i=ifv(j,l),ilv(j,l)
if (v(i,j).ge.0.) then
q=fld(i,ja )
else
q=fld(i,j )
end if
fly(i,j)=v(i,j)*q
q=fld(i,ja )+fld(i,j) ! 2nd order
if (ip(i,jb )+iv(i,ja).eq.2)
. q=1.125*q-.125*(fld(i,jb )+fld(i,jaa)) ! 4th order
3 van(i,j)=.5*q*v(i,j)-fly(i,j)
c$OMP END PARALLEL DO
c
c$OMP PARALLEL DO
do 22 j=1,jj
do 22 l=1,isp(j)
flx(ifp(j,l) ,j)=0.
flx(ilp(j,l)+1,j)=0.
uan(ifp(j,l) ,j)=0.
uan(ilp(j,l)+1,j)=0.
22 continue
c$OMP END PARALLEL DO
c
c$OMP PARALLEL DO PRIVATE(j,wrap)
do 33 i=1,ii1
wrap=jfv(i,1).eq.1 ! true if j=1 and j=jj are both water points
do 33 l=1,jsp(i)
j=jfp(i,l)
if (j.gt.1 .or. .not.wrap) then
fly(i,j)=0.
van(i,j)=0.
end if
j=mod(jlp(i,l),jj)+1
if (j.gt.1 .or. .not.wrap) then
fly(i,j)=0.
van(i,j)=0.
end if
33 continue
c$OMP END PARALLEL DO
c
c$OMP PARALLEL DO PRIVATE(ja,jb,ia,ib)
do 11 j=1,jj
do 11 l=1,isp(j)
do 11 i=ifp(j,l),ilp(j,l)
ia=max( 1,i-1)
if (ip(ia,j).eq.0) ia=i
ib=min(ii,i+1)
if (ip(ib,j).eq.0) ib=i
ja=mod(j-2+jj,jj)+1
if (ip(i,ja).eq.0) ja=j
jb=mod(j ,jj)+1
if (ip(i,jb).eq.0) jb=j
fmx(i,j)=max(fld(i,j),fld(ia,j),fld(ib,j),fld(i,ja),fld(i,jb))
11 fmn(i,j)=min(fld(i,j),fld(ia,j),fld(ib,j),fld(i,ja),fld(i,jb))
c$OMP END PARALLEL DO
c
cdiag i=itest
cdiag j=jtest
cdiag ja=mod(j+2+jj,jj)+1
cdiag jb=mod(j ,jj)+1
cdiag write (lp,101) 'advem (1)',i,j,fld(i-1,j),flx(i,j),fld(i,ja )
cdiag.,fly(i,j),fld(i,j),fly(i,jb ),fld(i,jb ),flx(i+1,j),fld(i+1,j)
101 format(a,2i5,f20.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3,1pe9.2,0pf9.3/
. 1pe39.2/0pf39.3)
c
c$OMP PARALLEL DO PRIVATE(jb,q,amount)
do 61 j=1,jj
jb=mod(j ,jj)+1
vlumj(j)=0.
clipj(j)=0.
do 61 l=1,isp(j)
do 61 i=ifp(j,l),ilp(j,l)
flxdiv(i,j)=(flx(i+1,j)-flx(i,j)+fly(i,jb )-fly(i,j))*dt
. *scali(i,j)
q=fld(i,j)*fco(i,j)-flxdiv(i,j)
amount=max(fmn(i,j)*fc(i,j),min(q,fmx(i,j)*fc(i,j)))
if (recovr) then
vlumj(j)=vlumj(j)+scal(i,j)*fc(i,j)
clipj(j)=clipj(j)+(q-amount)*scal(i,j)
end if
61 fld(i,j)=(fld(i,j)*onemu+amount)/(onemu+fc(i,j))
c$OMP END PARALLEL DO
c
if (iord.le.1) go to 100
c
c --- at each grid point, determine the ratio of the largest permissible
c --- pos. (neg.) change in -fld- to the sum of all incoming (outgoing) fluxes
c
c$OMP PARALLEL DO PRIVATE(jb)
do 12 j=1,jj
jb=mod(j ,jj)+1
do 12 l=1,isp(j)
do 12 i=ifp(j,l),ilp(j,l)
flp(i,j)=(fmx(i,j)-fld(i,j))*fc(i,j)
./((max(0.,uan(i,j))-min(0.,uan(i+1,j))
. +max(0.,van(i,j))-min(0.,van(i,jb ))+epsil)*dt*scali(i,j))
c
12 fln(i,j)=(fmn(i,j)-fld(i,j))*fc(i,j)
./((min(0.,uan(i,j))-max(0.,uan(i+1,j))
. +min(0.,van(i,j))-max(0.,van(i,jb ))-epsil)*dt*scali(i,j))
c$OMP END PARALLEL DO
c
c---- limit antidiffusive fluxes
c
call cpy_p(flp)
call cpy_p(fln)
c
c$OMP PARALLEL DO PRIVATE(ja,clip)
do 8 j=1,jj
ja=mod(j-2+jj,jj)+1
c
do 7 l=1,isu(j)
do 7 i=ifu(j,l),ilu(j,l)
if (uan(i,j).ge.0.) then
clip=min(1.,flp(i,j),fln(i-1,j))
else
clip=min(1.,fln(i,j),flp(i-1,j))
end if
7 flx(i,j)=uan(i,j)*clip
c
do 8 l=1,isv(j)
do 8 i=ifv(j,l),ilv(j,l)
if (van(i,j).ge.0.) then
clip=min(1.,flp(i,j),fln(i,ja ))
else
clip=min(1.,fln(i,j),flp(i,ja ))
end if
8 fly(i,j)=van(i,j)*clip
c$OMP END PARALLEL DO
c
cdiag i=itest
cdiag j=jtest
cdiag ja=mod(j+2+jj,jj)+1
cdiag jb=mod(j ,jj)+1
cdiag write (lp,101) 'advem (2)',i,j,fld(i-1,j),flx(i,j),fld(i,ja )
cdiag.,fly(i,j),fld(i,j),fly(i,jb ),fld(i,jb ),flx(i+1,j),fld(i+1,j)
c
c$OMP PARALLEL DO PRIVATE(jb,amount,q)
do 62 j=1,jj
jb=mod(j ,jj)+1
do 62 l=1,isp(j)
do 62 i=ifp(j,l),ilp(j,l)
flxdiv(i,j)=(flx(i+1,j)-flx(i,j)+fly(i,jb )-fly(i,j))*dt
. *scali(i,j)
q=fld(i,j)*fc(i,j)-flxdiv(i,j)
amount=max(fmn(i,j)*fc(i,j),min(q,fmx(i,j)*fc(i,j)))
if (recovr) clipj(j)=clipj(j)+(q-amount)*scal(i,j)
62 fld(i,j)=(fld(i,j)*onemu+amount)/(onemu+fc(i,j))
c$OMP END PARALLEL DO
c
100 continue
c
c --- revover 'clipped' amount and return to field
c
if (recovr) then
vlume=0.
clip=0.
c$OMP PARALLEL DO REDUCTION(+:vlume,clip)
do 19 j=1,jj
vlume=vlume+vlumj(j)
19 clip=clip+clipj(j)
c$OMP END PARALLEL DO
c
if (vlume.ne.0.) then
clip=clip/vlume
cdiag write (lp,'(a,1pe11.3)') 'tracer drift in advem:',-clip
c$OMP PARALLEL DO
do 13 j=1,jj
do 13 l=1,isp(j)
do 13 i=ifp(j,l),ilp(j,l)
13 fld(i,j)=fld(i,j)+clip
c$OMP END PARALLEL DO
end if
end if
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c --- optional code for checking conservation properties
ccc after=0.
ccc do 15 j=1,jj
ccc do 15 l=1,isp(j)
ccc do 15 i=ifp(j,l),ilp(j,l)
ccc 15 after=after+fld(i,j)*fc(i,j)*scal(i,j)
ccc write (lp,'(a,1p,3e14.6,e11.1)') 'advem conservation:',
ccc . bfore,after,after-bfore,(after-bfore)/bfore
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
return
end
c
c
c> Revision history:
c>
c> Mar. 2000 - removed 'cushn' and added logic to assure global conservation
c> Apr. 2000 - conversion to SI units
c> Apr. 2000 - changed i/j loop nesting to j/i
c> May 2000 - modified j-1,j+1 to accomodate both channel & closed basin b.c.
c> Sep. 2000 - fixed cyclicity problem in loop 33
c> Nov. 2004 - switched from mpdata to fct
c> Mar. 2006 - added bering strait exchange logic