<p><b>duda</b> 2012-01-19 16:55:17 -0700 (Thu, 19 Jan 2012)</p><p>BRANH COMMIT<br>
<br>
Simplifications to the computation of w from Sang-Hun.<br>
<br>
To do: Ultimately, we should be able to get rid of zf and zf3 altogether.<br>
<br>
<br>
M    src/core_nhyd_atmos/mpas_atm_time_integration.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F        2012-01-19 20:04:56 UTC (rev 1398)
+++ branches/atmos_physics/src/core_nhyd_atmos/mpas_atm_time_integration.F        2012-01-19 23:55:17 UTC (rev 1399)
@@ -674,14 +674,19 @@
       type (tend_type) :: tend
       type (diag_type) :: diag
       type (mesh_type) :: grid
-      integer :: iCell, iEdge, k, cell1, cell2
+      !SHP-w
+      integer :: iCell, iEdge, k, cell1, cell2, coef_3rd_order
       integer, dimension(:,:), pointer :: cellsOnEdge
       real (kind=RKIND), dimension(:,:,:), pointer :: zf, zf3
       real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, areaCell
       real (kind=RKIND) :: flux
+      !SHP-w
+      coef_3rd_order = config_coef_3rd_order
+      if(config_theta_adv_order /=3) coef_3rd_order = 0
 
-      zf =&gt; grid % zf % array
-      zf3 =&gt; grid % zf3 % array
+      !SHP-w
+      !zf =&gt; grid % zf % array
+      !zf3 =&gt; grid % zf3 % array
       fzm =&gt; grid % fzm % array
       fzp =&gt; grid % fzp % array
       dvEdge =&gt; grid % dvEdge % array
@@ -708,18 +713,27 @@
          cell1 = cellsOnEdge(1,iEdge)
          cell2 = cellsOnEdge(2,iEdge)
 
+         !do k = 2, grid%nVertLevels
+         !   flux = fzm(k) * tend % u % array(k,iEdge) + fzp(k) * tend % u % array(k-1,iEdge)
+         !   tend % w % array(k,cell2) = tend % w % array(k,cell2) + zf(k,2,iEdge)*flux
+         !   tend % w % array(k,cell1) = tend % w % array(k,cell1) - zf(k,1,iEdge)*flux
+!3rd order stencil
+         !   if (config_theta_adv_order == 3) then
+         !      tend % w % array(k,cell2) = tend % w % array(k,cell2) + sign(1.0_RKIND,tend % u % array(k,iEdge))  &amp;
+         !                                                      *config_coef_3rd_order*zf3(k,2,iEdge)*flux
+         !      tend % w % array(k,cell1) = tend % w % array(k,cell1) - sign(1.0_RKIND,tend % u % array(k,iEdge))  &amp;
+         !                                                      *config_coef_3rd_order*zf3(k,1,iEdge)*flux
+         !   end if
+ 
+         !SHP-w
          do k = 2, grid%nVertLevels
             flux = fzm(k) * tend % u % array(k,iEdge) + fzp(k) * tend % u % array(k-1,iEdge)
-            tend % w % array(k,cell2) = tend % w % array(k,cell2) + zf(k,2,iEdge)*flux
-            tend % w % array(k,cell1) = tend % w % array(k,cell1) - zf(k,1,iEdge)*flux
-!3rd order stencil
-            if (config_theta_adv_order == 3) then
-               tend % w % array(k,cell2) = tend % w % array(k,cell2) + sign(1.0_RKIND,tend % u % array(k,iEdge))  &amp;
-                                                               *config_coef_3rd_order*zf3(k,2,iEdge)*flux
-               tend % w % array(k,cell1) = tend % w % array(k,cell1) - sign(1.0_RKIND,tend % u % array(k,iEdge))  &amp;
-                                                               *config_coef_3rd_order*zf3(k,1,iEdge)*flux
-            end if
-               
+            tend % w % array(k,cell2) = tend % w % array(k,cell2)   &amp;
+                     + (grid % zb % array(k,2,iEdge) + coef_3rd_order*sign(1.0_RKIND,tend % u % array(k,iEdge))*grid %zb3 % array(k,2,iEdge))*flux   &amp;
+                     * (fzm(k) * grid % zz % array(k,cell2) + fzp(k) * grid % zz % array(k-1,cell2)) 
+            tend % w % array(k,cell1) = tend % w % array(k,cell1)   &amp;
+                     - (grid % zb % array(k,1,iEdge) + coef_3rd_order*sign(1.0_RKIND,tend % u % array(k,iEdge))*grid %zb3 % array(k,1,iEdge))*flux   &amp;
+                     * (fzm(k) * grid % zz % array(k,cell1) + fzp(k) * grid % zz % array(k-1,cell1)) 
          end do
 
       end do
@@ -2084,7 +2098,7 @@
 
                do k=2,nVertLevels
 
-                  kr = min(nVertLevels,k+ nint(.5-sign(.5,zx(k,iEdge)+zx(k+1,iEdge))))
+                  kr = min(nVertLevels,k+ nint(.5-sign(0.5_RKIND,zx(k,iEdge)+zx(k+1,iEdge))))
                   kl = min(nVertLevels,2*k+1-kr)
 
                   pr = pp(k,cell2)+.5*(zgrid(k   ,cell1)+zgrid(k +1,cell1)-zgrid(k ,cell2)-zgrid(k +1,cell2))   &amp;
@@ -2669,7 +2683,7 @@
          do k=2,nVertLevels
 
             tend_w(k,iCell) = tend_w(k,iCell)/areaCell(iCell) -rdzu(k)*(wdwz(k+1)-wdwz(k))    &amp;
-!SHP-w
+!SHP-buoy
                                   - cqw(k,iCell)*( rdzu(k)*(pp(k,iCell)-pp(k-1,iCell))        &amp;
                                   + gravity*  &amp;
                                    ( fzm(k)*(rb(k,iCell)*(qtot(k,iCell)) +         &amp;
@@ -3331,9 +3345,14 @@
       type (diag_type), intent(inout) :: diag
       type (mesh_type), intent(inout) :: grid
 
-      integer :: k,iCell,iEdge,i,iCell1,iCell2, cell1, cell2
+      !SHP-w
+      integer :: k,iCell,iEdge,i,iCell1,iCell2, cell1, cell2, coef_3rd_order
       real (kind=RKIND) :: p0, rcv, flux
 
+      !SHP-w
+      coef_3rd_order = config_coef_3rd_order
+      if(config_theta_adv_order /=3) coef_3rd_order = 0
+
       rcv = rgas / (cp-rgas)
       p0 = 1.e5  ! this should come from somewhere else...
 
@@ -3378,26 +3397,42 @@
          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.0_RKIND,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.0_RKIND,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
+  
+      !SHP-w
+      ! 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;
+            diag % rw % array(k,cell2) = diag % rw % array(k,cell2)   &amp;
+                          + (grid % zb % array(k,2,iEdge) + coef_3rd_order * sign(1.0_RKIND,flux) * 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) - grid % zb % array(k,1,iEdge)*flux   &amp;
+            diag % rw % array(k,cell1) = diag % rw % array(k,cell1)   &amp;
+                          - (grid % zb % array(k,1,iEdge) + coef_3rd_order * sign(1.0_RKIND,flux) * 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))
-!3rd order! stencil
-            if (config_theta_adv_order ==3) then
-               diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + sign(1.0_RKIND,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.0_RKIND,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 do
+      end do
 
 !  end WCS bug fix
 

</font>
</pre>