/[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.63 by jmc, Fri Dec 27 15:52:17 2013 UTC revision 1.73 by rpa, Wed Jun 3 13:39:22 2015 UTC
# Line 47  C !USES: =============================== Line 47  C !USES: ===============================
47  #include "GRID.h"  #include "GRID.h"
48  #include "SURFACE.h"  #include "SURFACE.h"
49  #include "GAD.h"  #include "GAD.h"
   
50  #ifdef ALLOW_AUTODIFF  #ifdef ALLOW_AUTODIFF
51  # include "AUTODIFF_PARAMS.h"  # include "AUTODIFF_PARAMS.h"
52  #endif /* ALLOW_AUTODIFF */  #endif /* ALLOW_AUTODIFF */
 #ifdef ALLOW_AUTODIFF_TAMC  
 # include "tamc.h"  
 # include "tamc_keys.h"  
 #endif /* ALLOW_AUTODIFF_TAMC */  
53    
54  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
55  C bi,bj            :: tile indices  C bi, bj           :: tile indices
56  C iMin,iMax        :: loop range for called routines  C iMin, iMax       :: for called routines, to get valid output "gTracer"
57  C jMin,jMax        :: loop range for called routines  C jMin, jMax       ::                      over this range of indices
58  C k                :: vertical index  C k                :: vertical index
59  C kM1              :: =k-1 for k>1, =1 for k=1  C kM1              :: =k-1 for k>1, =1 for k=1
60  C kUp              :: index into 2 1/2D array, toggles between 1|2  C kUp              :: index into 2 1/2D array, toggles between 1|2
61  C kDown            :: index into 2 1/2D array, toggles between 2|1  C kDown            :: index into 2 1/2D array, toggles between 2|1
62  C xA,yA            :: areas of X and Y face of tracer cells  C xA, yA           :: areas of X and Y face of tracer cells
63  C maskUp           :: 2-D array for mask at W points  C maskUp           :: 2-D array for mask at W points
64  C uFld,vFld,wFld   :: Local copy of velocity field (3 components)  C uFld, vFld, wFld :: Local copy of velocity field (3 components)
65  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points  C uTrans, vTrans   :: 2-D arrays of volume transports at U,V points
66  C rTrans           :: 2-D arrays of volume transports at W points  C rTrans           :: 2-D arrays of volume transports at W points
67  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
68  C diffKh           :: horizontal diffusion coefficient  C diffKh           :: horizontal diffusion coefficient
# Line 105  C myThid           :: thread number Line 100  C myThid           :: thread number
100        _RL diffKh, diffK4        _RL diffKh, diffK4
101        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL diffKr4(Nr)        _RL diffKr4(Nr)
103        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
105        _RL deltaTLev(Nr)        _RL deltaTLev(Nr)
106        INTEGER trIdentity        INTEGER trIdentity
107        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
# Line 121  C gTracer          :: tendency array Line 116  C gTracer          :: tendency array
116  C fZon             :: zonal flux  C fZon             :: zonal flux
117  C fMer             :: meridional flux  C fMer             :: meridional flux
118  C fVerT            :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
119        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
120        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
# Line 157  C locABT           :: local copy of (AB- Line 152  C locABT           :: local copy of (AB-
152  #endif  #endif
153  CEOP  CEOP
154    
155  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
156  C--   only the kUp part of fverT is set in this subroutine  C--   only the kUp part of fverT is set in this subroutine
157  C--   the kDown is still required  C--   the kDown is still required
158        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
# Line 189  C--   Make local copy of tracer array Line 184  C--   Make local copy of tracer array
184        IF ( applyAB_onTracer ) THEN        IF ( applyAB_onTracer ) THEN
185          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
186           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
187            localT(i,j)=TracerN(i,j,k,bi,bj)            localT(i,j)=TracerN(i,j,k)
188            locABT(i,j)= TracAB(i,j,k,bi,bj)            locABT(i,j)= TracAB(i,j,k)
189           ENDDO           ENDDO
190          ENDDO          ENDDO
191        ELSE        ELSE
192          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
193           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
194            localT(i,j)= TracAB(i,j,k,bi,bj)            localT(i,j)=TracerN(i,j,k)
195            locABT(i,j)= TracAB(i,j,k,bi,bj)            locABT(i,j)=TracerN(i,j,k)
196           ENDDO           ENDDO
197          ENDDO          ENDDO
198        ENDIF        ENDIF
# Line 262  cph with limiters in forward, but withou Line 257  cph with limiters in forward, but withou
257             CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),             CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
258       I              uTrans, uFld, maskLocW, locABT,       I              uTrans, uFld, maskLocW, locABT,
259       O              af, myThid )       O              af, myThid )
260  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF
261           ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN           ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
262             CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),             CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
263       I              uTrans, uFld, maskLocW, locABT,       I              uTrans, uFld, maskLocW, locABT,
# Line 291  C-      replace advective flux with 1st Line 286  C-      replace advective flux with 1st
286            diagName = 'ADVx'//diagSufx            diagName = 'ADVx'//diagSufx
287            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
288          ENDIF          ENDIF
289    #ifdef ALLOW_LAYERS
290            IF ( useLayers ) THEN
291              CALL LAYERS_FILL( af, trIdentity, 'AFX',
292         &                      k, 1, 2,bi,bj, myThid )
293            ENDIF
294    #endif /* ENDIF */
295  #endif  #endif
296        ENDIF        ENDIF
297    
# Line 313  C-    Add bi-harmonic diffusive flux in Line 314  C-    Add bi-harmonic diffusive flux in
314  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
315  C-    GM/Redi flux in X  C-    GM/Redi flux in X
316        IF ( trUseGMRedi ) THEN        IF ( trUseGMRedi ) THEN
317  C *note* should update GMREDI_XTRANSPORT to set df  *aja*          CALL GMREDI_XTRANSPORT(
318          IF ( applyAB_onTracer ) THEN       I         trIdentity, bi, bj, k, iMin, iMax+1, jMin, jMax,
319            CALL GMREDI_XTRANSPORT(       I         xA, TracerN,
      I         iMin,iMax,jMin,jMax,bi,bj,k,  
      I         xA,TracerN,trIdentity,  
      U         df,  
      I         myThid)  
         ELSE  
           CALL GMREDI_XTRANSPORT(  
      I         iMin,iMax,jMin,jMax,bi,bj,k,  
      I         xA,TracAB, trIdentity,  
320       U         df,       U         df,
321       I         myThid)       I         myThid )
         ENDIF  
