/[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.6 by adcroft, Tue Nov 5 19:58:21 2002 UTC revision 1.10 by jmc, Tue Feb 11 04:06:15 2003 UTC
# 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        phi_hyd,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 59  C  k                    :: vertical leve Line 59  C  k                    :: vertical leve
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  phi_hyd              :: hydrostatic pressure (perturbation)
62    C  dPhiHydX,Y           :: Gradient (X & Y dir.) of Hydrostatic Potential
63  C  KappaRU              :: vertical viscosity  C  KappaRU              :: vertical viscosity
64  C  KappaRV              :: vertical viscosity  C  KappaRV              :: vertical viscosity
65  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
66  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
67  C  myCurrentTime        :: current time  C  myTime               :: current time
68  C  myIter               :: current time-step number  C  myIter               :: current time-step number
69  C  myThid               :: thread number  C  myThid               :: thread number
70        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
71        INTEGER k,kUp,kDown        INTEGER k,kUp,kDown
72        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73          _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74          _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
76        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
77        _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
78        _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
79        _RL     myCurrentTime        _RL     myTime
80        INTEGER myIter        INTEGER myIter
81        INTEGER myThid        INTEGER myThid
82    
# Line 119  C     uDudxFac, AhDudxFac, etc ... indiv Line 122  C     uDudxFac, AhDudxFac, etc ... indiv
122        _RL  vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL  vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123        _RL  uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL  uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124        _RL  vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL  vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125          _RL  rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126          _RL  rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127  C     I,J,K - Loop counters  C     I,J,K - Loop counters
128  C     rVelMaskOverride - Factor for imposing special surface boundary conditions  C     rVelMaskOverride - Factor for imposing special surface boundary conditions
129  C                        ( set according to free-surface condition ).  C                        ( set according to free-surface condition ).
# Line 174  C     Initialise intermediate terms Line 179  C     Initialise intermediate terms
179          pF(i,j)   = 0.          pF(i,j)   = 0.
180          fZon(i,j) = 0.          fZon(i,j) = 0.
181          fMer(i,j) = 0.          fMer(i,j) = 0.
182            rTransU(i,j) = 0.
183            rTransV(i,j) = 0.
184         ENDDO         ENDDO
185        ENDDO        ENDDO
186    
# Line 250  C     Calculate velocity field "volume t Line 257  C     Calculate velocity field "volume t
257    
258        CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)        CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)
259    
260    C---  First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
261          IF (momAdvection.AND.k.EQ.1) THEN
262    
263    C-    Calculate vertical transports above U & V points (West & South face):
264           CALL MOM_CALC_RTRANS( k, bi, bj,
265         O                       rTransU, rTransV,
266         I                       myTime, myIter, myThid)
267    
268    C-    Free surface correction term (flux at k=1)
269           CALL MOM_U_ADV_WU(bi,bj,k,uVel,wVel,rTransU,af,myThid)
270           DO j=jMin,jMax
271            DO i=iMin,iMax
272             fVerU(i,j,kUp) = af(i,j)
273            ENDDO
274           ENDDO
275    
276           CALL MOM_V_ADV_WV(bi,bj,k,vVel,wVel,rTransV,af,myThid)
277           DO j=jMin,jMax
278            DO i=iMin,iMax
279             fVerV(i,j,kUp) = af(i,j)
280            ENDDO
281           ENDDO
282    
283    C---  endif momAdvection & k=1
284          ENDIF
285    
286    
287    C---  Calculate vertical transports (at k+1) below U & V points :
288          IF (momAdvection) THEN
289           CALL MOM_CALC_RTRANS( k+1, bi, bj,
290         O                       rTransU, rTransV,
291         I                       myTime, myIter, myThid)
292          ENDIF
293    
294    
295  C---- Zonal momentum equation starts here  C---- Zonal momentum equation starts here
296    
297  C     Bi-harmonic term del^2 U -> v4F  C     Bi-harmonic term del^2 U -> v4F
298        IF (momViscosity)        IF (momViscosity .AND. viscA4.NE.0. )
299       & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)       & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
300    
301  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 328  C     Laplacian and bi-harmonic term
328       & CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)       & CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
329    
330  C     Combine fluxes -> fMer  C     Combine fluxes -> fMer
331        DO j=jMin,jMax        DO j=jMin,jMax+1
332         DO i=iMin,iMax         DO i=iMin,iMax
333          fMer(i,j) = vDudyFac*aF(i,j) + AhDudyFac*vF(i,j)          fMer(i,j) = vDudyFac*aF(i,j) + AhDudyFac*vF(i,j)
334         ENDDO         ENDDO
# Line 294  C     Combine fluxes -> fMer Line 336  C     Combine fluxes -> fMer
336    
337  C--   Vertical flux (fVer is at upper face of "u" cell)  C--   Vertical flux (fVer is at upper face of "u" cell)
338    
 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  
