/[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.60 by jmc, Thu Nov 23 00:45:21 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 156  c       mT(i,j)    = 0. Line 176  c       mT(i,j)    = 0.
176          vort3(i,j) = 0.          vort3(i,j) = 0.
177          omega3(i,j)= 0.          omega3(i,j)= 0.
178          KE(i,j)    = 0.          KE(i,j)    = 0.
179    c       hDiv(i,j)  = 0.
180          viscAh_Z(i,j) = 0.          viscAh_Z(i,j) = 0.
181          viscAh_D(i,j) = 0.          viscAh_D(i,j) = 0.
182          viscA4_Z(i,j) = 0.          viscA4_Z(i,j) = 0.
183          viscA4_D(i,j) = 0.          viscA4_D(i,j) = 0.
184    
 #ifdef ALLOW_AUTODIFF_TAMC  
185          strain(i,j)  = 0. _d 0          strain(i,j)  = 0. _d 0
186          tension(i,j) = 0. _d 0          tension(i,j) = 0. _d 0
187    #ifdef ALLOW_AUTODIFF_TAMC
188          hFacZ(i,j)   = 0. _d 0          hFacZ(i,j)   = 0. _d 0
189  #endif  #endif
190         ENDDO         ENDDO
# Line 172  c       mT(i,j)    = 0. Line 193  c       mT(i,j)    = 0.
193  C--   Term by term tracer parmeters  C--   Term by term tracer parmeters
194  C     o U momentum equation  C     o U momentum equation
195        ArDudrFac    = vfFacMom*1.        ArDudrFac    = vfFacMom*1.
 c     mTFacU       = mtFacMom*1.  
196  C     o V momentum equation  C     o V momentum equation
197        ArDvdrFac    = vfFacMom*1.        ArDvdrFac    = vfFacMom*1.
 c     mTFacV       = mtFacMom*1.  
198    
199  C note: using standard stencil (no mask) results in under-estimating  C note: using standard stencil (no mask) results in under-estimating
200  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 232  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFa
232        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
233    
234        IF (momViscosity) THEN        IF (momViscosity) THEN
235  C--    For viscous term, compute horizontal divergence, tension & strain  C--    For viscous term, compute horizontal divergence, tension & strain
236  C      and mask relative vorticity (free-slip case):  C      and mask relative vorticity (free-slip case):
237    
238    #ifdef ALLOW_AUTODIFF_TAMC
239    CADJ STORE vort3(:,:) =
240    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
241    #endif
242    
243         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
244    
245         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 315  C      in terms of vorticity and diverge
315       I                       viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,       I                       viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
316       I                       harmonic,biharmonic,useVariableViscosity,       I                       harmonic,biharmonic,useVariableViscosity,
317       O                       guDiss,gvDiss,       O                       guDiss,gvDiss,
318       &                       myThid)               &                       myThid)
319         ENDIF         ENDIF
320  C--   if (momViscosity) end of block.  C--   if (momViscosity) end of block.
321        ENDIF        ENDIF
# Line 327  C--   Tendency is minus divergence of th Line 351  C--   Tendency is minus divergence of th
351         ENDDO         ENDDO
352        ENDIF        ENDIF
353    
354  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
355        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
356  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
357         CALL MOM_U_SIDEDRAG(         CALL MOM_U_SIDEDRAG(
# Line 392  C--   Tendency is minus divergence of th Line 416  C--   Tendency is minus divergence of th
416         ENDDO         ENDDO
417        ENDIF        ENDIF
418    
419  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
420        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
421  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
422         CALL MOM_V_SIDEDRAG(         CALL MOM_V_SIDEDRAG(
# Line 481  C- jmc: change it to keep the Coriolis t Line 505  C- jmc: change it to keep the Coriolis t
505           gV(i,j,k,bi,bj) = vCf(i,j)           gV(i,j,k,bi,bj) = vCf(i,j)
506          ENDDO          ENDDO
507         ENDDO         ENDDO
   
508         IF ( writeDiag ) THEN         IF ( writeDiag ) THEN
509           IF (snapshot_mdsio) THEN           IF (snapshot_mdsio) THEN
510             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 525  C- jmc: change it to keep the Coriolis t
525           CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)           CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
526         ENDIF         ENDIF
527  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
   
528        ELSE        ELSE
529         DO j=jMin,jMax         DO j=jMin,jMax
530          DO i=iMin,iMax          DO i=iMin,iMax
# Line 634  C--   Bernoulli term Line 656  C--   Bernoulli term
656  C--   end if momAdvection  C--   end if momAdvection
657        ENDIF        ENDIF
658    
659  C--   Metric terms for curvilinear grid systems  C--   3.D Coriolis term (horizontal momentum, Eastward component: -f'*w)
660  c     IF (usingSphericalPolarMTerms) THEN        IF ( use3dCoriolis ) THEN
661  C      o Spherical polar grid metric terms          CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
662  c      CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)          DO j=jMin,jMax
663  c      DO j=jMin,jMax           DO i=iMin,iMax
664  c       DO i=iMin,iMax            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
665  c        gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)           ENDDO
666  c       ENDDO          ENDDO
667  c      ENDDO         IF ( usingCurvilinearGrid ) THEN
668  c      CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)  C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
669  c      DO j=jMin,jMax          CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
670  c       DO i=iMin,iMax          DO j=jMin,jMax
671  c        gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)           DO i=iMin,iMax
672  c       ENDDO            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
673  c      ENDDO           ENDDO
674  c     ENDIF          ENDDO
675           ENDIF
676          ENDIF
677    
678    C--   Non-Hydrostatic (spherical) metric terms
679          IF ( useNHMTerms ) THEN
680           CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
681           DO j=jMin,jMax
682            DO i=iMin,iMax
683             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
684            ENDDO
685           ENDDO
686           CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
687           DO j=jMin,jMax
688            DO i=iMin,iMax
689             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
690            ENDDO
691           ENDDO
692          ENDIF
693    
694  C--   Set du/dt & dv/dt on boundaries to zero  C--   Set du/dt & dv/dt on boundaries to zero
695        DO j=jMin,jMax        DO j=jMin,jMax
Legend:
Removed from v.1.56  
changed lines
  Added in v.1.60
  ViewVC Help
Powered by ViewVC 1.1.22