/[MITgcm]/MITgcm/pkg/mom_fluxform/mom_fluxform.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_fluxform/mom_fluxform.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by adcroft, Fri Aug 17 18:40:30 2001 UTC revision 1.15 by jmc, Sat Oct 11 16:37:55 2003 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  CBOI
5    C !TITLE: pkg/mom\_advdiff
6    C !AUTHORS: adcroft@mit.edu
7    C !INTRODUCTION: Flux-form Momentum Equations Package
8    C
9    C Package "mom\_fluxform" provides methods for calculating explicit terms
10    C in the momentum equation cast in flux-form:
11    C \begin{eqnarray*}
12    C G^u & = & -\frac{1}{\rho} \partial_x \phi_h
13    C           -\nabla \cdot {\bf v} u
14    C           -fv
15    C           +\frac{1}{\rho} \nabla \cdot {\bf \tau}^x
16    C           + \mbox{metrics}
17    C \\
18    C G^v & = & -\frac{1}{\rho} \partial_y \phi_h
19    C           -\nabla \cdot {\bf v} v
20    C           +fu
21    C           +\frac{1}{\rho} \nabla \cdot {\bf \tau}^y
22    C           + \mbox{metrics}
23    C \end{eqnarray*}
24    C where ${\bf v}=(u,v,w)$ and $\tau$, the stress tensor, includes surface
25    C stresses as well as internal viscous stresses.
26    CEOI
27    
28    #include "MOM_FLUXFORM_OPTIONS.h"
29    
30    CBOP
31    C !ROUTINE: MOM_FLUXFORM
32    
33    C !INTERFACE: ==========================================================
34        SUBROUTINE MOM_FLUXFORM(        SUBROUTINE MOM_FLUXFORM(
35       I        bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
36       I        phi_hyd,KappaRU,KappaRV,       I        dPhihydX,dPhiHydY,KappaRU,KappaRV,
37       U        fVerU, fVerV,       U        fVerU, fVerV,
38       I        myCurrentTime, myIter, myThid)       I        myTime,myIter,myThid)
39  C     /==========================================================\  
40  C     | S/R MOM_FLUXFORM                                         |  C !DESCRIPTION:
41  C     | o Form the right hand-side of the momentum equation.     |  C Calculates all the horizontal accelerations except for the implicit surface
42  C     |==========================================================|  C pressure gradient and implciit vertical viscosity.
 C     | Terms are evaluated one layer at a time working from     |  
 C     | the bottom to the top. The vertically integrated         |  
 C     | barotropic flow tendency term is evluated by summing the |  
 C     | tendencies.                                              |  
 C     | Notes:                                                   |  
 C     | We have not sorted out an entirely satisfactory formula  |  
 C     | for the diffusion equation bc with lopping. The present  |  
 C     | form produces a diffusive flux that does not scale with  |  
 C     | open-area. Need to do something to solidfy this and to   |  
 C     | deal "properly" with thin walls.                         |  
 C     \==========================================================/  
       IMPLICIT NONE  
43    
44    C !USES: ===============================================================
45  C     == Global variables ==  C     == Global variables ==
46          IMPLICIT NONE
47  #include "SIZE.h"  #include "SIZE.h"
48  #include "DYNVARS.h"  #include "DYNVARS.h"
49  #include "FFIELDS.h"  #include "FFIELDS.h"
# Line 34  C     == Global variables == Line 52  C     == Global variables ==
52  #include "GRID.h"  #include "GRID.h"
53  #include "SURFACE.h"  #include "SURFACE.h"
54    
55  C     == Routine arguments ==  C !INPUT PARAMETERS: ===================================================
56  C     fZon    - Work array for flux of momentum in the east-west  C  bi,bj                :: tile indices
57  C               direction at the west face of a cell.  C  iMin,iMax,jMin,jMAx  :: loop ranges
58  C     fMer    - Work array for flux of momentum in the north-south  C  k                    :: vertical level
59  C               direction at the south face of a cell.  C  kUp                  :: =1 or 2 for consecutive k
60  C     fVerU   - Flux of momentum in the vertical  C  kDown                :: =2 or 1 for consecutive k
61  C     fVerV     direction out of the upper face of a cell K  C  dPhiHydX,Y           :: Gradient (X & Y dir.) of Hydrostatic Potential
62  C               ( flux into the cell above ).  C  KappaRU              :: vertical viscosity
63  C     phi_hyd - Hydrostatic pressure  C  KappaRV              :: vertical viscosity
64  C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation  C  fVerU                :: vertical flux of U, 2 1/2 dim for pipe-lining
65  C                                      results will be set.  C  fVerV                :: vertical flux of V, 2 1/2 dim for pipe-lining
66  C     kUp, kDown                     - Index for upper and lower layers.  C  myTime               :: current time
67  C     myThid - Instance number for this innvocation of CALC_MOM_RHS  C  myIter               :: current time-step number
68        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C  myThid               :: thread number
69          INTEGER bi,bj,iMin,iMax,jMin,jMax
70          INTEGER k,kUp,kDown
71          _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72          _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
74        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
75        _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
76        _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
77        INTEGER kUp,kDown        _RL     myTime
       _RL     myCurrentTime  