339  C     Mean flow component of vertical flux (at k+1) -> aF  C     Mean flow component of vertical flux (at k+1) -> aF
340        IF (momAdvection)        IF (momAdvection)
341       & 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)
342    
343  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
344        IF (momViscosity.AND..NOT.implicitViscosity)        IF (momViscosity.AND..NOT.implicitViscosity)
# Line 318  C     Combine fluxes Line 351  C     Combine fluxes
351         ENDDO         ENDDO
352        ENDDO        ENDDO
353    
 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  
   
354  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
355        DO j=jMin,jMax        DO j=jMin,jMax
356         DO i=iMin,iMax         DO i=iMin,iMax
# Line 343  C--   Tendency is minus divergence of th Line 366  C--   Tendency is minus divergence of th
366       &   +fMer(i,j+1)          - fMer(i  ,j)       &   +fMer(i,j+1)          - fMer(i  ,j)
367       &   +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac       &   +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac
368       &   )       &   )
369       & _PHM( +phxFac * pf(i,j) )       &  - phxFac*dPhiHydX(i,j)
370         ENDDO         ENDDO
371        ENDDO        ENDDO
372    
373    #ifdef NONLIN_FRSURF
374    C-- account for 3.D divergence of the flow in rStar coordinate:
375          IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
376           DO j=jMin,jMax
377            DO i=iMin,iMax
378             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
379         &     - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
380         &       *uVel(i,j,k,bi,bj)
381            ENDDO
382           ENDDO
383          ENDIF
384          IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
385           DO j=jMin,jMax
386            DO i=iMin,iMax
387             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
388         &     - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
389            ENDDO
390           ENDDO
391          ENDIF
392    #endif /* NONLIN_FRSURF */
393    
394  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
395        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
396  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
# Line 371  C--   Forcing term Line 415  C--   Forcing term
415        IF (momForcing)        IF (momForcing)
416       &  CALL EXTERNAL_FORCING_U(       &  CALL EXTERNAL_FORCING_U(
417       I     iMin,iMax,jMin,jMax,bi,bj,k,       I     iMin,iMax,jMin,jMax,bi,bj,k,
418       I     myCurrentTime,myThid)       I     myTime,myThid)
419    
420  C--   Metric terms for curvilinear grid systems  C--   Metric terms for curvilinear grid systems
421        IF (useNHMTerms) THEN        IF (useNHMTerms) THEN
# Line 403  C--   Set du/dt on boundaries to zero Line 447  C--   Set du/dt on boundaries to zero
447  C---- Meridional momentum equation starts here  C---- Meridional momentum equation starts here
448    
449  C     Bi-harmonic term del^2 V -> v4F  C     Bi-harmonic term del^2 V -> v4F
450        IF (momViscosity)        IF (momViscosity .AND. viscA4.NE.0. )
451       & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)       & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
452    
453  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 464  C     Laplacian and bi-harmonic terms ->
464    
465  C     Combine fluxes -> fZon  C     Combine fluxes -> fZon
466        DO j=jMin,jMax        DO j=jMin,jMax
467         DO i=iMin,iMax         DO i=iMin,iMax+1
468          fZon(i,j) = uDvdxFac*aF(i,j) + AhDvdxFac*vF(i,j)          fZon(i,j) = uDvdxFac*aF(i,j) + AhDvdxFac*vF(i,j)
469         ENDDO         ENDDO
470        ENDDO        ENDDO
# Line 444  C     Combine fluxes -> fMer Line 488  C     Combine fluxes -> fMer
488    
489  C--   Vertical flux (fVer is at upper face of "v" cell)  C--   Vertical flux (fVer is at upper face of "v" cell)
490    
 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  
