/[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.51 by jahn, Fri Apr 18 19:39:48 2008 UTC revision 1.52 by jahn, Wed Apr 23 18:32:20 2008 UTC
# Line 137  C locABT           :: local copy of (AB- Line 137  C locABT           :: local copy of (AB-
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  #ifdef GAD_SMOLARKIEWICZ_HACK
140        _RL outFlux, trac, fac        _RL outFlux, trac, fac, gTrFac
141  #endif  #endif
142  CEOP  CEOP
143    
# Line 590  C-    Set non local KPP transport term ( Line 590  C-    Set non local KPP transport term (
590  #endif  #endif
591    
592  #ifdef GAD_SMOLARKIEWICZ_HACK  #ifdef GAD_SMOLARKIEWICZ_HACK
593  coj   hack to make redi (and everything else in this s/r) positive  coj   Hack to make redi (and everything else in this s/r) positive
594  coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989)  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  coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
596  coj  coj
597  coj   apply to all tracers except temperature  coj   Apply to all tracers except temperature
598        IF (tracerIdentity.NE.GAD_TEMPERATURE) THEN        IF (tracerIdentity.NE.GAD_TEMPERATURE .AND.
599         &    tracerIdentity.NE.GAD_SALINITY) THEN
600         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
601          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
602  coj   add outgoing fluxes  coj   Add outgoing fluxes
603           outFlux=dTtracerLev(k)*           outFlux=dTtracerLev(k)*
604       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
605       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
# Line 612  coj   add outgoing fluxes Line 613  coj   add outgoing fluxes
613           ELSE           ELSE
614             trac=TracAB(i,j,k,bi,bj)             trac=TracAB(i,j,k,bi,bj)
615           ENDIF           ENDIF
616  coj   if they would reduce tracer by a fraction of more than  coj   If they would reduce tracer by a fraction of more than
617  coj   SmolarkiewiczMaxFrac, ...  coj   SmolarkiewiczMaxFrac, scale them down
618           IF (outFlux.GT.0. _d 0 .AND.           IF (outFlux.GT.0. _d 0 .AND.
619       &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN       &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
620  coj   ... scale them down  coj   If tracer is already negative, scale flux to zero
 coj   if trac is already negative, set flux to zero  
621             fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)             fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
622    
623             IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)             IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
624             IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)                   IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)      
625             IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)             IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
626             IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)             IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
627             IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)             IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
628       &      fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)       &       fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
629             IF (fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN  
630              fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)             IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
631              IF (k.LT.Nr) THEN  coj   Down flux is special: it has already been applied in lower layer,
632  coj   down flux has already been applied in lower layer - undo it  coj   so we have to readjust this.
633    coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!
634    coj   thus it has an extra factor dTtracerLev(k+1)
635                 gTrFac=dTtracerLev(k+1)
636    coj   Other factors that have been applied to gTracer since the last call:
637    #ifdef NONLIN_FRSURF
638                 IF (nonlinFreeSurf.GT.0) THEN
639                  IF (select_rStar.GT.0) THEN
640    #ifndef DISABLE_RSTAR_CODE
641                    gTrFac = gTrFac/rStarExpC(i,j,bi,bj)
642    #endif /* DISABLE_RSTAR_CODE */
643                  ENDIF
644                 ENDIF
645    #endif /* NONLIN_FRSURF */
646    coj   Now: undo down flux, ...
647               gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)               gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
648       &        -_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)       &        +gTrFac
649       &        *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)       &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
650       &        *recip_rhoFacC(k+1)       &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
651       &        *( (-(fac-1. _d 0)*fVerT(i,j,kDown))*rkSign )       &         *recip_rhoFacC(k+1)
652              ENDIF       &         *( -fVerT(i,j,kDown)*rkSign )
653    coj   ... scale ...
654                 fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
655    coj   ... and reapply
656                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
657         &        +gTrFac
658         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
659         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
660         &         *recip_rhoFacC(k+1)
661         &         *( fVerT(i,j,kDown)*rkSign )
662             ENDIF             ENDIF
663    
664           ENDIF           ENDIF
665          ENDDO          ENDDO
666         ENDDO         ENDDO
Legend:
Removed from v.1.51  
changed lines
  Added in v.1.52
  ViewVC Help
Powered by ViewVC 1.1.22