<p><b>gaw06e@fsu.edu</b> 2011-05-30 14:15:01 -0600 (Mon, 30 May 2011)</p><p>- modified compute_tend, compute_scalar_tend, compute_solve diagnostics for dvEdge approximation (where boundaryEdge = 1)<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F        2011-05-29 00:18:43 UTC (rev 859)
+++ branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F        2011-05-30 20:15:01 UTC (rev 860)
@@ -264,6 +264,7 @@
       integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
       integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
       integer, dimension(:,:), pointer :: boundaryCell
+      integer, dimension(:,:), pointer :: boundaryEdge
       real (kind=RKIND) :: r, u_diffusion
 
       real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
@@ -273,6 +274,7 @@
       real (kind=RKIND), dimension(:,:), pointer :: u_src
       real (kind=RKIND), parameter :: rho_ref = 1000.0
       real (kind=RKIND) :: ke_edge
+      real (kind=RKIND) :: thisdvEdge
 
       h           =&gt; s % h % array
       u           =&gt; s % u % array
@@ -314,7 +316,8 @@
       u_src =&gt; grid % u_src % array
       
       boundaryCell =&gt; grid % boundaryCell % array
-
+      boundaryEdge =&gt; grid % boundaryEdge % array
+      
       !
       ! Compute height tendency for each cell
       !
@@ -322,8 +325,14 @@
       do iEdge=1,nEdges
          cell1 = cellsOnEdge(1,iEdge)
          cell2 = cellsOnEdge(2,iEdge)
+         if (boundaryEdge(1,iEdge).eq.1) then
+            thisdvEdge = dvEdge(iEdge) / 2
+         else
+            thisdvEdge = dvEdge(iEdge)
+         end if
          do k=1,nVertLevels
-            flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+            !flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+            flux = u(k,iEdge) * thisdvEdge * h_edge(k,iEdge)
             tend_h(k,cell1) = tend_h(k,cell1) - flux
             tend_h(k,cell2) = tend_h(k,cell2) + flux
          end do
@@ -395,24 +404,29 @@
       end do
 #endif
 
-     ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
-     !                    only valid for visc == constant
-     if (config_h_mom_eddy_visc2 &gt; 0.0) then
-        do iEdge=1,grid % nEdgesSolve
-           cell1 = cellsOnEdge(1,iEdge)
-           cell2 = cellsOnEdge(2,iEdge)
-           vertex1 = verticesOnEdge(1,iEdge)
-           vertex2 = verticesOnEdge(2,iEdge)
+      ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
+      !                    only valid for visc == constant
+      if (config_h_mom_eddy_visc2 &gt; 0.0) then
+         do iEdge=1,grid % nEdgesSolve
+            cell1 = cellsOnEdge(1,iEdge)
+            cell2 = cellsOnEdge(2,iEdge)
+            vertex1 = verticesOnEdge(1,iEdge)
+            vertex2 = verticesOnEdge(2,iEdge)
+            if (boundaryEdge(1,iEdge).eq.1) then
+               thisdvEdge = dvEdge(iEdge) / 2
+            else
+               thisdvEdge = dvEdge(iEdge)
+            end if
+            do k=1,nVertLevels
+               u_diffusion =   ( divergence(k,cell2)  -  divergence(k,cell1) ) / dcEdge(iEdge) &amp;
+                  ! -(vorticity(k,vertex2)  - vorticity(k,vertex1) ) / dvEdge(iEdge)
+                  -(vorticity(k,vertex2)  - vorticity(k,vertex1) ) / thisdvEdge
+               u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
+               tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+            end do
+         end do
+      end if
 
-           do k=1,nVertLevels
-              u_diffusion =   ( divergence(k,cell2)  -  divergence(k,cell1) ) / dcEdge(iEdge) &amp;
-                   -(vorticity(k,vertex2)  - vorticity(k,vertex1) ) / dvEdge(iEdge)
-              u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
-              tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
-           end do
-        end do
-     end if
-
      !
      ! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="black">abla^4 u
      !   computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity
@@ -461,28 +475,33 @@
            end do
         end do
 
