<p><b>ringler@lanl.gov</b> 2010-08-06 11:44:48 -0600 (Fri, 06 Aug 2010)</p><p><br>
incremental update with minor code clean up<br>
</p><hr noshade><pre><font color="gray">Modified: branches/dissipation/src/core_ocean/module_time_integration.F
===================================================================
--- branches/dissipation/src/core_ocean/module_time_integration.F        2010-08-06 15:55:08 UTC (rev 467)
+++ branches/dissipation/src/core_ocean/module_time_integration.F        2010-08-06 17:44:48 UTC (rev 468)
@@ -461,7 +461,7 @@
 
                !
                ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
-               !                    only valid for h_mom_eddy_visc2 == constant
+               !                    strictly only valid for h_mom_eddy_visc2 == constant
                !
                u_diffusion =   ( divergence(k,cell2)  - divergence(k,cell1) ) / dcEdge(iEdge)  &amp;
                               -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
@@ -493,7 +493,7 @@
    
                !
                ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
-               !                    only valid for h_mom_eddy_visc4 == constant
+               !                    strictly only valid for h_mom_eddy_visc4 == constant
                !
                u_diffusion =   ( divergence(k,cell2)  - divergence(k,cell1) ) / dcEdge(iEdge)  &amp;
                               -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
@@ -505,9 +505,11 @@
 
          delsq_circulation(:,:) = 0.0
          do iEdge=1,nEdges
+            vertex1 = verticesOnEdge(1,iEdge)
+            vertex2 = verticesOnEdge(2,iEdge)
             do k=1,nVertLevels
-               delsq_circulation(k,verticesOnEdge(1,iEdge)) = delsq_circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * delsq_u(k,iEdge)
-               delsq_circulation(k,verticesOnEdge(2,iEdge)) = delsq_circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * delsq_u(k,iEdge)
+               delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) - dcEdge(iEdge) * delsq_u(k,iEdge)
+               delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) + dcEdge(iEdge) * delsq_u(k,iEdge)
             end do
          end do
          do iVertex=1,nVertices
@@ -562,7 +564,7 @@
       !
       ! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
       !
-       do iEdge=1,grid % nEdgesSolve
+      do iEdge=1,grid % nEdgesSolve
          cell1 = cellsOnEdge(1,iEdge)
          cell2 = cellsOnEdge(2,iEdge)
 
@@ -597,13 +599,6 @@
            - 1.0e-3*u(nVertLevels,iEdge) &amp;
              *sqrt(2.0*ke_edge(nVertLevels,iEdge))
 
-         ! mrp 100603 The following method is more straight forward, 
-         ! that the above method of computing ke_edge, but I have
-         ! not verified that v is working correctly yet.
-         !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge)  &amp;
-         !  - 1.0e-3*u(nVertLevels,iEdge) &amp;
-         !    *sqrt(u(nVertLevels,iEdge)**2 + v(nVertLevels,iEdge)**2)
-
          ! old bottom drag, just linear friction
          ! du/dt = u/tau where tau=100 days.
          !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge)  &amp;
@@ -672,7 +667,7 @@
 
       integer :: iCell, iEdge, k, iTracer, cell1, cell2, upwindCell,&amp;
         nEdges, nCells, nVertLevels
-      real (kind=RKIND) :: h_tracer_eddy_visc2, h_tracer_eddy_visc4, areaCell1, areaCell2, tracer_turb_flux
+      real (kind=RKIND) :: h_tracer_eddy_visc2, h_tracer_eddy_visc4, invAreaCell1, invAreaCell2, tracer_turb_flux
       real (kind=RKIND) :: flux, tracer_edge, r
       real (kind=RKIND) :: dist
       real (kind=RKIND), dimension(:), pointer :: &amp;
@@ -816,19 +811,25 @@
 
       if ( h_tracer_eddy_visc2 &gt; 0.0 ) then
 
+
+      !
+      ! mask
+      !
+
+
          do iEdge=1,grid % nEdges
             cell1 = grid % cellsOnEdge % array(1,iEdge)
             cell2 = grid % cellsOnEdge % array(2,iEdge)
-            areaCell1 = areaCell(cell1)
-            areaCell2 = areaCell(cell2)
+            invAreaCell1 = 1.0/areaCell(cell1)
+            invAreaCell2 = 1.0/areaCell(cell2)
 
             do k=1,grid % nVertLevels
               do iTracer=1,num_tracers
                  tracer_turb_flux = h_tracer_eddy_visc2*(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1))/dcEdge(iEdge)
                  flux = dvEdge (iEdge) * h_edge(k,iEdge) * tracer_turb_flux
-                 tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) + flux/areaCell1
-                 tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) - flux/areaCell2
-              enddo
+                 tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) + flux * invAreaCell1
+                 tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) - flux * invAreaCell2
+              end do
             end do
 
          end do
@@ -866,16 +867,16 @@
          do iEdge=1,grid % nEdges
             cell1 = grid % cellsOnEdge % array(1,iEdge)
             cell2 = grid % cellsOnEdge % array(2,iEdge)
-            areaCell1 = areaCell(cell1)
-            areaCell2 = areaCell(cell2)
+            invAreaCell1 = 1.0 / areaCell(cell1)
+            invAreaCell2 = 1.0 / areaCell(cell2)
 
             do k=1,grid % nVertLevels
             do iTracer=1,num_tracers
                tracer_turb_flux = h_tracer_eddy_visc4*(delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
                flux = dvEdge (iEdge) * tracer_turb_flux
 
-               tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell1
-               tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell2
+               tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux * invAreaCell1
+               tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux * invAreaCell2
 
             end do
             enddo

</font>
</pre>