/[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.12.6.1 by heimbach, Sun Mar 24 17:22:22 2002 UTC revision 1.22 by dimitri, Thu Sep 25 03:01:59 2003 UTC
# Line 11  C !INTERFACE: ========================== Line 11  C !INTERFACE: ==========================
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,maskUp,       I           xA,yA,uTrans,vTrans,rTrans,maskUp,
13       I           diffKh, diffK4, KappaRT, Tracer,       I           diffKh, diffK4, KappaRT, Tracer,
14       I           tracerIdentity, advectionScheme,       I           tracerIdentity, advectionScheme, calcAdvection,
15       U           fVerT, gTracer,       U           fVerT, gTracer,
16       I           myThid )       I           myThid )
17    
# Line 41  C !USES: =============================== Line 41  C !USES: ===============================
41  #include "PARAMS.h"  #include "PARAMS.h"
42  #include "GRID.h"  #include "GRID.h"
43  #include "DYNVARS.h"  #include "DYNVARS.h"
44    #include "SURFACE.h"
45  #include "GAD.h"  #include "GAD.h"
46    
47  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 61  C  diffKh               :: horizontal di Line 62  C  diffKh               :: horizontal di
62  C  diffK4               :: bi-harmonic diffusion coefficient  C  diffK4               :: bi-harmonic diffusion coefficient
63  C  KappaRT              :: 3-D array for vertical diffusion coefficient  C  KappaRT              :: 3-D array for vertical diffusion coefficient
64  C  Tracer               :: tracer field  C  Tracer               :: tracer field
65  C  tracerIdentity       :: identifier for the tracer (required only for KPP)  C  tracerIdentity       :: identifier for the tracer (required for KPP and GM)
66  C  advectionScheme      :: advection scheme to use  C  advectionScheme      :: advection scheme to use
67    C  calcAdvection        :: =False if Advec terms computed with multiDim scheme
68  C  myThid               :: thread number  C  myThid               :: thread number
69        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
70        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
# Line 77  C  myThid               :: thread number Line 79  C  myThid               :: thread number
79        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
80        INTEGER tracerIdentity        INTEGER tracerIdentity
81        INTEGER advectionScheme        INTEGER advectionScheme
82          LOGICAL calcAdvection
83        INTEGER myThid        INTEGER myThid
84    
85  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
# Line 128  C--   Make local copy of tracer array Line 131  C--   Make local copy of tracer array
131    
132  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
133  C     the tendency to zero.  C     the tendency to zero.
134        IF (.NOT. multiDimAdvection .OR.  C     <== now done earlier at the beginning of thermodynamics.
135       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  c     IF (calcAdvection) THEN
136       &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.  c      DO j=1-Oly,sNy+Oly
137       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN  c       DO i=1-Olx,sNx+Olx
138         DO j=1-Oly,sNy+Oly  c        gTracer(i,j,k,bi,bj)=0. _d 0
139          DO i=1-Olx,sNx+Olx  c       ENDDO
140           gTracer(i,j,k,bi,bj)=0. _d 0  c      ENDDO
141          ENDDO  c     ENDIF
        ENDDO  
       ENDIF  
142    
143  C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero  C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
144        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 154  C--   Initialize net flux in X direction Line 155  C--   Initialize net flux in X direction
155        ENDDO        ENDDO
156    
157  C-    Advective flux in X  C-    Advective flux in X
158        IF (.NOT. multiDimAdvection .OR.        IF (calcAdvection) THEN
      &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  
      &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.  
      &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN  
159        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
160         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
161        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 200  C-    GM/Redi flux in X Line 198  C-    GM/Redi flux in X
198  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*
199          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
200       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
201       I     xA,Tracer,       I     xA,Tracer,tracerIdentity,
202       U     df,       U     df,
203       I     myThid)       I     myThid)
204        ENDIF        ENDIF
# Line 229  C--   Initialize net flux in Y direction Line 227  C--   Initialize net flux in Y direction
227        ENDDO        ENDDO
228    
229  C-    Advective flux in Y  C-    Advective flux in Y
230        IF (.NOT. multiDimAdvection .OR.        IF (calcAdvection) THEN
      &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  
      &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.  
      &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN  
231        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
232         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
233        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 275  C-    GM/Redi flux in Y Line 270  C-    GM/Redi flux in Y
270  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*
271         CALL GMREDI_YTRANSPORT(         CALL GMREDI_YTRANSPORT(
272       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
273       I     yA,Tracer,       I     yA,Tracer,tracerIdentity,
274       U     df,       U     df,
275       I     myThid)       I     myThid)
276        ENDIF        ENDIF
# Line 296  C-    Bi-harmonic flux in Y Line 291  C-    Bi-harmonic flux in Y
291         ENDDO         ENDDO
292        ENDIF        ENDIF
293    
294    #ifdef NONLIN_FRSURF
295    C--   Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):
296          IF ( calcAdvection .AND. K.EQ.Nr .AND.
297         &     useRealFreshWaterFlux .AND.
298         &     buoyancyRelation .EQ. 'OCEANICP' ) THEN  
299           DO j=1-Oly,sNy+Oly
300            DO i=1-Olx,sNx+Olx
301             fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)
302         &     *rA(i,j,bi,bj)*maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)
303            ENDDO
304           ENDDO
305          ENDIF
306    #endif /* NONLIN_FRSURF */
307    
308    C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
309  C-    Advective flux in R  C-    Advective flux in R
310        IF (.NOT. multiDimAdvection .OR.        IF (calcAdvection) THEN
      &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  
      &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.  
      &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN  
311  C     Note: wVel needs to be masked  C     Note: wVel needs to be masked
312        IF (K.GE.2) THEN        IF (K.GE.2) THEN
313  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
# Line 358  C           boundary condition. Line 365  C           boundary condition.
365        ELSE        ELSE
366         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
367        ENDIF        ENDIF
 c     DO j=1-Oly,sNy+Oly  
 c      DO i=1-Olx,sNx+Olx  
 c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)  
 c      ENDDO  
 c     ENDDO  
