/[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.71 by jmc, Sun Feb 9 18:46:46 2014 UTC revision 1.72 by jmc, Tue Feb 11 20:24:06 2014 UTC
# Line 88  C     == Functions == Line 88  C     == Functions ==
88        EXTERNAL DIFFERENT_MULTIPLE        EXTERNAL DIFFERENT_MULTIPLE
89    
90  C     == Local variables ==  C     == Local variables ==
91    C     strainBC :: same as strain but account for no-slip BC
92    C     vort3BC  :: same as vort3  but account for no-slip BC
93        _RL      vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL      vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL      uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL      uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 103  C     == Local variables == Line 105  C     == Local variables ==
105        _RL zStar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL zStar   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _RL tension (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tension (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RL strain  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL strain  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108          _RL strainBC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109        _RL KE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KE      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _RL omega3  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL omega3  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111        _RL vort3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vort3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112          _RL vort3BC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL hDiv    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL hDiv    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 200  c       viscAh_D(i,j) = 0. Line 204  c       viscAh_D(i,j) = 0.
204  c       viscA4_Z(i,j) = 0.  c       viscA4_Z(i,j) = 0.
205  c       viscA4_D(i,j) = 0.  c       viscA4_D(i,j) = 0.
206          strain(i,j)  = 0. _d 0          strain(i,j)  = 0. _d 0
207            strainBC(i,j)= 0. _d 0
208          tension(i,j) = 0. _d 0          tension(i,j) = 0. _d 0
209  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
210          hFacZ(i,j)   = 0. _d 0          hFacZ(i,j)   = 0. _d 0
# Line 248  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFa Line 253  c     CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFa
253    
254        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
255    
256    C-    mask vort3 and account for no-slip / free-slip BC in vort3BC:
257          DO j=1-OLy,sNy+OLy
258           DO i=1-OLx,sNx+OLx
259             vort3BC(i,j) = vort3(i,j)
260             IF ( hFacZ(i,j).EQ.zeroRS ) THEN
261               vort3BC(i,j) = sideMaskFac*vort3BC(i,j)
262               vort3(i,j)   = 0.
263             ENDIF
264           ENDDO
265          ENDDO
266    
267        IF (momViscosity) THEN        IF (momViscosity) THEN
268  C--    For viscous term, compute horizontal divergence, tension & strain  C--    For viscous term, compute horizontal divergence, tension & strain
269  C      and mask relative vorticity (free-slip case):  C      and mask relative vorticity (free-slip case):
# Line 269  C      and mask relative vorticity (free Line 285  C      and mask relative vorticity (free
285         ENDIF         ENDIF
286  #endif /* NONLIN_FRSURF */  #endif /* NONLIN_FRSURF */
287    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE vort3(:,:) =  
 CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte  
 #endif  
   
288         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)         CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
289    
290         CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)         IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN
291            CALL MOM_CALC_TENSION( bi,bj,k,uFld,vFld,tension,myThid )
292         CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)          CALL MOM_CALC_STRAIN( bi,bj,k,uFld,vFld,hFacZ,strain,myThid )
293    C-    mask strain and account for no-slip / free-slip BC in strainBC:
294  C-     account for no-slip / free-slip BC:          DO j=1-OLy,sNy+OLy
295         DO j=1-OLy,sNy+OLy           DO i=1-OLx,sNx+OLx
296          DO i=1-OLx,sNx+OLx             strainBC(i,j) = strain(i,j)
297            IF ( hFacZ(i,j).EQ.0. ) THEN             IF ( hFacZ(i,j).EQ.zeroRS ) THEN
298              vort3(i,j)  = sideMaskFac*vort3(i,j)               strainBC(i,j) = sideMaskFac*strainBC(i,j)
299              strain(i,j) = sideMaskFac*strain(i,j)               strain(i,j)   = 0.
300            ENDIF             ENDIF
301             ENDDO
302          ENDDO          ENDDO
303         ENDDO         ENDIF
304    
305  C--    Calculate Lateral Viscosities  C--    Calculate Lateral Viscosities
306         DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
# Line 300  C--    Calculate Lateral Viscosities Line 312  C--    Calculate Lateral Viscosities
312          ENDDO          ENDDO
313         ENDDO         ENDDO
314         IF ( useVariableVisc ) THEN         IF ( useVariableVisc ) THEN
315    C-     uses vort3BC & strainBC which account for no-slip / free-slip BC
316           CALL MOM_CALC_VISC( bi, bj, k,           CALL MOM_CALC_VISC( bi, bj, k,
317       O            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,       O            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
318       I            hDiv, vort3, tension, strain, KE, hfacZ,       I            hDiv, vort3BC, tension, strainBC, KE, hfacZ,
319       I            myThid )       I            myThid )
320         ENDIF         ENDIF
321    
# Line 310  C      Calculate del^2 u and del^2 v for Line 323  C      Calculate del^2 u and del^2 v for
323         IF (useBiharmonicVisc) THEN         IF (useBiharmonicVisc) THEN
324           CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,           CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,
325       O                      del2u,del2v,       O                      del2u,del2v,
326       &                      myThid)       I                      myThid)
327           CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)           CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)
328           CALL MOM_CALC_RELVORT3(bi,bj,k,           CALL MOM_CALC_RELVORT3(bi,bj,k,
329       &                          del2u,del2v,hFacZ,zStar,myThid)       &                          del2u,del2v,hFacZ,zStar,myThid)
          IF ( writeDiag ) THEN  
            CALL WRITE_LOCAL_RL( 'del2u', 'I10', 1, del2u,  
      &                           bi,bj,k, myIter, myThid )  
            CALL WRITE_LOCAL_RL( 'del2v', 'I10', 1, del2v,  
      &                           bi,bj,k, myIter, myThid )  
            CALL WRITE_LOCAL_RL( 'dStar', 'I10', 1, dStar,  
      &                           bi,bj,k, myIter, myThid )  
            CALL WRITE_LOCAL_RL( 'zStar', 'I10', 1, zStar,  
      &                           bi,bj,k, myIter, myThid )  
          ENDIF  
        ENDIF  
   
 C-    Strain diagnostics:  
        IF ( writeDiag ) THEN  
         IF (snapshot_mdsio) THEN  
           CALL WRITE_LOCAL_RL('Ds','I10',1,strain,bi,bj,k,myIter,myThid)  
         ENDIF  
 #ifdef ALLOW_MNC  
         IF (useMNC .AND. snapshot_mnc) THEN  
           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strain,  
      &          offsets, myThid)  
         ENDIF  
 #endif /*  ALLOW_MNC  */  
        ENDIF  
 #ifdef ALLOW_DIAGNOSTICS  
        IF ( useDiagnostics ) THEN  
         CALL DIAGNOSTICS_FILL(strain, 'Strain  ',k,1,2,bi,bj,myThid)  
