/[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.8 by adcroft, Mon Sep 10 01:22:48 2001 UTC revision 1.12.6.1 by heimbach, Sun Mar 24 17:22:22 2002 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6    CBOP
7    C !ROUTINE: GAD_CALC_RHS
8    
9    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,maskUp,
# Line 10  C $Name$ Line 14  C $Name$
14       I           tracerIdentity, advectionScheme,       I           tracerIdentity, advectionScheme,
15       U           fVerT, gTracer,       U           fVerT, gTracer,
16       I           myThid )       I           myThid )
 C     /==========================================================\  
 C     | SUBROUTINE GAD_CALC_RHS                                  |  
 C     |==========================================================|  
 C     \==========================================================/  
       IMPLICIT NONE  
17    
18  C     == GLobal variables ==  C !DESCRIPTION:
19    C Calculates the tendancy of a tracer due to advection and diffusion.
20    C It calculates the fluxes in each direction indepentently and then
21    C sets the tendancy to the divergence of these fluxes. The advective
22    C fluxes are only calculated here when using the linear advection schemes
23    C otherwise only the diffusive and parameterized fluxes are calculated.
24    C
25    C Contributions to the flux are calculated and added:
26    C \begin{equation*}
27    C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
28    C \end{equation*}
29    C
30    C The tendancy is the divergence of the fluxes:
31    C \begin{equation*}
32    C G_\theta = G_\theta + \nabla \cdot {\bf F}
33    C \end{equation*}
34    C
35    C The tendancy is assumed to contain data on entry.
36    
37    C !USES: ===============================================================
38          IMPLICIT NONE
39  #include "SIZE.h"  #include "SIZE.h"
40  #include "EEPARAMS.h"  #include "EEPARAMS.h"
41  #include "PARAMS.h"  #include "PARAMS.h"
# Line 24  C     == GLobal variables == Line 43  C     == GLobal variables ==
43  #include "DYNVARS.h"  #include "DYNVARS.h"
44  #include "GAD.h"  #include "GAD.h"
45    
46  C     == Routine arguments ==  #ifdef ALLOW_AUTODIFF_TAMC
47        INTEGER k,kUp,kDown,kM1  #include "tamc.h"
48    #include "tamc_keys.h"
49    #endif /* ALLOW_AUTODIFF_TAMC */
50    
51    C !INPUT PARAMETERS: ===================================================
52    C  bi,bj                :: tile indices
53    C  iMin,iMax,jMin,jMax  :: loop range for called routines
54    C  kup                  :: index into 2 1/2D array, toggles between 1 and 2
55    C  kdown                :: index into 2 1/2D array, toggles between 2 and 1
56    C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr
57    C  xA,yA                :: areas of X and Y face of tracer cells
58    C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points
59    C  maskUp               :: 2-D array for mask at W points
60    C  diffKh               :: horizontal diffusion coefficient
61    C  diffK4               :: bi-harmonic diffusion coefficient
62    C  KappaRT              :: 3-D array for vertical diffusion coefficient
63    C  Tracer               :: tracer field
64    C  tracerIdentity       :: identifier for the tracer (required only for KPP)
65    C  advectionScheme      :: advection scheme to use
66    C  myThid               :: thread number
67        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
68          INTEGER k,kUp,kDown,kM1
69        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 38  C     == Routine arguments == Line 77  C     == Routine arguments ==
77        _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)
78        INTEGER tracerIdentity        INTEGER tracerIdentity
79        INTEGER advectionScheme        INTEGER advectionScheme
       _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
80        INTEGER myThid        INTEGER myThid
81    
82  C     == Local variables ==  C !OUTPUT PARAMETERS: ==================================================
83  C     I, J, K - Loop counters  C  gTracer              :: tendancy array
84    C  fVerT                :: 2 1/2D arrays for vertical advective flux
85          _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
86          _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
87    
88    C !LOCAL VARIABLES: ====================================================
89    C  i,j                  :: loop indices
90    C  df4                  :: used for storing del^2 T for bi-harmonic term
91    C  fZon                 :: zonal flux
92    C  fmer                 :: meridional flux
93    C  af                   :: advective flux
94    C  df                   :: diffusive flux
95    C  localT               :: local copy of tracer field
96        INTEGER i,j        INTEGER i,j
       LOGICAL TOP_LAYER  
       _RL afFacT, dfFacT  
97        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103    CEOP
104    
105  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
106  C--   only the kUp part of fverT is set in this subroutine  C--   only the kUp part of fverT is set in this subroutine
107  C--   the kDown is still required  C--   the kDown is still required
108        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
109  #endif  #endif
110    
111        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
112         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
113          fZon(i,j)      = 0.0          fZon(i,j)      = 0. _d 0
114          fMer(i,j)      = 0.0          fMer(i,j)      = 0. _d 0
115          fVerT(i,j,kUp) = 0.0          fVerT(i,j,kUp) = 0. _d 0
116            df(i,j)        = 0. _d 0
117            df4(i,j)       = 0. _d 0
118            localT(i,j)    = 0. _d 0
119         ENDDO         ENDDO
120        ENDDO        ENDDO
121    
       afFacT = 1. _d 0  
       dfFacT = 1. _d 0  
       TOP_LAYER = K .EQ. 1  
   
