/[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.7 by adcroft, Thu Nov 7 21:51:15 2002 UTC revision 1.14 by heimbach, Fri Oct 10 23:00:01 2003 UTC
# Line 25  C where ${\bf v}=(u,v,w)$ and $\tau$, th Line 25  C where ${\bf v}=(u,v,w)$ and $\tau$, th
25  C stresses as well as internal viscous stresses.  C stresses as well as internal viscous stresses.
26  CEOI  CEOI
27    
28  #include "CPP_OPTIONS.h"  #include "MOM_FLUXFORM_OPTIONS.h"
29    
30  CBOP  CBOP
31  C !ROUTINE: MOM_FLUXFORM  C !ROUTINE: MOM_FLUXFORM
# Line 33  C !ROUTINE: MOM_FLUXFORM Line 33  C !ROUTINE: MOM_FLUXFORM
33  C !INTERFACE: ==========================================================  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    
40  C !DESCRIPTION:  C !DESCRIPTION:
41  C Calculates all the horizontal accelerations except for the implicit surface  C Calculates all the horizontal accelerations except for the implicit surface
# Line 58  C  iMin,iMax,jMin,jMAx  :: loop ranges Line 58  C  iMin,iMax,jMin,jMAx  :: loop ranges
58  C  k                    :: vertical level  C  k                    :: vertical level
59  C  kUp                  :: =1 or 2 for consecutive k  C  kUp                  :: =1 or 2 for consecutive k
60  C  kDown                :: =2 or 1 for consecutive k  C  kDown                :: =2 or 1 for consecutive k
61  C  phi_hyd              :: hydrostatic pressure (perturbation)  C  dPhiHydX,Y           :: Gradient (X & Y dir.) of Hydrostatic Potential
62  C  KappaRU              :: vertical viscosity  C  KappaRU              :: vertical viscosity
63  C  KappaRV              :: vertical viscosity  C  KappaRV              :: vertical viscosity
64  C  fVerU                :: vertical flux of U, 2 1/2 dim for pipe-lining  C  fVerU                :: vertical flux of U, 2 1/2 dim for pipe-lining
65  C  fVerV                :: vertical flux of V, 2 1/2 dim for pipe-lining  C  fVerV                :: vertical flux of V, 2 1/2 dim for pipe-lining
66  C  myCurrentTime        :: current time  C  myTime               :: current time
67  C  myIter               :: current time-step number  C  myIter               :: current time-step number
68  C  myThid               :: thread number  C  myThid               :: thread number
69        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
70        INTEGER k,kUp,kDown        INTEGER k,kUp,kDown
71        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _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        _RL     myCurrentTime        _RL     myTime
78        INTEGER myIter        INTEGER myIter
79        INTEGER myThid        INTEGER myThid
80    
# Line 119  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
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 ).
# Line 174  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            fVerU(i,j,1) = 0. _d 0
183            fVerU(i,j,2) = 0. _d 0
184            fVerV(i,j,1) = 0. _d 0
185            fVerV(i,j,2) = 0. _d 0
186         ENDDO         ENDDO
187        ENDDO        ENDDO
188    
# Line 250  C     Calculate velocity field "volume t Line 259  C     Calculate velocity field "volume t
259    
260        CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)        CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)
261    
262    C---  First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
263          IF (momAdvection.AND.k.EQ.1) THEN
264    
265    C-    Calculate vertical transports above U & V points (West & South face):
266           CALL MOM_CALC_RTRANS( k, bi, bj,
267         O                       rTransU, rTransV,
268         I                       myTime, myIter, myThid)
269    
270    C-    Free surface correction term (flux at k=1)
271           CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,rTransU,af,myThid)
272           DO j=jMin,jMax
273            DO i=iMin,iMax
274             fVerU(i,j,kUp) = af(i,j)
275            ENDDO
276           ENDDO
277    
278           CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,rTransV,af,myThid)
279           DO j=jMin,jMax
280            DO i=iMin,iMax
281             fVerV(i,j,kUp) = af(i,j)
282            ENDDO
283           ENDDO
284    
285    C---  endif momAdvection & k=1
286          ENDIF
287    
288    
289    C---  Calculate vertical transports (at k+1) below U & V points :
290          IF (momAdvection) THEN
291           CALL MOM_CALC_RTRANS( k+1, bi, bj,
292         O                       rTransU, rTransV,
293         I                       myTime, myIter, myThid)
294          ENDIF
295    
296    
297  C---- Zonal momentum equation starts here  C---- Zonal momentum equation starts here
298    
299  C     Bi-harmonic term del^2 U -> v4F  C     Bi-harmonic term del^2 U -> v4F
300        IF (momViscosity)        IF (momViscosity .AND. viscA4.NE.0. )
301       & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)       & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
302    
303  C---  Calculate mean and eddy fluxes between cells for zonal flow.  C---  Calculate mean and eddy fluxes between cells for zonal flow.
# Line 286  C     Laplacian and bi-harmonic term Line 330  C     Laplacian and bi-harmonic term
330       & CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)       & CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
331    
332  C     Combine fluxes -> fMer  C     Combine fluxes -> fMer
333        DO j=jMin,jMax        DO j=jMin,jMax+1
334         DO i=iMin,iMax         DO i=iMin,iMax
335          fMer(i,j) = vDudyFac*aF(i,j) + AhDudyFac*vF(i,j)          fMer(i,j) = vDudyFac*aF(i,j) + AhDudyFac*vF(i,j)
336         ENDDO         ENDDO
# Line 294  C     Combine fluxes -> fMer Line 338  C     Combine fluxes -> fMer
338    
339  C--   Vertical flux (fVer is at upper face of "u" cell)  C--   Vertical flux (fVer is at upper face of "u" cell)
340    
 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  