322        ENDIF        ENDIF
323  #endif  #endif
324  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
# Line 343  C       excluding advective terms: Line 335  C       excluding advective terms:
335       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
336            diagName = 'DFxE'//diagSufx            diagName = 'DFxE'//diagSufx
337            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
338    #ifdef ALLOW_LAYERS
339              IF ( useLayers ) THEN
340               CALL LAYERS_FILL( df, trIdentity, 'DFX',
341         &                       k, 1, 2,bi,bj, myThid )
342              ENDIF
343    #endif /* ALLOW_LAYERS */
344        ENDIF        ENDIF
345  #endif  #endif
346    
# Line 399  cph with limiters in forward, but withou Line 397  cph with limiters in forward, but withou
397             CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),             CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
398       I              vTrans, vFld, maskLocS, locABT,       I              vTrans, vFld, maskLocS, locABT,
399       O              af, myThid )       O              af, myThid )
400  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF
401           ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN           ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
402             CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),             CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
403       I              vTrans, vFld, maskLocS, locABT,       I              vTrans, vFld, maskLocS, locABT,
# Line 428  C-      replace advective flux with 1st Line 426  C-      replace advective flux with 1st
426            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
427            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
428          ENDIF          ENDIF
429    #ifdef ALLOW_LAYERS
430            IF ( useLayers ) THEN
431              CALL LAYERS_FILL( af, trIdentity, 'AFY',
432         &                          k, 1, 2,bi,bj, myThid )
433            ENDIF
434    #endif /* ALLOW_LAYES */
435  #endif  #endif
436        ENDIF        ENDIF
437    
# Line 450  C-    Add bi-harmonic flux in Y Line 454  C-    Add bi-harmonic flux in Y
454  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
455  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
456        IF ( trUseGMRedi ) THEN        IF ( trUseGMRedi ) THEN
457  C *note* should update GMREDI_YTRANSPORT to set df  *aja*          CALL GMREDI_YTRANSPORT(
458          IF ( applyAB_onTracer ) THEN       I         trIdentity, bi, bj, k, iMin, iMax, jMin, jMax+1,
459            CALL GMREDI_YTRANSPORT(       I         yA, TracerN,
      I         iMin,iMax,jMin,jMax,bi,bj,k,  
      I         yA,TracerN,trIdentity,  
      U         df,  
      I         myThid)  
         ELSE  
           CALL GMREDI_YTRANSPORT(  
      I         iMin,iMax,jMin,jMax,bi,bj,k,  
      I         yA,TracAB, trIdentity,  
460       U         df,       U         df,
461       I         myThid)       I         myThid )
         ENDIF  