-        ! Divergence using </font>
<font color="red">abla^2 u
-        delsq_divergence(:,:) = 0.0
-        do iEdge=1,nEdges
-           cell1 = cellsOnEdge(1,iEdge)
-           cell2 = cellsOnEdge(2,iEdge)
-           do k=1,nVertLevels
-              delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &amp;
-                   + delsq_u(k,iEdge)*dvEdge(iEdge)
-              delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &amp;
-                   - delsq_u(k,iEdge)*dvEdge(iEdge)
-           end do
-        end do
-        do iCell = 1,nCells
-           if (boundaryCell(1,iCell).eq.1) then 
-              r = 1.0 / ( areaCell(iCell) * 2.0 )
+         ! Divergence using </font>
<font color="red">abla^2 u
+         delsq_divergence(:,:) = 0.0
+         do iEdge=1,nEdges
+            cell1 = cellsOnEdge(1,iEdge)
+            cell2 = cellsOnEdge(2,iEdge)
+            if (boundaryEdge(1,iEdge).eq.1) then
+               thisdvEdge = dvEdge(iEdge) / 2
+            else
+               thisdvEdge = dvEdge(iEdge)
+            end if
+            do k=1,nVertLevels
+               delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &amp;
+                     + delsq_u(k,iEdge)*thisdvEdge
+               delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &amp;
+                     - delsq_u(k,iEdge)*thisdvEdge
+            end do
+         end do
+         do iCell = 1,nCells
+            if (boundaryCell(1,iCell).eq.1) then 
+               r = 1.0 / ( areaCell(iCell) * 2.0 )
            else
-              r = 1.0 / areaCell(iCell)
-           end if
-           do k = 1,nVertLevels
-              delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
-           end do
-        end do
+               r = 1.0 / areaCell(iCell)
+            end if
+            do k = 1,nVertLevels
+               delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
+            end do
+         end do
 
         ! Compute - \kappa </font>
<font color="black">abla^4 u 
         ! as  </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="gray">abla^2 u) )
@@ -491,13 +510,17 @@
            cell2 = cellsOnEdge(2,iEdge)
            vertex1 = verticesOnEdge(1,iEdge)
            vertex2 = verticesOnEdge(2,iEdge)
-
+            if (boundaryEdge(1,iEdge).eq.1) then
+               thisdvEdge = dvEdge(iEdge) / 2
+            else
+               thisdvEdge = dvEdge(iEdge)
+            end if
            do k=1,nVertLevels
 
               u_diffusion = (  delsq_divergence(k,cell2) &amp;
                    - delsq_divergence(k,cell1) ) / dcEdge(iEdge)  &amp;
                    -(  delsq_vorticity(k,vertex2) &amp;
-                   - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+                   - delsq_vorticity(k,vertex1) ) / thisdvEdge
 
               u_diffusion = config_h_mom_eddy_visc4 * u_diffusion
               tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
@@ -569,6 +592,8 @@
       real (kind=RKIND) :: coef_3rd_order
       real (kind=RKIND), dimension(:,:), pointer :: u, h_edge
 
+      real (kind=RKIND) :: thisdvEdge
+
       u           =&gt; s % u % array
       h_edge      =&gt; s % h_edge % array
       dcEdge      =&gt; grid % dcEdge % array
@@ -594,6 +619,11 @@
       do iEdge=1,grid % nEdges
             cell1 = cellsOnEdge(1,iEdge)
             cell2 = cellsOnEdge(2,iEdge)
+            if (boundaryEdge(1,iEdge).eq.1) then
+               thisdvEdge = dvEdge(iEdge) / 2
+            else
+               thisdvEdge = dvEdge(iEdge)
+            end if
             if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
                if (boundaryCell(1,cell1).eq.1) then 
                   invAreaCell1 = 1.0 / ( areaCell(cell1) * 2.0 )
@@ -608,7 +638,7 @@
                do k=1,grid % nVertLevels
                   do iTracer=1,grid % nTracers
                      tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
-                     flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+                     flux = u(k,iEdge) * thisdvEdge * h_edge(k,iEdge) * tracer_edge
                      tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1
                      tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2
                   end do 
@@ -739,6 +769,11 @@
          do iEdge=1,grid % nEdges
             cell1 = cellsOnEdge(1,iEdge)
             cell2 = cellsOnEdge(2,iEdge)
+            if (boundaryEdge(1,iEdge).eq.1) then
+               thisdvEdge = dvEdge(iEdge) / 2
+            else
+               thisdvEdge = dvEdge(iEdge)
+            end if
             if (boundaryCell(1,cell1).eq.1) then 
                invAreaCell1 = 1.0 / ( areaCell(cell1) * 2.0 )
             else
@@ -756,7 +791,7 @@
                     *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
 
                  ! div(h \kappa_2 </font>
<font color="gray">abla \phi) at cell center
-                 flux = dvEdge(iEdge) * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
+                 flux = thisdvEdge * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
                  tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
                  tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) - flux * invAreaCell2
               end do
@@ -789,13 +824,17 @@
          do iEdge=1,grid % nEdges
             cell1 = cellsOnEdge(1,iEdge)
             cell2 = cellsOnEdge(2,iEdge)
-
+            if (boundaryEdge(1,iEdge).eq.1) then
+               thisdvEdge = dvEdge(iEdge) / 2
+            else
+               thisdvEdge = dvEdge(iEdge)
+            end if
             do k=1,grid % nVertLevels
               do iTracer=1, grid % nTracers
                  delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &amp;
-                    + dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+                    + thisdvEdge * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
                  delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &amp;
-                    - dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+                    - thisdvEdge * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
               end do
             end do
 