122  C--   Make local copy of tracer array  C--   Make local copy of tracer array
123        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
124         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 86  C     the tendency to zero. Line 134  C     the tendency to zero.
134       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
135         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
136          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
137           gTracer(i,j,k,bi,bj)=0.           gTracer(i,j,k,bi,bj)=0. _d 0
138          ENDDO          ENDDO
139         ENDDO         ENDDO
140        ENDIF        ENDIF
# Line 101  C--   Pre-calculate del^2 T if bi-harmon Line 149  C--   Pre-calculate del^2 T if bi-harmon
149  C--   Initialize net flux in X direction  C--   Initialize net flux in X direction
150        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
151         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
152          fZon(i,j) = 0.          fZon(i,j) = 0. _d 0
153         ENDDO         ENDDO
154        ENDDO        ENDDO
155    
# Line 141  C-    Diffusive flux in X Line 189  C-    Diffusive flux in X
189        ELSE        ELSE
190         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
191          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
192           df(i,j) = 0.           df(i,j) = 0. _d 0
193          ENDDO          ENDDO
194         ENDDO         ENDDO
195        ENDIF        ENDIF
# Line 176  C-    Bi-harmonic duffusive flux in X Line 224  C-    Bi-harmonic duffusive flux in X
224  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
225        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
226         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
227          fMer(i,j) = 0.          fMer(i,j) = 0. _d 0
228         ENDDO         ENDDO
229        ENDDO        ENDDO
230    
# Line 216  C-    Diffusive flux in Y Line 264  C-    Diffusive flux in Y
264        ELSE        ELSE
265         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
266          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
267           df(i,j) = 0.           df(i,j) = 0. _d 0
268          ENDDO          ENDDO
269         ENDDO         ENDDO
270        ENDIF        ENDIF
# Line 248  C-    Bi-harmonic flux in Y Line 296  C-    Bi-harmonic flux in Y
296         ENDDO         ENDDO
297        ENDIF        ENDIF
298    
 C--   Initialize net flux in R  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = 0.  
        ENDDO  
       ENDDO  
   
299  C-    Advective flux in R  C-    Advective flux in R
300        IF (.NOT. multiDimAdvection .OR.        IF (.NOT. multiDimAdvection .OR.
301       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
# Line 273  C-    Compute vertical advective flux in Line 314  C-    Compute vertical advective flux in
314         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
315          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
316         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
317  c       CALL GAD_DST3_ADV_R(          CALL GAD_DST3_ADV_R(
318  c    &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
         STOP 'GAD_CALC_RHS: GAD_DST3_ADV_R not coded yet'  
319         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
320  c       CALL GAD_DST3FL_ADV_R(          CALL GAD_DST3FL_ADV_R(
321  c    &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
         STOP 'GAD_CALC_RHS: GAD_DST3FL_ADV_R not coded yet'  
322         ELSE         ELSE
323          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
324         ENDIF         ENDIF
# Line 302  C-    Surface "correction" term at k=1 : Line 341  C-    Surface "correction" term at k=1 :
341  C-    add the advective flux to fVerT  C-    add the advective flux to fVerT
342        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
343         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
344          fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)          fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
345         ENDDO         ENDDO
346        ENDDO        ENDDO
347        ENDIF        ENDIF
# Line 313  C           boundary condition. Line 352  C           boundary condition.
352        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
353         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
354          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
355           df(i,j) = 0.           df(i,j) = 0. _d 0
356          ENDDO          ENDDO
357         ENDDO         ENDDO
358        ELSE        ELSE
# Line 321  C           boundary condition. Line 360  C           boundary condition.
360        ENDIF        ENDIF
361  c     DO j=1-Oly,sNy+Oly  c     DO j=1-Oly,sNy+Oly
362  c      DO i=1-Olx,sNx+Olx  c      DO i=1-Olx,sNx+Olx
363  c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
364  c      ENDDO  c      ENDDO
365  c     ENDDO  c     ENDDO
366    
# Line 336  C *note* should update GMREDI_RTRANSPORT Line 375  C *note* should update GMREDI_RTRANSPORT
375       I     myThid)       I     myThid)
376  c      DO j=1-Oly,sNy+Oly  c      DO j=1-Oly,sNy+Oly
377  c       DO i=1-Olx,sNx+Olx  c       DO i=1-Olx,sNx+Olx
378  c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
379  c       ENDDO  c       ENDDO
380  c      ENDDO  c      ENDDO
381        ENDIF        ENDIF
# Line 344  c      ENDDO Line 383  c      ENDDO
383    
384        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
385         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
386          fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)          fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
387         ENDDO         ENDDO
388        ENDDO        ENDDO
389    
# Line 353  C-    Add non local KPP transport term ( Line 392  C-    Add non local KPP transport term (
392        IF (useKPP) THEN        IF (useKPP) THEN
393         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
394          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
395           df(i,j) = 0.           df(i,j) = 0. _d 0
396          ENDDO          ENDDO
397         ENDDO         ENDDO
398         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
# Line 372  C *note* should update KPP_TRANSPORT_T t Line 411  C *note* should update KPP_TRANSPORT_T t
411         ENDIF         ENDIF
412         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
413          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
414           fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
415          ENDDO          ENDDO
416         ENDDO         ENDDO
417        ENDIF        ENDIF
418  #endif  #endif
419    
420  C--   Divergence of fluxes  C--   Divergence of fluxes
421        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly-1
422         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx-1
423          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
424       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
425       &    *recip_rA(i,j,bi,bj)       &    *recip_rA(i,j,bi,bj)
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.12.6.1
  ViewVC Help
Powered by ViewVC 1.1.22