78        INTEGER myIter        INTEGER myIter
79        INTEGER myThid        INTEGER myThid
       INTEGER bi,bj,iMin,iMax,jMin,jMax  
80    
81  C     == Local variables ==  C !OUTPUT PARAMETERS: ==================================================
82  C     ab15, ab05    - Weights for Adams-Bashforth time stepping scheme.  C None - updates gU() and gV() in common blocks
83  C     i,j,k         - Loop counters  
84    C !LOCAL VARIABLES: ====================================================
85    C  i,j                  :: loop indices
86    C  aF                   :: advective flux
87    C  vF                   :: viscous flux
88    C  v4F                  :: bi-harmonic viscous flux
89    C  vrF                  :: vertical viscous flux
90    C  cF                   :: Coriolis acceleration
91    C  mT                   :: Metric terms
92    C  pF                   :: Pressure gradient
93    C  fZon                 :: zonal fluxes
94    C  fMer                 :: meridional fluxes
95          INTEGER i,j
96          _RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97          _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98          _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99          _RL vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100          _RL cF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101          _RL mT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102          _RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103          _RL fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104          _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105  C     wMaskOverride - Land sea flag override for top layer.  C     wMaskOverride - Land sea flag override for top layer.
106  C     afFacMom      - Tracer parameters for turning terms  C     afFacMom      - Tracer parameters for turning terms
107  C     vfFacMom        on and off.  C     vfFacMom        on and off.
# Line 70  C     mTFacMom        pfFacMom - Pressur Line 111  C     mTFacMom        pfFacMom - Pressur
111  C                     cfFacMom - Coriolis terms  C                     cfFacMom - Coriolis terms
112  C                     foFacMom - Forcing  C                     foFacMom - Forcing
113  C                     mTFacMom - Metric term  C                     mTFacMom - Metric term
 C     vF            - Temporary holding viscous term (Laplacian)  
 C     v4F           - Temporary holding viscous term (Biharmonic)  
 C     cF            - Temporary holding coriolis term.  
 C     mT            - Temporary holding metric terms(s).  
 C     pF            - Temporary holding pressure|potential gradient terms.  
114  C     uDudxFac, AhDudxFac, etc ... individual term tracer parameters  C     uDudxFac, AhDudxFac, etc ... individual term tracer parameters
       _RL      aF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      vrF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      cF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      pF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL    fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL    fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
115        _RS    hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RS  r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS  r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RS      xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS      xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 93  C     uDudxFac, AhDudxFac, etc ... indiv Line 120  C     uDudxFac, AhDudxFac, etc ... indiv
120        _RL  vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL  vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL  uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL  uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL  vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL  vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123          _RL  rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124          _RL  rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125  C     I,J,K - Loop counters  C     I,J,K - Loop counters
       INTEGER i,j,k  
