/[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.20 by jmc, Mon Sep 22 22:13:11 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    #ifdef ALLOW_PTRACERS
47    #include "PTRACERS_OPTIONS.h"
48    #include "PTRACERS.h"
49    #endif
50    
51  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
52  #include "tamc.h"  #include "tamc.h"
# Line 61  C  diffKh               :: horizontal di Line 66  C  diffKh               :: horizontal di
66  C  diffK4               :: bi-harmonic diffusion coefficient  C  diffK4               :: bi-harmonic diffusion coefficient
67  C  KappaRT              :: 3-D array for vertical diffusion coefficient  C  KappaRT              :: 3-D array for vertical diffusion coefficient
68  C  Tracer               :: tracer field  C  Tracer               :: tracer field
69  C  tracerIdentity       :: identifier for the tracer (required only for KPP)  C  tracerIdentity       :: identifier for the tracer (required for KPP and GM)
70  C  advectionScheme      :: advection scheme to use  C  advectionScheme      :: advection scheme to use
71    C  calcAdvection        :: =False if Advec terms computed with multiDim scheme
72  C  myThid               :: thread number  C  myThid               :: thread number
73        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
74        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
# Line 77  C  myThid               :: thread number Line 83  C  myThid               :: thread number
83        _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)
84        INTEGER tracerIdentity        INTEGER tracerIdentity
85        INTEGER advectionScheme        INTEGER advectionScheme
86          LOGICAL calcAdvection
87        INTEGER myThid        INTEGER myThid
88    
89  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
# Line 128  C--   Make local copy of tracer array Line 135  C--   Make local copy of tracer array
135    
136  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
137  C     the tendency to zero.  C     the tendency to zero.
138        IF (.NOT. multiDimAdvection .OR.  C     <== now done earlier at the beginning of thermodynamics.
139       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  c     IF (calcAdvection) THEN
140       &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.  c      DO j=1-Oly,sNy+Oly
141       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN  c       DO i=1-Olx,sNx+Olx
142         DO j=1-Oly,sNy+Oly  c        gTracer(i,j,k,bi,bj)=0. _d 0
143          DO i=1-Olx,sNx+Olx  c       ENDDO
144           gTracer(i,j,k,bi,bj)=0. _d 0  c      ENDDO
145          ENDDO  c     ENDIF
        ENDDO  
       ENDIF  
146    
147  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
148        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 154  C--   Initialize net flux in X direction Line 159  C--   Initialize net flux in X direction
159        ENDDO        ENDDO
160    
161  C-    Advective flux in X  C-    Advective flux in X
162        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  
163        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
164         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
165        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 200  C-    GM/Redi flux in X Line 202  C-    GM/Redi flux in X
202  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*
203          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
204       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
205       I     xA,Tracer,       I     xA,Tracer,tracerIdentity,
206       U     df,       U     df,
207       I     myThid)       I     myThid)
208        ENDIF        ENDIF
# Line 229  C--   Initialize net flux in Y direction Line 231  C--   Initialize net flux in Y direction
231        ENDDO        ENDDO
232    
233  C-    Advective flux in Y  C-    Advective flux in Y
234        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  
235        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
236         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
237        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 275  C-    GM/Redi flux in Y Line 274  C-    GM/Redi flux in Y
274  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*
275         CALL GMREDI_YTRANSPORT(         CALL GMREDI_YTRANSPORT(
276       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
277       I     yA,Tracer,       I     yA,Tracer,tracerIdentity,
278       U     df,       U     df,
279       I     myThid)       I     myThid)
280        ENDIF        ENDIF
# Line 296  C-    Bi-harmonic flux in Y Line 295  C-    Bi-harmonic flux in Y
295         ENDDO         ENDDO
296        ENDIF        ENDIF
297    
298    #ifdef NONLIN_FRSURF
299    C--   Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):
300          IF ( calcAdvection .AND. K.EQ.Nr .AND.
301         &     useRealFreshWaterFlux .AND.
302         &     buoyancyRelation .EQ. 'OCEANICP' ) THEN  
303           DO j=1-Oly,sNy+Oly
304            DO i=1-Olx,sNx+Olx
305             fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)
306         &     *rA(i,j,bi,bj)*maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)
307            ENDDO
308           ENDDO
309          ENDIF
310    #endif /* NONLIN_FRSURF */
311    
312    C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
313  C-    Advective flux in R  C-    Advective flux in R
314        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  
315  C     Note: wVel needs to be masked  C     Note: wVel needs to be masked
316        IF (K.GE.2) THEN        IF (K.GE.2) THEN
317  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
# Line 358  C           boundary condition. Line 369  C           boundary condition.
369        ELSE        ELSE
370         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
371        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  
372    
373  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
374  C-    GM/Redi flux in R  C-    GM/Redi flux in R
# Line 370  C-    GM/Redi flux in R Line 376  C-    GM/Redi flux in R
376  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
377         CALL GMREDI_RTRANSPORT(         CALL GMREDI_RTRANSPORT(
378       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
379       I     Tracer,       I     Tracer,tracerIdentity,
380       U     df,       U     df,
381       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  
382        ENDIF        ENDIF
383  #endif  #endif
384    
# Line 406  C *note* should update KPP_TRANSPORT_T t Line 407  C *note* should update KPP_TRANSPORT_T t
407       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
408       I     KappaRT,       I     KappaRT,
409       U     df )       U     df )
410    #ifdef ALLOW_PTRACERS
411           ELSEIF (tracerIdentity .GE. GAD_TR1 .AND.
412         &         tracerIdentity .LE. (GAD_TR1+PTRACERS_numInUse-1)) THEN
413            CALL KPP_TRANSPORT_PTR(
414         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
415         I     tracerIdentity,KappaRT,
416         U     df )
417    #endif
418         ELSE         ELSE
419            PRINT*,'invalid tracer indentity: ', tracerIdentity
420          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
421         ENDIF         ENDIF
422         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 431  C--   Divergence of fluxes Line 441  C--   Divergence of fluxes
441         ENDDO         ENDDO
442        ENDDO        ENDDO
443    
444    #ifdef NONLIN_FRSURF
445    C-- account for 3.D divergence of the flow in rStar coordinate:
446          IF (calcAdvection .AND. select_rStar.GT.0) THEN
447           DO j=1-Oly,sNy+Oly-1
448            DO i=1-Olx,sNx+Olx-1
449             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
450         &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
451         &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
452            ENDDO
453           ENDDO
454          ENDIF
455          IF (calcAdvection .AND. select_rStar.LT.0) THEN
456           DO j=1-Oly,sNy+Oly-1
457            DO i=1-Olx,sNx+Olx-1
458             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
459         &     - rStarDhCDt(i,j,bi,bj)
460         &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
461            ENDDO
462           ENDDO
463          ENDIF
464    #endif /* NONLIN_FRSURF */
465          
466    
467        RETURN        RETURN
468        END        END
Legend:
Removed from v.1.12.6.1  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22