<p><b>duda</b> 2010-08-05 17:38:23 -0600 (Thu, 05 Aug 2010)</p><p>BRANCH COMMIT<br>
<br>
Add corrections to high-order advection setup for periodic planar domains.<br>
<br>
M    module_advection.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_advection.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_advection.F        2010-08-05 22:57:05 UTC (rev 464)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_advection.F        2010-08-05 23:38:23 UTC (rev 465)
@@ -58,6 +58,7 @@
       logical, parameter :: reset_poly = .true.
 
       real (kind=RKIND) :: rcell, cos2t, costsint, sin2t
+      real (kind=RKIND), dimension(grid%maxEdges) :: angle_2d
 
 !---
 
@@ -152,8 +153,18 @@
          else     ! On an x-y plane
 
             do i=1,n-1
-               xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
-               yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+
+               angle_2d(i) = grid%angleEdge%array(grid % EdgesOnCell % array(i,iCell))
+               iEdge = grid % EdgesOnCell % array(i,iCell)
+               if ( iCell .ne. grid % CellsOnEdge % array(1,iEdge)) &amp;
+                  angle_2d(i) = angle_2d(i) - pii
+
+!               xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
+!               yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+
+               xp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * cos(angle_2d(i))
+               yp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * sin(angle_2d(i))
+
             end do
 
          end if
@@ -271,9 +282,11 @@
                else
                   thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
                end if
-            else
-               xe(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (xv1 + xv2)
-               ye(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (yv1 + yv2)
+!            else
+!
+!               xe(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (xv1 + xv2)
+!               ye(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (yv1 + yv2)
+
             end if
   
          end do
@@ -312,13 +325,36 @@
                                             + 2.*sin2t*bmatrix(6,j)
                   end do
                end if
+
             else
-               do j=1,n
-                  deriv_two(j,1,iEdge) =   2.*xe(iEdge)*xe(iEdge)*bmatrix(4,j)  &amp;
-                                         + 2.*xe(iEdge)*ye(iEdge)*bmatrix(5,j)  &amp;
-                                         + 2.*ye(iEdge)*ye(iEdge)*bmatrix(6,j)
-                  deriv_two(j,2,iEdge) = deriv_two(j,1,iEdge)
-               end do
+
+               cos2t = cos(angle_2d(i))
+               sin2t = sin(angle_2d(i))
+               costsint = cos2t*sin2t
+               cos2t = cos2t**2
+               sin2t = sin2t**2
+
+!               do j=1,n
+!
+!                  deriv_two(j,1,iEdge) =   2.*xe(iEdge)*xe(iEdge)*bmatrix(4,j)  &amp;
+!                                         + 2.*xe(iEdge)*ye(iEdge)*bmatrix(5,j)  &amp;
+!                                         + 2.*ye(iEdge)*ye(iEdge)*bmatrix(6,j)
+!               end do
+
+               if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
+                  do j=1,n
+                     deriv_two(j,1,iEdge) =   2.*cos2t*bmatrix(4,j)  &amp;
+                                            + 2.*costsint*bmatrix(5,j)  &amp;
+                                            + 2.*sin2t*bmatrix(6,j)
+                  end do
+               else
+                  do j=1,n
+                     deriv_two(j,2,iEdge) =   2.*cos2t*bmatrix(4,j)  &amp;
+                                            + 2.*costsint*bmatrix(5,j)  &amp;
+                                            + 2.*sin2t*bmatrix(6,j)
+                  end do
+               end if
+
             end if
          end do
  
@@ -326,6 +362,25 @@
 
       if (debug) stop
 
+
+!      write(0,*) ' check for deriv2 coefficients, iEdge 4 '
+!
+!      iEdge = 4
+!      j = 1
+!      iCell = grid % cellsOnEdge % array(1,iEdge)
+!      write(0,*) ' j, icell, coef ',j,iCell,deriv_two(j,1,iEdge)
+!      do j=2,7
+!         write(0,*) ' j, icell, coef ',j,grid % CellsOnCell % array(j-1,iCell),deriv_two(j,1,iEdge)
+!      end do
+!
+!      j = 1
+!      iCell = grid % cellsOnEdge % array(2,iEdge)
+!      write(0,*) ' j, icell, coef ',j,iCell,deriv_two(j,2,iEdge)
+!      do j=2,7
+!         write(0,*) ' j, icell, coef ',j,grid % CellsOnCell % array(j-1,iCell),deriv_two(j,2,iEdge)
+!      end do
+!      stop
+
    end subroutine initialize_advection_rk
 
 
@@ -755,7 +810,7 @@
 
 !  check to see if we are reaching outside the halo
 
-        if (debug) write(0,*) ' points ', n
+         if (debug) write(0,*) ' points ', n
 
          do_the_cell = .true.
          do i=1,n
@@ -763,11 +818,11 @@
          end do
 
 
-         if ( .not. do_the_cell ) cycle
+         if (.not. do_the_cell) cycle
 
 
 !  compute poynomial fit for this cell if all needed neighbors exist
-         if ( grid % on_a_sphere ) then
+         if (grid % on_a_sphere) then
 
             xc(1) = grid % xCell % array(iCell)/a
             yc(1) = grid % yCell % array(iCell)/a

</font>
</pre>