462        ENDIF        ENDIF
463  #endif  #endif
464  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
# Line 480  C       excluding advective terms: Line 475  C       excluding advective terms:
475       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
476            diagName = 'DFyE'//diagSufx            diagName = 'DFyE'//diagSufx
477            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
478    #ifdef ALLOW_LAYERS
479              IF ( useLayers ) THEN
480               CALL LAYERS_FILL( df, trIdentity, 'DFY',
481         &                           k, 1, 2,bi,bj, myThid )
482              ENDIF
483    #endif /* ALLOW_LAYERS */
484        ENDIF        ENDIF
485  #endif  #endif
486    
# Line 493  C- a hack to prevent Water-Vapor vert.tr Line 494  C- a hack to prevent Water-Vapor vert.tr
494  #else  #else
495        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
496  #endif  #endif
497  C-    Compute vertical advective flux in the interior:         IF ( applyAB_onTracer ) THEN
498    C-    Compute vertical advective flux in the interior using TracAB:
499          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
500             CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )             CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
501          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
502       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
503             CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),             CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
504       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),       I              rTrans, wFld, TracAB,
505       O              af, myThid )       O              af, myThid )
506          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
507             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
508       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),       I              rTrans, wFld, TracAB,
509       O              af, myThid )       O              af, myThid )
510          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
511             CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )             CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
# Line 520  cph with limiters in forward, but withou Line 522  cph with limiters in forward, but withou
522          ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN          ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
523  #endif /* ALLOW_AUTODIFF */  #endif /* ALLOW_AUTODIFF */
524             CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),             CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
525       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),       I              rTrans, wFld, TracAB,
526       O              af, myThid )       O              af, myThid )
527          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
528             CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),             CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
529       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),       I              rTrans, wFld, TracAB,
530       O              af, myThid )       O              af, myThid )
531  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF
532          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
533             CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),             CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
534       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),       I              rTrans, wFld, TracAB,
535       O              af, myThid )       O              af, myThid )
536  #endif  #endif
537          ELSE          ELSE
538            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
539          ENDIF          ENDIF
540           ELSE
541    C-    Compute vertical advective flux in the interior using TracerN:
542            IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
543               CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
544            ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
545         &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
546               CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
547         I              rTrans, wFld, TracerN,
548         O              af, myThid )
549            ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
550               CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
551         I              rTrans, wFld, TracerN,
552         O              af, myThid )
553            ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
554               CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
555            ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
556               CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
557    #ifdef ALLOW_AUTODIFF
558            ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 .OR.
559         &         (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT .AND. inAdMode)
560         &        ) THEN
561    cph This block is to trick the adjoint:
562    cph If inAdExact=.FALSE., we want to use DST3
563    cph with limiters in forward, but without limiters in reverse.
564    #else /* ALLOW_AUTODIFF */
565            ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
566    #endif /* ALLOW_AUTODIFF */
567               CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
568         I              rTrans, wFld, TracerN,
569         O              af, myThid )
570            ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
571               CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
572         I              rTrans, wFld, TracerN,
573         O              af, myThid )
574    #ifndef ALLOW_AUTODIFF
575            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
576               CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
577         I              rTrans, wFld, TracerN,
578         O              af, myThid )
579    #endif
580            ELSE
581              STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
582            ENDIF
583           ENDIF
584  C-     add the advective flux to fVerT  C-     add the advective flux to fVerT
585          DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
586           DO i=1-OLx,sNx+OLx          DO i=1-OLx,sNx+OLx
587            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)*maskInC(i,j,bi,bj)            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)*maskInC(i,j,bi,bj)
          ENDDO  
