/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_calc_rhs.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.37 by jmc, Sat Oct 22 20:09:52 2005 UTC revision 1.38 by jmc, Sun Nov 6 22:14:02 2005 UTC
# Line 11  C !INTERFACE: ========================== Line 11  C !INTERFACE: ==========================
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I           xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
13       I           uVel, vVel, wVel,       I           uVel, vVel, wVel,
14       I           diffKh, diffK4, KappaR, Tracer,       I           diffKh, diffK4, KappaR, Tracer, TracAB,
15       I           tracerIdentity, advectionScheme, vertAdvecScheme,       I           tracerIdentity, advectionScheme, vertAdvecScheme,
16       I           calcAdvection, implicitAdvection,       I           calcAdvection, implicitAdvection,
17       U           fVerT, gTracer,       U           fVerT, gTracer,
18       I           myTime, myIter, myThid )       I           myTime, myIter, myThid )
19    
20  C !DESCRIPTION:  C !DESCRIPTION:
21  C Calculates the tendancy of a tracer due to advection and diffusion.  C Calculates the tendency of a tracer due to advection and diffusion.
22  C It calculates the fluxes in each direction indepentently and then  C It calculates the fluxes in each direction indepentently and then
23  C sets the tendancy to the divergence of these fluxes. The advective  C sets the tendency to the divergence of these fluxes. The advective
24  C fluxes are only calculated here when using the linear advection schemes  C fluxes are only calculated here when using the linear advection schemes
25  C otherwise only the diffusive and parameterized fluxes are calculated.  C otherwise only the diffusive and parameterized fluxes are calculated.
26  C  C
# Line 29  C \begin{equation*} Line 29  C \begin{equation*}
29  C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}  C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
30  C \end{equation*}  C \end{equation*}
31  C  C
32  C The tendancy is the divergence of the fluxes:  C The tendency is the divergence of the fluxes:
33  C \begin{equation*}  C \begin{equation*}
34  C G_\theta = G_\theta + \nabla \cdot {\bf F}  C G_\theta = G_\theta + \nabla \cdot {\bf F}
35  C \end{equation*}  C \end{equation*}
36  C  C
37  C The tendancy is assumed to contain data on entry.  C The tendency is assumed to contain data on entry.
38    
39  C !USES: ===============================================================  C !USES: ===============================================================
40        IMPLICIT NONE        IMPLICIT NONE
# Line 67  C diffKh           :: horizontal diffusi Line 67  C diffKh           :: horizontal diffusi
67  C diffK4           :: bi-harmonic diffusion coefficient  C diffK4           :: bi-harmonic diffusion coefficient
68  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
69  C Tracer           :: tracer field  C Tracer           :: tracer field
70    C TracAB           :: tracer field but extrapolated fwd in time (to nIter+1/2)
71    C                     if applying AB on tracer field (instead of on tendency)
72  C tracerIdentity   :: tracer identifier (required for KPP,GM)  C tracerIdentity   :: tracer identifier (required for KPP,GM)
73  C advectionScheme  :: advection scheme to use (Horizontal plane)  C advectionScheme  :: advection scheme to use (Horizontal plane)
74  C vertAdvecScheme  :: advection scheme to use (Vertical direction)  C vertAdvecScheme  :: advection scheme to use (Vertical direction)
# Line 90  C myThid           :: thread number Line 92  C myThid           :: thread number
92        _RL diffKh, diffK4        _RL diffKh, diffK4
93        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
95          _RL TracAB(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
96        INTEGER tracerIdentity        INTEGER tracerIdentity
97        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
98        LOGICAL calcAdvection        LOGICAL calcAdvection
# Line 98  C myThid           :: thread number Line 101  C myThid           :: thread number
101        INTEGER myIter, myThid        INTEGER myIter, myThid
102    
103  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
104  C gTracer          :: tendancy array  C gTracer          :: tendency array
105  C fVerT            :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
106        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
107        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
# Line 111  C fMer             :: meridional flux Line 114  C fMer             :: meridional flux
114  C af               :: advective flux  C af               :: advective flux
115  C df               :: diffusive flux  C df               :: diffusive flux
116  C localT           :: local copy of tracer field  C localT           :: local copy of tracer field
117    C locABT           :: local copy of (AB-extrapolated) tracer field
118  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
119        CHARACTER*8 diagName        CHARACTER*8 diagName
120        CHARACTER*4 GAD_DIAG_SUFX, diagSufx        CHARACTER*4 GAD_DIAG_SUFX, diagSufx
# Line 123  C localT           :: local copy of trac Line 127  C localT           :: local copy of trac
127        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130          _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131        _RL advFac, rAdvFac        _RL advFac, rAdvFac
132  CEOP  CEOP
133    
# Line 158  C--   Make local copy of tracer array Line 163  C--   Make local copy of tracer array
163        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
164         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
165          localT(i,j)=Tracer(i,j,k,bi,bj)          localT(i,j)=Tracer(i,j,k,bi,bj)
166            locABT(i,j)=TracAB(i,j,k,bi,bj)
167         ENDDO         ENDDO
168        ENDDO        ENDDO
169    
# Line 189  C--   Initialize net flux in X direction Line 195  C--   Initialize net flux in X direction
195  C-    Advective flux in X  C-    Advective flux in X
196        IF (calcAdvection) THEN        IF (calcAdvection) THEN
197          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
198            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
199          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
200       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
201            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
202       I            dTtracerLev(k), uTrans, uVel, localT,       I            dTtracerLev(k), uTrans, uVel, locABT,
203       O            af, myThid )       O            af, myThid )
204          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
205            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
206       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
207       O            af, myThid )       O            af, myThid )
208          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
209            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
210          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
211            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
212          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
213            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
214       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
215       O            af, myThid )       O            af, myThid )
216          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
217           IF ( inAdMode ) THEN           IF ( inAdMode ) THEN
# Line 213  cph This block is to trick the adjoint: Line 219  cph This block is to trick the adjoint:
219  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
220  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
221            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
222       I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,       I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
223       O           af, myThid )       O           af, myThid )
224           ELSE           ELSE
225            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
226       I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,       I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
227       O           af, myThid )       O           af, myThid )
228           ENDIF           ENDIF
229          ELSE          ELSE
# Line 255  C-    Add bi-harmonic diffusive flux in Line 261  C-    Add bi-harmonic diffusive flux in
261  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
262  C-    GM/Redi flux in X  C-    GM/Redi flux in X
263        IF (useGMRedi) THEN        IF (useGMRedi) THEN
264  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
265          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
266       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
267       I     xA,Tracer,tracerIdentity,       I     xA,Tracer,tracerIdentity,
# Line 289  C--   Initialize net flux in Y direction Line 295  C--   Initialize net flux in Y direction
295  C-    Advective flux in Y  C-    Advective flux in Y
296        IF (calcAdvection) THEN        IF (calcAdvection) THEN
297          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
298            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
299          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
300       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
301            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
302       I            dTtracerLev(k), vTrans, vVel, localT,       I            dTtracerLev(k), vTrans, vVel, locABT,
303       O            af, myThid )       O            af, myThid )
304          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
305            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
306       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
307       O            af, myThid )       O            af, myThid )
308          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
309            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
310          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
311            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
312          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
313            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
314       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
315       O            af, myThid )       O            af, myThid )
316          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
317           IF ( inAdMode ) THEN           IF ( inAdMode ) THEN
# Line 313  cph This block is to trick the adjoint: Line 319  cph This block is to trick the adjoint:
319  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
320  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
321            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
322       I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,       I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
323       O           af, myThid )       O           af, myThid )
324           ELSE           ELSE
325            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
326       I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,       I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
327       O           af, myThid )       O           af, myThid )
328           ENDIF           ENDIF
329          ELSE          ELSE
# Line 355  C-    Add bi-harmonic flux in Y Line 361  C-    Add bi-harmonic flux in Y
361  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
362  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
363        IF (useGMRedi) THEN        IF (useGMRedi) THEN
364  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
365         CALL GMREDI_YTRANSPORT(         CALL GMREDI_YTRANSPORT(
366       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
367       I     yA,Tracer,tracerIdentity,       I     yA,Tracer,tracerIdentity,
# Line 391  C- a hack to prevent Water-Vapor vert.tr Line 397  C- a hack to prevent Water-Vapor vert.tr
397  #endif  #endif
398  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
399          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
400            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
401          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
402       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
403            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
404       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),       I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
405       O         af, myThid )       O         af, myThid )
406          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
407            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
408       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),       I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
409       O         af, myThid )       O         af, myThid )
410          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
411            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
412          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
413            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
414          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
415            CALL GAD_DST3_ADV_R( bi,bj,k,            CALL GAD_DST3_ADV_R( bi,bj,k,
416       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),       I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
417       O         af, myThid )       O         af, myThid )
418          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
419  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
# Line 415  cph IF inAdExact=.FALSE., we want to use Line 421  cph IF inAdExact=.FALSE., we want to use
421  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
422            IF ( inAdMode ) THEN            IF ( inAdMode ) THEN
423             CALL GAD_DST3_ADV_R( bi,bj,k,             CALL GAD_DST3_ADV_R( bi,bj,k,
424       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),       I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
425       O         af, myThid )       O         af, myThid )
426            ELSE            ELSE
427             CALL GAD_DST3FL_ADV_R( bi,bj,k,             CALL GAD_DST3FL_ADV_R( bi,bj,k,
428       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),       I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
429       O         af, myThid )       O         af, myThid )
430            ENDIF            ENDIF
431          ELSE          ELSE
Legend:
Removed from v.1.37  
changed lines
  Added in v.1.38
  ViewVC Help
Powered by ViewVC 1.1.22