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

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

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

revision 1.10 by jmc, Wed Jun 22 00:27:47 2005 UTC revision 1.23 by rpa, Wed Jun 3 13:39:22 2015 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  CBOP  CBOP
7  C     !ROUTINE: GAD_IMPLICIT_R  C     !ROUTINE: GAD_IMPLICIT_R
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE GAD_IMPLICIT_R(        SUBROUTINE GAD_IMPLICIT_R(
10       I      implicitAdvection, advectionScheme, tracerIdentity,       I      implicitAdvection, advectionScheme, tracerIdentity,
11       I      kappaRX, wVel, tracer,       I      deltaTLev,
12         I      kappaRX, recip_hFac, wFld, tracer,
13       U      gTracer,       U      gTracer,
14       I      bi, bj, myTime, myIter, myThid )       I      bi, bj, myTime, myIter, myThid )
15  C     !DESCRIPTION:  C     !DESCRIPTION:
# Line 21  C     == Global data == Line 22  C     == Global data ==
22  #include "EEPARAMS.h"  #include "EEPARAMS.h"
23  #include "PARAMS.h"  #include "PARAMS.h"
24  #include "GRID.h"  #include "GRID.h"
25    #include "SURFACE.h"
26  #include "GAD.h"  #include "GAD.h"
27    
28  C !INPUT/OUTPUT PARAMETERS:  C !INPUT/OUTPUT PARAMETERS:
# Line 29  C implicitAdvection :: if True, treat ve Line 31  C implicitAdvection :: if True, treat ve
31  C advectionScheme   :: advection scheme to use  C advectionScheme   :: advection scheme to use
32  C tracerIdentity    :: Identifier for the tracer  C tracerIdentity    :: Identifier for the tracer
33  C kappaRX           :: 3-D array for vertical diffusion coefficient  C kappaRX           :: 3-D array for vertical diffusion coefficient
34  C wVel              :: vertical component of the velcity field  C recip_hFac        :: inverse of cell open-depth factor
35    C wFld              :: Advection velocity field, vertical component
36  C tracer            :: tracer field at current time step  C tracer            :: tracer field at current time step
37  C gTracer           :: future tracer field  C gTracer           :: future tracer field
38  C bi,bj             :: tile indices  C bi,bj             :: tile indices
# Line 39  C myThid            :: thread number Line 42  C myThid            :: thread number
42        LOGICAL implicitAdvection        LOGICAL implicitAdvection
43        INTEGER advectionScheme        INTEGER advectionScheme
44        INTEGER tracerIdentity        INTEGER tracerIdentity
45        _RL kappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL     deltaTLev(Nr)
46        _RL wVel   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL kappaRX   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
47        _RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
48        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL wFld      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
49          _RL tracer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
50          _RL gTracer   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51        INTEGER bi, bj        INTEGER bi, bj
52        _RL myTime        _RL myTime
53        INTEGER myIter, myThid        INTEGER myIter, myThid
54    
55    #ifdef ALLOW_DIAGNOSTICS
56    C !FUNCTIONS:
57          CHARACTER*4 GAD_DIAG_SUFX
58          EXTERNAL    GAD_DIAG_SUFX
59          LOGICAL     DIAGNOSTICS_IS_ON
60          EXTERNAL    DIAGNOSTICS_IS_ON
61    #endif
62    
63  C !LOCAL VARIABLES:  C !LOCAL VARIABLES:
64  C == Local variables ==  C == Local variables ==
65  C iMin,iMax,jMin,jMax  :: computational domain  C iMin,iMax,jMin,jMax :: computational domain
66  C i,j,k  :: loop indices  C i,j,k     :: loop indices
67  C a5d    :: 2nd  lower diagonal of the pentadiagonal matrix  C a5d       :: 2nd  lower diagonal of the pentadiagonal matrix
68  C b5d    :: 1rst lower diagonal of the pentadiagonal matrix  C b5d       :: 1rst lower diagonal of the pentadiagonal matrix
69  C c5d    :: main diagonal       of the pentadiagonal matrix  C c5d       :: main diagonal       of the pentadiagonal matrix
70  C d5d    :: 1rst upper diagonal of the pentadiagonal matrix  C d5d       :: 1rst upper diagonal of the pentadiagonal matrix
71  C e5d    :: 2nd  upper diagonal of the pentadiagonal matrix  C e5d       :: 2nd  upper diagonal of the pentadiagonal matrix
72  C rTrans    :: vertical volume transport at inteface k  C rTrans    :: vertical volume transport at interface k
 C rTransKp1 :: vertical volume transport at inteface k+1  
