/[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.13 by heimbach, Sun Mar 24 02:12:50 2002 UTC revision 1.24 by edhill, Mon Mar 29 03:33:51 2004 UTC
# Line 9  C !ROUTINE: GAD_CALC_RHS Line 9  C !ROUTINE: GAD_CALC_RHS
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_CALC_RHS(        SUBROUTINE GAD_CALC_RHS(
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,rTransKp1,maskUp,
13         I           uVel, vVel, wVel,
14       I           diffKh, diffK4, KappaRT, Tracer,       I           diffKh, diffK4, KappaRT, Tracer,
15       I           tracerIdentity, advectionScheme,       I           tracerIdentity, advectionScheme,
16         I           calcAdvection, implicitAdvection,
17       U           fVerT, gTracer,       U           fVerT, gTracer,
18       I           myThid )       I           myThid )
19    
# Line 40  C !USES: =============================== Line 42  C !USES: ===============================
42  #include "EEPARAMS.h"  #include "EEPARAMS.h"
43  #include "PARAMS.h"  #include "PARAMS.h"
44  #include "GRID.h"  #include "GRID.h"
45  #include "DYNVARS.h"  #include "SURFACE.h"
46  #include "GAD.h"  #include "GAD.h"
47    
48  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 49  C !USES: =============================== Line 51  C !USES: ===============================
51  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
52    
53  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
54  C  bi,bj                :: tile indices  C bi,bj            :: tile indices
55  C  iMin,iMax,jMin,jMax  :: loop range for called routines  C iMin,iMax        :: loop range for called routines
56  C  kup                  :: index into 2 1/2D array, toggles between 1 and 2  C jMin,jMax        :: loop range for called routines
57  C  kdown                :: index into 2 1/2D array, toggles between 2 and 1  C kup              :: index into 2 1/2D array, toggles between 1|2
58  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C kdown            :: index into 2 1/2D array, toggles between 2|1
59  C  xA,yA                :: areas of X and Y face of tracer cells  C kp1              :: =k+1 for k<Nr, =Nr for k=Nr
60  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C xA,yA            :: areas of X and Y face of tracer cells
61  C  maskUp               :: 2-D array for mask at W points  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
62  C  diffKh               :: horizontal diffusion coefficient  C rTrans           :: 2-D arrays of volume transports at W points
63  C  diffK4               :: bi-harmonic diffusion coefficient  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
64  C  KappaRT              :: 3-D array for vertical diffusion coefficient  C maskUp           :: 2-D array for mask at W points
65  C  Tracer               :: tracer field  C uVel,vVel,wVel   :: 3 components of the velcity field (3-D array)
66  C  tracerIdentity       :: identifier for the tracer (required only for KPP)  C diffKh           :: horizontal diffusion coefficient
67  C  advectionScheme      :: advection scheme to use  C diffK4           :: bi-harmonic diffusion coefficient
68  C  myThid               :: thread number  C KappaRT          :: 3-D array for vertical diffusion coefficient
69    C Tracer           :: tracer field
70    C tracerIdentity   :: tracer identifier (required for KPP,GM)
71    C advectionScheme  :: advection scheme to use
72    C calcAdvection    :: =False if Advec computed with multiDim scheme
73    C implicitAdvection:: =True if vertical Advec computed implicitly
74    C myThid           :: thread number
75        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
76        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
77        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 71  C  myThid               :: thread number Line 79  C  myThid               :: thread number
79        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84          _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
85          _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
86          _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
87        _RL diffKh, diffK4        _RL diffKh, diffK4
88        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _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)
90        INTEGER tracerIdentity        INTEGER tracerIdentity
91        INTEGER advectionScheme        INTEGER advectionScheme
92          LOGICAL calcAdvection
93          LOGICAL implicitAdvection
94        INTEGER myThid        INTEGER myThid
95    
96  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
97  C  gTracer              :: tendancy array  C gTracer          :: tendancy array
98  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
99        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
100        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
101    
102  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
103  C  i,j                  :: loop indices  C i,j              :: loop indices
104  C  df4                  :: used for storing del^2 T for bi-harmonic term  C df4              :: used for storing del^2 T for bi-harmonic term
105  C  fZon                 :: zonal flux  C fZon             :: zonal flux
106  C  fmer                 :: meridional flux  C fmer             :: meridional flux
107  C  af                   :: advective flux  C af               :: advective flux
108  C  df                   :: diffusive flux  C df               :: diffusive flux
109  C  localT               :: local copy of tracer field  C localT           :: local copy of tracer field
110        INTEGER i,j        INTEGER i,j
111        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 100  C  localT               :: local copy of Line 114  C  localT               :: local copy of
114        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117          _RL advFac, rAdvFac
118  CEOP  CEOP
119    
120  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 108  C--   the kDown is still required Line 123  C--   the kDown is still required
123        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
124  #endif  #endif
125    
126          advFac  = 0. _d 0
127          IF (calcAdvection) advFac = 1. _d 0
128          rAdvFac = rkFac*advFac
129          IF (implicitAdvection) rAdvFac = 0. _d 0
130    
131        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
132         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
133          fZon(i,j)      = 0. _d 0          fZon(i,j)      = 0. _d 0
# Line 115  C--   the kDown is still required Line 135  C--   the kDown is still required
135          fVerT(i,j,kUp) = 0. _d 0          fVerT(i,j,kUp) = 0. _d 0
136          df(i,j)        = 0. _d 0          df(i,j)        = 0. _d 0
137          df4(i,j)       = 0. _d 0          df4(i,j)       = 0. _d 0
         localT(i,j)    = 0. _d 0  