330         ENDIF         ENDIF
 #endif /* ALLOW_DIAGNOSTICS */  
331    
332  C---   Calculate dissipation terms for U and V equations  C---   Calculate dissipation terms for U and V equations
333    
334  C      in terms of tension and strain  C-     in terms of tension and strain
335         IF (useStrainTensionVisc) THEN         IF (useStrainTensionVisc) THEN
336  C        mask strain as if free-slip since side-drag is computed separately  C      use masked strain as if free-slip since side-drag is computed separately
          DO j=1-OLy,sNy+OLy  
           DO i=1-OLx,sNx+OLx  
             IF ( hFacZ(i,j).EQ.0. ) strain(i,j) = 0. _d 0  
           ENDDO  
          ENDDO  
337           CALL MOM_HDISSIP( bi, bj, k,           CALL MOM_HDISSIP( bi, bj, k,
338       I            hDiv, vort3, tension, strain, KE, hFacZ,       I            tension, strain, hFacZ,
339       I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,       I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
340       I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,       I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
341       O            guDiss, gvDiss,       O            guDiss, gvDiss,
342       I            myThid )       I            myThid )
343         ELSE         ELSE
344  C      in terms of vorticity and divergence  C-     in terms of vorticity and divergence
345           CALL MOM_VI_HDISSIP( bi, bj, k,           CALL MOM_VI_HDISSIP( bi, bj, k,
346       I            hDiv, vort3, tension, strain, KE, hFacZ,dStar,zStar,       I            hDiv, vort3, dStar, zStar, hFacZ,
347       I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,       I            viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
348       I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,       I            useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
349       O            guDiss, gvDiss,       O            guDiss, gvDiss,
350       &            myThid )       I            myThid )
351         ENDIF         ENDIF
352    
353  C---  Other dissipation terms in Zonal momentum equation  C---  Other dissipation terms in Zonal momentum equation
# Line 494  C--   if (momViscosity) end of block. Line 474  C--   if (momViscosity) end of block.
474  C-    Return to standard hfacZ (min-4) and mask vort3 accordingly:  C-    Return to standard hfacZ (min-4) and mask vort3 accordingly:
475  c     CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)  c     CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)
476    
 C-    Vorticity diagnostics:  
       IF ( writeDiag ) THEN  
         IF (snapshot_mdsio) THEN  
           CALL WRITE_LOCAL_RL('Z3','I10',1,vort3, bi,bj,k,myIter,myThid)  
         ENDIF  
 #ifdef ALLOW_MNC  
         IF (useMNC .AND. snapshot_mnc) THEN  
           CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3,  
      &          offsets, myThid)  
         ENDIF  
 #endif /*  ALLOW_MNC  */  
       ENDIF  
 #ifdef ALLOW_DIAGNOSTICS  
       IF ( useDiagnostics ) THEN  
         CALL DIAGNOSTICS_FILL(vort3,  'momVort3',k,1,2,bi,bj,myThid)  
       ENDIF  
 #endif /* ALLOW_DIAGNOSTICS */  
   