73  C localTijk :: local copy of tracer (for Non-Lin Adv.Scheme)  C localTijk :: local copy of tracer (for Non-Lin Adv.Scheme)
74  C diagonalNumber :: number of non-zero diagonals in the matrix  C diagonalNumber :: number of non-zero diagonals in the matrix
75  C errCode   :: > 0 if singular matrix  C errCode   :: > 0 if singular matrix
76        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
77          PARAMETER( iMin = 1, iMax = sNx )
78          PARAMETER( jMin = 1, jMax = sNy )
79        INTEGER i,j,k        INTEGER i,j,k
80        INTEGER diagonalNumber, errCode        INTEGER diagonalNumber, errCode
81        _RL a5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL a5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
82        _RL b5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL b5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
83        _RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL c5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
84        _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL d5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
85        _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL e5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86        _RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rTrans   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87        _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
88  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
89        CHARACTER*8 diagName        CHARACTER*8 diagName
90        CHARACTER*4 GAD_DIAG_SUFX, diagSufx        CHARACTER*4 diagSufx
91        EXTERNAL    GAD_DIAG_SUFX        LOGICAL     diagDif, diagAdv
92        LOGICAL     DIAGNOSTICS_IS_ON        INTEGER km1, km2, kp1, kp2
93        EXTERNAL    DIAGNOSTICS_IS_ON        _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL df (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95          _RL div(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96          _RL flx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97  #endif  #endif
98  CEOP  CEOP
99    
# Line 86  C--   no need to solve anything with onl Line 101  C--   no need to solve anything with onl
101        IF (Nr.GT.1) THEN        IF (Nr.GT.1) THEN
102    
103  C--   Initialise  C--   Initialise
       iMin = 1  
       jMin = 1  
       iMax = sNx  
       jMax = sNy  
104        DO k=1,Nr        DO k=1,Nr
105         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
106          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
107           a5d(i,j,k) = 0. _d 0           a5d(i,j,k) = 0. _d 0
108           b5d(i,j,k) = 0. _d 0           b5d(i,j,k) = 0. _d 0
109           c5d(i,j,k) = 1. _d 0           c5d(i,j,k) = 1. _d 0
# Line 104  C--   Initialise Line 115  C--   Initialise
115        diagonalNumber = 1        diagonalNumber = 1
116    
117  C--   Non-Linear Advection scheme: keep a local copy of tracer field  C--   Non-Linear Advection scheme: keep a local copy of tracer field
118        IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.        IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
119       &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
120          IF ( multiDimAdvection ) THEN          IF ( multiDimAdvection ) THEN
121           DO k=1,Nr           DO k=1,Nr
122            DO j=1-Oly,sNy+Oly            DO j=1-OLy,sNy+OLy
123             DO i=1-Olx,sNx+Olx             DO i=1-OLx,sNx+OLx
124              localTijk(i,j,k) = gTracer(i,j,k,bi,bj)              localTijk(i,j,k) = gTracer(i,j,k)
125             ENDDO             ENDDO
126            ENDDO            ENDDO
127           ENDDO           ENDDO
128          ELSE          ELSE
129           DO k=1,Nr           DO k=1,Nr
130            DO j=1-Oly,sNy+Oly            DO j=1-OLy,sNy+OLy
131             DO i=1-Olx,sNx+Olx             DO i=1-OLx,sNx+OLx
132              localTijk(i,j,k) = tracer(i,j,k,bi,bj)              localTijk(i,j,k) = tracer(i,j,k)
133             ENDDO             ENDDO
134            ENDDO            ENDDO
135           ENDDO           ENDDO
# Line 132  C-     1rst lower diagonal : Line 143  C-     1rst lower diagonal :
143         DO k=2,Nr         DO k=2,Nr
144          DO j=jMin,jMax          DO j=jMin,jMax
145           DO i=iMin,iMax           DO i=iMin,iMax
146             b5d(i,j,k) = -dTtracerLev(k)*maskC(i,j,k-1,bi,bj)             b5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k-1,bi,bj)
147       &                  *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                  *recip_hFac(i,j,k)*recip_drF(k)
148       &                  *kappaRX(i,j, k )*recip_drC( k )       &                  *kappaRX(i,j, k )*recip_drC( k )
149           ENDDO           ENDDO
150          ENDDO          ENDDO
# Line 142  C-     1rst upper diagonal : Line 153  C-     1rst upper diagonal :
153         DO k=1,Nr-1         DO k=1,Nr-1
154          DO j=jMin,jMax          DO j=jMin,jMax
155           DO i=iMin,iMax           DO i=iMin,iMax
156             d5d(i,j,k) = -dTtracerLev(k)*maskC(i,j,k+1,bi,bj)             d5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k+1,bi,bj)
157       &                 *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                 *recip_hFac(i,j,k)*recip_drF(k)
158       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)
159           ENDDO           ENDDO
160          ENDDO          ENDDO
# Line 165  C--   end if implicitDiffusion Line 176  C--   end if implicitDiffusion
176         DO k=Nr,1,-1         DO k=Nr,1,-1
177    
178  C--    Compute transport  C--    Compute transport
         IF (k.EQ.Nr) THEN  
          DO j=1-Oly,sNy+Oly  
           DO i=1-Olx,sNx+Olx  
             rTransKp1(i,j) = 0.  
           ENDDO  
          ENDDO  
         ELSE  
          DO j=1-Oly,sNy+Oly  
           DO i=1-Olx,sNx+Olx  
             rTransKp1(i,j) = rTrans(i,j)  
           ENDDO  
          ENDDO  
         ENDIF  
   
