/[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.52 by jahn, Wed Apr 23 18:32:20 2008 UTC revision 1.53 by jahn, Fri Jun 26 23:10:09 2009 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,       I           trUseGMRedi, trUseKPP,
19       U           fVerT, gTracer,       U           fVerT, gTracer,
# Line 99  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
# Line 218  C-    Advective flux in X Line 220  C-    Advective flux in X
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, .TRUE.,            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, .TRUE., 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 229  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, .TRUE., 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 237  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, .TRUE., 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, .TRUE., 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          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
251            CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
252       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
253       O            af, myThid )       O            af, myThid )
254          ELSE          ELSE
# Line 331  C-    Advective flux in Y Line 333  C-    Advective flux in Y
333          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
334       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
335            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
336       I            dTtracerLev(k), vTrans, vFld, locABT,       I            deltaTLev(k), vTrans, vFld, locABT,
337       O            af, myThid )       O            af, myThid )
338          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
339            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
340       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
341       O            af, myThid )       O            af, myThid )
342          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
# Line 342  C-    Advective flux in Y Line 344  C-    Advective flux in Y
344          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
345            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
346          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
347            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
348       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
349       O            af, myThid )       O            af, myThid )
350          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
# Line 350  C-    Advective flux in Y Line 352  C-    Advective flux in Y
352  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
353  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
354  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
355            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
356       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
357       O           af, myThid )       O           af, myThid )
358           ELSE           ELSE
359            CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
360       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
361       O           af, myThid )       O           af, myThid )
362           ENDIF           ENDIF
363          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
364            CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
365       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
366       O            af, myThid )       O            af, myThid )
367          ELSE          ELSE
# Line 446  C-    Compute vertical advective flux in Line 448  C-    Compute vertical advective flux in
448          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
449       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
450            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
451       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),
452       O         af, myThid )       O         af, myThid )
453          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
454            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
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_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
458            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 458  C-    Compute vertical advective flux in Line 460  C-    Compute vertical advective flux in
460            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
461          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
462            CALL GAD_DST3_ADV_R( bi,bj,k,            CALL GAD_DST3_ADV_R( bi,bj,k,
463       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),
464       O         af, myThid )       O         af, myThid )
465          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
466  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
# Line 466  cph IF inAdExact=.FALSE., we want to use Line 468  cph IF inAdExact=.FALSE., we want to use
468  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
469            IF ( inAdMode ) THEN            IF ( inAdMode ) THEN
470             CALL GAD_DST3_ADV_R( bi,bj,k,             CALL GAD_DST3_ADV_R( bi,bj,k,
471       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),
472       O         af, myThid )       O         af, myThid )
473            ELSE            ELSE
474             CALL GAD_DST3FL_ADV_R( bi,bj,k,             CALL GAD_DST3FL_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            ENDIF            ENDIF
478          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
479             CALL GAD_OS7MP_ADV_R( bi,bj,k,             CALL GAD_OS7MP_ADV_R( bi,bj,k,
480       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),
481       O         af, myThid )       O         af, myThid )
482          ELSE          ELSE
483            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
# Line 600  coj   Apply to all tracers except temper Line 602  coj   Apply to all tracers except temper
602         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
603          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
604  coj   Add outgoing fluxes  coj   Add outgoing fluxes
605           outFlux=dTtracerLev(k)*           outFlux=deltaTLev(k)*
606       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
607       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
608       &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))       &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
# Line 631  coj   If tracer is already negative, sca Line 633  coj   If tracer is already negative, sca
633  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,
634  coj   so we have to readjust this.  coj   so we have to readjust this.
635  coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!  coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!
636  coj   thus it has an extra factor dTtracerLev(k+1)  coj   thus it has an extra factor deltaTLev(k+1)
637               gTrFac=dTtracerLev(k+1)               gTrFac=deltaTLev(k+1)
638  coj   Other factors that have been applied to gTracer since the last call:  coj   Other factors that have been applied to gTracer since the last call:
639  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
640               IF (nonlinFreeSurf.GT.0) THEN               IF (nonlinFreeSurf.GT.0) THEN
Legend:
Removed from v.1.52  
changed lines
  Added in v.1.53
  ViewVC Help
Powered by ViewVC 1.1.22