588          ENDDO          ENDDO
589           ENDDO
590  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
591          IF ( useDiagnostics ) THEN         IF ( useDiagnostics ) THEN
592            diagName = 'ADVr'//diagSufx            diagName = 'ADVr'//diagSufx
593            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
594  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
595  C        does it only if k=1 (never the case here)  C        does it only if k=1 (never the case here)
596            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
597          ENDIF  #ifdef ALLOW_LAYERS
598              IF ( useLayers ) THEN
599                CALL LAYERS_FILL(af,trIdentity,'AFR',k,1,2,bi,bj,myThid)
600              ENDIF
601    #endif /* ALLOW_LAYERS */
602           ENDIF
603  #endif  #endif
604        ENDIF        ENDIF
605    
# Line 562  C           boundary condition. Line 613  C           boundary condition.
613          ENDDO          ENDDO
614         ENDDO         ENDDO
615        ELSE        ELSE
616         IF ( applyAB_onTracer ) THEN          CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
          CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)  
        ELSE  
          CALL GAD_DIFF_R(bi,bj,k,KappaR,TracAB, df,myThid)  
        ENDIF  
617        ENDIF        ENDIF
618    
619        IF ( trUseDiffKr4 ) THEN        IF ( trUseDiffKr4 ) THEN
620         IF ( applyAB_onTracer ) THEN          CALL GAD_BIHARM_R( bi,bj,k, diffKr4, TracerN, df, myThid )
          CALL GAD_BIHARM_R( bi,bj,k, diffKr4, TracerN, df, myThid )  
        ELSE  
          CALL GAD_BIHARM_R( bi,bj,k, diffKr4, TracAB,  df, myThid )  
        ENDIF  
621        ENDIF        ENDIF
622    
623  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
624  C-    GM/Redi flux in R  C-    GM/Redi flux in R
625        IF ( trUseGMRedi ) THEN        IF ( trUseGMRedi ) THEN
626  C *note* should update GMREDI_RTRANSPORT to set df  *aja*          CALL GMREDI_RTRANSPORT(
627          IF ( applyAB_onTracer ) THEN       I         trIdentity, bi, bj, k, iMin, iMax, jMin, jMax,
628            CALL GMREDI_RTRANSPORT(       I         TracerN,
      I         iMin,iMax,jMin,jMax,bi,bj,k,  
      I         TracerN,trIdentity,  
      U         df,  
      I         myThid)  
         ELSE  
           CALL GMREDI_RTRANSPORT(  
      I         iMin,iMax,jMin,jMax,bi,bj,k,  
      I         TracAB, trIdentity,  
629       U         df,       U         df,
630       I         myThid)       I         myThid )
         ENDIF  
