/[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.22 by dimitri, Thu Sep 25 03:01:59 2003 UTC revision 1.23 by jmc, Wed Jan 7 21:35:00 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, calcAdvection,       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"
 #include "DYNVARS.h"  
45  #include "SURFACE.h"  #include "SURFACE.h"
46  #include "GAD.h"  #include "GAD.h"
47    
# Line 57  C  kdown                :: index into 2 Line 58  C  kdown                :: index into 2
58  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr
59  C  xA,yA                :: areas of X and Y face of tracer cells  C  xA,yA                :: areas of X and Y face of tracer cells
60  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points
61    C  rTransKp1            :: 2-D array of volume transport at W pt, interface k+1
62  C  maskUp               :: 2-D array for mask at W points  C  maskUp               :: 2-D array for mask at W points
63    C  uVel, vVel, wVel     :: 3 components of the velcity field (3-D array)
64  C  diffKh               :: horizontal diffusion coefficient  C  diffKh               :: horizontal diffusion coefficient
65  C  diffK4               :: bi-harmonic diffusion coefficient  C  diffK4               :: bi-harmonic diffusion coefficient
66  C  KappaRT              :: 3-D array for vertical diffusion coefficient  C  KappaRT              :: 3-D array for vertical diffusion coefficient
# Line 65  C  Tracer               :: tracer field Line 68  C  Tracer               :: tracer field
68  C  tracerIdentity       :: identifier for the tracer (required for KPP and GM)  C  tracerIdentity       :: identifier for the tracer (required for KPP and GM)
69  C  advectionScheme      :: advection scheme to use  C  advectionScheme      :: advection scheme to use
70  C  calcAdvection        :: =False if Advec terms computed with multiDim scheme  C  calcAdvection        :: =False if Advec terms computed with multiDim scheme
71    C  implicitAdvection    :: =True if vertical Advec term is computed implicitly
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 73  C  myThid               :: thread number Line 77  C  myThid               :: thread number
77        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
83          _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
84          _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
85        _RL diffKh, diffK4        _RL diffKh, diffK4
86        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87        _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)
88        INTEGER tracerIdentity        INTEGER tracerIdentity
89        INTEGER advectionScheme        INTEGER advectionScheme
90        LOGICAL calcAdvection        LOGICAL calcAdvection
91          LOGICAL implicitAdvection
92        INTEGER myThid        INTEGER myThid
93    
94  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
# Line 103  C  localT               :: local copy of Line 112  C  localT               :: local copy of
112        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115          _RL advFac, rAdvFac
116  CEOP  CEOP
117    
118  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 111  C--   the kDown is still required Line 121  C--   the kDown is still required
121        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
122  #endif  #endif
123    
124          advFac  = 0. _d 0
125          IF (calcAdvection) advFac = 1. _d 0
126          rAdvFac = rkFac*advFac
127          IF (implicitAdvection) rAdvFac = 0. _d 0
128    
129        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
130         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
131          fZon(i,j)      = 0. _d 0          fZon(i,j)      = 0. _d 0
# Line 118  C--   the kDown is still required Line 133  C--   the kDown is still required
133          fVerT(i,j,kUp) = 0. _d 0          fVerT(i,j,kUp) = 0. _d 0
134          df(i,j)        = 0. _d 0          df(i,j)        = 0. _d 0
135          df4(i,j)       = 0. _d 0          df4(i,j)       = 0. _d 0
         localT(i,j)    = 0. _d 0  
136         ENDDO         ENDDO
137        ENDDO        ENDDO
138    
# Line 291  C-    Bi-harmonic flux in Y Line 305  C-    Bi-harmonic flux in Y
305         ENDDO         ENDDO
306        ENDIF        ENDIF
307    
 #ifdef NONLIN_FRSURF  
 C--   Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):  
       IF ( calcAdvection .AND. K.EQ.Nr .AND.  
      &     useRealFreshWaterFlux .AND.  
      &     buoyancyRelation .EQ. 'OCEANICP' ) THEN    
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)  
      &     *rA(i,j,bi,bj)*maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ENDIF  
 #endif /* NONLIN_FRSURF */  
   
308  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):  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 (calcAdvection) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
 C     Note: wVel needs to be masked  
       IF (K.GE.2) THEN  
311  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
312         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
313          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 329  C-    Compute vertical advective flux in Line 327  C-    Compute vertical advective flux in
327         ELSE         ELSE
328          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
329         ENDIF         ENDIF
330  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 :  
331         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
332          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
333           af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)           fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
334          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)  
335         ENDDO         ENDDO
       ENDDO  
336        ENDIF        ENDIF
337    
338  C-    Diffusive flux in R  C-    Diffusive flux in R
# Line 426  C--   Divergence of fluxes Line 408  C--   Divergence of fluxes
408        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
409         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
410          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
411       &   -_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)
412       &    *recip_rA(i,j,bi,bj)       &   *( (fZon(i+1,j)-fZon(i,j))
413       &    *(       &     +(fMer(i,j+1)-fMer(i,j))
414       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fVerT(i,j,kUp)-fVerT(i,j,kDown))*rkFac
415       &    +( fMer(i,j+1)-fMer(i,j) )       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
416       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &                   +(vTrans(i,j+1)-vTrans(i,j))
417         &                   +(rTrans(i,j)-rTransKp1(i,j))*rAdvFac
418         &                  )*advFac
419       &    )       &    )
420         ENDDO         ENDDO
421        ENDDO        ENDDO
422    
 #ifdef NONLIN_FRSURF  
 C-- account for 3.D divergence of the flow in rStar coordinate:  
       IF (calcAdvection .AND. select_rStar.GT.0) THEN  
        DO j=1-Oly,sNy+Oly-1  
         DO i=1-Olx,sNx+Olx-1  
          gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)  
      &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf  
      &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ENDIF  
       IF (calcAdvection .AND. select_rStar.LT.0) THEN  
        DO j=1-Oly,sNy+Oly-1  
         DO i=1-Olx,sNx+Olx-1  
          gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)  
      &     - rStarDhCDt(i,j,bi,bj)  
      &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ENDIF  
 #endif /* NONLIN_FRSURF */  
         
   
423        RETURN        RETURN
424        END        END
Legend:
Removed from v.1.22  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22