126  C     rVelMaskOverride - Factor for imposing special surface boundary conditions  C     rVelMaskOverride - Factor for imposing special surface boundary conditions
127  C                        ( set according to free-surface condition ).  C                        ( set according to free-surface condition ).
128  C     hFacROpen        - Lopped cell factos used tohold fraction of open  C     hFacROpen        - Lopped cell factos used tohold fraction of open
# Line 124  C     xxxFac - On-off tracer parameters Line 152  C     xxxFac - On-off tracer parameters
152        _RL  phyFac        _RL  phyFac
153        _RL  vForcFac        _RL  vForcFac
154        _RL  mtFacV        _RL  mtFacV
 C     ab05, ab15 - Adams-Bashforth time-stepping weights.  
       _RL  ab05, ab15  
155        INTEGER km1,kp1        INTEGER km1,kp1
156        _RL wVelBottomOverride        _RL wVelBottomOverride
157        LOGICAL bottomDragTerms        LOGICAL bottomDragTerms
158        _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159    CEOP
160    
161        km1=MAX(1,k-1)        km1=MAX(1,k-1)
162        kp1=MIN(Nr,k+1)        kp1=MIN(Nr,k+1)
# Line 150  C     Initialise intermediate terms Line 177  C     Initialise intermediate terms
177          pF(i,j)   = 0.          pF(i,j)   = 0.
178          fZon(i,j) = 0.          fZon(i,j) = 0.
179          fMer(i,j) = 0.          fMer(i,j) = 0.
180            rTransU(i,j) = 0.
181            rTransV(i,j) = 0.
182    #ifdef ALLOW_AUTODIFF_TAMC
183    C- jmc: this is wrong, but at least with #ifdef/endif TAMC, it does not break
184    C       the forward code ; (same thing in mom_vectinv)
185            fVerU(i,j,1) = 0. _d 0
186            fVerU(i,j,2) = 0. _d 0
187            fVerV(i,j,1) = 0. _d 0
188            fVerV(i,j,2) = 0. _d 0
189    #endif
190         ENDDO         ENDDO
191        ENDDO        ENDDO
192    
# Line 194  C-- with stagger time stepping, grad Phi Line 231  C-- with stagger time stepping, grad Phi
231          phyFac = 0.          phyFac = 0.
232        ENDIF        ENDIF
233    
 C--   Adams-Bashforth weighting factors  
       ab15   =  1.5 _d 0 + abEps  
       ab05   = -0.5 _d 0 - abEps  
     
234  C--   Calculate open water fraction at vorticity points  C--   Calculate open water fraction at vorticity points
235        CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)        CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
236    
# Line 230  C     Calculate velocity field "volume t Line 263  C     Calculate velocity field "volume t
263    
264        CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)        CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)
265    
266    C---  First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
267          IF (momAdvection.AND.k.EQ.1) THEN
268    
269    C-    Calculate vertical transports above U & V points (West & South face):
270           CALL MOM_CALC_RTRANS( k, bi, bj,
271         O                       rTransU, rTransV,
272         I                       myTime, myIter, myThid)
273    
274    C-    Free surface correction term (flux at k=1)
275           CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,rTransU,af,myThid)
276           DO j=jMin,jMax
277            DO i=iMin,iMax
278             fVerU(i,j,kUp) = af(i,j)
279            ENDDO
280           ENDDO
281    
282           CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,rTransV,af,myThid)
283           DO j=jMin,jMax
284            DO i=iMin,iMax
285             fVerV(i,j,kUp) = af(i,j)
286            ENDDO
287           ENDDO
288    
289    C---  endif momAdvection & k=1
290          ENDIF
291    
292    
293    C---  Calculate vertical transports (at k+1) below U & V points :
294          IF (momAdvection) THEN
295           CALL MOM_CALC_RTRANS( k+1, bi, bj,
296         O                       rTransU, rTransV,
297         I                       myTime, myIter, myThid)
298          ENDIF
299    
300    
301  C---- Zonal momentum equation starts here  C---- Zonal momentum equation starts here
302    
303  C     Bi-harmonic term del^2 U -> v4F  C     Bi-harmonic term del^2 U -> v4F
304        IF (momViscosity)        IF (momViscosity .AND. viscA4.NE.0. )
305       & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)       & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
306    
307  C---  Calculate mean and eddy fluxes between cells for zonal flow.  C---  Calculate mean and eddy fluxes between cells for zonal flow.
# Line 266  C     Laplacian and bi-harmonic term Line 334  C     Laplacian and bi-harmonic term
334       & CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)       & CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
335    
336  C     Combine fluxes -> fMer  C     Combine fluxes -> fMer
337        DO j=jMin,jMax        DO j=jMin,jMax+1
338         DO i=iMin,iMax         DO i=iMin,iMax
339          fMer(i,j) = vDudyFac*aF(i,j) + AhDudyFac*vF(i,j)          fMer(i,j) = vDudyFac*aF(i,j) + AhDudyFac*vF(i,j)
340         ENDDO         ENDDO
# Line 274  C     Combine fluxes -> fMer Line 342  C     Combine fluxes -> fMer
342    
343  C--   Vertical flux (fVer is at upper face of "u" cell)  C--   Vertical flux (fVer is at upper face of "u" cell)
344    
 C--   Free surface correction term (flux at k=1)  
       IF (momAdvection.AND.k.EQ.1) THEN  
        CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,af,myThid)  
        DO j=jMin,jMax  
         DO i=iMin,iMax  
          fVerU(i,j,kUp) = af(i,j)  
         ENDDO  
        ENDDO  
       ENDIF  