138         ENDDO         ENDDO
139        ENDDO        ENDDO
140    
# Line 128  C--   Make local copy of tracer array Line 147  C--   Make local copy of tracer array
147    
148  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
149  C     the tendency to zero.  C     the tendency to zero.
150        IF (.NOT. multiDimAdvection .OR.  C     <== now done earlier at the beginning of thermodynamics.
151       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  c     IF (calcAdvection) THEN
152       &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.  c      DO j=1-Oly,sNy+Oly
153       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN  c       DO i=1-Olx,sNx+Olx
154         DO j=1-Oly,sNy+Oly  c        gTracer(i,j,k,bi,bj)=0. _d 0
155          DO i=1-Olx,sNx+Olx  c       ENDDO
156           gTracer(i,j,k,bi,bj)=0. _d 0  c      ENDDO
157          ENDDO  c     ENDIF
        ENDDO  
       ENDIF  
158    
159  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
160        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 154  C--   Initialize net flux in X direction Line 171  C--   Initialize net flux in X direction
171        ENDDO        ENDDO
172    
173  C-    Advective flux in X  C-    Advective flux in X
174        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  
175        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
176         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
177        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 200  C-    GM/Redi flux in X Line 214  C-    GM/Redi flux in X
214  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*
215          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
216       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
217       I     xA,Tracer,       I     xA,Tracer,tracerIdentity,
218       U     df,       U     df,
219       I     myThid)       I     myThid)
220        ENDIF        ENDIF
# Line 229  C--   Initialize net flux in Y direction Line 243  C--   Initialize net flux in Y direction
243        ENDDO        ENDDO
244    
245  C-    Advective flux in Y  C-    Advective flux in Y
246        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  
247        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
248         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
249        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 275  C-    GM/Redi flux in Y Line 286  C-    GM/Redi flux in Y
286  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*
287         CALL GMREDI_YTRANSPORT(         CALL GMREDI_YTRANSPORT(
288       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
289       I     yA,Tracer,       I     yA,Tracer,tracerIdentity,
290       U     df,       U     df,
291       I     myThid)       I     myThid)
292        ENDIF        ENDIF
# Line 296  C-    Bi-harmonic flux in Y Line 307  C-    Bi-harmonic flux in Y
307         ENDDO         ENDDO
308        ENDIF        ENDIF
309    
310    C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
311  C-    Advective flux in R  C-    Advective flux in R
312        IF (.NOT. multiDimAdvection .OR.        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
      &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  
      &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.  
      &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN  
 C     Note: wVel needs to be masked  
       IF (K.GE.2) THEN  
313  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
314         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
315          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
# Line 322  C-    Compute vertical advective flux in Line 329  C-    Compute vertical advective flux in
329         ELSE         ELSE
330          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
331         ENDIF         ENDIF
332  C-    Surface "correction" term at k>1 :  C-     add the advective flux to fVerT
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = af(i,j)  
      &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  
      &             rTrans(i,j)*Tracer(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ELSE  
 C-    Surface "correction" term at k=1 :  
333         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
334          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
335           af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)           fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
336          ENDDO          ENDDO
        ENDDO  
       ENDIF  
 C-    add the advective flux to fVerT  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)  
337         ENDDO         ENDDO
       ENDDO  
338        ENDIF        ENDIF
339    
340  C-    Diffusive flux in R  C-    Diffusive flux in R
# Line 358  C           boundary condition. Line 349  C           boundary condition.
349        ELSE        ELSE
350         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
351        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  
352    
353  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
354  C-    GM/Redi flux in R  C-    GM/Redi flux in R
# Line 370  C-    GM/Redi flux in R Line 356  C-    GM/Redi flux in R
356  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
357         CALL GMREDI_RTRANSPORT(         CALL GMREDI_RTRANSPORT(
358       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
359       I     Tracer,       I     Tracer,tracerIdentity,
360       U     df,       U     df,
361       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  
362        ENDIF        ENDIF
363  #endif  #endif
364    
# Line 406  C *note* should update KPP_TRANSPORT_T t Line 387  C *note* should update KPP_TRANSPORT_T t
387       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
388       I     KappaRT,       I     KappaRT,
389       U     df )       U     df )
390    #ifdef ALLOW_PTRACERS
391           ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
392            CALL KPP_TRANSPORT_PTR(
393         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
394         I     tracerIdentity-GAD_TR1+1,KappaRT,
395         U     df )
396    #endif
397         ELSE         ELSE
398            PRINT*,'invalid tracer indentity: ', tracerIdentity
399          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
400         ENDIF         ENDIF
401         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 421  C--   Divergence of fluxes Line 410  C--   Divergence of fluxes
410        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
411         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
412          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
413       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
414       &    *recip_rA(i,j,bi,bj)       &   *( (fZon(i+1,j)-fZon(i,j))
415       &    *(       &     +(fMer(i,j+1)-fMer(i,j))
416       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fVerT(i,j,kUp)-fVerT(i,j,kDown))*rkFac
417       &    +( fMer(i,j+1)-fMer(i,j) )       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
418       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &                   +(vTrans(i,j+1)-vTrans(i,j))
419         &                   +(rTrans(i,j)-rTransKp1(i,j))*rAdvFac
420         &                  )*advFac
421       &    )       &    )
422         ENDDO         ENDDO
423        ENDDO        ENDDO
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22