477  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
478    
479  C---  Prepare for Advection & Coriolis terms:  C---  Prepare for Advection & Coriolis terms:
480  C-    Mask relative vorticity and calculate absolute vorticity  C-    calculate absolute vorticity
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
          IF ( hFacZ(i,j).EQ.0. ) vort3(i,j) = 0.  
        ENDDO  
       ENDDO  
481        IF (useAbsVorticity)        IF (useAbsVorticity)
482       &  CALL MOM_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)       &  CALL MOM_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)
483    
# Line 753  C--   Set du/dt & dv/dt on boundaries to Line 710  C--   Set du/dt & dv/dt on boundaries to
710  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
711    
712        IF ( writeDiag ) THEN        IF ( writeDiag ) THEN
713            IF (useBiharmonicVisc) THEN
714             CALL WRITE_LOCAL_RL( 'del2u', 'I10', 1, del2u,
715         &                        bi,bj,k, myIter, myThid )
716             CALL WRITE_LOCAL_RL( 'del2v', 'I10', 1, del2v,
717         &                        bi,bj,k, myIter, myThid )
718             CALL WRITE_LOCAL_RL( 'dStar', 'I10', 1, dStar,
719         &                        bi,bj,k, myIter, myThid )
720             CALL WRITE_LOCAL_RL( 'zStar', 'I10', 1, zStar,
721         &                        bi,bj,k, myIter, myThid )
722            ENDIF
723          IF (snapshot_mdsio) THEN          IF (snapshot_mdsio) THEN
724           CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)
725             CALL WRITE_LOCAL_RL('Z3','I10',1,vort3BC,bi,bj,k,myIter,myThid)
726           CALL WRITE_LOCAL_RL('KE','I10',1,KE,     bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('KE','I10',1,KE,     bi,bj,k,myIter,myThid)
727           CALL WRITE_LOCAL_RL('D', 'I10',1,hDiv,   bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('D', 'I10',1,hDiv,   bi,bj,k,myIter,myThid)
728           CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)
729             CALL WRITE_LOCAL_RL( 'Ds', 'I10', 1, strainBC,
730         &                        bi,bj,k, myIter, myThid )
731           CALL WRITE_LOCAL_RL('Du','I10',1,guDiss, bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('Du','I10',1,guDiss, bi,bj,k,myIter,myThid)
732           CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)           CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)
733          ENDIF          ENDIF
# Line 765  C--   Set du/dt & dv/dt on boundaries to Line 735  C--   Set du/dt & dv/dt on boundaries to
735          IF (useMNC .AND. snapshot_mnc) THEN          IF (useMNC .AND. snapshot_mnc) THEN
736            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,
737       &          offsets, myThid)       &          offsets, myThid)
738              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3BC,
739         &          offsets, myThid)
740            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,
741       &          offsets, myThid)       &          offsets, myThid)
742            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,
743       &          offsets, myThid)       &          offsets, myThid)
744            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,
745       &          offsets, myThid)       &          offsets, myThid)
746              CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strainBC,
747         &          offsets, myThid)
748            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,
749       &          offsets, myThid)       &          offsets, myThid)
750            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,            CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,
# Line 781  C--   Set du/dt & dv/dt on boundaries to Line 755  C--   Set du/dt & dv/dt on boundaries to
755    
756  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
757        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
758            CALL DIAGNOSTICS_FILL(vort3BC,'momVort3',k,1,2,bi,bj,myThid)
759          CALL DIAGNOSTICS_FILL(KE,     'momKE   ',k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(KE,     'momKE   ',k,1,2,bi,bj,myThid)
760         IF (momViscosity) THEN         IF (momViscosity) THEN
761          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(tension,'Tension ',k,1,2,bi,bj,myThid)  
762          CALL DIAGNOSTICS_FILL(guDiss, 'Um_Diss ',k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(guDiss, 'Um_Diss ',k,1,2,bi,bj,myThid)
763          CALL DIAGNOSTICS_FILL(gvDiss, 'Vm_Diss ',k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(gvDiss, 'Vm_Diss ',k,1,2,bi,bj,myThid)
764         ENDIF         ENDIF
765           IF ( useVariableVisc .OR. useStrainTensionVisc ) THEN
766            CALL DIAGNOSTICS_FILL(tension, 'Tension ',k,1,2,bi,bj,myThid)
767            CALL DIAGNOSTICS_FILL(strainBC,'Strain  ',k,1,2,bi,bj,myThid)
768           ENDIF
769          CALL DIAGNOSTICS_FILL(gU(1-OLx,1-OLy,k,bi,bj),          CALL DIAGNOSTICS_FILL(gU(1-OLx,1-OLy,k,bi,bj),
770       &                                'Um_Advec',k,1,2,bi,bj,myThid)       &                                'Um_Advec',k,1,2,bi,bj,myThid)
771          CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),          CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),
Legend:
Removed from v.1.71  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.22