--- MITgcm/pkg/generic_advdiff/gad_advection.F 2001/09/10 01:22:48 1.1 +++ MITgcm/pkg/generic_advdiff/gad_advection.F 2002/03/05 15:17:45 1.10 @@ -1,37 +1,113 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.1 2001/09/10 01:22:48 adcroft Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.10 2002/03/05 15:17:45 jmc Exp $ C $Name: $ +CBOI +C !TITLE: pkg/generic\_advdiff +C !AUTHORS: adcroft@mit.edu +C !INTRODUCTION: Generic Advection Diffusion Package +C +C Package "generic\_advdiff" provides a common set of routines for calculating +C advective/diffusive fluxes for tracers (cell centered quantities on a C-grid). +C +C Many different advection schemes are available: the standard centered +C second order, centered fourth order and upwind biased third order schemes +C are known as linear methods and require some stable time-stepping method +C such as Adams-Bashforth. Alternatives such as flux-limited schemes are +C stable in the forward sense and are best combined with the multi-dimensional +C method provided in gad\_advection. +C +C There are two high-level routines: +C \begin{itemize} +C \item{GAD\_CALC\_RHS} calculates all fluxes at time level "n" and is used +C for the standard linear schemes. This must be used in conjuction with +C Adams-Bashforth time-stepping. Diffusive and parameterized fluxes are +C always calculated here. +C \item{GAD\_ADVECTION} calculates just the advective fluxes using the +C non-linear schemes and can not be used in conjuction with Adams-Bashforth +C time-stepping. +C \end{itemize} +CEOI + #include "GAD_OPTIONS.h" +CBOP +C !ROUTINE: GAD_ADVECTION + +C !INTERFACE: ========================================================== SUBROUTINE GAD_ADVECTION(bi,bj,advectionScheme,tracerIdentity, U Tracer,Gtracer, I myTime,myIter,myThid) -C /==========================================================\ -C | SUBROUTINE GAD_ADVECTION | -C | o Solves the pure advection tracer equation. | -C |==========================================================| -C \==========================================================/ - IMPLICIT NONE -C == Global variables === +C !DESCRIPTION: +C Calculates the tendancy of a tracer due to advection. +C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection} +C and can only be used for the non-linear advection schemes such as the +C direct-space-time method and flux-limiters. +C +C The algorithm is as follows: +C \begin{itemize} +C \item{$\theta^{(n+1/3)} = \theta^{(n)} +C - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$} +C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)} +C - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$} +C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)} +C - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$} +C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$} +C \end{itemize} +C +C The tendancy (output) is over-written by this routine. + +C !USES: =============================================================== + IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "GRID.h" #include "GAD.h" - -C == Routine arguments == +#ifdef ALLOW_AUTODIFF_TAMC +# include "tamc.h" +# include "tamc_keys.h" +#endif + +C !INPUT PARAMETERS: =================================================== +C bi,bj :: tile indices +C advectionScheme :: advection scheme to use +C tracerIdentity :: identifier for the tracer (required only for OBCS) +C Tracer :: tracer field +C myTime :: current time +C myIter :: iteration number +C myThid :: thread number INTEGER bi,bj INTEGER advectionScheme INTEGER tracerIdentity _RL Tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) - _RL Gtracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) _RL myTime INTEGER myIter INTEGER myThid -C == Local variables +C !OUTPUT PARAMETERS: ================================================== +C gTracer :: tendancy array + _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) + +C !LOCAL VARIABLES: ==================================================== +C maskUp :: 2-D array for mask at W points +C iMin,iMax,jMin,jMax :: loop range for called routines +C i,j,k :: loop indices +C kup :: index into 2 1/2D array, toggles between 1 and 2 +C kdown :: index into 2 1/2D array, toggles between 2 and 1 +C kp1 :: =k+1 for k maxpass. check tamc.h' + ENDIF +#endif /* ALLOW_AUTODIFF_TAMC */ + + IF (nipass.EQ.3) THEN + calc_fluxes_X=.FALSE. + calc_fluxes_Y=.FALSE. + IF (ipass.EQ.1 .AND. (bi.EQ.1 .OR. bi.EQ.2) ) THEN + calc_fluxes_X=.TRUE. + ELSEIF (ipass.EQ.1 .AND. (bi.EQ.4 .OR. bi.EQ.5) ) THEN + calc_fluxes_Y=.TRUE. + ELSEIF (ipass.EQ.2 .AND. (bi.EQ.1 .OR. bi.EQ.6) ) THEN + calc_fluxes_Y=.TRUE. + ELSEIF (ipass.EQ.2 .AND. (bi.EQ.3 .OR. bi.EQ.4) ) THEN + calc_fluxes_X=.TRUE. + ELSEIF (ipass.EQ.3 .AND. (bi.EQ.2 .OR. bi.EQ.3) ) THEN + calc_fluxes_Y=.TRUE. + ELSEIF (ipass.EQ.3 .AND. (bi.EQ.5 .OR. bi.EQ.6) ) THEN + calc_fluxes_X=.TRUE. + ENDIF + ELSE + calc_fluxes_X=.TRUE. + calc_fluxes_Y=.TRUE. + ENDIF + +C-- X direction + IF (calc_fluxes_X) THEN + +C-- Internal exchange for calculations in X + IF (useCubedSphereExchange) THEN + DO j=1,Oly + DO i=1,Olx + localTij( 1-i , 1-j )=localTij( 1-j , i ) + localTij( 1-i ,sNy+j)=localTij( 1-j , sNy+1-i ) + localTij(sNx+i, 1-j )=localTij(sNx+j, i ) + localTij(sNx+i,sNy+j)=localTij(sNx+j, sNy+1-i ) + ENDDO + ENDDO + ENDIF + C- Advective flux in X DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx af(i,j) = 0. ENDDO ENDDO + +#ifdef ALLOW_AUTODIFF_TAMC +#ifndef DISABLE_MULTIDIM_ADVECTION +CADJ STORE localTij(:,:) = comlev1_bibj_pass, key=passkey, byte=isbyte +#endif +#endif /* ALLOW_AUTODIFF_TAMC */ + IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN CALL GAD_FLUXLIMIT_ADV_X( & bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid) @@ -100,8 +256,10 @@ CALL GAD_DST3FL_ADV_X( & bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid) ELSE - STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim' + write(0,*) advectionScheme + STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim' ENDIF + DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx-1 localTij(i,j)=localTij(i,j)-deltaTtracer* @@ -124,12 +282,37 @@ END IF #endif /* ALLOW_OBCS */ +C-- End of X direction + ENDIF + +C-- Y direction + IF (calc_fluxes_Y) THEN + +C-- Internal exchange for calculations in Y + IF (useCubedSphereExchange) THEN + DO j=1,Oly + DO i=1,Olx + localTij( 1-i , 1-j )=localTij( j , 1-i ) + localTij( 1-i ,sNy+j)=localTij( j ,sNy+i) + localTij(sNx+i, 1-j )=localTij(sNx+1-j, 1-i ) + localTij(sNx+i,sNy+j)=localTij(sNx+1-j,sNy+i) + ENDDO + ENDDO + ENDIF + C- Advective flux in Y DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx af(i,j) = 0. ENDDO ENDDO + +#ifdef ALLOW_AUTODIFF_TAMC +#ifndef DISABLE_MULTIDIM_ADVECTION +CADJ STORE localTij(:,:) = comlev1_bibj_pass, key=passkey, byte=isbyte +#endif +#endif /* ALLOW_AUTODIFF_TAMC */ + IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN CALL GAD_FLUXLIMIT_ADV_Y( & bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid) @@ -142,6 +325,7 @@ ELSE STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim' ENDIF + DO j=1-Oly,sNy+Oly-1 DO i=1-Olx,sNx+Olx localTij(i,j)=localTij(i,j)-deltaTtracer* @@ -152,6 +336,7 @@ & ) ENDDO ENDDO + #ifdef ALLOW_OBCS C-- Apply open boundary conditions IF (useOBCS) THEN @@ -162,18 +347,27 @@ END IF END IF #endif /* ALLOW_OBCS */ - DO j=1-Oly,sNy+Oly-1 + +C-- End of Y direction + ENDIF + + DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx localTijk(i,j,k)=localTij(i,j) ENDDO ENDDO +C-- End of ipass loop + ENDDO C-- End of K loop for horizontal fluxes ENDDO C-- Start of k loop for vertical flux DO k=Nr,1,-1 +#ifdef ALLOW_AUTODIFF_TAMC + kkey = (ikey-1)*Nr + k +#endif /* ALLOW_AUTODIFF_TAMC */ C-- kup Cycles through 1,2 to point to w-layer above C-- kDown Cycles through 2,1 to point to w-layer below @@ -193,6 +387,11 @@ ENDDO ENDDO +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE localTijk(:,:,k) +CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ + C Note: wVel needs to be masked IF (K.GE.2) THEN C- Compute vertical advective flux in the interior: @@ -200,13 +399,11 @@ CALL GAD_FLUXLIMIT_ADV_R( & bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid) ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN -c CALL GAD_DST3_ADV_R( -c & bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid) - STOP 'GAD_ADVECTION: adv. scheme not avail. yet' + CALL GAD_DST3_ADV_R( + & bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid) ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN -c CALL GAD_DST3FL_ADV_R( -c & bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid) - STOP 'GAD_ADVECTION: adv. scheme not avail. yet' + CALL GAD_DST3FL_ADV_R( + & bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid) ELSE STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim' ENDIF @@ -215,14 +412,16 @@ DO i=1-Olx,sNx+Olx af(i,j) = af(i,j) & + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))* - & rTrans(i,j)*localTijk(i,j,k) + & rTrans(i,j)*tracer(i,j,k,bi,bj) +c & rTrans(i,j)*localTijk(i,j,k) ENDDO ENDDO ELSE C- Surface "correction" term at k=1 : DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx - af(i,j) = rTrans(i,j)*localTijk(i,j,k) + af(i,j) = rTrans(i,j)*tracer(i,j,k,bi,bj) +c af(i,j) = rTrans(i,j)*localTijk(i,j,k) ENDDO ENDDO ENDIF