179          IF (k.EQ.1) THEN          IF (k.EQ.1) THEN
180           DO j=1-Oly,sNy+Oly           DO j=1-OLy,sNy+OLy
181            DO i=1-Olx,sNx+Olx            DO i=1-OLx,sNx+OLx
182              rTrans(i,j) = 0.              rTrans(i,j) = 0. _d 0
183            ENDDO            ENDDO
184           ENDDO           ENDDO
185          ELSE          ELSE
186           DO j=1-Oly,sNy+Oly           DO j=1-OLy,sNy+OLy
187            DO i=1-Olx,sNx+Olx            DO i=1-OLx,sNx+OLx
188              rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)              rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj)
189       &                   *maskC(i,j,k-1,bi,bj)       &                               *maskC(i,j,k-1,bi,bj)
190            ENDDO            ENDDO
191           ENDDO           ENDDO
 #ifdef ALLOW_GMREDI  
 C--   Residual transp = Bolus transp + Eulerian transp  
           IF (useGMRedi)  
      &      CALL GMREDI_CALC_WFLOW(  
      &                    rTrans, bi, bj, k, myThid)  
 #endif /* ALLOW_GMREDI */  
192          ENDIF          ENDIF
         DO j=jMin,jMax  
           DO i=iMin,iMax  
 c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)  
            gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)  
      &      + dTtracerLev(1)*recip_rA(i,j,bi,bj)  
      &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  
      &       *tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))*rkSign  
           ENDDO  
         ENDDO  
