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

Diff of /MITgcm/pkg/mom_vecinv/mom_vecinv.F

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

revision 1.72 by jmc, Tue Feb 11 20:24:06 2014 UTC revision 1.77 by jmc, Thu Sep 10 18:08:51 2015 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "MOM_VECINV_OPTIONS.h"  #include "MOM_VECINV_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF
6    # include "AUTODIFF_OPTIONS.h"
7    #endif
8  #ifdef ALLOW_MOM_COMMON  #ifdef ALLOW_MOM_COMMON
9  # include "MOM_COMMON_OPTIONS.h"  # include "MOM_COMMON_OPTIONS.h"
10  #endif  #endif
11    
12        SUBROUTINE MOM_VECINV(        SUBROUTINE MOM_VECINV(
13       I        bi,bj,k,iMin,iMax,jMin,jMax,       I        bi,bj,k,iMin,iMax,jMin,jMax,
14       I        KappaRU, KappaRV,       I        kappaRU, kappaRV,
15       I        fVerUkm, fVerVkm,       I        fVerUkm, fVerVkm,
16       O        fVerUkp, fVerVkp,       O        fVerUkp, fVerVkp,
17       O        guDiss, gvDiss,       O        guDiss, gvDiss,
# Line 56  C     bi,bj   :: current tile indices Line 59  C     bi,bj   :: current tile indices
59  C     k       :: current vertical level  C     k       :: current vertical level
60  C     iMin,iMax,jMin,jMax :: loop ranges  C     iMin,iMax,jMin,jMax :: loop ranges
61  C     fVerU   :: Flux of momentum in the vertical direction, out of the upper  C     fVerU   :: Flux of momentum in the vertical direction, out of the upper
62  C     fVerV   :: face of a cell K ( flux into the cell above ).  C     fVerV   :: face of a cell k ( flux into the cell above ).
63  C     fVerUkm :: vertical viscous flux of U, interface above (k-1/2)  C     fVerUkm :: vertical viscous flux of U, interface above (k-1/2)
64  C     fVerVkm :: vertical viscous flux of V, interface above (k-1/2)  C     fVerVkm :: vertical viscous flux of V, interface above (k-1/2)
65  C     fVerUkp :: vertical viscous flux of U, interface below (k+1/2)  C     fVerUkp :: vertical viscous flux of U, interface below (k+1/2)
# Line 69  C     myIter  :: current time-step numbe Line 72  C     myIter  :: current time-step numbe
72  C     myThid  :: my Thread Id number  C     myThid  :: my Thread Id number
73        INTEGER bi,bj,k        INTEGER bi,bj,k
74        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
75        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
76        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
77        _RL fVerUkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerUkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RL fVerVkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerVkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79        _RL fVerUkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerUkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 132  C     xxxFac :: On-off tracer parameters Line 135  C     xxxFac :: On-off tracer parameters
135        CHARACTER*(1) pf        CHARACTER*(1) pf
136  #endif  #endif
137    
138  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
139  C--   only the kDown part of fverU/V is set in this subroutine  C--   only the kDown part of fverU/V is set in this subroutine
140  C--   the kUp is still required  C--   the kUp is still required
141  C--   In the case of mom_fluxform Kup is set as well  C--   In the case of mom_fluxform kUp is set as well
142  C--   (at least in part)  C--   (at least in part)
143        fVerUkm(1,1) = fVerUkm(1,1)        fVerUkm(1,1) = fVerUkm(1,1)
144        fVerVkm(1,1) = fVerVkm(1,1)        fVerVkm(1,1) = fVerVkm(1,1)
# Line 206  c       viscA4_D(i,j) = 0. Line 209  c       viscA4_D(i,j) = 0.
209          strain(i,j)  = 0. _d 0          strain(i,j)  = 0. _d 0
210          strainBC(i,j)= 0. _d 0          strainBC(i,j)= 0. _d 0
211          tension(i,j) = 0. _d 0          tension(i,j) = 0. _d 0
212  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
213          hFacZ(i,j)   = 0. _d 0          hFacZ(i,j)   = 0. _d 0
214  #endif  #endif
215         ENDDO         ENDDO
# Line 227  C       vorticity at a no-slip boundary Line 230  C       vorticity at a no-slip boundary
230        ENDIF        ENDIF
231    
232        IF (     no_slip_bottom        IF (     no_slip_bottom
233       &    .OR. bottomDragQuadratic.NE.0.       &    .OR. selectBotDragQuadr.GE.0
234       &    .OR. bottomDragLinear.NE.0.) THEN       &    .OR. bottomDragLinear.NE.0.) THEN
235         bottomDragTerms=.TRUE.         bottomDragTerms=.TRUE.
236        ELSE        ELSE
# Line 355  C---  Other dissipation terms in Zonal m Line 358  C---  Other dissipation terms in Zonal m
358  C--   Vertical flux (fVer is at upper face of "u" cell)  C--   Vertical flux (fVer is at upper face of "u" cell)
359  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
360        IF ( .NOT.implicitViscosity ) THEN        IF ( .NOT.implicitViscosity ) THEN
361         CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,vrF,myThid)         CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,kappaRU,vrF,myThid)
362  C     Combine fluxes  C     Combine fluxes
363         DO j=jMin,jMax         DO j=jMin,jMax
364          DO i=iMin,iMax          DO i=iMin,iMax
# Line 363  C     Combine fluxes Line 366  C     Combine fluxes
366          ENDDO          ENDDO
367         ENDDO         ENDDO
368  C--   Tendency is minus divergence of the fluxes  C--   Tendency is minus divergence of the fluxes
369    C     vert.visc.flx is scaled by deepFac2F (deep-atmos) and rhoFac (anelastic)
370         DO j=jMin,jMax         DO j=jMin,jMax
371          DO i=iMin,iMax          DO i=iMin,iMax
372           guDiss(i,j) = guDiss(i,j)           guDiss(i,j) = guDiss(i,j)
373       &   -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
374       &   *recip_rAw(i,j,bi,bj)       &   *recip_rAw(i,j,bi,bj)
375       &   *( fVerUkp(i,j) - fVerUkm(i,j) )*rkSign       &   *( fVerUkp(i,j) - fVerUkm(i,j) )*rkSign
376         &   *recip_deepFac2C(k)*recip_rhoFacC(k)
377          ENDDO          ENDDO
378         ENDDO         ENDDO
379        ENDIF        ENDIF
# Line 391  C-     No-slip BCs impose a drag at wall Line 396  C-     No-slip BCs impose a drag at wall
396    
397  C-    No-slip BCs impose a drag at bottom  C-    No-slip BCs impose a drag at bottom
398        IF ( bottomDragTerms ) THEN        IF ( bottomDragTerms ) THEN
399         CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)         CALL MOM_U_BOTTOMDRAG( bi, bj, k,
400         I            uFld, vFld, KE, kappaRU,
401         O            vF,
402         I            myThid )
403         DO j=jMin,jMax         DO j=jMin,jMax
404          DO i=iMin,iMax          DO i=iMin,iMax
405           guDiss(i,j) = guDiss(i,j)+vF(i,j)           guDiss(i,j) = guDiss(i,j)+vF(i,j)
# Line 399  C-    No-slip BCs impose a drag at botto Line 407  C-    No-slip BCs impose a drag at botto
407         ENDDO         ENDDO
408        ENDIF        ENDIF
409  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
410        IF ( useShelfIce.AND.bottomDragTerms ) THEN        IF ( useShelfIce ) THEN
411         CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)         CALL SHELFICE_U_DRAG( bi, bj, k,
412         I            uFld, vFld, KE, kappaRU,
413         O            vF,
414         I            myThid )
415         DO j=jMin,jMax         DO j=jMin,jMax
416          DO i=iMin,iMax          DO i=iMin,iMax
417           guDiss(i,j) = guDiss(i,j) + vF(i,j)           guDiss(i,j) = guDiss(i,j) + vF(i,j)
# Line 414  C---  Other dissipation terms in Meridio Line 425  C---  Other dissipation terms in Meridio
425  C--   Vertical flux (fVer is at upper face of "v" cell)  C--   Vertical flux (fVer is at upper face of "v" cell)
426  C     Eddy component of vertical flux (interior component only) -> vrF  C     Eddy component of vertical flux (interior component only) -> vrF
427        IF ( .NOT.implicitViscosity ) THEN        IF ( .NOT.implicitViscosity ) THEN
428         CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,vrF,myThid)         CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,kappaRV,vrF,myThid)
429  C     Combine fluxes -> fVerV  C     Combine fluxes -> fVerV
430         DO j=jMin,jMax         DO j=jMin,jMax
431          DO i=iMin,iMax          DO i=iMin,iMax
# Line 422  C     Combine fluxes -> fVerV Line 433  C     Combine fluxes -> fVerV
433          ENDDO          ENDDO
434         ENDDO         ENDDO
435  C--   Tendency is minus divergence of the fluxes  C--   Tendency is minus divergence of the fluxes
436    C     vert.visc.flx is scaled by deepFac2F (deep-atmos) and rhoFac (anelastic)
437         DO j=jMin,jMax         DO j=jMin,jMax
438          DO i=iMin,iMax          DO i=iMin,iMax
439           gvDiss(i,j) = gvDiss(i,j)           gvDiss(i,j) = gvDiss(i,j)
440       &   -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
441       &   *recip_rAs(i,j,bi,bj)       &   *recip_rAs(i,j,bi,bj)
442       &   *( fVerVkp(i,j) - fVerVkm(i,j) )*rkSign       &   *( fVerVkp(i,j) - fVerVkm(i,j) )*rkSign
443         &   *recip_deepFac2C(k)*recip_rhoFacC(k)
444          ENDDO          ENDDO
445         ENDDO         ENDDO
446        ENDIF        ENDIF
# Line 450  C-     No-slip BCs impose a drag at wall Line 463  C-     No-slip BCs impose a drag at wall
463    
464  C-    No-slip BCs impose a drag at bottom  C-    No-slip BCs impose a drag at bottom
465        IF ( bottomDragTerms ) THEN        IF ( bottomDragTerms ) THEN
466         CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)         CALL MOM_V_BOTTOMDRAG( bi, bj, k,
467         I            uFld, vFld, KE, kappaRV,
468         O            vF,
469         I            myThid )
470         DO j=jMin,jMax         DO j=jMin,jMax
471          DO i=iMin,iMax          DO i=iMin,iMax
472           gvDiss(i,j) = gvDiss(i,j)+vF(i,j)           gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
# Line 458  C-    No-slip BCs impose a drag at botto Line 474  C-    No-slip BCs impose a drag at botto
474         ENDDO         ENDDO
475        ENDIF        ENDIF
476  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
477        IF  (useShelfIce.AND.bottomDragTerms ) THEN        IF ( useShelfIce ) THEN
478           CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)         CALL SHELFICE_V_DRAG( bi, bj, k,
479           DO j=jMin,jMax       I            uFld, vFld, KE, kappaRV,
480            DO i=iMin,iMax       O            vF,
481             gvDiss(i,j) = gvDiss(i,j) + vF(i,j)       I            myThid )
482            ENDDO         DO j=jMin,jMax
483           ENDDO          DO i=iMin,iMax
484          ENDIF           gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
485            ENDDO
486           ENDDO
487          ENDIF
488  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
489    
490  C--   if (momViscosity) end of block.  C--   if (momViscosity) end of block.
# Line 489  C- jmc: change it to keep the Coriolis t Line 508  C- jmc: change it to keep the Coriolis t
508       &     .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )       &     .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )
509       &   ) THEN       &   ) THEN
510         IF (useAbsVorticity) THEN         IF (useAbsVorticity) THEN
511          CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,          CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,omega3,hFacZ,r_hFacZ,
512       &                         uCf,myThid)       &                         uCf,myThid)
513          CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,          CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,omega3,hFacZ,r_hFacZ,
514       &                         vCf,myThid)       &                         vCf,myThid)
515         ELSE         ELSE
516          CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,          CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,
# Line 542  C--   Horizontal advection of relative ( Line 561  C--   Horizontal advection of relative (
561          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,vort3, r_hFacZ,          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,vort3, r_hFacZ,
562       &                         uCf,myThid)       &                         uCf,myThid)
563         ELSEIF ( useAbsVorticity ) THEN         ELSEIF ( useAbsVorticity ) THEN
564          CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,          CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,omega3,hFacZ,r_hFacZ,
565       &                         uCf,myThid)       &                         uCf,myThid)
566         ELSE         ELSE
567          CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,vort3, hFacZ,r_hFacZ,          CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,vort3, hFacZ,r_hFacZ,
# Line 555  C--   Horizontal advection of relative ( Line 574  C--   Horizontal advection of relative (
574         ENDDO         ENDDO
575         IF ( (highOrderVorticity.OR.upwindVorticity)         IF ( (highOrderVorticity.OR.upwindVorticity)
576       &     .AND.useAbsVorticity ) THEN       &     .AND.useAbsVorticity ) THEN
577          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,omega3,r_hFacZ,          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,uFld,omega3,r_hFacZ,
578       &                         vCf,myThid)       &                         vCf,myThid)
579         ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN         ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
580          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3, r_hFacZ,          CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,uFld,vort3, r_hFacZ,
581       &                         vCf,myThid)       &                         vCf,myThid)
582         ELSEIF ( useAbsVorticity ) THEN         ELSEIF ( useAbsVorticity ) THEN
583          CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,          CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,omega3,hFacZ,r_hFacZ,
584       &                         vCf,myThid)       &                         vCf,myThid)
585         ELSE         ELSE
586          CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,vort3, hFacZ,r_hFacZ,          CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,vort3, hFacZ,r_hFacZ,
# Line 605  C--   Horizontal advection of relative ( Line 624  C--   Horizontal advection of relative (
624    
625  C--   Vertical shear terms (-w*du/dr & -w*dv/dr)  C--   Vertical shear terms (-w*du/dr & -w*dv/dr)
626         IF ( .NOT. momImplVertAdv ) THEN         IF ( .NOT. momImplVertAdv ) THEN
627          CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)          CALL MOM_VI_U_VERTSHEAR(bi,bj,k,uVel,wVel,uCf,myThid)
628          DO j=jMin,jMax          DO j=jMin,jMax
629           DO i=iMin,iMax           DO i=iMin,iMax
630            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
631           ENDDO           ENDDO
632          ENDDO          ENDDO
633          CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)          CALL MOM_VI_V_VERTSHEAR(bi,bj,k,vVel,wVel,vCf,myThid)
634          DO j=jMin,jMax          DO j=jMin,jMax
635           DO i=iMin,iMax           DO i=iMin,iMax
636            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
# Line 626  C--   Vertical shear terms (-w*du/dr & - Line 645  C--   Vertical shear terms (-w*du/dr & -
645         ENDIF         ENDIF
646    
647  C--   Bernoulli term  C--   Bernoulli term
648         CALL MOM_VI_U_GRAD_KE(bi,bj,K,KE,uCf,myThid)         CALL MOM_VI_U_GRAD_KE(bi,bj,k,KE,uCf,myThid)
649         DO j=jMin,jMax         DO j=jMin,jMax
650          DO i=iMin,iMax          DO i=iMin,iMax
651           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
652          ENDDO          ENDDO
653         ENDDO         ENDDO
654         CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)         CALL MOM_VI_V_GRAD_KE(bi,bj,k,KE,vCf,myThid)
655         DO j=jMin,jMax         DO j=jMin,jMax
656          DO i=iMin,iMax          DO i=iMin,iMax
657           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
# Line 759  C--   Set du/dt & dv/dt on boundaries to Line 778  C--   Set du/dt & dv/dt on boundaries to
778          CALL DIAGNOSTICS_FILL(KE,     'momKE   ',k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(KE,     'momKE   ',k,1,2,bi,bj,myThid)
779         IF (momViscosity) THEN         IF (momViscosity) THEN
780          CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)
         CALL DIAGNOSTICS_FILL(guDiss, 'Um_Diss ',k,1,2,bi,bj,myThid)  
         CALL DIAGNOSTICS_FILL(gvDiss, 'Vm_Diss ',k,1,2,bi,bj,myThid)  
781         ENDIF         ENDIF
782         IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN         IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN
783          CALL DIAGNOSTICS_FILL(tension, 'Tension ',k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(tension, 'Tension ',k,1,2,bi,bj,myThid)
Legend:
Removed from v.1.72  
changed lines
  Added in v.1.77
  ViewVC Help
Powered by ViewVC 1.1.22