491  C     o Mean flow component of vertical flux  C     o Mean flow component of vertical flux
492        IF (momAdvection)        IF (momAdvection)
493       & 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)
494    
495  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
496        IF (momViscosity.AND..NOT.implicitViscosity)        IF (momViscosity.AND..NOT.implicitViscosity)
# Line 468  C     Combine fluxes -> fVerV Line 503  C     Combine fluxes -> fVerV
503         ENDDO         ENDDO
504        ENDDO        ENDDO
505    
 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  
   
506  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
507        DO j=jMin,jMax        DO j=jMin,jMax
508         DO i=iMin,iMax         DO i=iMin,iMax
# Line 493  C--   Tendency is minus divergence of th Line 518  C--   Tendency is minus divergence of th
518       &   +fMer(i,j  )          - fMer(i,j-1)       &   +fMer(i,j  )          - fMer(i,j-1)
519       &   +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac       &   +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac
520       &   )       &   )
521       & _PHM( +phyFac*pf(i,j) )       &  - phyFac*dPhiHydY(i,j)
522         ENDDO         ENDDO
523        ENDDO        ENDDO
524    
525    #ifdef NONLIN_FRSURF
526    C-- account for 3.D divergence of the flow in rStar coordinate:
527          IF ( momAdvection .AND. select_rStar.GT.0 ) THEN
528           DO j=jMin,jMax
529            DO i=iMin,iMax
530             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
531         &     - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
532         &       *vVel(i,j,k,bi,bj)
533            ENDDO
534           ENDDO
535          ENDIF
536          IF ( momAdvection .AND. select_rStar.LT.0 ) THEN
537           DO j=jMin,jMax
538            DO i=iMin,iMax
539             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
540         &     - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
541            ENDDO
542           ENDDO
543          ENDIF
544    #endif /* NONLIN_FRSURF */
545    
546  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
547        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
548  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
# Line 521  C--   Forcing term Line 567  C--   Forcing term
567        IF (momForcing)        IF (momForcing)
568       & CALL EXTERNAL_FORCING_V(       & CALL EXTERNAL_FORCING_V(
569       I     iMin,iMax,jMin,jMax,bi,bj,k,       I     iMin,iMax,jMin,jMax,bi,bj,k,
570       I     myCurrentTime,myThid)       I     myTime,myThid)
571    
572  C--   Metric terms for curvilinear grid systems  C--   Metric terms for curvilinear grid systems
573        IF (useNHMTerms) THEN        IF (useNHMTerms) THEN
# Line 552  C--   Set dv/dt on boundaries to zero Line 598  C--   Set dv/dt on boundaries to zero
598  C--   Coriolis term  C--   Coriolis term
599  C     Note. As coded here, coriolis will not work with "thin walls"  C     Note. As coded here, coriolis will not work with "thin walls"
600  #ifdef INCLUDE_CD_CODE  #ifdef INCLUDE_CD_CODE
601        CALL MOM_CDSCHEME(bi,bj,k,phi_hyd,myThid)        CALL MOM_CDSCHEME(bi,bj,k,phi_hyd,dPhiHydX,dPhiHydY,myThid)
602  #else  #else
603        CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)        CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
604        DO j=jMin,jMax        DO j=jMin,jMax
# Line 567  C     Note. As coded here, coriolis will Line 613  C     Note. As coded here, coriolis will
613         ENDDO         ENDDO
614        ENDDO        ENDDO
615  #endif /* INCLUDE_CD_CODE */  #endif /* INCLUDE_CD_CODE */
616        IF (nonHydrostatic) THEN        IF (nonHydrostatic.OR.quasiHydrostatic) THEN
617         CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)         CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
618         DO j=jMin,jMax         DO j=jMin,jMax
619          DO i=iMin,iMax          DO i=iMin,iMax
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22