631        ENDIF        ENDIF
632  #endif  #endif
633    
# Line 610  C       Explicit terms only & excluding Line 644  C       Explicit terms only & excluding
644       &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN       &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
645            diagName = 'DFrE'//diagSufx            diagName = 'DFrE'//diagSufx
646            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
647    #ifdef ALLOW_LAYERS
648              IF ( useLayers ) THEN
649               CALL LAYERS_FILL(df,trIdentity,'DFR',k,1,2,bi,bj,myThid)
650              ENDIF
651    #endif /* ALLOW_LAYERS */
652        ENDIF        ENDIF
653  #endif  #endif
654    
# Line 658  C-    Diagnostics of Non-Local Tracer (v Line 697  C-    Diagnostics of Non-Local Tracer (v
697  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
698  C        does it only if k=1 (never the case here)  C        does it only if k=1 (never the case here)
699           IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)           IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
700    #ifdef ALLOW_LAYERS
701             IF ( useLayers ) THEN
702              CALL LAYERS_FILL(df,trIdentity,'DFR',k,1,2,bi,bj,myThid)
703             ENDIF
704    #endif /* ALLOW_LAYERS */
705         ENDIF         ENDIF
706  #endif  #endif
707        ENDIF        ENDIF
# Line 682  coj   Add outgoing fluxes Line 726  coj   Add outgoing fluxes
726       &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)       &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
727       &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)       &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
728       &     )       &     )
729           IF ( applyAB_onTracer ) THEN           trac = localT(i,j)
            trac=TracerN(i,j,k,bi,bj)  
          ELSE  
            trac=TracAB(i,j,k,bi,bj)  
          ENDIF  
730  coj   If they would reduce tracer by a fraction of more than  coj   If they would reduce tracer by a fraction of more than
731  coj   SmolarkiewiczMaxFrac, scale them down  coj   SmolarkiewiczMaxFrac, scale them down
732           IF (outFlux.GT.0. _d 0 .AND.           IF (outFlux.GT.0. _d 0 .AND.
# Line 704  coj   If tracer is already negative, sca Line 744  coj   If tracer is already negative, sca
744             IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN             IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
745  coj   Down flux is special: it has already been applied in lower layer,  coj   Down flux is special: it has already been applied in lower layer,
746  coj   so we have to readjust this.  coj   so we have to readjust this.
747  coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!  coj   undo down flux, ...
748  coj   thus it has an extra factor deltaTLev(k+1)               gTracer(i,j,k+1) = gTracer(i,j,k+1)
749               gTrFac=deltaTLev(k+1)       &        +_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
 coj   Other factors that have been applied to gTracer since the last call:  
 #ifdef NONLIN_FRSURF  
              IF (nonlinFreeSurf.GT.0) THEN  
               IF (select_rStar.GT.0) THEN  
 #ifndef DISABLE_RSTAR_CODE  
                 gTrFac = gTrFac/rStarExpC(i,j,bi,bj)  
 #endif /* DISABLE_RSTAR_CODE */  
               ENDIF  
              ENDIF  
 #endif /* NONLIN_FRSURF */  
 coj   Now: undo down flux, ...  
              gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)  
      &        +gTrFac  
      &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)  
750       &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)       &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
751       &         *recip_rhoFacC(k+1)       &         *recip_rhoFacC(k+1)
752       &         *( -fVerT(i,j,kDown)*rkSign )       &         *( -fVerT(i,j,kDown)*rkSign )
753  coj   ... scale ...  coj   ... scale ...
754               fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)               fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
755  coj   ... and reapply  coj   ... and reapply
756               gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)               gTracer(i,j,k+1) = gTracer(i,j,k+1)
757       &        +gTrFac       &        +_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
      &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)  
758       &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)       &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
759       &         *recip_rhoFacC(k+1)       &         *recip_rhoFacC(k+1)
760       &         *( fVerT(i,j,kDown)*rkSign )       &         *( fVerT(i,j,kDown)*rkSign )
# Line 746  C     Anelastic: scale vertical fluxes b Line 771  C     Anelastic: scale vertical fluxes b
771  C     for Stevens OBC: keep only vertical diffusive contribution on boundaries  C     for Stevens OBC: keep only vertical diffusive contribution on boundaries
772        DO j=1-OLy,sNy+OLy-1        DO j=1-OLy,sNy+OLy-1
773         DO i=1-OLx,sNx+OLx-1         DO i=1-OLx,sNx+OLx-1
774          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k) = gTracer(i,j,k)
775       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
776       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
777       &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)       &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
Legend:
Removed from v.1.63  
changed lines
  Added in v.1.73
  ViewVC Help
Powered by ViewVC 1.1.22