--- MITgcm/pkg/generic_advdiff/gad_advection.F 2004/03/27 03:51:51 1.19 +++ MITgcm/pkg/generic_advdiff/gad_advection.F 2007/10/19 19:18:01 1.53 @@ -1,34 +1,8 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.19 2004/03/27 03:51:51 edhill Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.53 2007/10/19 19:18:01 heimbach 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" +#undef MULTIDIM_OLD_VERSION C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP @@ -36,16 +10,17 @@ C !INTERFACE: ========================================================== SUBROUTINE GAD_ADVECTION( - I implicitAdvection, advectionScheme, tracerIdentity, - I uVel, vVel, wVel, tracer, - O gTracer, - I bi,bj, myTime,myIter,myThid) + I implicitAdvection, advectionScheme, vertAdvecScheme, + I tracerIdentity, + I uVel, vVel, wVel, tracer, + O gTracer, + I bi,bj, myTime,myIter,myThid) C !DESCRIPTION: -C Calculates the tendancy of a tracer due to advection. +C Calculates the tendency 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 direct-space-time method and flux-limiters. C C The algorithm is as follows: C \begin{itemize} @@ -58,7 +33,7 @@ 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 The tendency (output) is over-written by this routine. C !USES: =============================================================== IMPLICIT NONE @@ -70,22 +45,30 @@ #ifdef ALLOW_AUTODIFF_TAMC # include "tamc.h" # include "tamc_keys.h" +# ifdef ALLOW_PTRACERS +# include "PTRACERS_SIZE.h" +# endif #endif +#ifdef ALLOW_EXCH2 +#include "W2_EXCH2_TOPOLOGY.h" +#include "W2_EXCH2_PARAMS.h" +#endif /* ALLOW_EXCH2 */ C !INPUT PARAMETERS: =================================================== -C implicitAdvection :: vertical advection treated implicitly (later on) -C advectionScheme :: advection scheme to use -C tracerIdentity :: identifier for the tracer (required only for OBCS) -C uVel :: velocity, zonal component -C vVel :: velocity, meridional component -C wVel :: velocity, vertical component -C tracer :: tracer field -C bi,bj :: tile indices -C myTime :: current time -C myIter :: iteration number -C myThid :: thread number +C implicitAdvection :: implicit vertical advection (later on) +C advectionScheme :: advection scheme to use (Horizontal plane) +C vertAdvecScheme :: advection scheme to use (vertical direction) +C tracerIdentity :: tracer identifier (required only for OBCS) +C uVel :: velocity, zonal component +C vVel :: velocity, meridional component +C wVel :: velocity, vertical component +C tracer :: tracer field +C bi,bj :: tile indices +C myTime :: current time +C myIter :: iteration number +C myThid :: thread number LOGICAL implicitAdvection - INTEGER advectionScheme + INTEGER advectionScheme, vertAdvecScheme INTEGER tracerIdentity _RL uVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) _RL vVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) @@ -97,48 +80,83 @@ INTEGER myThid C !OUTPUT PARAMETERS: ================================================== -C gTracer :: tendancy array +C gTracer :: tendency 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 maxcube. 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. + interiorOnly = .FALSE. + overlapOnly = .FALSE. + IF (useCubedSphereExchange) THEN +#ifdef MULTIDIM_OLD_VERSION +C- CubedSphere : pass 3 times, with full update of local tracer field + IF (ipass.EQ.1) THEN + calc_fluxes_X = nCFace.EQ.1 .OR. nCFace.EQ.2 + calc_fluxes_Y = nCFace.EQ.4 .OR. nCFace.EQ.5 + ELSEIF (ipass.EQ.2) THEN + calc_fluxes_X = nCFace.EQ.3 .OR. nCFace.EQ.4 + calc_fluxes_Y = nCFace.EQ.6 .OR. nCFace.EQ.1 +#else /* MULTIDIM_OLD_VERSION */ +C- CubedSphere : pass 3 times, with partial update of local tracer field + IF (ipass.EQ.1) THEN + overlapOnly = MOD(nCFace,3).EQ.0 + interiorOnly = MOD(nCFace,3).NE.0 + calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2 + calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5 + ELSEIF (ipass.EQ.2) THEN + overlapOnly = MOD(nCFace,3).EQ.2 + calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4 + calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1 +#endif /* MULTIDIM_OLD_VERSION */ + ELSE + calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6 + calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3 ENDIF ELSE - calc_fluxes_X=.TRUE. - calc_fluxes_Y=.TRUE. +C- not CubedSphere + calc_fluxes_X = MOD(ipass,2).EQ.1 + calc_fluxes_Y = .NOT.calc_fluxes_X 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 ) +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +C-- X direction +C- Advective flux in X + DO j=1-Oly,sNy+Oly + DO i=1-Olx,sNx+Olx + af(i,j) = 0. + ENDDO ENDDO - ENDDO - ENDIF +C +#ifdef ALLOW_AUTODIFF_TAMC +# ifndef DISABLE_MULTIDIM_ADVECTION +CADJ STORE localTij(:,:) = +CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte +CADJ STORE af(:,:) = +CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ +C + IF (calc_fluxes_X) THEN -C- Advective flux in X - DO j=1-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx - af(i,j) = 0. - ENDDO - ENDDO +C- Do not compute fluxes if +C a) needed in overlap only +C and b) the overlap of myTile are not cube-face Edges + IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN + +cph-exch2#ifndef ALLOW_AUTODIFF_TAMC +C- Internal exchange for calculations in X +#ifdef MULTIDIM_OLD_VERSION + IF ( useCubedSphereExchange ) THEN +#else + IF ( useCubedSphereExchange .AND. + & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN +#endif + CALL FILL_CS_CORNER_TR_RL( .TRUE., .FALSE., + & localTij, bi,bj, myThid ) + ENDIF +cph-exch2#endif #ifdef ALLOW_AUTODIFF_TAMC -#ifndef DISABLE_MULTIDIM_ADVECTION +# ifndef DISABLE_MULTIDIM_ADVECTION CADJ STORE localTij(:,:) = CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte -#endif +# 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) - ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN - CALL GAD_DST3_ADV_X( - & bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid) - ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN - CALL GAD_DST3FL_ADV_X( - & bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid) - ELSE - STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim' - ENDIF + IF ( advectionScheme.EQ.ENUM_UPWIND_1RST + & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN + CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE., + I dTtracerLev(k),uTrans,uFld,localTij, + O af, myThid ) + ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN + CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k), + I uTrans, uFld, maskLocW, localTij, + O af, myThid ) + ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN + CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k), + I uTrans, uFld, maskLocW, localTij, + O af, myThid ) + ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN + CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k), + I uTrans, uFld, maskLocW, localTij, + O af, myThid ) +#ifndef ALLOW_AUTODIFF_TAMC + ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN + CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k), + I uTrans, uFld, maskLocW, localTij, + O af, myThid ) +#endif + ELSE + 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* - & _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) - & *recip_rA(i,j,bi,bj) - & *( af(i+1,j)-af(i,j) - & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j)) - & ) - ENDDO - ENDDO +C- Advective flux in X : done + ENDIF + +cph-exch2#ifndef ALLOW_AUTODIFF_TAMC +C- Internal exchange for next calculations in Y + IF ( overlapOnly .AND. ipass.EQ.1 ) THEN + CALL FILL_CS_CORNER_TR_RL(.FALSE., .FALSE., + & localTij, bi,bj, myThid ) + ENDIF +cph-exch2#endif + +C- Update the local tracer field where needed: + +C update in overlap-Only + IF ( overlapOnly ) THEN + iMinUpd = 1-Olx+1 + iMaxUpd = sNx+Olx-1 +C- notes: these 2 lines below have no real effect (because recip_hFac=0 +C in corner region) but safer to keep them. + IF ( W_edge ) iMinUpd = 1 + IF ( E_edge ) iMaxUpd = sNx + + IF ( S_edge ) THEN + DO j=1-Oly,0 + DO i=iMinUpd,iMaxUpd + localTij(i,j) = localTij(i,j) + & -dTtracerLev(k)*recip_rhoFacC(k) + & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) + & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) + & *( af(i+1,j)-af(i,j) + & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j)) + & ) + ENDDO + ENDDO + ENDIF + IF ( N_edge ) THEN + DO j=sNy+1,sNy+Oly + DO i=iMinUpd,iMaxUpd + localTij(i,j) = localTij(i,j) + & -dTtracerLev(k)*recip_rhoFacC(k) + & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) + & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) + & *( af(i+1,j)-af(i,j) + & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j)) + & ) + ENDDO + ENDDO + ENDIF + + ELSE +C do not only update the overlap + jMinUpd = 1-Oly + jMaxUpd = sNy+Oly + IF ( interiorOnly .AND. S_edge ) jMinUpd = 1 + IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy + DO j=jMinUpd,jMaxUpd + DO i=1-Olx+1,sNx+Olx-1 + localTij(i,j) = localTij(i,j) + & -dTtracerLev(k)*recip_rhoFacC(k) + & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) + & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) + & *( af(i+1,j)-af(i,j) + & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j)) + & ) + ENDDO + ENDDO +C- keep advective flux (for diagnostics) + DO j=1-Oly,sNy+Oly + DO i=1-Olx,sNx+Olx + afx(i,j) = af(i,j) + ENDDO + ENDDO #ifdef ALLOW_OBCS -C-- Apply open boundary conditions - IF (useOBCS) THEN - IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN - CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid ) - ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN - CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid ) - END IF - END IF +C- Apply open boundary conditions + IF ( useOBCS ) THEN + IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN + CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid ) + ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN + CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid ) +#ifdef ALLOW_PTRACERS + ELSEIF (tracerIdentity.GE.GAD_TR1) THEN + CALL OBCS_APPLY_PTRACER( bi, bj, k, + & tracerIdentity-GAD_TR1+1, localTij, myThid ) +#endif /* ALLOW_PTRACERS */ + ENDIF + ENDIF #endif /* ALLOW_OBCS */ +C- end if/else update overlap-Only + ENDIF + C-- End of X direction ENDIF +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Y direction +cph-test +C- Advective flux in Y + DO j=1-Oly,sNy+Oly + DO i=1-Olx,sNx+Olx + af(i,j) = 0. + ENDDO + ENDDO +C +#ifdef ALLOW_AUTODIFF_TAMC +# ifndef DISABLE_MULTIDIM_ADVECTION +CADJ STORE localTij(:,:) = +CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte +CADJ STORE af(:,:) = +CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ +C 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- Do not compute fluxes if +C a) needed in overlap only +C and b) the overlap of myTile are not cube-face edges + IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN + +cph-exch2#ifndef ALLOW_AUTODIFF_TAMC +C- Internal exchange for calculations in Y +#ifdef MULTIDIM_OLD_VERSION + IF ( useCubedSphereExchange ) THEN +#else + IF ( useCubedSphereExchange .AND. + & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN +#endif + CALL FILL_CS_CORNER_TR_RL(.FALSE., .FALSE., + & localTij, bi,bj, myThid ) + ENDIF +cph-exch2#endif -C- Advective flux in Y - DO j=1-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx - af(i,j) = 0. - ENDDO - ENDDO +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 @@ -347,41 +543,125 @@ #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) - ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN - CALL GAD_DST3_ADV_Y( - & bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid) - ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN - CALL GAD_DST3FL_ADV_Y( - & bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid) - ELSE - STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim' - ENDIF + IF ( advectionScheme.EQ.ENUM_UPWIND_1RST + & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN + CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE., + I dTtracerLev(k),vTrans,vFld,localTij, + O af, myThid ) + ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN + CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k), + I vTrans, vFld, maskLocS, localTij, + O af, myThid ) + ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN + CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k), + I vTrans, vFld, maskLocS, localTij, + O af, myThid ) + ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN + CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k), + I vTrans, vFld, maskLocS, localTij, + O af, myThid ) +#ifndef ALLOW_AUTODIFF_TAMC + ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN + CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k), + I vTrans, vFld, maskLocS, localTij, + O af, myThid ) +#endif + 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* - & _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) - & *recip_rA(i,j,bi,bj) - & *( af(i,j+1)-af(i,j) - & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j)) - & ) - ENDDO - ENDDO +C- Advective flux in Y : done + ENDIF + +cph-exch2#ifndef ALLOW_AUTODIFF_TAMC +C- Internal exchange for next calculations in X + IF ( overlapOnly .AND. ipass.EQ.1 ) THEN + CALL FILL_CS_CORNER_TR_RL( .TRUE., .FALSE., + & localTij, bi,bj, myThid ) + ENDIF +cph-exch2#endif + +C- Update the local tracer field where needed: + +C update in overlap-Only + IF ( overlapOnly ) THEN + jMinUpd = 1-Oly+1 + jMaxUpd = sNy+Oly-1 +C- notes: these 2 lines below have no real effect (because recip_hFac=0 +C in corner region) but safer to keep them. + IF ( S_edge ) jMinUpd = 1 + IF ( N_edge ) jMaxUpd = sNy + + IF ( W_edge ) THEN + DO j=jMinUpd,jMaxUpd + DO i=1-Olx,0 + localTij(i,j) = localTij(i,j) + & -dTtracerLev(k)*recip_rhoFacC(k) + & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) + & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) + & *( af(i,j+1)-af(i,j) + & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j)) + & ) + ENDDO + ENDDO + ENDIF + IF ( E_edge ) THEN + DO j=jMinUpd,jMaxUpd + DO i=sNx+1,sNx+Olx + localTij(i,j) = localTij(i,j) + & -dTtracerLev(k)*recip_rhoFacC(k) + & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) + & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) + & *( af(i,j+1)-af(i,j) + & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j)) + & ) + ENDDO + ENDDO + ENDIF + + ELSE +C do not only update the overlap + iMinUpd = 1-Olx + iMaxUpd = sNx+Olx + IF ( interiorOnly .AND. W_edge ) iMinUpd = 1 + IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx + DO j=1-Oly+1,sNy+Oly-1 + DO i=iMinUpd,iMaxUpd + localTij(i,j) = localTij(i,j) + & -dTtracerLev(k)*recip_rhoFacC(k) + & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) + & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) + & *( af(i,j+1)-af(i,j) + & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j)) + & ) + ENDDO + ENDDO +C- keep advective flux (for diagnostics) + DO j=1-Oly,sNy+Oly + DO i=1-Olx,sNx+Olx + afy(i,j) = af(i,j) + ENDDO + ENDDO #ifdef ALLOW_OBCS -C-- Apply open boundary conditions - IF (useOBCS) THEN - IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN - CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid ) - ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN - CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid ) - END IF - END IF +C- Apply open boundary conditions + IF (useOBCS) THEN + IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN + CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid ) + ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN + CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid ) +#ifdef ALLOW_PTRACERS + ELSEIF (tracerIdentity.GE.GAD_TR1) THEN + CALL OBCS_APPLY_PTRACER( bi, bj, k, + & tracerIdentity-GAD_TR1+1, localTij, myThid ) +#endif /* ALLOW_PTRACERS */ + ENDIF + ENDIF #endif /* ALLOW_OBCS */ +C end if/else update overlap-Only + ENDIF + C-- End of Y direction ENDIF @@ -393,122 +673,194 @@ DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx gTracer(i,j,k,bi,bj)= - & (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer + & (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k) ENDDO ENDDO ELSE C- horizontal advection done; store intermediate result in 3D array: - DO j=1-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx - localTijk(i,j,k)=localTij(i,j) + DO j=1-Oly,sNy+Oly + DO i=1-Olx,sNx+Olx + localTijk(i,j,k)=localTij(i,j) + ENDDO ENDDO - ENDDO ENDIF +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + diagName = 'ADVx'//diagSufx + CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid) + diagName = 'ADVy'//diagSufx + CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid) + ENDIF +#endif + +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevB + & .AND. tracerIdentity.EQ.GAD_TEMPERATURE + & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0 + & .AND. nPx.EQ.1 .AND. nPy.EQ.1 + & .AND. useCubedSphereExchange ) THEN + CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION', + & afx,afy, k, standardMessageUnit,bi,bj,myThid ) + ENDIF +#endif /* ALLOW_DEBUG */ + C-- End of K loop for horizontal fluxes ENDDO +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + IF ( .NOT.implicitAdvection ) THEN C-- Start of k loop for vertical flux DO k=Nr,1,-1 #ifdef ALLOW_AUTODIFF_TAMC - kkey = (igadkey-1)*Nr + k + kkey = (igadkey-1)*Nr + (Nr-k+1) #endif /* ALLOW_AUTODIFF_TAMC */ -C-- kup Cycles through 1,2 to point to w-layer above +C-- kUp Cycles through 1,2 to point to w-layer above C-- kDown Cycles through 2,1 to point to w-layer below - kup = 1+MOD(k+1,2) + kUp = 1+MOD(k+1,2) kDown= 1+MOD(k,2) c kp1=min(Nr,k+1) kp1Msk=1. if (k.EQ.Nr) kp1Msk=0. +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE rtrans(:,:) = +CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte +cphCADJ STORE wfld(:,:) = +cphCADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte +#endif + C-- Compute Vertical transport - IF (k.EQ.1) THEN +#ifdef ALLOW_AIM +C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr + IF ( k.EQ.1 .OR. + & (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr) + & ) THEN +#else + IF ( k.EQ.1 ) THEN +#endif + +#ifdef ALLOW_AUTODIFF_TAMC +cphmultiCADJ STORE wfld(:,:) = +cphmultiCADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ C- Surface interface : DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx - rTransKp1(i,j) = rTrans(i,j) + rTransKp1(i,j) = kp1Msk*rTrans(i,j) + wFld(i,j) = 0. rTrans(i,j) = 0. fVerT(i,j,kUp) = 0. - af(i,j) = 0. ENDDO ENDDO ELSE -C- Interior interface : +#ifdef ALLOW_AUTODIFF_TAMC +cphmultiCADJ STORE wfld(:,:) = +cphmultiCADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ + +C- Interior interface : DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx rTransKp1(i,j) = kp1Msk*rTrans(i,j) + wFld(i,j) = wVel(i,j,k,bi,bj) rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj) + & *deepFac2F(k)*rhoFacF(k) & *maskC(i,j,k-1,bi,bj) - af(i,j) = 0. + fVerT(i,j,kUp) = 0. ENDDO ENDDO #ifdef ALLOW_GMREDI C-- Residual transp = Bolus transp + Eulerian transp - IF (useGMRedi) - & CALL GMREDI_CALC_WFLOW( - & rTrans, bi, bj, k, myThid) + IF (useGMRedi) + & CALL GMREDI_CALC_WFLOW( + U wFld, rTrans, + I k, bi, bj, myThid ) #endif /* ALLOW_GMREDI */ #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE localTijk(:,:,k) -CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte -CADJ STORE rTrans(:,:) -CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte +cphmultiCADJ STORE localTijk(:,:,k) +cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte +cphmultiCADJ STORE rTrans(:,:) +cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C- Compute vertical advective flux in the interior: - IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN - CALL GAD_FLUXLIMIT_ADV_R( - & bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid) - ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN - CALL GAD_DST3_ADV_R( - & bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid) - ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN - CALL GAD_DST3FL_ADV_R( - & bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid) + IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST + & .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN + CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme, + I dTtracerLev(k),rTrans,wFld,localTijk, + O fVerT(1-Olx,1-Oly,kUp), myThid ) + ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN + CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k), + I rTrans, wFld, localTijk, + O fVerT(1-Olx,1-Oly,kUp), myThid ) + ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN + CALL GAD_DST3_ADV_R( bi,bj,k, dTtracerLev(k), + I rTrans, wFld, localTijk, + O fVerT(1-Olx,1-Oly,kUp), myThid ) + ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN + CALL GAD_DST3FL_ADV_R( bi,bj,k, dTtracerLev(k), + I rTrans, wFld, localTijk, + O fVerT(1-Olx,1-Oly,kUp), myThid ) +#ifndef ALLOW_AUTODIFF_TAMC + ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN + CALL GAD_OS7MP_ADV_R( bi,bj,k, dTtracerLev(k), + I rTrans, wFld, localTijk, + O fVerT(1-Olx,1-Oly,kUp), myThid ) +#endif ELSE STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim' ENDIF -C- add the advective flux to fVerT - DO j=1-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx - fVerT(i,j,kUp) = af(i,j) - ENDDO - ENDDO C- end Surface/Interior if bloc ENDIF #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE rTrans(:,:) -CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte -CADJ STORE rTranskp1(:,:) +cphmultiCADJ STORE rTrans(:,:) +cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte +cphmultiCADJ STORE rTranskp1(:,:) +cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte +cph --- following storing of fVerT is critical for correct +cph --- gradient with multiDimAdvection +cph --- Without it, kDown component is not properly recomputed +cph --- This is a TAF bug (and no warning available) +CADJ STORE fVerT(:,:,:) CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C-- Divergence of vertical fluxes DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx - localTij(i,j)=localTijk(i,j,k)-deltaTtracer* - & _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) - & *recip_rA(i,j,bi,bj) - & *( fVerT(i,j,kUp)-fVerT(i,j,kDown) - & -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j)) - & )*rkFac + localTij(i,j) = localTijk(i,j,k) + & -dTtracerLev(k)*recip_rhoFacC(k) + & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k) + & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) + & *( fVerT(i,j,kDown)-fVerT(i,j,kUp) + & -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j)) + & )*rkSign gTracer(i,j,k,bi,bj)= - & (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer + & (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k) ENDDO ENDDO - + +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + diagName = 'ADVr'//diagSufx + CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp), + & diagName, k,1, 2,bi,bj, myThid) + ENDIF +#endif + C-- End of K loop for vertical flux ENDDO C-- end of if not.implicitAdvection block - ENDIF + ENDIF RETURN END