341  C     Mean flow component of vertical flux (at k+1) -> aF  C     Mean flow component of vertical flux (at k+1) -> aF
342        IF (momAdvection)        IF (momAdvection)
343       & 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)
344    
345  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
346        IF (momViscosity.AND..NOT.implicitViscosity)        IF (momViscosity.AND..NOT.implicitViscosity)
# Line 318  C     Combine fluxes Line 353  C     Combine fluxes
353         ENDDO         ENDDO
354        ENDDO        ENDDO
355    
 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  
   
356  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
357        DO j=jMin,jMax        DO j=jMin,jMax
358         DO i=iMin,iMax         DO i=iMin,iMax
# Line 343  C--   Tendency is minus divergence of th Line 368  C--   Tendency is minus divergence of th
368       &   +fMer(i,j+1)          - fMer(i  ,j)       &   +fMer(i,j+1)          - fMer(i  ,j)
369       &   +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac       &   +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac
370       &   )       &   )
371       & _PHM( +phxFac * pf(i,j) )       &  - phxFac*dPhiHydX(i,j)
372         ENDDO         ENDDO
373        ENDDO        ENDDO
374    
375    #ifdef NONLIN_FRSURF
376    C-- account for 3.D divergence of the flow in rStar coordinate:
377          IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
378           DO j=jMin,jMax
379            DO i=iMin,iMax
380             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
381         &     - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
382         &       *uVel(i,j,k,bi,bj)
383            ENDDO
384           ENDDO
385          ENDIF
386          IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
387           DO j=jMin,jMax
388            DO i=iMin,iMax
389             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
390         &     - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
391            ENDDO
392           ENDDO
393          ENDIF
394    #endif /* NONLIN_FRSURF */
395    
396  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
397        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
398  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
# Line 367  C-    No-slip BCs impose a drag at botto Line 413  C-    No-slip BCs impose a drag at botto
413         ENDDO         ENDDO
414        ENDIF        ENDIF
415    
416  C--   Forcing term  C--   Forcing term (moved to timestep.F)
417        IF (momForcing)  c     IF (momForcing)
418       &  CALL EXTERNAL_FORCING_U(  c    &  CALL EXTERNAL_FORCING_U(
419       I     iMin,iMax,jMin,jMax,bi,bj,k,  c    I     iMin,iMax,jMin,jMax,bi,bj,k,
420       I     myCurrentTime,myThid)  c    I     myTime,myThid)
421    
422  C--   Metric terms for curvilinear grid systems  C--   Metric terms for curvilinear grid systems
423        IF (useNHMTerms) THEN        IF (useNHMTerms) THEN
# Line 403  C--   Set du/dt on boundaries to zero Line 449  C--   Set du/dt on boundaries to zero
449  C---- Meridional momentum equation starts here  C---- Meridional momentum equation starts here
450    
451  C     Bi-harmonic term del^2 V -> v4F  C     Bi-harmonic term del^2 V -> v4F
452        IF (momViscosity)        IF (momViscosity .AND. viscA4.NE.0. )
453       & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)       & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
454    
455  C---  Calculate mean and eddy fluxes between cells for meridional flow.  C---  Calculate mean and eddy fluxes between cells for meridional flow.
# Line 420  C     Laplacian and bi-harmonic terms -> Line 466  C     Laplacian and bi-harmonic terms ->
466    
467  C     Combine fluxes -> fZon  C     Combine fluxes -> fZon
468        DO j=jMin,jMax        DO j=jMin,jMax
469         DO i=iMin,iMax         DO i=iMin,iMax+1
470          fZon(i,j) = uDvdxFac*aF(i,j) + AhDvdxFac*vF(i,j)          fZon(i,j) = uDvdxFac*aF(i,j) + AhDvdxFac*vF(i,j)
471         ENDDO         ENDDO
472        ENDDO        ENDDO
# Line 444  C     Combine fluxes -> fMer Line 490  C     Combine fluxes -> fMer
490    
491  C--   Vertical flux (fVer is at upper face of "v" cell)  C--   Vertical flux (fVer is at upper face of "v" cell)
492    
 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  