193    
194  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
195  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
196          IF ( K.GE.2 .AND.          IF ( k.GE.2 .AND.
197       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
198       &              ) THEN       &              ) THEN
199  #else  #else
200          IF ( K.GE.2 ) THEN          IF ( k.GE.2 ) THEN
201  #endif  #endif
202    
203           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
204            diagonalNumber = 3            diagonalNumber = 3
205            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
206       I                        dTtracerLev, rTrans,       I                        deltaTLev, rTrans, recip_hFac,
207       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
208       I                        myThid)       I                        myThid )
209           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
210         &       .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
211              diagonalNumber = 3
212              CALL GAD_DST2U1_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
213         I                        advectionScheme, deltaTLev,
214         I                        rTrans, recip_hFac,
215         U                        b5d, c5d, d5d,
216         I                        myThid )
217             ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
218            diagonalNumber = 3            diagonalNumber = 3
219            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
220       I                        dTtracerLev, rTrans, localTijk,       I                        deltaTLev, rTrans, recip_hFac, localTijk,
221       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
222       I                        myThid)       I                        myThid )
223           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.           ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD
224       &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       &       .OR. advectionScheme.EQ.ENUM_CENTERED_4TH
225         &       .OR. advectionScheme.EQ.ENUM_DST3 ) THEN
226            diagonalNumber = 5            diagonalNumber = 5
227            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
228       I                        advectionScheme, dTtracerLev, rTrans,       I                        advectionScheme, deltaTLev,
229         I                        rTrans, recip_hFac,
230         U                        a5d, b5d, c5d, d5d, e5d,
231         I                        myThid )
232             ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
233              diagonalNumber = 5
234              CALL GAD_DST3FL_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
235         I                        deltaTLev, rTrans, recip_hFac, localTijk,
236       U                        a5d, b5d, c5d, d5d, e5d,       U                        a5d, b5d, c5d, d5d, e5d,
237       I                        myThid)       I                        myThid )
238           ELSE           ELSE
239            STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'            STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
240           ENDIF           ENDIF
# Line 275  C--   Solve penta-diagonal system : Line 273  C--   Solve penta-diagonal system :
273    
274  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
275  C--   Set diagnostic suffix for the current tracer  C--   Set diagnostic suffix for the current tracer
276        IF ( useDiagnostics .AND. implicitDiffusion ) THEN        IF ( useDiagnostics ) THEN
277          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
278          diagName = 'DFrI'//diagSufx          diagName = 'DFrI'//diagSufx
279          IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN          diagDif = implicitDiffusion
280           DO k= 1,Nr          IF ( diagDif ) diagDif = DIAGNOSTICS_IS_ON(diagName,myThid)
281            IF ( k.EQ.1 ) THEN          diagName = 'ADVr'//diagSufx
282  C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0          diagAdv = implicitAdvection
283  C         otherwise counter is not incremented !!          IF ( diagAdv ) diagAdv = DIAGNOSTICS_IS_ON(diagName,myThid)
284    
285            IF ( diagDif .OR. diagAdv ) THEN
286             DO j=1-OLy,sNy+OLy
287              DO i=1-OLx,sNx+OLx
288                flx(i,j) = 0. _d 0
289              ENDDO
290             ENDDO
291             DO k= Nr,1,-1
292              IF ( implicitDiffusion .AND. k.GE.2 ) THEN
293                DO j=jMin,jMax
294                 DO i=iMin,iMax
295                   df(i,j) =
296    cc#ifdef ALLOW_AUTODIFF_OPENAD
297    cc     &             -rA(i,j,bi,bj)%v
298    cc#else
299         &             -rA(i,j,bi,bj)
300    cc#endif
301         &            * KappaRX(i,j,k)*recip_drC(k)*rkSign
302         &            * (gTracer(i,j,k) - gTracer(i,j,k-1))
303         &            * maskC(i,j,k,bi,bj)
304         &            * maskC(i,j,k-1,bi,bj)
305                 ENDDO
306                ENDDO
307              ELSE
308              DO j=1-OLy,sNy+OLy              DO j=1-OLy,sNy+OLy
309               DO i=1-OLx,sNx+OLx               DO i=1-OLx,sNx+OLx
310                 df(i,j) = 0. _d 0                 df(i,j) = 0. _d 0
311               ENDDO               ENDDO
312              ENDDO              ENDDO
313            ELSE            ENDIF
314              DO j=1,sNy  C-  Note: Needs to explicitly increment counter (call DIAGNOSTICS_COUNT)
315               DO i=1,sNx  C         since skipping k=1 DIAGNOSTICS_FILL call.
316                 df(i,j) =            IF ( diagDif .AND. k.GE.2 ) THEN
317       &              rA(i,j,bi,bj)             diagName = 'DFrI'//diagSufx
318       &            * KappaRX(i,j,k)*recip_drC(k)             CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
319       &            * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))             IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
320    #ifdef ALLOW_LAYERS
321               IF ( useLayers ) THEN
322                 CALL LAYERS_FILL( df, tracerIdentity, 'DFR',
323         &                           k, 1, 2,bi,bj, myThid )
324               ENDIF
325    #endif /* ALLOW_LAYERS */
326              ENDIF
327              IF ( diagAdv ) THEN
328               km1=MAX(1,k-1)
329               km2=MAX(1,k-2)
330               kp1=MIN(Nr,k+1)
331               kp2=MIN(Nr,k+2)
332    C--   Flux_divergence*deltaT = Tr^n - Tr^n+1 = [A-I](Tr^n+1)
333    C                            = deltaT*rkSign*[ Flx_k+1 - Flx_k ]/dz
334               DO j=jMin,jMax
335                DO i=iMin,iMax
336                  div(i,j) = gTracer(i,j,k)*( c5d(i,j,k) - 1. _d 0 )
337         &                 + gTracer(i,j,km1)*b5d(i,j,k)
338         &                 + gTracer(i,j,kp1)*d5d(i,j,k)
339                ENDDO
340               ENDDO
341               IF ( diagonalNumber .EQ. 5 ) THEN
342                DO j=jMin,jMax
343                 DO i=iMin,iMax
344                  div(i,j) = div(i,j)
345         &                 + gTracer(i,j,km2)*a5d(i,j,k)
346         &                 + gTracer(i,j,kp2)*e5d(i,j,k)
347               ENDDO               ENDDO
348              ENDDO              ENDDO
349               ENDIF
350    #ifdef NONLIN_FRSURF
351               IF ( nonlinFreeSurf.GT.0 ) THEN
352    C--    use future hFac to stay consistent with solver matrix
353                IF ( select_rStar.GT.0 ) THEN
354                 DO j=jMin,jMax
355                  DO i=iMin,iMax
356                   div(i,j) = div(i,j)*h0FacC(i,j,k,bi,bj)*drF(k)
357         &                            *rStarFacC(i,j,bi,bj)
358                  ENDDO
359                 ENDDO
360                ELSEIF ( selectSigmaCoord.NE.0 ) THEN
361                 DO j=jMin,jMax
362                  DO i=iMin,iMax
363                   div(i,j) = div(i,j)*(
364         &               _hFacC(i,j,k,bi,bj)*drF(k)
365         &              + dBHybSigF(k)*dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
366         &                             )
367                  ENDDO
368                 ENDDO
369                ELSE
370                 DO j=jMin,jMax
371                  DO i=iMin,iMax
372                   IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
373                    div(i,j) = div(i,j)*hFac_surfC(i,j,bi,bj)*drF(k)
374                   ELSE
375                    div(i,j) = div(i,j)*_hFacC(i,j,k,bi,bj)*drF(k)
376                   ENDIF
377                  ENDDO
378                 ENDDO
379                ENDIF
380               ELSE
381    #else /* NONLIN_FRSURF */
382               IF ( .TRUE. ) THEN
383    #endif /* NONLIN_FRSURF */
384    C--    use current hFac (consistent with solver matrix)
385                 DO j=jMin,jMax
386                  DO i=iMin,iMax
387                   div(i,j) = div(i,j)*_hFacC(i,j,k,bi,bj)*drF(k)
388                  ENDDO
389                 ENDDO
390               ENDIF
391               DO j=jMin,jMax
392                DO i=iMin,iMax
393                  flx(i,j) = flx(i,j)
394    cc#ifdef ALLOW_AUTODIFF_OPENAD
395    cc     &                - rkSign*div(i,j)*rA(i,j,bi,bj)%v/deltaTLev(k)
396    cc#else
397         &                - rkSign*div(i,j)*rA(i,j,bi,bj)/deltaTLev(k)
398    cc#endif
399                  af(i,j) = flx(i,j) - df(i,j)
400                ENDDO
401               ENDDO
402               diagName = 'ADVr'//diagSufx
403               CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
404    #ifdef ALLOW_LAYERS
405               IF ( useLayers ) THEN
406                 CALL LAYERS_FILL(af,tracerIdentity,'AFR',
407         &                            k,1,2,bi,bj,myThid)
408               ENDIF
409    #endif /* ALLOW_LAYERS */
410            ENDIF            ENDIF
           CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)  
411           ENDDO           ENDDO
412          ENDIF          ENDIF
413        ENDIF        ENDIF
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22