@@ -819,6 +858,11 @@
          do iEdge=1,grid % nEdges
             cell1 = grid % cellsOnEdge % array(1,iEdge)
             cell2 = grid % cellsOnEdge % array(2,iEdge)
+            if (boundaryEdge(1,iEdge).eq.1) then
+               thisdvEdge = dvEdge(iEdge) / 2
+            else
+               thisdvEdge = dvEdge(iEdge)
+            end if
             if (boundaryCell(1,cell1).eq.1) then 
                invAreaCell1 = 1.0 / ( areaCell(cell1) * 2.0 )
             else
@@ -835,7 +879,7 @@
             do k=1,grid % nVertLevels
             do iTracer=1,grid % nTracers
                tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) / dcEdge(iEdge)
-               flux = dvEdge(iEdge) * tracer_turb_flux
+               flux = thisdvEdge * tracer_turb_flux
                tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1 * boundaryMask(k,iEdge)
                tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2 * boundaryMask(k,iEdge)
             end do
@@ -882,6 +926,7 @@
       real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
       real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
       real (kind=RKIND) :: coef_3rd_order
+      real (kind=RKIND) :: thisdvEdge
 
       h           =&gt; s % h % array
       u           =&gt; s % u % array
@@ -1094,22 +1139,32 @@
       do iEdge=1,nEdges
          cell1 = cellsOnEdge(1,iEdge)
          cell2 = cellsOnEdge(2,iEdge)
+         if (boundaryEdge(1,iEdge).eq.1) then
+            thisdvEdge = dvEdge(iEdge) / 2
+         else
+            thisdvEdge = dvEdge(iEdge)
+         end if
          if (cell1 &lt;= nCells) then
             do k=1,nVertLevels
-              divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+              divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*thisdvEdge
             enddo
          endif
          if(cell2 &lt;= nCells) then
             do k=1,nVertLevels
-              divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+              divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*thisdvEdge
             enddo
          end if
       end do
       do iCell = 1,nCells
-        r = 1.0 / areaCell(iCell)
-        do k = 1,nVertLevels
-           divergence(k,iCell) = divergence(k,iCell) * r
-        enddo
+         if (boundaryCell(1,iCell).eq.1) then 
+            r = 1.0 / ( areaCell(iCell) * 2.0 )
+         else
+            r = 1.0 / areaCell(iCell)
+         end if
+        !r = 1.0 / areaCell(iCell)
+         do k = 1,nVertLevels
+            divergence(k,iCell) = divergence(k,iCell) * r
+         enddo
       enddo
 
       !
@@ -1117,14 +1172,24 @@
       !
       ke(:,:) = 0.0
       do iCell=1,nCells
+         if (boundaryCell(1,iCell).eq.1) then 
+            r = 1.0 / ( areaCell(iCell) * 2.0 )
+         else
+            r = 1.0 / areaCell(iCell)
+         end if
          do i=1,nEdgesOnCell(iCell)
             iEdge = edgesOnCell(i,iCell)
+            if (boundaryEdge(1,iEdge).eq.1) then
+               thisdvEdge = dvEdge(iEdge) / 2
+            else
+               thisdvEdge = dvEdge(iEdge)
+            end if
             do k=1,nVertLevels
-               ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+               ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * thisdvEdge * u(k,iEdge)**2.0
             end do
          end do
          do k=1,nVertLevels
-            ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
+            ke(k,iCell) = ke(k,iCell) * r
          end do
       end do
 
@@ -1179,9 +1244,14 @@
       !   ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
       !
       do iEdge = 1,nEdges
+         if (boundaryEdge(1,iEdge).eq.1) then
+            thisdvEdge = dvEdge(iEdge) / 2
+         else
+            thisdvEdge = dvEdge(iEdge)
+         end if
          do k = 1,nVertLevels
            gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &amp;
-                              dvEdge(iEdge)
+                              thisdvEdge
          enddo
       enddo
 
@@ -1219,10 +1289,15 @@
       do iVertex = 1, nVertices
        do i=1,grid % vertexDegree
          iCell = cellsOnVertex(i,iVertex)
+         if (boundaryCell(1,iCell).eq.1) then 
+            r = 1.0 / ( areaCell(iCell) * 2.0 )
+         else
+            r = 1.0 / areaCell(iCell)
+         end if
          if (iCell &lt;= nCells) then
            do k = 1,nVertLevels
-             pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
-             vorticity_cell(k,iCell) = vorticity_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) / areaCell(iCell)
+             pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) * r
+             vorticity_cell(k,iCell) = vorticity_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) * r
            enddo
          endif
        enddo
@@ -1285,7 +1360,7 @@
    !
    ! Input: grid - grid metadata
    !
-   ! Output: tend_u set to zero at boundaryEdge == 1 locations
+   ! Output: tend_u set to zero at boundaryEdge == 2 locations
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 

</font>
</pre>