345  C     Mean flow component of vertical flux (at k+1) -> aF  C     Mean flow component of vertical flux (at k+1) -> aF
346        IF (momAdvection)        IF (momAdvection)
347       & CALL MOM_U_ADV_WU(bi,bj,k+1,uVel,wVel,af,myThid)       & CALL MOM_U_ADV_WU(bi,bj,k+1,uVel,wVel,rTransU,af,myThid)
348    
349  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
350        IF (momViscosity.AND..NOT.implicitViscosity)        IF (momViscosity.AND..NOT.implicitViscosity)
# Line 298  C     Combine fluxes Line 357  C     Combine fluxes
357         ENDDO         ENDDO
358        ENDDO        ENDDO
359    
 C---  Hydrostatic term ( -1/rhoConst . dphi/dx )  
       IF (momPressureForcing) THEN  
        DO j=jMin,jMax  
         DO i=iMin,iMax  
          pf(i,j) = - _recip_dxC(i,j,bi,bj)  
      &    *(phi_hyd(i,j,k)-phi_hyd(i-1,j,k))  
         ENDDO  
        ENDDO  
       ENDIF  
   
360  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
361        DO j=jMin,jMax        DO j=jMin,jMax
362         DO i=iMin,iMax         DO i=iMin,iMax
# Line 323  C--   Tendency is minus divergence of th Line 372  C--   Tendency is minus divergence of th
372       &   +fMer(i,j+1)          - fMer(i  ,j)       &   +fMer(i,j+1)          - fMer(i  ,j)
373       &   +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac       &   +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac
374       &   )       &   )
375       & _PHM( +phxFac * pf(i,j) )       &  - phxFac*dPhiHydX(i,j)
376         ENDDO         ENDDO
377        ENDDO        ENDDO
378    
379    #ifdef NONLIN_FRSURF
380    C-- account for 3.D divergence of the flow in rStar coordinate:
381          IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
382           DO j=jMin,jMax
383            DO i=iMin,iMax
384             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
385         &     - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
386         &       *uVel(i,j,k,bi,bj)
387            ENDDO
388           ENDDO
389          ENDIF
390          IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
391           DO j=jMin,jMax
392            DO i=iMin,iMax
393             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
394         &     - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
395            ENDDO
396           ENDDO
397          ENDIF
398    #endif /* NONLIN_FRSURF */
399    
400  C-- No-slip and drag BCs appear as body forces in cell abutting topography  C-- No-slip and drag BCs appear as body forces in cell abutting topography
401        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
402  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
# Line 347  C-    No-slip BCs impose a drag at botto Line 417  C-    No-slip BCs impose a drag at botto
417         ENDDO         ENDDO
418        ENDIF        ENDIF
419    
420  C--   Forcing term  C--   Forcing term (moved to timestep.F)
421        IF (momForcing)  c     IF (momForcing)
422       &  CALL EXTERNAL_FORCING_U(  c    &  CALL EXTERNAL_FORCING_U(
423       I     iMin,iMax,jMin,jMax,bi,bj,k,  c    I     iMin,iMax,jMin,jMax,bi,bj,k,
424       I     myCurrentTime,myThid)  c    I     myTime,myThid)
425    
426  C--   Metric terms for curvilinear grid systems  C--   Metric terms for curvilinear grid systems
427        IF (usingSphericalPolarMTerms) THEN        IF (useNHMTerms) THEN
428  C      o Spherical polar grid metric terms  C      o Non-hydrosatic metric terms
429         CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)         CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
430         DO j=jMin,jMax         DO j=jMin,jMax
431          DO i=iMin,iMax          DO i=iMin,iMax
432           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
433          ENDDO          ENDDO
434         ENDDO         ENDDO
435          ENDIF
436          IF (usingSphericalPolarMTerms) THEN
437         CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)         CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)
438         DO j=jMin,jMax         DO j=jMin,jMax
439          DO i=iMin,iMax          DO i=iMin,iMax
# Line 381  C--   Set du/dt on boundaries to zero Line 453  C--   Set du/dt on boundaries to zero
453  C---- Meridional momentum equation starts here  C---- Meridional momentum equation starts here
454    
455  C     Bi-harmonic term del^2 V -> v4F  C     Bi-harmonic term del^2 V -> v4F
456        IF (momViscosity)        IF (momViscosity .AND. viscA4.NE.0. )
457       & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)       & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
458    
459  C---  Calculate mean and eddy fluxes between cells for meridional flow.  C---  Calculate mean and eddy fluxes between cells for meridional flow.
# Line 398  C     Laplacian and bi-harmonic terms -> Line 470  C     Laplacian and bi-harmonic terms ->
470    
471  C     Combine fluxes -> fZon  C     Combine fluxes -> fZon
472        DO j=jMin,jMax        DO j=jMin,jMax
473         DO i=iMin,iMax         DO i=iMin,iMax+1
474          fZon(i,j) = uDvdxFac*aF(i,j) + AhDvdxFac*vF(i,j)          fZon(i,j) = uDvdxFac*aF(i,j) + AhDvdxFac*vF(i,j)
475         ENDDO         ENDDO
476        ENDDO        ENDDO
# Line 422  C     Combine fluxes -> fMer Line 494  C     Combine fluxes -> fMer
494    
495  C--   Vertical flux (fVer is at upper face of "v" cell)  C--   Vertical flux (fVer is at upper face of "v" cell)
496    
 C--   Free surface correction term (flux at k=1)  
       IF (momAdvection.AND.k.EQ.1) THEN  
        CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,af,myThid)  
        DO j=jMin,jMax  
         DO i=iMin,iMax  
          fVerV(i,j,kUp) = af(i,j)  
         ENDDO  
        ENDDO  
       ENDIF  
