/[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.55 by heimbach, Thu Dec 8 15:44:34 2005 UTC revision 1.62 by jmc, Tue Apr 1 01:27:33 2008 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 135  C--   (at least in part) Line 156  C--   (at least in part)
156            offsets(i) = 0            offsets(i) = 0
157          ENDDO          ENDDO
158          offsets(3) = k          offsets(3) = k
159  C       write(*,*) 'offsets = ',(offsets(i),i=1,9)  c       write(*,*) 'offsets = ',(offsets(i),i=1,9)
160        ENDIF        ENDIF
161  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
162    
163  C     Initialise intermediate terms  C--   Initialise intermediate terms
164        DO J=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
165         DO I=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
166          vF(i,j)    = 0.          vF(i,j)    = 0.
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-    need to initialise hDiv for MOM_VI_DEL2UV(call FILL_CS_CORNER_TR_RL)
180            hDiv(i,j)  = 0.
181          viscAh_Z(i,j) = 0.          viscAh_Z(i,j) = 0.
182          viscAh_D(i,j) = 0.          viscAh_D(i,j) = 0.
183          viscA4_Z(i,j) = 0.          viscA4_Z(i,j) = 0.
184          viscA4_D(i,j) = 0.          viscA4_D(i,j) = 0.
185    
 #ifdef ALLOW_AUTODIFF_TAMC  
186          strain(i,j)  = 0. _d 0          strain(i,j)  = 0. _d 0
187          tension(i,j) = 0. _d 0          tension(i,j) = 0. _d 0
188    #ifdef ALLOW_AUTODIFF_TAMC
189          hFacZ(i,j)   = 0. _d 0          hFacZ(i,j)   = 0. _d 0
190  #endif  #endif
191         ENDDO         ENDDO
# Line 172  c       mT(i,j)    = 0. Line 194  c       mT(i,j)    = 0.
194  C--   Term by term tracer parmeters  C--   Term by term tracer parmeters
195  C     o U momentum equation  C     o U momentum equation
196        ArDudrFac    = vfFacMom*1.        ArDudrFac    = vfFacMom*1.
 c     mTFacU       = mtFacMom*1.  
197  C     o V momentum equation  C     o V momentum equation
198        ArDvdrFac    = vfFacMom*1.        ArDvdrFac    = vfFacMom*1.
 c     mTFacV       = mtFacMom*1.  
199    
200  C note: using standard stencil (no mask) results in under-estimating  C note: using standard stencil (no mask) results in under-estimating
201  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 233  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFa
233        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
234    
235        IF (momViscosity) THEN        IF (momViscosity) THEN
236  C--    For viscous term, compute horizontal divergence, tension & strain  C--    For viscous term, compute horizontal divergence, tension & strain
237  C      and mask relative vorticity (free-slip case):  C      and mask relative vorticity (free-slip case):
238    
239    #ifdef ALLOW_AUTODIFF_TAMC
240    CADJ STORE vort3(:,:) =
241    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
242    #endif
243    
244         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
245    
246         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 316  C      in terms of vorticity and diverge
316       I                       viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,       I                       viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
317       I                       harmonic,biharmonic,useVariableViscosity,       I                       harmonic,biharmonic,useVariableViscosity,
318       O                       guDiss,gvDiss,       O                       guDiss,gvDiss,
319       &                       myThid)               &                       myThid)
320         ENDIF         ENDIF
321  C--   if (momViscosity) end of block.  C--   if (momViscosity) end of block.
322        ENDIF        ENDIF
# Line 327  C--   Tendency is minus divergence of th Line 352  C--   Tendency is minus divergence of th
352         ENDDO         ENDDO
353        ENDIF        ENDIF
354    
355  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
356        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
357  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
358         CALL MOM_U_SIDEDRAG(         CALL MOM_U_SIDEDRAG(
# Line 352  C-    No-slip BCs impose a drag at botto Line 377  C-    No-slip BCs impose a drag at botto
377          ENDDO          ENDDO
378         ENDDO         ENDDO
379        ENDIF        ENDIF
380    #ifdef ALLOW_SHELFICE
381          IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN
382           CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
383           DO j=jMin,jMax
384            DO i=iMin,iMax
385             guDiss(i,j) = guDiss(i,j) + vF(i,j)
386            ENDDO
387           ENDDO
388          ENDIF
389    #endif /* ALLOW_SHELFICE */
390    
391    
392  C---  Other dissipation terms in Meridional momentum equation  C---  Other dissipation terms in Meridional momentum equation
393    
# Line 381  C--   Tendency is minus divergence of th Line 417  C--   Tendency is minus divergence of th
417         ENDDO         ENDDO
418        ENDIF        ENDIF
419    
420  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
421        IF (momViscosity.AND.no_slip_sides) THEN        IF (momViscosity.AND.no_slip_sides) THEN
422  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
423         CALL MOM_V_SIDEDRAG(         CALL MOM_V_SIDEDRAG(
# Line 406  C-    No-slip BCs impose a drag at botto Line 442  C-    No-slip BCs impose a drag at botto
442          ENDDO          ENDDO
443         ENDDO         ENDDO
444        ENDIF        ENDIF
445    #ifdef ALLOW_SHELFICE
446          IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN
447             CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid)
448             DO j=jMin,jMax
449              DO i=iMin,iMax
450               gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
451              ENDDO
452             ENDDO
453            ENDIF
454    #endif /* ALLOW_SHELFICE */
455    
456    
457  C-    Vorticity diagnostics:  C-    Vorticity diagnostics:
458        IF ( writeDiag ) THEN        IF ( writeDiag ) THEN
# Line 459  C- jmc: change it to keep the Coriolis t Line 506  C- jmc: change it to keep the Coriolis t
506           gV(i,j,k,bi,bj) = vCf(i,j)           gV(i,j,k,bi,bj) = vCf(i,j)
507          ENDDO          ENDDO
508         ENDDO         ENDDO
   
509         IF ( writeDiag ) THEN         IF ( writeDiag ) THEN
510           IF (snapshot_mdsio) THEN           IF (snapshot_mdsio) THEN
511             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 480  C- jmc: change it to keep the Coriolis t Line 526  C- jmc: change it to keep the Coriolis t
526           CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)           CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
527         ENDIF         ENDIF
528  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
   
529        ELSE        ELSE
530         DO j=jMin,jMax         DO j=jMin,jMax
531          DO i=iMin,iMax          DO i=iMin,iMax
# Line 492  C- jmc: change it to keep the Coriolis t Line 537  C- jmc: change it to keep the Coriolis t
537    
538        IF (momAdvection) THEN        IF (momAdvection) THEN
539  C--   Horizontal advection of relative (or absolute) vorticity  C--   Horizontal advection of relative (or absolute) vorticity
540         IF (highOrderVorticity.AND.useAbsVorticity) THEN         IF ( (highOrderVorticity.OR.upwindVorticity)
541         &     .AND.useAbsVorticity ) THEN
542          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,omega3,r_hFacZ,          CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,omega3,r_hFacZ,
543       &                         uCf,myThid)       &                         uCf,myThid)
544         ELSEIF (highOrderVorticity) THEN         ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
545          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,
546       &                         uCf,myThid)       &                         uCf,myThid)
547         ELSEIF (useAbsVorticity) THEN         ELSEIF ( useAbsVorticity ) THEN
548          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,
549       &                         uCf,myThid)       &                         uCf,myThid)
550         ELSE         ELSE
# Line 510  C--   Horizontal advection of relative ( Line 556  C--   Horizontal advection of relative (
556           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)
557          ENDDO          ENDDO
558         ENDDO         ENDDO
559         IF (highOrderVorticity.AND.useAbsVorticity) THEN         IF ( (highOrderVorticity.OR.upwindVorticity)
560         &     .AND.useAbsVorticity ) THEN
561          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,
562       &                         vCf,myThid)       &                         vCf,myThid)
563         ELSEIF (highOrderVorticity) THEN         ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
564          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,
565       &                         vCf,myThid)       &                         vCf,myThid)
566         ELSEIF (useAbsVorticity) THEN         ELSEIF ( useAbsVorticity ) THEN
567          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,
568       &                         vCf,myThid)       &                         vCf,myThid)
569         ELSE         ELSE
# Line 612  C--   Bernoulli term Line 659  C--   Bernoulli term
659  C--   end if momAdvection  C--   end if momAdvection
660        ENDIF        ENDIF
661    
662  C--   Metric terms for curvilinear grid systems  C--   3.D Coriolis term (horizontal momentum, Eastward component: -f'*w)
663  c     IF (usingSphericalPolarMTerms) THEN        IF ( use3dCoriolis ) THEN
664  C      o Spherical polar grid metric terms          CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
665  c      CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)          DO j=jMin,jMax
666  c      DO j=jMin,jMax           DO i=iMin,iMax
667  c       DO i=iMin,iMax            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
668  c        gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)           ENDDO
669  c       ENDDO          ENDDO
670  c      ENDDO         IF ( usingCurvilinearGrid ) THEN
671  c      CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)  C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
672  c      DO j=jMin,jMax          CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
673  c       DO i=iMin,iMax          DO j=jMin,jMax
674  c        gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)           DO i=iMin,iMax
675  c       ENDDO            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
676  c      ENDDO           ENDDO
677  c     ENDIF          ENDDO
678           ENDIF
679          ENDIF
680    
681    C--   Non-Hydrostatic (spherical) metric terms
682          IF ( useNHMTerms ) THEN
683           CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
684           DO j=jMin,jMax
685            DO i=iMin,iMax
686             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
687            ENDDO
688           ENDDO
689           CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
690           DO j=jMin,jMax
691            DO i=iMin,iMax
692             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
693            ENDDO
694           ENDDO
695          ENDIF
696    
697  C--   Set du/dt & dv/dt on boundaries to zero  C--   Set du/dt & dv/dt on boundaries to zero
698        DO j=jMin,jMax        DO j=jMin,jMax
Legend:
Removed from v.1.55  
changed lines
  Added in v.1.62
  ViewVC Help
Powered by ViewVC 1.1.22