SUBROUTINE ADVT(FB,F,FCLIM,DTI2,FF,NItera) C********************************************************************** C THIS SUBROUTINE CALCULATES THE 3D ADVECTION OF PASSIVE SCALAR * C TRACERS BY A FIRST-ORDER UPSTREAM SCHEME. TO REDUCE IMPLICIT * C DIFFUSION THE SMOLARKIEWICZ ITERATIVE UPSTREAM SCHEME WITH A * C SPECIAL DEFINED ANTIDIFFUSIVE VELOCITY IS USED. * C * C Gianmaria Sannino' Vincenzo Artale'' * C * C ' Inter-University Computing Consortium (CASPUR) * C Universita' "La Sapienza", Rome, ITALY * C * C '' Italian National Agency for New Technology and Environment * C (ENEA C.R. Casaccia), Rome, ITALY * C * C Input: * C FB -----> As Subroutine ADVT in POM * C F ------> As Subroutine ADVT in POM * C FCLIM --> As Subroutine ADVT in POM * C DTI2 ---> As Subroutine ADVT in POM * C NItera -> Number of iteration * C Output: * C FF -----> As Subroutine ADVT in POM source code * C * C Smolarkiewicz, P.K. "A Fully Multidimensional Positive Definite * C Advection Transport Algorithm with Small Implicit Diffusion" * C Journal of Computational Physics 54, 325-362 (1984) * C * C Notes for POM users: * C 1. replaces subroutine ADVT * C 2. suggested number of iterations: 1 to 4 (1 is standard upstream; * C more iterations are less diffusive; 3 itr. adds 50% to POM) * C 3. can run with no additional diff. (TPRNI=0.) * C 4. the shock switch option may not work for most cases, to ignore * C it set SW=Constant (should be 1, but smaller value gives * C smoother solution with less overshoot for NItera>1 ). * C * C (Try on your own risk, no advice or support will be provided) T.E. * C********************************************************************** C INCLUDE 'comblk97.h' C DIMENSION FB(IM,JM,KB),F(IM,JM,KB),FF(IM,JM,KB) DIMENSION FCLIM(IM,JM,KB),FBMem(IM,JM,KB),Eta(IM,JM) DIMENSION XMassFlux(IM,JM,KB),YMassFlux(IM,JM,KB),ZWFlux(IM,JM,KB) DIMENSION XFlux(IM,JM,KB),YFlux(IM,JM,KB),ZFlux(IM,JM,KB) C C----------------------------------------------------------------------- C HORIZONTAL MASS FLUXES: DY*U*D & DX*V*D C----------------------------------------------------------------------- DO 100 K=1,KB DO 100 J=1,JM DO 100 I=1,IM XMassFlux(I,J,K)=0. YMassFlux(I,J,K)=0. 100 CONTINUE C DO 110 K=1,KBM1 DO 111 J=2,JMM1 DO 111 I=2,IM XMassFlux(I,J,K)=0.25*(DY(I-1,J)+DY(I,J))*(DT(I-1,J)+DT(I,J))* & U(I,J,K) 111 CONTINUE DO 112 J=2,JM DO 112 I=2,IMM1 YMassFlux(I,J,K)=0.25*(DX(I,J-1)+DX(I,J))*(DT(I,J-1)+DT(I,J))* & V(I,J,K) 112 CONTINUE 110 CONTINUE C DO 120 J=1,JM DO 120 I=1,IM FB(I,J,KB)=FB(I,J,KBM1) 120 Eta(I,J)=ETB(I,J) C DO 130 K=1,KB DO 130 J=1,JM DO 130 I=1,IM ZWFlux(I,J,K)=W(I,J,K) 130 FBMem(I,J,K)=FB(I,J,K) C----------------------------------------------------------------------- C START C SMOLARKIEWICZ SCHEME C----------------------------------------------------------------------- DO 999 ITERA=1,NItera C----------------------------------------------------------------------- C UPWIND ADVECTION SCHEME C----------------------------------------------------------------------- DO 140 K=1,KBM1 DO 140 J=2,JM DO 140 I=2,IM XFlux(I,J,K)=0.5*( & (XMassFlux(I,J,K)+ABS(XMassFlux(I,J,K))) & *FBMem(I-1,J,K)+ & (XMassFlux(I,J,K)-ABS(XMassFlux(I,J,K))) & *FBMem(I,J,K)) C YFlux(I,J,K)=0.5*( & (YMassFlux(I,J,K)+ABS(YMassFlux(I,J,K))) & *FBMem(I,J-1,K)+ & (YMassFlux(I,J,K)-ABS(YMassFlux(I,J,K))) & *FBMem(I,J,K)) 140 CONTINUE C DO 150 J=2,JMM1 DO 150 I=2,IMM1 ZFlux(I,J,1)=0.0 150 ZFlux(I,J,KB)=0.0 C DO 160 K=2,KBM1 DO 160 J=2,JMM1 DO 160 I=2,IMM1 ZFlux(I,J,K)=0.5*( & (ZWFlux(I,J,K)+ABS(ZWFlux(I,J,K))) & *FBMem(I,J,K)+ & (ZWFlux(I,J,K)-ABS(ZWFlux(I,J,K))) & *FBMem(I,J,K-1)) 160 ZFlux(I,J,K)=ZFlux(I,J,K)*ART(I,J) C----------------------------------------------------------------------- C ADD NET ADVECTIVE FLUXES; THEN STEP FORWARD IN TIME C----------------------------------------------------------------------- DO 170 K=1,KBM1 DO 170 J=2,JMM1 DO 170 I=2,IMM1 FF(I,J,K)= XFlux(I+1,J,K)-XFlux(I,J,K) & +YFlux(I,J+1,K)-YFlux(I,J,K) & +(ZFlux(I,J,K)-ZFlux(I,J,K+1))/DZ(K) 170 FF(I,J,K)=(FBMem(I,J,K)*(H(I,J)+Eta(I,J))*ART(I,J)- & DTI2*FF(I,J,K))/((H(I,J)+ETF(I,J))*ART(I,J)) C----------------------------------------------------------------------- C ANTIDIFFUSION VELOCITY C----------------------------------------------------------------------- CALL SMOL_ADIF(XMassFlux,YMassFlux,ZWFlux,FB,FF,DTI2) C DO 180 J=1,JM DO 180 I=1,IM Eta(I,J)=ETF(I,J) DO 180 K=1,KB FBMem(I,J,K)=FF(I,J,K) 180 CONTINUE C----------------------------------------------------------------------- C END C SMOLARKIEWICZ SCHEME C----------------------------------------------------------------------- 999 CONTINUE C----------------------------------------------------------------------- C ADD HORIZONTAL DIFFUSIVE FLUXES C----------------------------------------------------------------------- DO 190 K=1,KB DO 190 J=1,JM DO 190 I=1,IM 190 FB(I,J,K)=FB(I,J,K)-FCLIM(I,J,K) C DO 200 K=1,KBM1 DO 200 J=2,JM DO 200 I=2,IM XMassFlux(I,J,K)=0.5*(AAM(I,J,K)+AAM(I-1,J,K)) YMassFlux(I,J,K)=0.5*(AAM(I,J,K)+AAM(I,J-1,K)) 200 CONTINUE C DO 210 K=1,KBM1 DO 210 J=2,JM DO 210 I=2,IM XFlux(I,J,K)= & -XMassFlux(I,J,K)*(H(I,J)+H(I-1,J))*TPRNI & *(FB(I,J,K)-FB(I-1,J,K))*DUM(I,J)/(DX(I,J)+DX(I-1,J)) & *0.5*(DY(I,J)+DY(I-1,J)) YFlux(I,J,K)= & -YMassFlux(I,J,K)*(H(I,J)+H(I,J-1))*TPRNI & *(FB(I,J,K)-FB(I,J-1,K))*DVM(I,J)/(DY(I,J)+DY(I,J-1)) & *0.5*(DX(I,J)+DX(I,J-1)) 210 CONTINUE C DO 220 K=1,KB DO 220 J=1,JM DO 220 I=1,IM 220 FB(I,J,K)=FB(I,J,K)+FCLIM(I,J,K) C C----- ADD NET HORIZONTAL FLUXES; THEN STEP FORWARD IN TIME -------- DO 230 K=1,KBM1 DO 230 J=2,JMM1 DO 230 I=2,IMM1 230 FF(I,J,K)=FF(I,J,K)-DTI2*( & +XFlux(I+1,J,K)-XFlux(I,J,K) & +YFlux(I,J+1,K)-YFlux(I,J,K)) & /((H(I,J)+ETF(I,J))*ART(I,J)) C-------------------------------------------------------------------- RETURN END SUBROUTINE SMOL_ADIF(XFlux,YFlux,ZFlux,FB,FF,DTI2) C********************************************************************** C THIS SUBROUTINE CALCULATES THE ANTIDIFFUSIVE VELOCITY USED TO * C REDUCE THE NUMERICAL DIFFUSION ASSOCIATED WITH THE UPSTREAM * C DIFFERENCING SCHEME. A SWITCH OPTION IS ALSO IMPLEMENTED FOR * C THE SHOCK-TYPE INITIAL CONDITIONS * C********************************************************************** C INCLUDE 'comblk97.h' C DIMENSION FF(IM,JM,KB),FB(IM,JM,KB) DIMENSION XFlux(IM,JM,KB),YFlux(IM,JM,KB),ZFlux(IM,JM,KB) DIMENSION SW_U(IM,JM,KB),SW_V(IM,JM,KB),SW_Z(IM,JM,KB) REAL Mol PARAMETER (VALUE_MAX=1.E20,VALUE_MIN=1.E-9,EPSILON=1.0E-14) PARAMETER (ES=.15E-2) C C------------------------------------------------------------------------ C TEMPERATURE AND SALINITY MASK C------------------------------------------------------------------------ DO 100 K=1,KB DO 100 I=1,IM DO 100 J=1,JM FF(I,J,K)=FF(I,J,K)*FSM(I,J) 100 CONTINUE c----------------------------------------------------------------------- c SWITCH OPTION c----------------------------------------------------------------------- DO 102 K=1,KBM1 DO 102 J=2,JMM1 DO 102 I=2,IM A2=(FF(I,J,K)+FF(I-1,J,K))/(FF(I-1,J,K)+FF(I,J,K)+EPSILON) A1=(abs(FF(I-1,J,K)-FF(I,J,K)))/(FF(I-1,J,K)+FF(I,J,K)+EPSILON) SW_U(I,J,K)=.5*(1.+sign(1.,-1.*((A1-ES)*A2))) 102 CONTINUE C DO 104 K=1,KBM1 DO 104 J=2,JM DO 104 I=2,IMM1 A1=(abs(-FF(I,J,K)+FF(I,J-1,K)))/(FF(I,J-1,K)+FF(I,J,K)+EPSILON) A2=(FF(I,J,K)+FF(I,J-1,K))/(FF(I,J-1,K)+FF(I,J,K)+EPSILON) SW_V(I,J,K)=.5*(1.+sign(1.,-1.*((A1-ES)*A2))) 104 CONTINUE C DO 105 K=2,KBM1 DO 105 J=2,JMM1 DO 105 I=2,IMM1 A1=(abs(FF(I,J,K-1)-FF(I,J,K)))/(FF(I,J,K)+FF(I,J,K-1)+EPSILON) A2=(FF(I,J,K-1)+FF(I,J,K))/(FF(I,J,K)+FF(I,J,K-1)+EPSILON) SW_Z(I,J,K)=.5*(1.+sign(1.,-1.*((A1-ES)*A2))) 105 CONTINUE C C----------------------------------------------------------------------- C RECALCULATE MASS FLUXES WITH ANTIDIFFUSION VELOCITY C----------------------------------------------------------------------- DO 110 K=1,KBM1 DO 110 J=2,JMM1 DO 110 I=2,IM IF(FF(I,J,K).LT.VALUE_MIN.OR.FF(I-1,J,K).LT.VALUE_MIN) THEN XFlux(I,J,K)=0.0 ELSE SW=min(SW_U(I-1,J,K),SW_U(I,J,K)) UDx=ABS(XFlux(I,J,K)) U2Dt=( 2.*DTI2*XFlux(I,J,K)*XFlux(I,J,K)/ & (ARU(I,J)*(DT(I-1,J)+DT(I,J))) ) Mol=(FF(I,J,K)-FF(I-1,J,K))/(FF(I-1,J,K)+FF(I,J,K)+EPSILON) XFlux(I,J,K)= (UDx - U2Dt)*Mol*SW Abs_1=ABS(UDx) Abs_2=ABS(U2Dt) IF(Abs_1 .LT. Abs_2) XFlux(I,J,K)=0.0 END IF 110 CONTINUE C DO 120 K=1,KBM1 DO 120 J=2,JM DO 120 I=2,IMM1 IF(FF(I,J,K).LT.VALUE_MIN.OR.FF(I,J-1,K).LT.VALUE_MIN) THEN YFlux(I,J,K)=0.0 ELSE SW=min(SW_V(I,J-1,K),SW_V(I,J,K)) VDy=ABS(YFlux(I,J,K)) V2Dt=( 2.*DTI2*YFlux(I,J,K)*YFlux(I,J,K)/ & (ARV(I,J)*(DT(I,J-1)+DT(I,J))) ) Mol=(FF(I,J,K)-FF(I,J-1,K))/(FF(I,J-1,K)+FF(I,J,K)+EPSILON) YFlux(I,J,K)= (VDy - V2Dt)*Mol*SW Abs_1=ABS(VDy) Abs_2=ABS(V2Dt) IF(Abs_1 .LT. Abs_2)YFlux(I,J,K)=0.0 END IF 120 CONTINUE C DO 130 K=2,KBM1 DO 130 J=2,JMM1 DO 130 I=2,IMM1 IF(FF(I,J,K).LT.VALUE_MIN.OR.FF(I,J,K-1).LT.VALUE_MIN) THEN ZFlux(I,J,K)=0.0 ELSE SW=min(SW_Z(I,J,K-1),SW_Z(I,J,K)) WDz=ABS(ZFlux(I,J,K)) W2Dt=( DTI2*ZFlux(I,J,K)*ZFlux(I,J,K)/ & (DZZ(K-1)*DT(I,J)) ) Mol=(FF(I,J,K-1)-FF(I,J,K))/(FF(I,J,K)+FF(I,J,K-1)+EPSILON) ZFlux(I,J,K)= (WDz - W2Dt)*Mol*SW Abs_1=ABS(WDz) Abs_2=ABS(W2Dt) IF(Abs_1 .LT. Abs_2)ZFlux(I,J,K)=0.0 END IF 130 CONTINUE C RETURN END