497  C     o Mean flow component of vertical flux  C     o Mean flow component of vertical flux
498        IF (momAdvection)        IF (momAdvection)
499       & CALL MOM_V_ADV_WV(bi,bj,k+1,vVel,wVel,af,myThid)       & CALL MOM_V_ADV_WV(bi,bj,k+1,vVel,wVel,rTransV,af,myThid)
500    
501  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
502        IF (momViscosity.AND..NOT.implicitViscosity)        IF (momViscosity.AND..NOT.implicitViscosity)
# Line 446  C     Combine fluxes -> fVerV Line 509  C     Combine fluxes -> fVerV
509         ENDDO         ENDDO
510        ENDDO        ENDDO
511    
 C---  Hydorstatic term (-1/rhoConst . dphi/dy )  
       IF (momPressureForcing) THEN  
        DO j=jMin,jMax  
         DO i=iMin,iMax  
          pF(i,j) = -_recip_dyC(i,j,bi,bj)  
      &    *(phi_hyd(i,j,k)-phi_hyd(i,j-1,k))  
         ENDDO  
        ENDDO  
       ENDIF  
   
512  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
513        DO j=jMin,jMax        DO j=jMin,jMax
514         DO i=iMin,iMax         DO i=iMin,iMax
# Line 471  C--   Tendency is minus divergence of th Line 524  C--   Tendency is minus divergence of th
524       &   +fMer(i,j  )          - fMer(i,j-1)       &   +fMer(i,j  )          - fMer(i,j-1)
525       &   +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac       &   +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac
526       &   )       &   )
527       & _PHM( +phyFac*pf(i,j) )       &  - phyFac*dPhiHydY(i,j)
528         ENDDO         ENDDO
529        ENDDO        ENDDO
530    
531    #ifdef NONLIN_FRSURF
532    C-- account for 3.D divergence of the flow in rStar coordinate:
533          IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
534           DO j=jMin,jMax
535            DO i=iMin,iMax
536             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
537         &     - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
538         &       *vVel(i,j,k,bi,bj)
539            ENDDO
540           ENDDO
541          ENDIF
542          IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
543           DO j=jMin,jMax
544            DO i=iMin,iMax
545             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
546         &     - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
547            ENDDO
548           ENDDO
549          ENDIF
550    #endif /* NONLIN_FRSURF */
551    
552  C-- No-slip and drag BCs appear as body forces in cell abutting topography  C-- No-slip and drag BCs appear as body forces in cell abutting topography
553        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
554  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
# Line 495  C-    No-slip BCs impose a drag at botto Line 569  C-    No-slip BCs impose a drag at botto
569         ENDDO         ENDDO
570        ENDIF        ENDIF
571    
572  C--   Forcing term  C--   Forcing term (moved to timestep.F)
573        IF (momForcing)  c     IF (momForcing)
574       & CALL EXTERNAL_FORCING_V(  c    & CALL EXTERNAL_FORCING_V(
575       I     iMin,iMax,jMin,jMax,bi,bj,k,  c    I     iMin,iMax,jMin,jMax,bi,bj,k,
576       I     myCurrentTime,myThid)  c    I     myTime,myThid)
577    
578  C--   Metric terms for curvilinear grid systems  C--   Metric terms for curvilinear grid systems
579        IF (usingSphericalPolarMTerms) THEN        IF (useNHMTerms) THEN
580  C      o Spherical polar grid metric terms  C      o Spherical polar grid metric terms
581         CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)         CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
582         DO j=jMin,jMax         DO j=jMin,jMax
# Line 510  C      o Spherical polar grid metric ter Line 584  C      o Spherical polar grid metric ter
584           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
585          ENDDO          ENDDO
586         ENDDO         ENDDO
587          ENDIF
588          IF (usingSphericalPolarMTerms) THEN
589         CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)         CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
590         DO j=jMin,jMax         DO j=jMin,jMax
591          DO i=iMin,iMax          DO i=iMin,iMax
# Line 527  C--   Set dv/dt on boundaries to zero Line 603  C--   Set dv/dt on boundaries to zero
603    
604  C--   Coriolis term  C--   Coriolis term
605  C     Note. As coded here, coriolis will not work with "thin walls"  C     Note. As coded here, coriolis will not work with "thin walls"
606  #ifdef INCLUDE_CD_CODE  c     IF (useCDscheme) THEN
607        CALL MOM_CDSCHEME(bi,bj,k,phi_hyd,myThid)  c       CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)
608  #else  c     ELSE
609        CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)        IF (.NOT.useCDscheme) THEN
610        DO j=jMin,jMax          CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
611         DO i=iMin,iMax          DO j=jMin,jMax
612          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)           DO i=iMin,iMax
613         ENDDO            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
614        ENDDO           ENDDO
615        CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)          ENDDO
616        DO j=jMin,jMax          CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
617         DO i=iMin,iMax          DO j=jMin,jMax
618          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)           DO i=iMin,iMax
619              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
620             ENDDO
621            ENDDO
622          ENDIF
623    
624          IF (nonHydrostatic.OR.quasiHydrostatic) THEN
625           CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
626           DO j=jMin,jMax
627            DO i=iMin,iMax
628             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
629            ENDDO
630         ENDDO         ENDDO
631        ENDDO        ENDIF
 #endif /* INCLUDE_CD_CODE */  
632    
633        RETURN        RETURN
634        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22