/[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.42 by jmc, Sun Nov 19 22:04:09 2006 UTC revision 1.55 by heimbach, Thu Oct 21 03:32:43 2010 UTC
# Line 12  C !INTERFACE: ========================== Line 12  C !INTERFACE: ==========================
12       I           xA, yA, maskUp, uFld, vFld, wFld,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           uTrans, vTrans, rTrans, rTransKp1,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           diffKh, diffK4, KappaR, TracerN, TracAB,       I           diffKh, diffK4, KappaR, TracerN, TracAB,
15       I           tracerIdentity, advectionScheme, vertAdvecScheme,       I           deltaTLev, tracerIdentity,
16         I           advectionScheme, vertAdvecScheme,
17       I           calcAdvection, implicitAdvection, applyAB_onTracer,       I           calcAdvection, implicitAdvection, applyAB_onTracer,
18         I           trUseGMRedi, trUseKPP,
19       U           fVerT, gTracer,       U           fVerT, gTracer,
20       I           myTime, myIter, myThid )       I           myTime, myIter, myThid )
21    
# Line 77  C vertAdvecScheme  :: advection scheme t Line 79  C vertAdvecScheme  :: advection scheme t
79  C calcAdvection    :: =False if Advec computed with multiDim scheme  C calcAdvection    :: =False if Advec computed with multiDim scheme
80  C implicitAdvection:: =True if vertical Advec computed implicitly  C implicitAdvection:: =True if vertical Advec computed implicitly
81  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
82    C trUseGMRedi      :: true if this tracer uses GM-Redi
83    C trUseKPP         :: true if this tracer uses KPP
84  C myTime           :: current time  C myTime           :: current time
85  C myIter           :: iteration number  C myIter           :: iteration number
86  C myThid           :: thread number  C myThid           :: thread number
# Line 96  C myThid           :: thread number Line 100  C myThid           :: thread number
100        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
102        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
103          _RL deltaTLev(Nr)
104        INTEGER tracerIdentity        INTEGER tracerIdentity
105        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
106        LOGICAL calcAdvection        LOGICAL calcAdvection
107        LOGICAL implicitAdvection, applyAB_onTracer        LOGICAL implicitAdvection, applyAB_onTracer
108          LOGICAL trUseGMRedi, trUseKPP
109        _RL     myTime        _RL     myTime
110        INTEGER myIter, myThid        INTEGER myIter, myThid
111    
# Line 132  C locABT           :: local copy of (AB- Line 138  C locABT           :: local copy of (AB-
138        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140        _RL advFac, rAdvFac        _RL advFac, rAdvFac
141    #ifdef GAD_SMOLARKIEWICZ_HACK
142          _RL outFlux, trac, fac, gTrFac
143    #endif
144  CEOP  CEOP
145    
146  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 210  C-    Advective flux in X Line 219  C-    Advective flux in X
219            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
220          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
221       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
222            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
223       I            dTtracerLev(k), uTrans, uFld, locABT,       I            deltaTLev(k), uTrans, uFld, locABT,
224       O            af, myThid )       O            af, myThid )
225          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
226            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
227       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
228       O            af, myThid )       O            af, myThid )
229          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
# Line 222  C-    Advective flux in X Line 231  C-    Advective flux in X
231          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
232            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
233          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
234            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
235       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
236       O            af, myThid )       O            af, myThid )
237          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
# Line 230  C-    Advective flux in X Line 239  C-    Advective flux in X
239  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
240  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
241  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
242            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
243       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
244       O           af, myThid )       O           af, myThid )
245           ELSE           ELSE
246            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
247       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
248       O           af, myThid )       O           af, myThid )
249           ENDIF           ENDIF
250    #ifndef ALLOW_AUTODIFF_TAMC
251            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
252              CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
253         I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
254         O            af, myThid )
255    #endif
256          ELSE          ELSE
257           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
258          ENDIF          ENDIF
# Line 272  C-    Add bi-harmonic diffusive flux in Line 287  C-    Add bi-harmonic diffusive flux in
287    
288  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
289  C-    GM/Redi flux in X  C-    GM/Redi flux in X
290        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
291  C *note* should update GMREDI_XTRANSPORT to set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
292          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
293            CALL GMREDI_XTRANSPORT(            CALL GMREDI_XTRANSPORT(
# Line 289  C *note* should update GMREDI_XTRANSPORT Line 304  C *note* should update GMREDI_XTRANSPORT
304          ENDIF          ENDIF
305        ENDIF        ENDIF
306  #endif  #endif
307    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
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          fZon(i,j) = fZon(i,j) + df(i,j)          fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
311         ENDDO         ENDDO
312        ENDDO        ENDDO
313    
# Line 299  C *note* should update GMREDI_XTRANSPORT Line 315  C *note* should update GMREDI_XTRANSPORT
315  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
316  C       excluding advective terms:  C       excluding advective terms:
317        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
318       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
319            diagName = 'DFxE'//diagSufx            diagName = 'DFxE'//diagSufx
320            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
321        ENDIF        ENDIF
# Line 318  C-    Advective flux in Y Line 334  C-    Advective flux in Y
334            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
335          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
336       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
337            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
338       I            dTtracerLev(k), vTrans, vFld, locABT,       I            deltaTLev(k), vTrans, vFld, locABT,
339       O            af, myThid )       O            af, myThid )
340          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
341            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
342       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
343       O            af, myThid )       O            af, myThid )
344          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
# Line 330  C-    Advective flux in Y Line 346  C-    Advective flux in Y
346          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
347            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
348          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
349            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
350       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
351       O            af, myThid )       O            af, myThid )
352          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
# Line 338  C-    Advective flux in Y Line 354  C-    Advective flux in Y
354  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
355  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
356  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
357            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
358       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
359       O           af, myThid )       O           af, myThid )
360           ELSE           ELSE
361            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
362       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
363       O           af, myThid )       O           af, myThid )
364           ENDIF           ENDIF
365    #ifndef ALLOW_AUTODIFF_TAMC
366            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
367              CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
368         I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
369         O            af, myThid )
370    #endif
371          ELSE          ELSE
372            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
373          ENDIF          ENDIF
# Line 380  C-    Add bi-harmonic flux in Y Line 402  C-    Add bi-harmonic flux in Y
402    
403  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
404  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
405        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
406  C *note* should update GMREDI_YTRANSPORT to set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
407          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
408            CALL GMREDI_YTRANSPORT(            CALL GMREDI_YTRANSPORT(
# Line 397  C *note* should update GMREDI_YTRANSPORT Line 419  C *note* should update GMREDI_YTRANSPORT
419          ENDIF          ENDIF
420        ENDIF        ENDIF
421  #endif  #endif
422    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
423        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
424         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
425          fMer(i,j) = fMer(i,j) + df(i,j)          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
426         ENDDO         ENDDO
427        ENDDO        ENDDO
428    
# Line 407  C *note* should update GMREDI_YTRANSPORT Line 430  C *note* should update GMREDI_YTRANSPORT
430  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
431  C       excluding advective terms:  C       excluding advective terms:
432        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
433       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
434            diagName = 'DFyE'//diagSufx            diagName = 'DFyE'//diagSufx
435            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
436        ENDIF        ENDIF
# Line 429  C-    Compute vertical advective flux in Line 452  C-    Compute vertical advective flux in
452          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
453       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
454            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
455       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
456       O         af, myThid )       O         af, myThid )
457          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
458            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
459       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
460       O         af, myThid )       O         af, myThid )
461          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
462            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 441  C-    Compute vertical advective flux in Line 464  C-    Compute vertical advective flux in
464            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
465          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
466            CALL GAD_DST3_ADV_R( bi,bj,k,            CALL GAD_DST3_ADV_R( bi,bj,k,
467       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
468       O         af, myThid )       O         af, myThid )
469          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
470  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
# Line 449  cph IF inAdExact=.FALSE., we want to use Line 472  cph IF inAdExact=.FALSE., we want to use
472  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
473            IF ( inAdMode ) THEN            IF ( inAdMode ) THEN
474             CALL GAD_DST3_ADV_R( bi,bj,k,             CALL GAD_DST3_ADV_R( bi,bj,k,
475       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
476       O         af, myThid )       O         af, myThid )
477            ELSE            ELSE
478             CALL GAD_DST3FL_ADV_R( bi,bj,k,             CALL GAD_DST3FL_ADV_R( bi,bj,k,
479       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
480       O         af, myThid )       O         af, myThid )
481            ENDIF            ENDIF
482    #ifndef ALLOW_AUTODIFF_TAMC
483            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
484               CALL GAD_OS7MP_ADV_R( bi,bj,k,
485         I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
486         O         af, myThid )
487    #endif
488          ELSE          ELSE
489            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
490          ENDIF          ENDIF
# Line 495  C           boundary condition. Line 524  C           boundary condition.
524    
525  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
526  C-    GM/Redi flux in R  C-    GM/Redi flux in R
527        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
528  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
529          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
530            CALL GMREDI_RTRANSPORT(            CALL GMREDI_RTRANSPORT(
# Line 523  C *note* should update GMREDI_RTRANSPORT Line 552  C *note* should update GMREDI_RTRANSPORT
552  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
553  C       Explicit terms only & excluding advective terms:  C       Explicit terms only & excluding advective terms:
554        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
555       &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN       &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
556            diagName = 'DFrE'//diagSufx            diagName = 'DFrE'//diagSufx
557            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
558        ENDIF        ENDIF
# Line 531  C       Explicit terms only & excluding Line 560  C       Explicit terms only & excluding
560    
561  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
562  C-    Set non local KPP transport term (ghat):  C-    Set non local KPP transport term (ghat):
563        IF ( useKPP .AND. k.GE.2 ) THEN        IF ( trUseKPP .AND. k.GE.2 ) THEN
564         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
565          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
566           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
# Line 539  C-    Set non local KPP transport term ( Line 568  C-    Set non local KPP transport term (
568         ENDDO         ENDDO
569         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
570          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
571       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
572       O     df )       O           df,
573         I           myTime, myIter, myThid )
574         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
575          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
576       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
577       O     df )       O           df,
578         I           myTime, myIter, myThid )
579  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
580         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
581          CALL KPP_TRANSPORT_PTR(          CALL KPP_TRANSPORT_PTR(
582       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
583       I     tracerIdentity-GAD_TR1+1,       I           tracerIdentity-GAD_TR1+1,
584       O     df )       O           df,
585         I           myTime, myIter, myThid )
586  #endif  #endif
587         ELSE         ELSE
588          PRINT*,'invalid tracer indentity: ', tracerIdentity          WRITE(errorMessageUnit,*)
589          STOP 'GAD_CALC_RHS: Ooops'       &    'tracer identity =', tracerIdentity, ' is not valid => STOP'
590            STOP 'ABNORMAL END: S/R GAD_CALC_RHS: invalid tracer identity'
591         ENDIF         ENDIF
592         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
593          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
594           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
595         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
596            ENDDO
597           ENDDO
598    #ifdef ALLOW_DIAGNOSTICS
599    C-    Diagnostics of Non-Local Tracer (vertical) flux
600           IF ( useDiagnostics ) THEN
601             diagName = 'KPPg'//diagSufx
602             CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
603    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
604    C        does it only if k=1 (never the case here)
605             IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
606           ENDIF
607    #endif
608          ENDIF
609    #endif /* ALLOW_KPP */
610    
611    #ifdef GAD_SMOLARKIEWICZ_HACK
612    coj   Hack to make redi (and everything else in this s/r) positive
613    coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989).
614    coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
615    coj
616    coj   Apply to all tracers except temperature
617          IF (tracerIdentity.NE.GAD_TEMPERATURE .AND.
618         &    tracerIdentity.NE.GAD_SALINITY) THEN
619           DO j=1-Oly,sNy+Oly-1
620            DO i=1-Olx,sNx+Olx-1
621    coj   Add outgoing fluxes
622             outFlux=deltaTLev(k)*
623         &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
624         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
625         &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
626         &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
627         &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
628         &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
629         &     )
630             IF ( applyAB_onTracer ) THEN
631               trac=TracerN(i,j,k,bi,bj)
632             ELSE
633               trac=TracAB(i,j,k,bi,bj)
634             ENDIF
635    coj   If they would reduce tracer by a fraction of more than
636    coj   SmolarkiewiczMaxFrac, scale them down
637             IF (outFlux.GT.0. _d 0 .AND.
638         &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
639    coj   If tracer is already negative, scale flux to zero
640               fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
641    
642               IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
643               IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)
644               IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
645               IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
646               IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
647         &       fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
648    
649               IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
650    coj   Down flux is special: it has already been applied in lower layer,
651    coj   so we have to readjust this.
652    coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!
653    coj   thus it has an extra factor deltaTLev(k+1)
654                 gTrFac=deltaTLev(k+1)
655    coj   Other factors that have been applied to gTracer since the last call:
656    #ifdef NONLIN_FRSURF
657                 IF (nonlinFreeSurf.GT.0) THEN
658                  IF (select_rStar.GT.0) THEN
659    #ifndef DISABLE_RSTAR_CODE
660                    gTrFac = gTrFac/rStarExpC(i,j,bi,bj)
661    #endif /* DISABLE_RSTAR_CODE */
662                  ENDIF
663                 ENDIF
664    #endif /* NONLIN_FRSURF */
665    coj   Now: undo down flux, ...
666                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
667         &        +gTrFac
668         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
669         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
670         &         *recip_rhoFacC(k+1)
671         &         *( -fVerT(i,j,kDown)*rkSign )
672    coj   ... scale ...
673                 fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
674    coj   ... and reapply
675                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
676         &        +gTrFac
677         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
678         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
679         &         *recip_rhoFacC(k+1)
680         &         *( fVerT(i,j,kDown)*rkSign )
681               ENDIF
682    
683             ENDIF
684          ENDDO          ENDDO
685         ENDDO         ENDDO
686        ENDIF        ENDIF
687  #endif  #endif
688    
689  C--   Divergence of fluxes  C--   Divergence of fluxes
690    C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
691        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
692         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
693          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
694       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
695         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
696       &   *( (fZon(i+1,j)-fZon(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))
697       &     +(fMer(i,j+1)-fMer(i,j))       &     +(fMer(i,j+1)-fMer(i,j))
698       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
Legend:
Removed from v.1.42  
changed lines
  Added in v.1.55
  ViewVC Help
Powered by ViewVC 1.1.22