493  C     o Mean flow component of vertical flux  C     o Mean flow component of vertical flux
494        IF (momAdvection)        IF (momAdvection)
495       & 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)
496    
497  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
498        IF (momViscosity.AND..NOT.implicitViscosity)        IF (momViscosity.AND..NOT.implicitViscosity)
# Line 468  C     Combine fluxes -> fVerV Line 505  C     Combine fluxes -> fVerV
505         ENDDO         ENDDO
506        ENDDO        ENDDO
507    
 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  
   
508  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
509        DO j=jMin,jMax        DO j=jMin,jMax
510         DO i=iMin,iMax         DO i=iMin,iMax
# Line 493  C--   Tendency is minus divergence of th Line 520  C--   Tendency is minus divergence of th
520       &   +fMer(i,j  )          - fMer(i,j-1)       &   +fMer(i,j  )          - fMer(i,j-1)
521       &   +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac       &   +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac
522       &   )       &   )
523       & _PHM( +phyFac*pf(i,j) )       &  - phyFac*dPhiHydY(i,j)
524         ENDDO         ENDDO
525        ENDDO        ENDDO
526    
527    #ifdef NONLIN_FRSURF
528    C-- account for 3.D divergence of the flow in rStar coordinate:
529          IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
530           DO j=jMin,jMax
531            DO i=iMin,iMax
532             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
533         &     - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
534         &       *vVel(i,j,k,bi,bj)
535            ENDDO
536           ENDDO
537          ENDIF
538          IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
539           DO j=jMin,jMax
540            DO i=iMin,iMax
541             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
542         &     - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
543            ENDDO
544           ENDDO
545          ENDIF
546    #endif /* NONLIN_FRSURF */
547    
548  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
549        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
550  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
# Line 517  C-    No-slip BCs impose a drag at botto Line 565  C-    No-slip BCs impose a drag at botto
565         ENDDO         ENDDO
566        ENDIF        ENDIF
567    
568  C--   Forcing term  C--   Forcing term (moved to timestep.F)
569        IF (momForcing)  c     IF (momForcing)
570       & CALL EXTERNAL_FORCING_V(  c    & CALL EXTERNAL_FORCING_V(
571       I     iMin,iMax,jMin,jMax,bi,bj,k,  c    I     iMin,iMax,jMin,jMax,bi,bj,k,
572       I     myCurrentTime,myThid)  c    I     myTime,myThid)
573    
574  C--   Metric terms for curvilinear grid systems  C--   Metric terms for curvilinear grid systems
575        IF (useNHMTerms) THEN        IF (useNHMTerms) THEN
# Line 551  C--   Set dv/dt on boundaries to zero Line 599  C--   Set dv/dt on boundaries to zero
599    
600  C--   Coriolis term  C--   Coriolis term
601  C     Note. As coded here, coriolis will not work with "thin walls"  C     Note. As coded here, coriolis will not work with "thin walls"
602  #ifdef INCLUDE_CD_CODE  c     IF (useCDscheme) THEN
603        CALL MOM_CDSCHEME(bi,bj,k,phi_hyd,myThid)  c       CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)
604  #else  c     ELSE
605        CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)        IF (.NOT.useCDscheme) THEN
606        DO j=jMin,jMax          CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
607         DO i=iMin,iMax          DO j=jMin,jMax
608          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)           DO i=iMin,iMax
609         ENDDO            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
610        ENDDO           ENDDO
611        CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)          ENDDO
612        DO j=jMin,jMax          CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
613         DO i=iMin,iMax          DO j=jMin,jMax
614          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)           DO i=iMin,iMax
615         ENDDO            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
616        ENDDO           ENDDO
617  #endif /* INCLUDE_CD_CODE */          ENDDO
618          ENDIF
619    
620        IF (nonHydrostatic.OR.quasiHydrostatic) THEN        IF (nonHydrostatic.OR.quasiHydrostatic) THEN
621         CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)         CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
622         DO j=jMin,jMax         DO j=jMin,jMax
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22