368    
369  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
370  C-    GM/Redi flux in R  C-    GM/Redi flux in R
# Line 370  C-    GM/Redi flux in R Line 372  C-    GM/Redi flux in R
372  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
373         CALL GMREDI_RTRANSPORT(         CALL GMREDI_RTRANSPORT(
374       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
375       I     Tracer,       I     Tracer,tracerIdentity,
376       U     df,       U     df,
377       I     myThid)       I     myThid)
 c      DO j=1-Oly,sNy+Oly  
 c       DO i=1-Olx,sNx+Olx  
 c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)  
 c       ENDDO  
 c      ENDDO  
378        ENDIF        ENDIF
379  #endif  #endif
380    
# Line 406  C *note* should update KPP_TRANSPORT_T t Line 403  C *note* should update KPP_TRANSPORT_T t
403       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
404       I     KappaRT,       I     KappaRT,
405       U     df )       U     df )
406    #ifdef ALLOW_PTRACERS
407           ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
408            CALL KPP_TRANSPORT_PTR(
409         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
410         I     tracerIdentity-GAD_TR1+1,KappaRT,
411         U     df )
412    #endif
413         ELSE         ELSE
414            PRINT*,'invalid tracer indentity: ', tracerIdentity
415          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
416         ENDIF         ENDIF
417         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 431  C--   Divergence of fluxes Line 436  C--   Divergence of fluxes
436         ENDDO         ENDDO
437        ENDDO        ENDDO
438    
439    #ifdef NONLIN_FRSURF
440    C-- account for 3.D divergence of the flow in rStar coordinate:
441          IF (calcAdvection .AND. select_rStar.GT.0) THEN
442           DO j=1-Oly,sNy+Oly-1
443            DO i=1-Olx,sNx+Olx-1
444             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
445         &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
446         &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
447            ENDDO
448           ENDDO
449          ENDIF
450          IF (calcAdvection .AND. select_rStar.LT.0) THEN
451           DO j=1-Oly,sNy+Oly-1
452            DO i=1-Olx,sNx+Olx-1
453             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
454         &     - rStarDhCDt(i,j,bi,bj)
455         &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
456            ENDDO
457           ENDDO
458          ENDIF
459    #endif /* NONLIN_FRSURF */
460          
461    
462        RETURN        RETURN
463        END        END
Legend:
Removed from v.1.12.6.1  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22