/[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.51 by jahn, Fri Apr 18 19:39:48 2008 UTC
# Line 14  C !INTERFACE: ========================== Line 14  C !INTERFACE: ==========================
14       I           diffKh, diffK4, KappaR, TracerN, TracAB,       I           diffKh, diffK4, KappaR, TracerN, TracAB,
15       I           tracerIdentity, advectionScheme, vertAdvecScheme,       I           tracerIdentity, advectionScheme, vertAdvecScheme,
16       I           calcAdvection, implicitAdvection, applyAB_onTracer,       I           calcAdvection, implicitAdvection, applyAB_onTracer,
17         I           trUseGMRedi, trUseKPP,
18       U           fVerT, gTracer,       U           fVerT, gTracer,
19       I           myTime, myIter, myThid )       I           myTime, myIter, myThid )
20    
# Line 77  C vertAdvecScheme  :: advection scheme t Line 78  C vertAdvecScheme  :: advection scheme t
78  C calcAdvection    :: =False if Advec computed with multiDim scheme  C calcAdvection    :: =False if Advec computed with multiDim scheme
79  C implicitAdvection:: =True if vertical Advec computed implicitly  C implicitAdvection:: =True if vertical Advec computed implicitly
80  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
81    C trUseGMRedi      :: true if this tracer uses GM-Redi
82    C trUseKPP         :: true if this tracer uses KPP
83  C myTime           :: current time  C myTime           :: current time
84  C myIter           :: iteration number  C myIter           :: iteration number
85  C myThid           :: thread number  C myThid           :: thread number
# Line 100  C myThid           :: thread number Line 103  C myThid           :: thread number
103        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
104        LOGICAL calcAdvection        LOGICAL calcAdvection
105        LOGICAL implicitAdvection, applyAB_onTracer        LOGICAL implicitAdvection, applyAB_onTracer
106          LOGICAL trUseGMRedi, trUseKPP
107        _RL     myTime        _RL     myTime
108        INTEGER myIter, myThid        INTEGER myIter, myThid
109    
# Line 132  C locABT           :: local copy of (AB- Line 136  C locABT           :: local copy of (AB-
136        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138        _RL advFac, rAdvFac        _RL advFac, rAdvFac
139    #ifdef GAD_SMOLARKIEWICZ_HACK
140          _RL outFlux, trac, fac
141    #endif
142  CEOP  CEOP
143    
144  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 210  C-    Advective flux in X Line 217  C-    Advective flux in X
217            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
218          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
219       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
220            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
221       I            dTtracerLev(k), uTrans, uFld, locABT,       I            dTtracerLev(k), uTrans, uFld, locABT,
222       O            af, myThid )       O            af, myThid )
223          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
224            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
225       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
226       O            af, myThid )       O            af, myThid )
227          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
# Line 222  C-    Advective flux in X Line 229  C-    Advective flux in X
229          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
230            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
231          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
232            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
233       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
234       O            af, myThid )       O            af, myThid )
235          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
# Line 230  C-    Advective flux in X Line 237  C-    Advective flux in X
237  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
238  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
239  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
240            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
241       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
242       O           af, myThid )       O           af, myThid )
243           ELSE           ELSE
244            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
245       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
246       O           af, myThid )       O           af, myThid )
247           ENDIF           ENDIF
248            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
249              CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
250         I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
251         O            af, myThid )
252          ELSE          ELSE
253           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
254          ENDIF          ENDIF
# Line 272  C-    Add bi-harmonic diffusive flux in Line 283  C-    Add bi-harmonic diffusive flux in
283    
284  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
285  C-    GM/Redi flux in X  C-    GM/Redi flux in X
286        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
287  C *note* should update GMREDI_XTRANSPORT to set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
288          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
289            CALL GMREDI_XTRANSPORT(            CALL GMREDI_XTRANSPORT(
# Line 289  C *note* should update GMREDI_XTRANSPORT Line 300  C *note* should update GMREDI_XTRANSPORT
300          ENDIF          ENDIF
301        ENDIF        ENDIF
302  #endif  #endif
303    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
304        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
305         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
306          fZon(i,j) = fZon(i,j) + df(i,j)          fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
307         ENDDO         ENDDO
308        ENDDO        ENDDO
309    
# Line 299  C *note* should update GMREDI_XTRANSPORT Line 311  C *note* should update GMREDI_XTRANSPORT
311  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
312  C       excluding advective terms:  C       excluding advective terms:
313        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
314       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
315            diagName = 'DFxE'//diagSufx            diagName = 'DFxE'//diagSufx
316            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
317        ENDIF        ENDIF
# Line 318  C-    Advective flux in Y Line 330  C-    Advective flux in Y
330            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
331          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
332       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
333            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
334       I            dTtracerLev(k), vTrans, vFld, locABT,       I            dTtracerLev(k), vTrans, vFld, locABT,
335       O            af, myThid )       O            af, myThid )
336          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
337            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
338       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
339       O            af, myThid )       O            af, myThid )
340          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
# Line 330  C-    Advective flux in Y Line 342  C-    Advective flux in Y
342          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
343            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
344          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
345            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
346       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
347       O            af, myThid )       O            af, myThid )
348          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
# Line 338  C-    Advective flux in Y Line 350  C-    Advective flux in Y
350  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
351  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
352  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
353            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
354       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
355       O           af, myThid )       O           af, myThid )
356           ELSE           ELSE
357            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(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           ENDIF           ENDIF
361            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
362              CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
363         I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
364         O            af, myThid )
365          ELSE          ELSE
366            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
367          ENDIF          ENDIF
# Line 380  C-    Add bi-harmonic flux in Y Line 396  C-    Add bi-harmonic flux in Y
396    
397  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
398  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
399        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
400  C *note* should update GMREDI_YTRANSPORT to set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
401          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
402            CALL GMREDI_YTRANSPORT(            CALL GMREDI_YTRANSPORT(
# Line 397  C *note* should update GMREDI_YTRANSPORT Line 413  C *note* should update GMREDI_YTRANSPORT
413          ENDIF          ENDIF
414        ENDIF        ENDIF
415  #endif  #endif
416    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
417        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
418         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
419          fMer(i,j) = fMer(i,j) + df(i,j)          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
420         ENDDO         ENDDO
421        ENDDO        ENDDO
422    
# Line 407  C *note* should update GMREDI_YTRANSPORT Line 424  C *note* should update GMREDI_YTRANSPORT
424  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
425  C       excluding advective terms:  C       excluding advective terms:
426        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
427       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
428            diagName = 'DFyE'//diagSufx            diagName = 'DFyE'//diagSufx
429            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
430        ENDIF        ENDIF
# Line 456  cph with limiters in forward, but withou Line 473  cph with limiters in forward, but withou
473       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
474       O         af, myThid )       O         af, myThid )
475            ENDIF            ENDIF
476            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
477               CALL GAD_OS7MP_ADV_R( bi,bj,k,
478         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
479         O         af, myThid )
480          ELSE          ELSE
481            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
482          ENDIF          ENDIF
# Line 495  C           boundary condition. Line 516  C           boundary condition.
516    
517  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
518  C-    GM/Redi flux in R  C-    GM/Redi flux in R
519        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
520  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
521          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
522            CALL GMREDI_RTRANSPORT(            CALL GMREDI_RTRANSPORT(
# Line 523  C *note* should update GMREDI_RTRANSPORT Line 544  C *note* should update GMREDI_RTRANSPORT
544  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
545  C       Explicit terms only & excluding advective terms:  C       Explicit terms only & excluding advective terms:
546        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
547       &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN       &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
548            diagName = 'DFrE'//diagSufx            diagName = 'DFrE'//diagSufx
549            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
550        ENDIF        ENDIF
# Line 531  C       Explicit terms only & excluding Line 552  C       Explicit terms only & excluding
552    
553  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
554  C-    Set non local KPP transport term (ghat):  C-    Set non local KPP transport term (ghat):
555        IF ( useKPP .AND. k.GE.2 ) THEN        IF ( trUseKPP .AND. k.GE.2 ) THEN
556         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
557          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
558           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
# Line 539  C-    Set non local KPP transport term ( Line 560  C-    Set non local KPP transport term (
560         ENDDO         ENDDO
561         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
562          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
563       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
564       O     df )       O           df,
565         I           myTime, myIter, myThid )
566         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
567          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
568       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
569       O     df )       O           df,
570         I           myTime, myIter, myThid )
571  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
572         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
573          CALL KPP_TRANSPORT_PTR(          CALL KPP_TRANSPORT_PTR(
574       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
575       I     tracerIdentity-GAD_TR1+1,       I           tracerIdentity-GAD_TR1+1,
576       O     df )       O           df,
577         I           myTime, myIter, myThid )
578  #endif  #endif
579         ELSE         ELSE
580          PRINT*,'invalid tracer indentity: ', tracerIdentity          PRINT*,'invalid tracer indentity: ', tracerIdentity
# Line 558  C-    Set non local KPP transport term ( Line 582  C-    Set non local KPP transport term (
582         ENDIF         ENDIF
583         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
584          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
585           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
586         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
587            ENDDO
588           ENDDO
589          ENDIF
590    #endif
591    
592    #ifdef GAD_SMOLARKIEWICZ_HACK
593    coj   hack to make redi (and everything else in this s/r) positive
594    coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989)
595    coj   only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
596    coj
597    coj   apply to all tracers except temperature
598          IF (tracerIdentity.NE.GAD_TEMPERATURE) THEN
599           DO j=1-Oly,sNy+Oly-1
600            DO i=1-Olx,sNx+Olx-1
601    coj   add outgoing fluxes
602             outFlux=dTtracerLev(k)*
603         &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
604         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
605         &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
606         &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
607         &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
608         &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
609         &     )
610             IF ( applyAB_onTracer ) THEN
611               trac=TracerN(i,j,k,bi,bj)
612             ELSE
613               trac=TracAB(i,j,k,bi,bj)
614             ENDIF
615    coj   if they would reduce tracer by a fraction of more than
616    coj   SmolarkiewiczMaxFrac, ...
617             IF (outFlux.GT.0. _d 0 .AND.
618         &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
619    coj   ... scale them down
620    coj   if trac is already negative, set flux to zero
621               fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
622               IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
623               IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)      
624               IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
625               IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
626               IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
627         &      fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
628               IF (fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
629                fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
630                IF (k.LT.Nr) THEN
631    coj   down flux has already been applied in lower layer - undo it
632                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
633         &        -_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
634         &        *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
635         &        *recip_rhoFacC(k+1)
636         &        *( (-(fac-1. _d 0)*fVerT(i,j,kDown))*rkSign )
637                ENDIF
638               ENDIF
639             ENDIF
640          ENDDO          ENDDO
641         ENDDO         ENDDO
642        ENDIF        ENDIF
643  #endif  #endif
644    
645  C--   Divergence of fluxes  C--   Divergence of fluxes
646    C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
647        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
648         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
649          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
650       &   -_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)
651         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
652       &   *( (fZon(i+1,j)-fZon(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))
653       &     +(fMer(i,j+1)-fMer(i,j))       &     +(fMer(i,j+1)-fMer(i,j))
654       &     +(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.51
  ViewVC Help
Powered by ViewVC 1.1.22