<p><b>duda</b> 2011-09-19 14:47:24 -0600 (Mon, 19 Sep 2011)</p><p>BRANCH COMMIT<br>
<br>
Bug fixes to omega computation from Bill, uncovered by Soyoung's testing.<br>
<br>
<br>
M    src/core_nhyd_atmos/module_time_integration.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2011-09-19 19:44:04 UTC (rev 1010)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2011-09-19 20:47:24 UTC (rev 1011)
@@ -1180,9 +1180,9 @@
           w(1,cell1) = w(1,cell1) + zb(1,1,iEdge)*flux/(cf1*rho_zz(1,cell1)+cf2*rho_zz(2,cell1)+cf3*rho_zz(3,cell1))
 !3rd order stencil
           if (config_theta_adv_order == 3) then
-             w(1,cell2) = w(1,cell2) - sign(1.,ru(1,iEdge))*config_coef_3rd_order      &amp;
+             w(1,cell2) = w(1,cell2) - sign(1.,flux)*config_coef_3rd_order      &amp;
                                       *zb3(1,2,iEdge)*flux/(cf1*rho_zz(1,cell2)+cf2*rho_zz(2,cell2)+cf3*rho_zz(3,cell2))
-             w(1,cell1) = w(1,cell1) + sign(1.,ru(1,iEdge))*config_coef_3rd_order      &amp;
+             w(1,cell1) = w(1,cell1) + sign(1.,flux)*config_coef_3rd_order      &amp;
                                       *zb3(1,1,iEdge)*flux/(cf1*rho_zz(1,cell1)+cf2*rho_zz(2,cell1)+cf3*rho_zz(3,cell1))
           end if
 
@@ -1192,9 +1192,9 @@
             w(k,cell1) = w(k,cell1) + zb(k,1,iEdge)*flux/(fzm(k)*rho_zz(k,cell1)+fzp(k)*rho_zz(k-1,cell1))
 !3rd order! stencil
             if (config_theta_adv_order ==3) then
-               w(k,cell2) = w(k,cell2) - sign(1.,ru(k,iEdge))*config_coef_3rd_order    &amp;
+               w(k,cell2) = w(k,cell2) - sign(1.,flux)*config_coef_3rd_order    &amp;
                                         *zb3(k,2,iEdge)*flux/(fzm(k)*rho_zz(k,cell2)+fzp(k)*rho_zz(k-1,cell2)) 
-               w(k,cell1) = w(k,cell1) + sign(1.,ru(k,iEdge))*config_coef_3rd_order    &amp;
+               w(k,cell1) = w(k,cell1) + sign(1.,flux)*config_coef_3rd_order    &amp;
                                         *zb3(k,1,iEdge)*flux/(fzm(k)*rho_zz(k,cell1)+fzp(k)*rho_zz(k-1,cell1)) 
             end if
           enddo
@@ -3485,8 +3485,8 @@
       type (diag_type), intent(inout) :: diag
       type (mesh_type), intent(inout) :: grid
 
-      integer :: k,iCell,iEdge,i,iCell1,iCell2
-      real (kind=RKIND) :: p0, rcv
+      integer :: k,iCell,iEdge,i,iCell1,iCell2, cell1, cell2
+      real (kind=RKIND) :: p0, rcv, flux
 
       rcv = rgas / (cp-rgas)
       p0 = 1.e5  ! this should come from somewhere else...
@@ -3506,16 +3506,55 @@
          end do
       end do
 
-      ! Compute w from rho_zz and rw
+!      ! Compute w from rho_zz and rw
+!      do iCell=1,grid%nCells
+!         diag % rw % array(1,iCell) = 0.
+!         diag % rw % array(grid%nVertLevels+1,iCell) = 0.
+!         do k=2,grid%nVertLevels
+!            diag % rw % array(k,iCell) = state % w % array(k,iCell)     &amp;
+!                          * (grid % fzp % array(k) * state % rho_zz % array(k-1,iCell) + grid % fzm % array(k) * state % rho_zz % array(k,iCell))
+!         end do
+!      end do
+
+
+! WCS bug fix 20110916
+
+      ! Compute rw (i.e. rho_zz * omega) from rho_zz, w, and ru.
+      ! We are reversing the procedure we use in subroutine recover_large_step_variables.
+      ! first, the piece that depends on w.
       do iCell=1,grid%nCells
          diag % rw % array(1,iCell) = 0.
          diag % rw % array(grid%nVertLevels+1,iCell) = 0.
          do k=2,grid%nVertLevels
             diag % rw % array(k,iCell) = state % w % array(k,iCell)     &amp;
-                          * (grid % fzp % array(k) * state % rho_zz % array(k-1,iCell) + grid % fzm % array(k) * state % rho_zz % array(k,iCell))
+                          * (grid % fzp % array(k) * state % rho_zz % array(k-1,iCell) + grid % fzm % array(k) * state % rho_zz % array(k,iCell)) &amp;
+                          * (grid % fzp % array(k) * grid % zz % array(k-1,iCell) + grid % fzm % array(k) * grid % zz % array(k,iCell))
          end do
       end do
+      ! next, the piece that depends on ru
+      do iEdge=1,grid%nEdges
+        cell1 = grid % CellsOnEdge % array(1,iEdge)
+        cell2 = grid % CellsOnEdge % array(2,iEdge)
+          do k = 2, grid % nVertLevels
+            flux = (grid % fzm % array(k) * diag % ru % array(k,iEdge)+grid % fzp % array(k) * diag % ru % array(k-1,iEdge))
+            diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + grid % zb % array(k,2,iEdge)*flux   &amp;
+                          * (grid % fzp % array(k) * grid % zz % array(k-1,cell2) + grid % fzm % array(k) * grid % zz % array(k,cell2))
+            diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - grid % zb % array(k,1,iEdge)*flux   &amp;
+                          * (grid % fzp % array(k) * grid % zz % array(k-1,cell1) + grid % fzm % array(k) * grid % zz % array(k,cell1))
+!3rd order! stencil
+            if (config_theta_adv_order ==3) then
+               diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + sign(1.,flux)*config_coef_3rd_order    &amp;
+                                        * grid % zb3 % array(k,2,iEdge)*flux                                    &amp;
+                          * (grid % fzp % array(k) * grid % zz % array(k-1,cell2) + grid % fzm % array(k) * grid % zz % array(k,cell2))
+               diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - sign(1.,flux)*config_coef_3rd_order    &amp;
+                                        * grid % zb3 % array(k,1,iEdge)*flux                                    &amp;
+                          * (grid % fzp % array(k) * grid % zz % array(k-1,cell1) + grid % fzm % array(k) * grid % zz % array(k,cell1))
+            end if
+          enddo
+      enddo
 
+!  end WCS bug fix
+
       do iCell = 1, grid % nCells
          do k=1,grid % nVertLevels
             diag % rho_p % array(k,iCell) = state % rho_zz % array(k,iCell) - diag % rho_base % array(k,iCell)

</font>
</pre>