SUBROUTINE ADVTL94(FB,F,FF) C THIS SUBROUTINE INTEGRATES CONSERVATIVE CONSTITUENT EQUATIONS C PASSIVE AND ACTIVE. c It uses a conservative upstream advection scheme developed in c Lin, S.J., W.C. Chao, Y.C. Sud, and G.K. Walker, "A Class of the van Leer c type transport schemes and its application to the moisture transport in a c general circulation model", Mon. Wea. Rev., 122, 1575-1593, 1994. c include 'comblk.h' include 'comblk.o' DIMENSION FB(IM,JM,KB),F(IM,JM,KB),FF(IM,JM,KB) DIMENSION XFLUX(IM,JM,KB),YFLUX(IM,JM,KB),ZFLUX(IM,JM,KB) real delfiavg(im,jm,kb),delfi(im,jm,kb) real fimin(im,jm,kb),fimax(im,jm,kb) real co,dti EQUIVALENCE (fimin,A),(fimax,C) C FF=0. F(:,:,KB)=F(:,:,KBM1) FB(:,:,KB)=FB(:,:,KBM1) XFLUX=0. YFLUX=0. ZFLUX=0. delfiavg=0. delfi=0. fimin=0. fimax=0. c solving equation du/dt+vc*du/dx=0 c to satisfy CFL-criterion: vc*dtx<1.0 c let dt(i,j)=h(i,j)+et(i,j) C ------------ ADVECTIVE FLUXES ----------- c ----------x-direction:----------- do i=2,im-1 do j=1,jm do k=1,kb delfiavg(i,j,k)=0.5*(f(i+1,j,k)*dt(i+1,j)-f(i-1,j,k)*dt(i-1,j)) fimin(i,j,k)=min((f(i-1,j,k)*dt(i-1,j)),(f(i,j,k)*dt(i,j)), * (f(i+1,j,k)*dt(i+1,j))) fimax(i,j,k)=max((f(i-1,j,k)*dt(i-1,j)),(f(i,j,k)*dt(i,j)), * (f(i+1,j,k)*dt(i+1,j))) delfi(i,j,k)=sign(1.0,delfiavg(i,j,k))*min(abs(delfiavg(i,j,k)) * ,2*dim((f(i,j,k)*dt(i,j)),fimin(i,j,k)),2*dim(fimax(i,j,k) * ,(f(i,j,k)*dt(i,j)))) enddo enddo enddo c set border values: do j=1,jm do k=1,kb delfi(1,j,k)=delfi(2,j,k) delfi(im,j,k)=delfi(imm1,j,k) enddo enddo do 220 i=2,im do 220 j=1,jm do 220 k=1,kb-1 c wall at i: if (u(i,j,k).ge.0.) then co=0.5*(1.-u(i,j,k)*dti/dx(i-1,j)) xflux(i,j,k)=u(i,j,k)* * (f(i-1,j,k)*dt(i-1,j)+delfi(i-1,j,k)*co) else co=0.5*(1.+u(i,j,k)*dti/dx(i,j)) xflux(i,j,k)=u(i,j,k)* * (f(i,j,k)*dt(i,j)-delfi(i,j,k)*co) endif 220 enddo c ----------- y-direction:-------------------- delfiavg=0. delfi=0. fimin=0. fimax=0. do i=1,im do j=2,jmm1 do k=1,kb delfiavg(i,j,k)=0.5*(f(i,j+1,k)*dt(i,j+1)-f(i,j-1,k)*dt(i,j-1)) fimin(i,j,k)=min((f(i,j-1,k)*dt(i,j-1)),(f(i,j,k)*dt(i,j)), * (f(i,j+1,k)*dt(i,j+1))) fimax(i,j,k)=max((f(i,j-1,k)*dt(i,j-1)),(f(i,j,k)*dt(i,j)), * (f(i,j+1,k)*dt(i,j+1))) delfi(i,j,k)=sign(1.0,delfiavg(i,j,k))*min(abs(delfiavg(i,j,k)) * ,2*dim((f(i,j,k)*dt(i,j)),fimin(i,j,k)),2*dim(fimax(i,j,k) * ,(f(i,j,k)*dt(i,j)))) enddo enddo enddo do i=1,im do k=1,kb delfi(i,1,k)=delfi(i,2,k) delfi(i,jm,k)=delfi(i,jmm1,k) enddo enddo c c wall at j: do 221 i=1,im do 221 j=2,jm do 221 k=1,kb-1 if (v(i,j,k).ge.0.) then co=0.5*(1.-v(i,j,k)*dti/dy(i,j-1)) yflux(i,j,k)=v(i,j,k)* * ((f(i,j-1,k)*dt(i,j-1))+delfi(i,j-1,k)*co) else co=0.5*(1.+v(i,j,k)*dti/dy(i,j)) yflux(i,j,k)=v(i,j,k)* * ((f(i,j,k)*dt(i,j))-delfi(i,j,k)*co) endif 221 enddo C--------------ADD HORIZONTAL DIFFUSION------------------------- c since we are only time stepping once, we must either use F c instead of FB in this calculation, or multiply by 2 since we c are only jumping one timestep. Lin claims he runs without c any diffusion at all. do 121 k=1,kbm1 do 121 j=2,jm do 121 i=2,im c xflux(i,j,k)=xflux(i,j,k) c 1 -0.5*(AAM(I,J,K)+AAM(I-1,J,K))*(H(I,J)+H(I-1,J)) c 2 *(FB(I,J,K)-FB(I-1,J,K))*DUM(I,J)/(TPRNU*(DX(I,J)+DX(I-1,J))) c yflux(i,j,k)=yflux(i,j,k) c 1 -0.5*(AAM(I,J,K)+AAM(I,J-1,K))*(H(I,J)+H(I,J-1)) c 2 *(FB(I,J,K)-FB(I,J-1,K))*DVM(I,J)/(TPRNU*(DY(I,J)+DY(I,J-1))) c xflux(i,j,k)=0.5*(dy(i,j)+dy(i-1,j))*xflux(i,j,k) yflux(i,j,k)=0.5*(dx(i,j)+dx(i,j-1))*yflux(i,j,k) 121 enddo C--------------VERTICAL ADVECTION------------------ delfiavg=0. delfi=0. fimin=0. fimax=0. c z-direction: do i=1,im do j=1,jm do k=2,kb-1 delfiavg(i,j,k)=0.5*(f(i,j,k-1)-f(i,j,k+1)) fimin(i,j,k)=min(f(i,j,k-1),f(i,j,k),f(i,j,k+1)) fimax(i,j,k)=max(f(i,j,k-1),f(i,j,k),f(i,j,k+1)) delfi(i,j,k)=sign(1.0,delfiavg(i,j,k))*min(abs(delfiavg(i,j,k)) * ,2*dim(f(i,j,k),fimin(i,j,k)),2*dim(fimax(i,j,k) * ,f(i,j,k))) enddo enddo enddo do i=1,im do j=1,jm delfi(i,j,1)=delfi(i,j,2) delfi(i,j,kb)=delfi(i,j,kbm1) enddo enddo c wall at k: k=1 do 505 j=1,jm do 505 i=1,im if (w(i,j,k).ge.0.) then co=0.5*(1.-w(i,j,k)*dti/(dz(k)*dt(i,j))) zflux(i,j,k)=w(i,j,k)* * (f(i,j,k)+delfi(i,j,k)*co) else zflux(i,j,k)=0.0 endif 505 continue DO 160 K=2,KBM1 DO 160 J=1,JM DO 160 I=1,IM if (w(i,j,k).ge.0.) then co=0.5*(1.-w(i,j,k)*dti/(dz(k)*dt(i,j))) zflux(i,j,k)=w(i,j,k)* * (f(i,j,k)+delfi(i,j,k)*co) else co=0.5*(1.+w(i,j,k)*dti/(dz(k-1)*dt(i,j))) zflux(i,j,k)=w(i,j,k)* * (f(i,j,k-1)-delfi(i,j,k-1)*co) endif 160 continue c wall at kb: k=kb DO 506 J=1,JM DO 506 I=1,IM if (w(i,j,k).ge.0.) then zflux(i,j,k)=0.0 else co=0.5*(1.+w(i,j,k)*dti/(dz(k-1)*dt(i,j))) zflux(i,j,k)=w(i,j,k)* * (f(i,j,k-1)-delfi(i,j,k-1)*co) endif 506 continue 668 format(16(2x,e10.4)) DO 161 K=1,KBM1 DO 161 J=1,JMM1 DO 161 I=1,IMM1 161 ff(i,j,k)=dzr(k)*(zflux(i,j,k)-zflux(i,j,k+1)) C-----------ADD NET HORIZONTAL FLUXES------------------ c divide by dxdy: DO 180 K=1,KBM1 DO 180 J=2,JMM1 DO 180 I=2,IMM1 180 FF(I,J,K)=FF(I,J,K) 1 +(XFLUX(I+1,J,K)-XFLUX(I,J,K) 2 + YFLUX(I,J+1,K)-YFLUX(I,J,K))/ART(I,J) DO 200 K=1,KBM1 DO 200 J=2,JMM1 DO 200 I=2,IMM1 C ---------- STEP FORWARD --------- 200 FF(I,J,K)=(f(I,J,K)*(H(I,J)+ET(I,J))-DTI*FF(I,J,K)) 1 /(H(I,J)+ETF(I,J)) RETURN END