/[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.56 by mlosch, Tue Feb 7 11:46:17 2006 UTC revision 1.59 by heimbach, Tue Jul 18 03:23:30 2006 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "MOM_VECINV_OPTIONS.h"  #include "MOM_VECINV_OPTIONS.h"
5    
6        SUBROUTINE MOM_VECINV(        SUBROUTINE MOM_VECINV(
7       I        bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
8       I        KappaRU, KappaRV,       I        KappaRU, KappaRV,
9       U        fVerU, fVerV,       U        fVerU, fVerV,
# Line 38  C     == Global variables == Line 38  C     == Global variables ==
38  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
39  #include "TIMEAVE_STATV.h"  #include "TIMEAVE_STATV.h"
40  #endif  #endif
41    #ifdef ALLOW_AUTODIFF_TAMC
42    # include "tamc.h"
43    # include "tamc_keys.h"
44    #endif
45    
46  C     == Routine arguments ==  C     == Routine arguments ==
47  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
# Line 71  C     == Local variables == Line 75  C     == Local variables ==
75        _RL      vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76        _RL      uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77        _RL      vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 c     _RL      mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
78        _RS hFacZ   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS hFacZ   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79        _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 94  C     i,j,k  :: Loop counters Line 97  C     i,j,k  :: Loop counters
97        INTEGER i,j,k        INTEGER i,j,k
98  C     xxxFac - On-off tracer parameters used for switching terms off.  C     xxxFac - On-off tracer parameters used for switching terms off.
99        _RL  ArDudrFac        _RL  ArDudrFac
 c     _RL  mtFacU  
100        _RL  ArDvdrFac        _RL  ArDvdrFac
 c     _RL  mtFacV  
101        _RL  sideMaskFac        _RL  sideMaskFac
102        LOGICAL bottomDragTerms        LOGICAL bottomDragTerms
103        LOGICAL writeDiag        LOGICAL writeDiag
104        LOGICAL harmonic,biharmonic,useVariableViscosity        LOGICAL harmonic,biharmonic,useVariableViscosity
105    #ifdef ALLOW_AUTODIFF_TAMC
106          INTEGER imomkey
107    #endif
108    
109  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
110        INTEGER offsets(9)        INTEGER offsets(9)
# Line 116  C--   (at least in part) Line 120  C--   (at least in part)
120        fVerV(1,1,kUp) = fVerV(1,1,kUp)        fVerV(1,1,kUp) = fVerV(1,1,kUp)
121  #endif  #endif
122    
123    #ifdef ALLOW_AUTODIFF_TAMC
124              act0 = k - 1
125              max0 = Nr
126              act1 = bi - myBxLo(myThid)
127              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
128              act2 = bj - myByLo(myThid)
129              max2 = myByHi(myThid) - myByLo(myThid) + 1
130              act3 = myThid - 1
131              max3 = nTx*nTy
132              act4 = ikey_dynamics - 1
133              imomkey = (act0 + 1)
134         &                    + act1*max0
135         &                    + act2*max0*max1
136         &                    + act3*max0*max1*max2
137         &                    + act4*max0*max1*max2*max3
138    #endif /* ALLOW_AUTODIFF_TAMC */
139    
140        writeDiag = DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock)        writeDiag = DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock)
141    
142  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
# Line 146  C     Initialise intermediate terms Line 167  C     Initialise intermediate terms
167          vrF(i,j)   = 0.          vrF(i,j)   = 0.
168          uCf(i,j)   = 0.          uCf(i,j)   = 0.
169          vCf(i,j)   = 0.          vCf(i,j)   = 0.
 c       mT(i,j)    = 0.  
170          del2u(i,j) = 0.          del2u(i,j) = 0.
171          del2v(i,j) = 0.          del2v(i,j) = 0.
172          dStar(i,j) = 0.          dStar(i,j) = 0.
# Line 172  c       mT(i,j)    = 0. Line 192  c       mT(i,j)    = 0.
192  C--   Term by term tracer parmeters  C--   Term by term tracer parmeters
193  C     o U momentum equation  C     o U momentum equation
194        ArDudrFac    = vfFacMom*1.        ArDudrFac    = vfFacMom*1.
 c     mTFacU       = mtFacMom*1.  
195  C     o V momentum equation  C     o V momentum equation
196        ArDvdrFac    = vfFacMom*1.        ArDvdrFac    = vfFacMom*1.
 c     mTFacV       = mtFacMom*1.  
197    
198  C note: using standard stencil (no mask) results in under-estimating  C note: using standard stencil (no mask) results in under-estimating
199  C       vorticity at a no-slip boundary by a factor of 2 = sideDragFactor  C       vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
# Line 213  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFa Line 231  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFa
231        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
232    
233        IF (momViscosity) THEN        IF (momViscosity) THEN
234  C--    For viscous term, compute horizontal divergence, tension & strain  C--    For viscous term, compute horizontal divergence, tension & strain
235  C      and mask relative vorticity (free-slip case):  C      and mask relative vorticity (free-slip case):
236    
237    #ifdef ALLOW_AUTODIFF_TAMC
238    CADJ STORE vort3(:,:) =
239    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
240    #endif
241    
242         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
243    
244         CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)         CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
# Line 291  C      in terms of vorticity and diverge Line 314  C      in terms of vorticity and diverge
314       I                       viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,       I                       viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
315       I                       harmonic,biharmonic,useVariableViscosity,       I                       harmonic,biharmonic,useVariableViscosity,
316       O                       guDiss,gvDiss,       O                       guDiss,gvDiss,
317       &                       myThid)               &                       myThid)
318         ENDIF         ENDIF
319  C--   if (momViscosity) end of block.  C--   if (momViscosity) end of block.
320        ENDIF        ENDIF
# Line 327  C--   Tendency is minus divergence of th Line 350  C--   Tendency is minus divergence of th
350         ENDDO         ENDDO
351        ENDIF        ENDIF
352    
353  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
354        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
355  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
356         CALL MOM_U_SIDEDRAG(         CALL MOM_U_SIDEDRAG(
# Line 392  C--   Tendency is minus divergence of th Line 415  C--   Tendency is minus divergence of th
415         ENDDO         ENDDO
416        ENDIF        ENDIF
417    
418  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
419        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
420  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
421         CALL MOM_V_SIDEDRAG(         CALL MOM_V_SIDEDRAG(
# Line 481  C- jmc: change it to keep the Coriolis t Line 504  C- jmc: change it to keep the Coriolis t
504           gV(i,j,k,bi,bj) = vCf(i,j)           gV(i,j,k,bi,bj) = vCf(i,j)
505          ENDDO          ENDDO
506         ENDDO         ENDDO
   
507         IF ( writeDiag ) THEN         IF ( writeDiag ) THEN
508           IF (snapshot_mdsio) THEN           IF (snapshot_mdsio) THEN
509             CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)             CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)
# Line 502  C- jmc: change it to keep the Coriolis t Line 524  C- jmc: change it to keep the Coriolis t
524           CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)           CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
525         ENDIF         ENDIF
526  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
   
527        ELSE        ELSE
528         DO j=jMin,jMax         DO j=jMin,jMax
529          DO i=iMin,iMax          DO i=iMin,iMax
# Line 634  C--   Bernoulli term Line 655  C--   Bernoulli term
655  C--   end if momAdvection  C--   end if momAdvection
656        ENDIF        ENDIF
657    
658  C--   Metric terms for curvilinear grid systems  C--   3.D Coriolis term (horizontal momentum, Eastward component: -f'*w)
659  c     IF (usingSphericalPolarMTerms) THEN        IF ( use3dCoriolis ) THEN
660  C      o Spherical polar grid metric terms          CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
661  c      CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)          DO j=jMin,jMax
662  c      DO j=jMin,jMax           DO i=iMin,iMax
663  c       DO i=iMin,iMax            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
664  c        gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)           ENDDO
665  c       ENDDO          ENDDO
666  c      ENDDO         IF ( usingCurvilinearGrid ) THEN
667  c      CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)  C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
668  c      DO j=jMin,jMax          CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
669  c       DO i=iMin,iMax          DO j=jMin,jMax
670  c        gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)           DO i=iMin,iMax
671  c       ENDDO            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
672  c      ENDDO           ENDDO
673  c     ENDIF          ENDDO
674           ENDIF
675          ENDIF
676    
677    C--   Non-Hydrostatic (spherical) metric terms
678          IF ( useNHMTerms ) THEN
679           CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
680           DO j=jMin,jMax
681            DO i=iMin,iMax
682             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
683            ENDDO
684           ENDDO
685           CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
686           DO j=jMin,jMax
687            DO i=iMin,iMax
688             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
689            ENDDO
690           ENDDO
691          ENDIF
692    
693  C--   Set du/dt & dv/dt on boundaries to zero  C--   Set du/dt & dv/dt on boundaries to zero
694        DO j=jMin,jMax        DO j=jMin,jMax
Legend:
Removed from v.1.56  
changed lines
  Added in v.1.59
  ViewVC Help
Powered by ViewVC 1.1.22