<p><b>duda</b> 2010-10-15 15:53:14 -0600 (Fri, 15 Oct 2010)</p><p>BRANCH COMMIT<br>
<br>
Update atmos_physics branch with respect to atmos_nonhydrostatic.<br>
<br>
<br>
M    src/core_physics/module_driver_microphysics.F<br>
M    src/core_physics/module_physics_driver.F<br>
M    src/core_physics/module_physics_interface_nhyd.F<br>
M    src/driver/module_subdriver.F<br>
M    src/core_nhyd_atmos/mpas_interface.F<br>
M    src/core_nhyd_atmos/module_test_cases.F<br>
M    src/core_nhyd_atmos/Registry<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/Registry
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/Registry        2010-10-15 20:55:10 UTC (rev 561)
+++ branches/atmos_physics/src/core_nhyd_atmos/Registry        2010-10-15 21:53:14 UTC (rev 562)
@@ -96,9 +96,9 @@
 var persistent real    areaCell ( nCells ) 0 iro areaCell mesh - -
 var persistent real    areaTriangle ( nVertices ) 0 iro areaTriangle mesh - -
 
-var persistent real    edgeNormalVectors ( R3 nEdges ) 0 o edgeNormalVectors mesh - -
-var persistent real    localVerticalUnitVectors ( R3 nCells ) 0 o localVerticalUnitVectors mesh - -
-var persistent real    cellTangentPlane ( R3 TWO nEdges ) 0 o cellTangentPlane mesh - -
+var persistent real    edgeNormalVectors ( R3 nEdges ) 0 ro edgeNormalVectors mesh - -
+var persistent real    localVerticalUnitVectors ( R3 nCells ) 0 ro localVerticalUnitVectors mesh - -
+var persistent real    cellTangentPlane ( R3 TWO nEdges ) 0 ro cellTangentPlane mesh - -
 
 var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
 var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
@@ -136,14 +136,14 @@
 # coefficients for the vertical tridiagonal solve
 # Note:  these could be local but...
 
-var persistent real    cofrz ( nVertLevels ) 0 - cofrz mesh - -
-var persistent real    cofwr ( nVertLevels nCells ) 0 - cofwr mesh - -
-var persistent real    cofwz ( nVertLevels nCells ) 0 - cofwz mesh - -
-var persistent real    coftz ( nVertLevelsP1 nCells ) 0 - coftz mesh - -
-var persistent real    cofwt ( nVertLevels nCells ) 0 - cofwt mesh - -
-var persistent real    a_tri ( nVertLevels nCells ) 0 - a_tri mesh - -
-var persistent real    alpha_tri ( nVertLevels nCells ) 0 - alpha_tri mesh - -
-var persistent real    gamma_tri ( nVertLevels nCells ) 0 - gamma_tri mesh - -
+var persistent real    cofrz ( nVertLevels Time ) 1 r cofrz diag - -
+var persistent real    cofwr ( nVertLevels nCells Time ) 1 r cofwr diag - -
+var persistent real    cofwz ( nVertLevels nCells Time ) 1 r cofwz diag - -
+var persistent real    coftz ( nVertLevelsP1 nCells Time ) 1 r coftz diag - -
+var persistent real    cofwt ( nVertLevels nCells Time ) 1 r cofwt diag - -
+var persistent real    a_tri ( nVertLevels nCells Time ) 1 r a_tri diag - -
+var persistent real    alpha_tri ( nVertLevels nCells Time ) 1 r alpha_tri diag - -
+var persistent real    gamma_tri ( nVertLevels nCells Time ) 1 r gamma_tri diag - -
 
 #  W-Rayleigh-damping coefficient
 
@@ -153,9 +153,7 @@
 var persistent real    u ( nVertLevels nEdges Time ) 2 iro u state - -
 var persistent real    w ( nVertLevelsP1 nCells Time ) 2 iro w state - -
 var persistent real    rho ( nVertLevels nCells Time ) 2 iro rho state - -
-var persistent real    rho_p ( nVertLevels nCells Time ) 2 iro rho_p state - -
 var persistent real    theta ( nVertLevels nCells Time ) 2 iro theta state - -
-var persistent real    rh ( nVertLevels nCells Time ) 2 iro rh state - -
 var persistent real    qv ( nVertLevels nCells Time ) 2 iro qv state scalars moist
 var persistent real    qc ( nVertLevels nCells Time ) 2 iro qc state scalars moist
 var persistent real    qr ( nVertLevels nCells Time ) 2 iro qr state scalars moist
@@ -166,13 +164,14 @@
 var persistent real    qni ( nVertLevels nCells Time ) 2 iro qni state scalars number
 
 # Tendency variables
-var persistent real    tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
-var persistent real    tend_w ( nVertLevelsP1 nCells Time ) 1 - w tend - -
-var persistent real    tend_rho ( nVertLevels nCells Time ) 1 - rho tend - -
-var persistent real    tend_theta ( nVertLevels nCells Time ) 1 - theta tend - -
-var persistent real    tend_qv ( nVertLevels nCells Time ) 1 - qv tend scalars moist
-var persistent real    tend_qc ( nVertLevels nCells Time ) 1 - qc tend scalars moist
-var persistent real    tend_qr ( nVertLevels nCells Time ) 1 - qr tend scalars moist
+var persistent real    tend_u ( nVertLevels nEdges Time ) 1 r u tend - -
+var persistent real    tend_w ( nVertLevelsP1 nCells Time ) 1 r w tend - -
+var persistent real    tend_rho ( nVertLevels nCells Time ) 1 r rho tend - -
+var persistent real    tend_theta ( nVertLevels nCells Time ) 1 r theta tend - -
+var persistent real    tend_qv ( nVertLevels nCells Time ) 1 r qv tend scalars moist
+var persistent real    tend_qc ( nVertLevels nCells Time ) 1 r qc tend scalars moist
+var persistent real    tend_qr ( nVertLevels nCells Time ) 1 r qr tend scalars moist
+var persistent real    rt_diabatic_tend ( nVertLevels nCells Time ) 1 r rt_diabatic_tend tend - -
 
 # state variables diagnosed from prognostic state
 var persistent real    pressure ( nVertLevels nCells Time ) 2 ro pressure state - -
@@ -182,76 +181,77 @@
 var persistent real    qv_init ( nVertLevels ) 0 iro qv_init mesh - -
 
 # Diagnostic fields: only written to output
-var persistent real    v ( nVertLevels nEdges Time ) 1 o v diag - -
-var persistent real    divergence ( nVertLevels nCells Time ) 1 o divergence diag - -
-var persistent real    vorticity ( nVertLevels nVertices Time ) 1 o vorticity diag - -
-var persistent real    pv_edge ( nVertLevels nEdges Time ) 1 o pv_edge diag - -
-var persistent real    rho_edge ( nVertLevels nEdges Time ) 1 o rho_edge diag - -
-var persistent real    ke ( nVertLevels nCells Time ) 1 o ke diag - -
-var persistent real    pv_vertex ( nVertLevels nVertices Time ) 1 o pv_vertex diag - -
-var persistent real    pv_cell ( nVertLevels nCells Time ) 1 o pv_cell diag - -
-var persistent real    uReconstructX ( nVertLevels nCells Time ) 1 o uReconstructX diag - -
-var persistent real    uReconstructY ( nVertLevels nCells Time ) 1 o uReconstructY diag - -
-var persistent real    uReconstructZ ( nVertLevels nCells Time ) 1 o uReconstructZ diag - -
-var persistent real    uReconstructZonal ( nVertLevels nCells Time ) 1 o uReconstructZonal diag - -
-var persistent real    uReconstructMeridional ( nVertLevels nCells Time ) 1 o uReconstructMeridional diag - -
+var persistent real    rh ( nVertLevels nCells Time ) 1 iro rh diag - -
+var persistent real    v ( nVertLevels nEdges Time ) 1 ro v diag - -
+var persistent real    divergence ( nVertLevels nCells Time ) 1 ro divergence diag - -
+var persistent real    vorticity ( nVertLevels nVertices Time ) 1 ro vorticity diag - -
+var persistent real    pv_edge ( nVertLevels nEdges Time ) 1 ro pv_edge diag - -
+var persistent real    rho_edge ( nVertLevels nEdges Time ) 1 ro rho_edge diag - -
+var persistent real    ke ( nVertLevels nCells Time ) 1 ro ke diag - -
+var persistent real    pv_vertex ( nVertLevels nVertices Time ) 1 ro pv_vertex diag - -
+var persistent real    pv_cell ( nVertLevels nCells Time ) 1 ro pv_cell diag - -
+var persistent real    uReconstructX ( nVertLevels nCells Time ) 1 ro uReconstructX diag - -
+var persistent real    uReconstructY ( nVertLevels nCells Time ) 1 ro uReconstructY diag - -
+var persistent real    uReconstructZ ( nVertLevels nCells Time ) 1 ro uReconstructZ diag - -
+var persistent real    uReconstructZonal ( nVertLevels nCells Time ) 1 ro uReconstructZonal diag - -
+var persistent real    uReconstructMeridional ( nVertLevels nCells Time ) 1 ro uReconstructMeridional diag - -
 
 # Other diagnostic variables: neither read nor written to any files
-var persistent real    rv ( nVertLevels nEdges Time ) 1 - rv diag - -
-var persistent real    circulation ( nVertLevels nVertices Time ) 1 - circulation diag - -
-var persistent real    gradPVt ( nVertLevels nEdges Time ) 1 - gradPVt diag - -
-var persistent real    gradPVn ( nVertLevels nEdges Time ) 1 - gradPVn diag - -
-var persistent real    h_divergence ( nVertLevels nCells ) 0 o h_divergence mesh - -
+var persistent real    rv ( nVertLevels nEdges Time ) 1 r rv diag - -
+var persistent real    circulation ( nVertLevels nVertices Time ) 1 r circulation diag - -
+var persistent real    gradPVt ( nVertLevels nEdges Time ) 1 r gradPVt diag - -
+var persistent real    gradPVn ( nVertLevels nEdges Time ) 1 r gradPVn diag - -
+var persistent real    h_divergence ( nVertLevels nCells Time ) 1 ro h_divergence diag - -
 
-var persistent real    exner ( nVertLevels nCells ) 0 - exner mesh - -
-var persistent real    exner_base ( nVertLevels nCells ) 0 or exner_base mesh - -
-var persistent real    rtheta_base ( nVertLevels nCells ) 0 or rtheta_base mesh - -
-var persistent real    pressure_base ( nVertLevels nCells ) 0 or pressure_base mesh - -
-var persistent real    rho_base ( nVertLevels nCells ) 0 or rho_base mesh - -
-var persistent real    theta_base ( nVertLevels nCells ) 0 or theta_base mesh - -
+var persistent real    exner ( nVertLevels nCells Time ) 1 r exner diag - -
+var persistent real    exner_base ( nVertLevels nCells Time ) 1 or exner_base diag - -
+var persistent real    rtheta_base ( nVertLevels nCells Time ) 1 or rtheta_base diag - -
+var persistent real    pressure_base ( nVertLevels nCells Time ) 1 or pressure_base diag - -
+var persistent real    rho_base ( nVertLevels nCells Time ) 1 or rho_base diag - -
+var persistent real    theta_base ( nVertLevels nCells Time ) 1 or theta_base diag - -
 
 
-var persistent real    ruAvg ( nVertLevels nEdges ) 0 - ruAvg mesh - -
-var persistent real    wwAvg ( nVertLevelsP1 nCells ) 0 - wwAvg mesh - -
-var persistent real    qtot ( nVertLevels nCells ) 0 - qtot mesh - -
-var persistent real    cqu  ( nVertLevels nEdges ) 0 - cqu mesh - -
-var persistent real    cqw  ( nVertLevels nCells ) 0 - cqw mesh - -
-var persistent real    rt_diabatic_tend  ( nVertLevels nCells ) 0 - rt_diabatic_tend mesh - -
+var persistent real    ruAvg ( nVertLevels nEdges Time ) 1 r ruAvg diag - -
+var persistent real    wwAvg ( nVertLevelsP1 nCells Time ) 1 r wwAvg diag - -
+var persistent real    qtot ( nVertLevels nCells Time ) 1 r qtot diag - -
+var persistent real    cqu  ( nVertLevels nEdges Time ) 1 r cqu diag - -
+var persistent real    cqw  ( nVertLevels nCells Time ) 1 r cqw diag - -
 
 #  coupled variables needed by the solver, but not output...
 
-var persistent real    ru ( nVertLevels nEdges ) 0 - ru mesh - -
-var persistent real    ru_p ( nVertLevels nEdges ) 0 - ru_p mesh - -
-var persistent real    ru_save ( nVertLevels nEdges ) 0 - ru_save mesh - -
+var persistent real    ru ( nVertLevels nEdges Time ) 1 r ru diag - -
+var persistent real    ru_p ( nVertLevels nEdges Time ) 1 r ru_p diag - -
+var persistent real    ru_save ( nVertLevels nEdges Time ) 1 r ru_save diag - -
 
 
-var persistent real    rw ( nVertLevelsP1 nCells ) 0 - rw mesh - -
-var persistent real    rw_p ( nVertLevelsP1 nCells ) 0 - rw_p mesh - -
-var persistent real    rw_save ( nVertLevelsP1 nCells ) 0 - rw_save mesh - -
+var persistent real    rw ( nVertLevelsP1 nCells Time ) 1 r rw diag - -
+var persistent real    rw_p ( nVertLevelsP1 nCells Time ) 1 r rw_p diag - -
+var persistent real    rw_save ( nVertLevelsP1 nCells Time ) 1 r rw_save diag - -
 
-var persistent real    rtheta_p ( nVertLevels nCells ) 0 - rtheta_p mesh - -
-var persistent real    rtheta_pp ( nVertLevels nCells ) 0 - rtheta_pp mesh - -
-var persistent real    rtheta_p_save ( nVertLevels nCells ) 0 - rtheta_p_save mesh - -
-var persistent real    rtheta_pp_old ( nVertLevels nCells ) 0 - rtheta_pp_old mesh - -
+var persistent real    rtheta_p ( nVertLevels nCells Time ) 1 r rtheta_p diag - -
+var persistent real    rtheta_pp ( nVertLevels nCells Time ) 1 r rtheta_pp diag - -
+var persistent real    rtheta_p_save ( nVertLevels nCells Time ) 1 r rtheta_p_save diag - -
+var persistent real    rtheta_pp_old ( nVertLevels nCells Time ) 1 r rtheta_pp_old diag - -
 
-var persistent real    rho_pp ( nVertLevels nCells ) 0 - rho_pp mesh - -
-var persistent real    rho_p_save ( nVertLevels nCells ) 0 - rho_p_save mesh - -
+var persistent real    rho_p ( nVertLevels nCells Time ) 1 iro rho_p diag - -
+var persistent real    rho_pp ( nVertLevels nCells Time ) 1 r rho_pp diag - -
+var persistent real    rho_p_save ( nVertLevels nCells Time ) 1 r rho_p_save diag - -
 
-var persistent real    qv_old ( nVertLevels nCells ) 0 - rqv mesh scalars_old moist_old
-var persistent real    qc_old ( nVertLevels nCells ) 0 - rqc mesh scalars_old moist_old
-var persistent real    qr_old ( nVertLevels nCells ) 0 - rqr mesh scalars_old moist_old
+var persistent real    qv_old ( nVertLevels nCells Time ) 1 r rqv diag scalars_old moist_old
+var persistent real    qc_old ( nVertLevels nCells Time ) 1 r rqc diag scalars_old moist_old
+var persistent real    qr_old ( nVertLevels nCells Time ) 1 r rqr diag scalars_old moist_old
 
 # Space needed for advection
-var persistent real    deriv_two ( FIFTEEN TWO nEdges ) 0 o deriv_two mesh - -
-var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
+var persistent real    deriv_two ( FIFTEEN TWO nEdges ) 0 ro deriv_two mesh - -
+var persistent integer advCells ( TWENTYONE nCells ) 0 r advCells mesh - -
 
 # Space needed for deformation calculation weights
-var persistent real    defc_a ( maxEdges nCells ) 0 - defc_a mesh - -
-var persistent real    defc_b ( maxEdges nCells ) 0 - defc_b mesh - -
-var persistent real    kdiff ( nVertLevels nCells Time ) 2 - kdiff diag - -
+var persistent real    defc_a ( maxEdges nCells ) 0 r defc_a mesh - -
+var persistent real    defc_b ( maxEdges nCells ) 0 r defc_b mesh - -
+var persistent real    kdiff ( nVertLevels nCells Time ) 2 r kdiff diag - -
 
 # Arrays required for reconstruction of velocity field
-var persistent real    coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
+var persistent real    coeffs_reconstruct ( R3 maxEdges nCells ) 0 r coeffs_reconstruct mesh - -
 
 #==================================================================================================
 # DECLARATIONS OF ALL PHYSICS VARIABLES (will need to be moved to a Physics Registry shared by the

Modified: branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F        2010-10-15 20:55:10 UTC (rev 561)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F        2010-10-15 21:53:14 UTC (rev 562)
@@ -208,20 +208,20 @@
       hx =&gt; grid % hx % array
       dss =&gt; grid % dss % array
 
-      pb =&gt; grid % exner_base % array
-      rb =&gt; grid % rho_base % array
-      tb =&gt; grid % theta_base % array
-      rtb =&gt; grid % rtheta_base % array
-      p =&gt; grid % exner % array
+      pb =&gt; diag % exner_base % array
+      rb =&gt; diag % rho_base % array
+      tb =&gt; diag % theta_base % array
+      rtb =&gt; diag % rtheta_base % array
+      p =&gt; diag % exner % array
 
-      ppb =&gt; grid % pressure_base % array
+      ppb =&gt; diag % pressure_base % array
       pp =&gt; state % pressure % array
 
       rho =&gt; state % rho % array
-      rr =&gt; state % rho_p % array
+      rr =&gt; diag % rho_p % array
       t =&gt; state % theta % array      
-      rt =&gt; grid % rtheta_p % array
-      rh =&gt; state % rh   % array
+      rt =&gt; diag % rtheta_p % array
+      rh =&gt; diag % rh % array
 
       scalars =&gt; state % scalars % array
 
@@ -669,7 +669,7 @@
          cell2 = grid % CellsOnEdge % array(2,i)
          if(cell1 &lt;= nCellsSolve .or. cell2 &lt;= nCellsSolve) then
            do k=1,nz1
-             grid % ru % array (k,iEdge)  = 0.5*(rho(k,cell1)+rho(k,cell2))*state % u % array (k,iEdge)
+             diag % ru % array (k,iEdge)  = 0.5*(rho(k,cell1)+rho(k,cell2))*state % u % array (k,iEdge)
            end do
          end if
 
@@ -750,7 +750,7 @@
        end do
 
       ! for including terrain
-      grid % rw % array = 0.
+      diag % rw % array = 0.
       state % w % array = 0.
       do iEdge = 1,grid % nEdges
 
@@ -759,15 +759,15 @@
 
          if (cell1 &lt;= nCellsSolve .or. cell2 &lt;= nCellsSolve ) then
          do k = 2, grid%nVertLevels
-            flux =  (fzm(k)*grid % ru % array(k,iEdge)+fzp(k)*grid % ru % array(k-1,iEdge))
-            grid % rw % array(k,cell2) = grid % rw % array(k,cell2) + zf(k,2,iEdge)*flux
-            grid % rw % array(k,cell1) = grid % rw % array(k,cell1) - zf(k,1,iEdge)*flux
+            flux =  (fzm(k)*diag % ru % array(k,iEdge)+fzp(k)*diag % ru % array(k-1,iEdge))
+            diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + zf(k,2,iEdge)*flux
+            diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - zf(k,1,iEdge)*flux
 
             if (config_theta_adv_order ==3) then 
-               grid % rw % array(k,cell2) = grid % rw % array(k,cell2)    &amp;
-                                            - sign(1.,grid % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
-               grid % rw % array(k,cell1) = grid % rw % array(k,cell1)    &amp;
-                                            + sign(1.,grid % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+               diag % rw % array(k,cell2) = diag % rw % array(k,cell2)    &amp;
+                                            - sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+               diag % rw % array(k,cell1) = diag % rw % array(k,cell1)    &amp;
+                                            + sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
             end if
 
          end do
@@ -941,22 +941,22 @@
       hx =&gt; grid % hx % array
       dss =&gt; grid % dss % array
 
-      ppb =&gt; grid % pressure_base % array
-      pb =&gt; grid % exner_base % array
-      rb =&gt; grid % rho_base % array
-      tb =&gt; grid % theta_base % array
-      rtb =&gt; grid % rtheta_base % array
-      p =&gt; grid % exner % array
-      cqw =&gt; grid % cqw % array
+      ppb =&gt; diag % pressure_base % array
+      pb =&gt; diag % exner_base % array
+      rb =&gt; diag % rho_base % array
+      tb =&gt; diag % theta_base % array
+      rtb =&gt; diag % rtheta_base % array
+      p =&gt; diag % exner % array
+      cqw =&gt; diag % cqw % array
 
       rho =&gt; state % rho % array
 
       pp =&gt; state % pressure % array
-      rr =&gt; state % rho_p % array
+      rr =&gt; diag % rho_p % array
       t =&gt; state % theta % array      
-      rt =&gt; grid % rtheta_p % array
+      rt =&gt; diag % rtheta_p % array
       u =&gt; state % u % array
-      ru =&gt; grid % ru % array
+      ru =&gt; diag % ru % array
 
       scalars =&gt; state % scalars % array
 
@@ -1338,7 +1338,7 @@
       !  we are assuming w and rw are zero for this initialization
       !  i.e., no terrain
       !
-       grid % rw % array = 0.
+       diag % rw % array = 0.
        state % w % array = 0.
 
        grid % zf % array = 0.
@@ -1484,22 +1484,22 @@
       xCell =&gt; grid % xCell % array
       yCell =&gt; grid % yCell % array
 
-      ppb =&gt; grid % pressure_base % array
-      pb =&gt; grid % exner_base % array
-      rb =&gt; grid % rho_base % array
-      tb =&gt; grid % theta_base % array
-      rtb =&gt; grid % rtheta_base % array
-      p =&gt; grid % exner % array
-      cqw =&gt; grid % cqw % array
+      ppb =&gt; diag % pressure_base % array
+      pb =&gt; diag % exner_base % array
+      rb =&gt; diag % rho_base % array
+      tb =&gt; diag % theta_base % array
+      rtb =&gt; diag % rtheta_base % array
+      p =&gt; diag % exner % array
+      cqw =&gt; diag % cqw % array
 
       rho =&gt; state % rho % array
 
       pp =&gt; state % pressure % array
-      rr =&gt; state % rho_p % array
+      rr =&gt; diag % rho_p % array
       t =&gt; state % theta % array      
-      rt =&gt; grid % rtheta_p % array
+      rt =&gt; diag % rtheta_p % array
       u =&gt; state % u % array
-      ru =&gt; grid % ru % array
+      ru =&gt; diag % ru % array
 
       scalars =&gt; state % scalars % array
 
@@ -1914,7 +1914,7 @@
 
 !     for including terrain
       state % w % array(:,:) = 0.0
-      grid % rw % array(:,:) = 0.0
+      diag % rw % array(:,:) = 0.0
 
 !
 !     calculation of omega, rw = zx * ru + zz * rw
@@ -1928,13 +1928,13 @@
          if (cell1 &lt;= nCellsSolve .or. cell2 &lt;= nCellsSolve ) then
          do k = 2, grid%nVertLevels
             flux =  (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))  
-            grid % rw % array(k,cell2) = grid % rw % array(k,cell2) + zf(k,2,iEdge)*flux 
-            grid % rw % array(k,cell1) = grid % rw % array(k,cell1) - zf(k,1,iEdge)*flux 
+            diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + zf(k,2,iEdge)*flux 
+            diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - zf(k,1,iEdge)*flux 
 
             if (config_theta_adv_order ==3) then
-               grid % rw % array(k,cell2) = grid % rw % array(k,cell2)    &amp;
+               diag % rw % array(k,cell2) = diag % rw % array(k,cell2)    &amp;
                                             - sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
-               grid % rw % array(k,cell1) = grid % rw % array(k,cell1)    &amp;
+               diag % rw % array(k,cell1) = diag % rw % array(k,cell1)    &amp;
                                             + sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
             end if
 

Modified: branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2010-10-15 20:55:10 UTC (rev 561)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2010-10-15 21:53:14 UTC (rev 562)
@@ -133,13 +133,27 @@
                                           block % mesh % nVertLevels, block % mesh % nEdges, &amp;
                                           block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
 ! rtheta_p
-         call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_p % array(:,:), &amp;
+         call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rtheta_p % array(:,:), &amp;
                                           block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                           block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
 ! ru
-         call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % ru % array(:,:), &amp;
+         call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % ru % array(:,:), &amp;
                                           block % mesh % nVertLevels, block % mesh % nEdges, &amp;
                                           block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+         call dmpar_exch_halo_field2dReal(domain % dminfo, &amp;
+                                          block % diag % cqu % array, &amp;
+                                          block % mesh % nVertLevels, block % mesh % nEdges,  &amp;
+                                          block % parinfo % EdgesToSend, block % parinfo % EdgesToRecv)
+         call dmpar_exch_halo_field2dReal(domain % dminfo, &amp;
+                                          block % diag % cqw % array, &amp;
+                                          block % mesh % nVertLevels, block % mesh % nCells,  &amp;
+                                          block % parinfo % CellsToSend, block % parinfo % CellsToRecv)
+         call dmpar_exch_halo_field2dReal(domain % dminfo, &amp;
+                                          block % diag % rw % array, &amp;
+                                          block % mesh % nVertLevels+1, block % mesh % nCells,  &amp;
+                                          block % parinfo % CellsToSend, block % parinfo % CellsToRecv)
+
          block =&gt; block % next
       end do
 
@@ -147,7 +161,7 @@
       do while (associated(block))
          ! We are setting values in the halo here, so no communications are needed.
          ! Alternatively, we could just set owned cells and edge values and communicate after this block loop.
-         call rk_integration_setup( block % state % time_levs(2) % state, block % state % time_levs(1) % state, block % mesh )
+         call rk_integration_setup( block % state % time_levs(2) % state, block % state % time_levs(1) % state, block % diag )
          block =&gt; block % next
       end do
 
@@ -164,7 +178,7 @@
            ! The coefficients are set for owned cells (cqw) and for all edges of owned cells, 
            ! thus no communications should be needed after this call.  
            ! We could consider combining this and the next block loop.
-           call compute_moist_coefficients( block % state % time_levs(2) % state, block % mesh )
+           call compute_moist_coefficients( block % state % time_levs(2) % state, block % diag, block % mesh )
            block =&gt; block % next
         end do
 
@@ -213,8 +227,8 @@
         block =&gt; domain % blocklist
           do while (associated(block))
             call set_smlstep_pert_variables( block % state % time_levs(1) % state, block % state % time_levs(2) % state,  &amp;
-                                             block % tend, block % mesh               )
-            call compute_vert_imp_coefs( block % state % time_levs(2) % state, block % mesh, rk_sub_timestep(rk_step) )
+                                             block % tend, block % diag, block % mesh                                     )
+            call compute_vert_imp_coefs( block % state % time_levs(2) % state, block % mesh, block % diag, rk_sub_timestep(rk_step) )
             block =&gt; block % next
         end do
 
@@ -224,8 +238,8 @@
       
            block =&gt; domain % blocklist
            do while (associated(block))
-              call advance_acoustic_step( block % state % time_levs(2) % state,  block % tend,  &amp;
-                                          block % mesh, rk_sub_timestep(rk_step)                          )
+              call advance_acoustic_step( block % state % time_levs(2) % state, block % diag, block % tend, &amp;
+                                          block % mesh, rk_sub_timestep(rk_step)                            )
               block =&gt; block % next
            end do
 
@@ -234,13 +248,13 @@
            !  will need communications here for rtheta_pp
             block =&gt; domain % blocklist
             do while (associated(block))
-               call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_pp % array(:,:), &amp;
+               call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rtheta_pp % array(:,:), &amp;
                                                 block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                                 block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-               call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rw_p % array(:,:), &amp;
+               call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rw_p % array(:,:), &amp;
                                                 block % mesh % nVertLevels+1, block % mesh % nCells, &amp;
                                                 block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-               call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % ru_p % array(:,:), &amp;
+               call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % ru_p % array(:,:), &amp;
                                                 block % mesh % nVertLevels, block % mesh % nEdges, &amp;
                                                 block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
                block =&gt; block % next
@@ -251,7 +265,7 @@
         !  will need communications here for rho_pp
          block =&gt; domain % blocklist
          do while (associated(block))
-            call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rho_pp % array(:,:), &amp;
+            call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rho_pp % array(:,:), &amp;
                                              block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                              block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
             block =&gt; block % next
@@ -259,9 +273,9 @@
 
         block =&gt; domain % blocklist
         do while (associated(block))
-           call recover_large_step_variables( block % state % time_levs(2) % state,     &amp;
-                                              block % mesh, rk_timestep(rk_step),       &amp;
-                                              number_sub_steps(rk_step), rk_step  )
+           call recover_large_step_variables( block % state % time_levs(2) % state,                     &amp;
+                                              block % diag, block % tend, block % mesh,                 &amp;
+                                              rk_timestep(rk_step), number_sub_steps(rk_step), rk_step  )
           
            block =&gt; block % next
         end do
@@ -338,7 +352,7 @@
             call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % rho % array(:,:), &amp;
                                              block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                              block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-            call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % rho_p % array(:,:), &amp;
+            call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rho_p % array(:,:), &amp;
                                              block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                              block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
             call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % w % array(:,:), &amp;
@@ -394,7 +408,7 @@
 
          !call microphysics schemes:
          if(config_microp_scheme .ne. 'off') &amp;
-            call microphysics_driver ( block % state % time_levs(2) % state, block % mesh, itimestep )
+            call microphysics_driver ( block % state % time_levs(2) % state, block % diag, block % tend, block % mesh, itimestep )
 
          block =&gt; block % next
       end do
@@ -410,10 +424,10 @@
 
       block =&gt; domain % blocklist
       do while (associated(block))
-         call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_p % array(:,:), &amp;
+         call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rtheta_p % array(:,:), &amp;
                                           block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                           block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-         call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % exner % array(:,:), &amp;
+         call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % exner % array(:,:), &amp;
                                           block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                           block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
          call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % theta % array(:,:), &amp;
@@ -422,7 +436,7 @@
          call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % pressure % array(:,:), &amp;
                                           block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                           block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-         call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rt_diabatic_tend % array(:,:), &amp;
+         call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % rt_diabatic_tend % array(:,:), &amp;
                                           block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                           block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
          call dmpar_exch_halo_field3dReal(domain % dminfo, block % state % time_levs(2) % state % scalars % array(:,:,:), &amp;
@@ -475,22 +489,21 @@
 
 !---
 
-   subroutine rk_integration_setup( s_old, s_new, grid )
+   subroutine rk_integration_setup( s_old, s_new, diag )
 
      implicit none
      type (state_type) :: s_new, s_old
-     type (mesh_type) :: grid
+     type (diag_type) :: diag
      integer :: iCell, k
 
-     grid % ru_save % array = grid % ru % array
-     grid % rw_save % array = grid % rw % array
-     grid % rtheta_p_save % array = grid % rtheta_p % array
-     grid % rho_p_save % array = s_new % rho_p % array
+     diag % ru_save % array = diag % ru % array
+     diag % rw_save % array = diag % rw % array
+     diag % rtheta_p_save % array = diag % rtheta_p % array
+     diag % rho_p_save % array = diag % rho_p % array
 
      s_old % u % array = s_new % u % array
      s_old % w % array = s_new % w % array
      s_old % theta % array = s_new % theta % array
-     s_old % rho_p % array = s_new % rho_p % array
      s_old % rho % array = s_new % rho % array
      s_old % pressure % array = s_new % pressure % array
      s_old % scalars % array = s_new % scalars % array
@@ -499,11 +512,12 @@
 
 !-----
 
-   subroutine compute_moist_coefficients( state, grid )
+   subroutine compute_moist_coefficients( state, diag, grid )
 
-     implicit none
-     type (state_type) :: state
-     type (mesh_type) :: grid
+      implicit none
+      type (state_type) :: state
+      type (diag_type) :: diag
+      type (mesh_type) :: grid
 
       integer :: iEdge, iCell, k, cell1, cell2, iq
       integer :: nCells, nEdges, nVertLevels, nCellsSolve
@@ -520,7 +534,7 @@
             do iq = state % moist_start, state % moist_end
               qtot = qtot + 0.5 * (state % scalars % array (iq, k, iCell) + state % scalars % array (iq, k-1, iCell))
             end do
-            grid % cqw % array(k,iCell) = 1./(1.+qtot)
+            diag % cqw % array(k,iCell) = 1./(1.+qtot)
           end do
         end do
 
@@ -533,7 +547,7 @@
               do iq = state % moist_start, state % moist_end
                  qtot = qtot + 0.5 * ( state % scalars % array (iq, k, cell1) + state % scalars % array (iq, k, cell2) )
               end do
-              grid % cqu % array(k,iEdge) = 1./( 1. + qtot)
+              diag % cqu % array(k,iEdge) = 1./( 1. + qtot)
             end do
           end if
         end do
@@ -542,20 +556,21 @@
 
 !---
 
-   subroutine compute_vert_imp_coefs(s, grid, dts)
+   subroutine compute_vert_imp_coefs(s, grid, diag, dts)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Compute coefficients for vertically implicit gravity-wave/acoustic computations
    !
    ! Input: s - current model state
    !        grid - grid metadata
    !
-   ! Output: tend - cofrz, cofwr, cofwz, coftz, cofwt, a, alpha and gamma
+   ! Output: diag - cofrz, cofwr, cofwz, coftz, cofwt, a, alpha and gamma
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       implicit none
 
       type (state_type), intent(in) :: s
-      type (mesh_type), intent(inout) :: grid
+      type (mesh_type), intent(in) :: grid
+      type (diag_type), intent(inout) :: diag
       real (kind=RKIND), intent(in) :: dts
 
       integer :: i, k, iq
@@ -582,22 +597,22 @@
       fzm =&gt; grid % fzm % array
       fzp =&gt; grid % fzp % array
       zz =&gt; grid % zz % array
-      cqw =&gt; grid % cqw % array
+      cqw =&gt; diag % cqw % array
 
-      p =&gt; grid % exner % array
-      pb =&gt; grid % exner_base % array
-      rt =&gt; grid % rtheta_p % array
-      rtb =&gt; grid % rtheta_base % array
-      rb =&gt; grid % rho_base % array
+      p =&gt; diag % exner % array
+      pb =&gt; diag % exner_base % array
+      rt =&gt; diag % rtheta_p % array
+      rtb =&gt; diag % rtheta_base % array
+      rb =&gt; diag % rho_base % array
 
-      alpha_tri =&gt; grid % alpha_tri % array
-      gamma_tri =&gt; grid % gamma_tri % array
-      a_tri =&gt; grid % a_tri % array
-      cofwr =&gt; grid % cofwr % array      
-      cofwz =&gt; grid % cofwz % array      
-      coftz =&gt; grid % coftz % array      
-      cofwt =&gt; grid % cofwt % array      
-      cofrz =&gt; grid % cofrz % array      
+      alpha_tri =&gt; diag % alpha_tri % array
+      gamma_tri =&gt; diag % gamma_tri % array
+      a_tri =&gt; diag % a_tri % array
+      cofwr =&gt; diag % cofwr % array      
+      cofwz =&gt; diag % cofwz % array      
+      coftz =&gt; diag % coftz % array      
+      cofwt =&gt; diag % cofwt % array      
+      cofrz =&gt; diag % cofrz % array      
 
       t =&gt; s % theta % array
 
@@ -663,11 +678,12 @@
 
 !------------------------
 
-   subroutine set_smlstep_pert_variables( s_old, s_new, tend, grid )
+   subroutine set_smlstep_pert_variables( s_old, s_new, tend, diag, grid )
 
       implicit none
       type (state_type) :: s_new, s_old
       type (tend_type) :: tend
+      type (diag_type) :: diag
       type (mesh_type) :: grid
       integer :: iCell, iEdge, k, cell1, cell2
       integer, dimension(:,:), pointer :: cellsOnEdge
@@ -683,12 +699,12 @@
       areaCell =&gt; grid % areaCell % array
       cellsOnEdge =&gt; grid % cellsOnEdge % array
 
-      grid % rho_pp % array = grid % rho_p_save % array - s_new % rho_p % array
+      diag % rho_pp % array = diag % rho_p_save % array - diag % rho_p % array
 
-      grid % ru_p % array = grid % ru_save % array - grid % ru % array
-      grid % rtheta_pp % array = grid % rtheta_p_save % array - grid % rtheta_p % array
-      grid % rtheta_pp_old % array = grid % rtheta_pp % array
-      grid % rw_p % array = grid % rw_save % array - grid % rw % array
+      diag % ru_p % array = diag % ru_save % array - diag % ru % array
+      diag % rtheta_pp % array = diag % rtheta_p_save % array - diag % rtheta_p % array
+      diag % rtheta_pp_old % array = diag % rtheta_pp % array
+      diag % rw_p % array = diag % rw_save % array - diag % rw % array
 
       do iCell = 1, grid % nCellsSolve
       do k = 2, grid % nVertLevels
@@ -719,18 +735,19 @@
 
       end do
 
-      grid % ruAvg % array = 0.
-      grid % wwAvg % array = 0.
+      diag % ruAvg % array = 0.
+      diag % wwAvg % array = 0.
 
    end subroutine set_smlstep_pert_variables
 
 !-------------------------------
 
-   subroutine advance_acoustic_step( s, tend, grid, dts )
+   subroutine advance_acoustic_step( s, diag, tend, grid, dts )
 
       implicit none
 
       type (state_type) :: s
+      type (diag_type) :: diag
       type (tend_type) :: tend
       type (mesh_type) :: grid
       real (kind=RKIND), intent(in) :: dts
@@ -767,24 +784,24 @@
       theta =&gt; s % theta % array
       w =&gt; s % w % array
 
-      rtheta_pp =&gt; grid % rtheta_pp % array
-      rtheta_pp_old =&gt; grid % rtheta_pp_old % array
-      h_divergence =&gt; grid % h_divergence % array
-      ru_p =&gt; grid % ru_p % array
-      rw_p =&gt; grid % rw_p % array
-      exner =&gt; grid % exner % array
-      cqu =&gt; grid % cqu % array
-      ruAvg =&gt; grid % ruAvg % array
-      wwAvg =&gt; grid % wwAvg % array
-      rho_pp =&gt; grid % rho_pp % array
-      cofwt =&gt; grid % cofwt % array
-      coftz =&gt; grid % coftz % array
-      cofrz =&gt; grid % cofrz % array
-      cofwr =&gt; grid % cofwr % array
-      cofwz =&gt; grid % cofwz % array
-      a_tri =&gt; grid % a_tri % array
-      alpha_tri =&gt; grid % alpha_tri % array
-      gamma_tri =&gt; grid % gamma_tri % array
+      rtheta_pp =&gt; diag % rtheta_pp % array
+      rtheta_pp_old =&gt; diag % rtheta_pp_old % array
+      h_divergence =&gt; diag % h_divergence % array
+      ru_p =&gt; diag % ru_p % array
+      rw_p =&gt; diag % rw_p % array
+      exner =&gt; diag % exner % array
+      cqu =&gt; diag % cqu % array
+      ruAvg =&gt; diag % ruAvg % array
+      wwAvg =&gt; diag % wwAvg % array
+      rho_pp =&gt; diag % rho_pp % array
+      cofwt =&gt; diag % cofwt % array
+      coftz =&gt; diag % coftz % array
+      cofrz =&gt; diag % cofrz % array
+      cofwr =&gt; diag % cofwr % array
+      cofwz =&gt; diag % cofwz % array
+      a_tri =&gt; diag % a_tri % array
+      alpha_tri =&gt; diag % alpha_tri % array
+      gamma_tri =&gt; diag % gamma_tri % array
       dss =&gt; grid % dss % array
 
       tend_ru =&gt; tend % u % array
@@ -941,10 +958,12 @@
 
 !------------------------
 
-      subroutine recover_large_step_variables( s, grid, dt, ns, rk_step )
+      subroutine recover_large_step_variables( s, diag, tend, grid, dt, ns, rk_step )
 
       implicit none
       type (state_type) :: s
+      type (diag_type) :: diag
+      type (tend_type) :: tend
       type (mesh_type) :: grid
       integer, intent(in) :: ns, rk_step
       real (kind=RKIND), intent(in) :: dt
@@ -967,33 +986,33 @@
 
 !---
 
-       wwAvg =&gt; grid % wwAvg % array
-       rw_save =&gt; grid % rw_save % array
-       rw =&gt; grid % rw % array
-       rw_p =&gt; grid % rw_p % array
+       wwAvg =&gt; diag % wwAvg % array
+       rw_save =&gt; diag % rw_save % array
+       rw =&gt; diag % rw % array
+       rw_p =&gt; diag % rw_p % array
        w =&gt; s % w % array
 
-       rtheta_p =&gt; grid % rtheta_p % array
-       rtheta_p_save =&gt; grid % rtheta_p_save % array
-       rtheta_pp =&gt; grid % rtheta_pp % array
-       rtheta_base =&gt; grid % rtheta_base % array
-       rt_diabatic_tend =&gt; grid % rt_diabatic_tend % array
+       rtheta_p =&gt; diag % rtheta_p % array
+       rtheta_p_save =&gt; diag % rtheta_p_save % array
+       rtheta_pp =&gt; diag % rtheta_pp % array
+       rtheta_base =&gt; diag % rtheta_base % array
+       rt_diabatic_tend =&gt; tend % rt_diabatic_tend % array
        theta =&gt; s % theta % array
 
        rho =&gt; s % rho % array
-       rho_p =&gt; s % rho_p % array
-       rho_p_save =&gt; grid % rho_p_save % array
-       rho_pp =&gt; grid % rho_pp % array
-       rho_base =&gt; grid % rho_base % array
+       rho_p =&gt; diag % rho_p % array
+       rho_p_save =&gt; diag % rho_p_save % array
+       rho_pp =&gt; diag % rho_pp % array
+       rho_base =&gt; diag % rho_base % array
 
-       ruAvg =&gt; grid % ruAvg % array
-       ru_save =&gt; grid % ru_save % array
-       ru_p =&gt; grid % ru_p % array
-       ru =&gt; grid % ru % array
+       ruAvg =&gt; diag % ruAvg % array
+       ru_save =&gt; diag % ru_save % array
+       ru_p =&gt; diag % ru_p % array
+       ru =&gt; diag % ru % array
        u =&gt; s % u % array
 
-       exner =&gt; grid % exner % array
-       exner_base =&gt; grid % exner_base % array
+       exner =&gt; diag % exner % array
+       exner_base =&gt; diag % exner_base % array
 
        pressure_p =&gt; s % pressure % array
 
@@ -1105,8 +1124,8 @@
         cell1 = CellsOnEdge(1,iEdge)
         cell2 = CellsOnEdge(2,iEdge)
 
-!        if( cell1 &lt;= nCellsSolve .or. cell2 &lt;= nCellsSolve ) then   ! If using this test, more halo exchanges will be needed
-
+        if( cell1 &lt;= nCellsSolve .or. cell2 &lt;= nCellsSolve ) then   ! If using this test, more halo exchanges will be needed
+                                                                    ! Keep for now to avoid division by zero, e.g., rho(:,cell2)
           do k = 1, nVertLevels
             ruAvg(k,iEdge) = ru(k,iEdge) + (ruAvg(k,iEdge) / float(ns))
 
@@ -1140,7 +1159,7 @@
             end if
           enddo
 
-!        end if
+        end if
 
       enddo
 
@@ -1204,7 +1223,7 @@
       kdiff       =&gt; diag % kdiff % array
       deriv_two   =&gt; grid % deriv_two % array
 !****      uhAvg       =&gt; grid % uhAvg % array
-      uhAvg       =&gt; grid % ruAvg % array
+      uhAvg       =&gt; diag % ruAvg % array
       dvEdge      =&gt; grid % dvEdge % array
       dcEdge      =&gt; grid % dcEdge % array
       cellsOnEdge =&gt; grid % cellsOnEdge % array
@@ -1213,7 +1232,7 @@
 !****      h_new       =&gt; s_new % h % array
       h_old       =&gt; s_old % rho % array
       h_new       =&gt; s_new % rho % array
-      wwAvg       =&gt; grid % wwAvg % array
+      wwAvg       =&gt; diag % wwAvg % array
       areaCell    =&gt; grid % areaCell % array
 
 !****      fnm         =&gt; grid % fnm % array
@@ -1245,7 +1264,7 @@
          do iEdge=1,grid%nEdges
             cell1 = cellsOnEdge(1,iEdge)
             cell2 = cellsOnEdge(2,iEdge)
-            if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+            if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
                do k=1,grid % nVertLevels
                   do iScalar=1,s_old % num_scalars
                      scalar_edge = 0.5 * (scalar_new(iScalar,k,cell1) + scalar_new(iScalar,k,cell2))
@@ -1262,7 +1281,7 @@
          do iEdge=1,grid%nEdges
             cell1 = cellsOnEdge(1,iEdge)
             cell2 = cellsOnEdge(2,iEdge)
-            if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+            if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
   
                do k=1,grid % nVertLevels
    
@@ -1270,12 +1289,12 @@
                      d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iScalar,k,cell1)
                      d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iScalar,k,cell2)
                      do i=1, grid % nEdgesOnCell % array (cell1)
-                        if ( grid % CellsOnCell % array (i,cell1) &gt; 0) &amp;
+                        if ( grid % CellsOnCell % array (i,cell1) &lt;= grid%nCells) &amp;
                         d2fdx2_cell1 = d2fdx2_cell1 + &amp;
                                        deriv_two(i+1,1,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell1))
                      end do
                      do i=1, grid % nEdgesOnCell % array (cell2)
-                        if ( grid % CellsOnCell % array (i,cell2) &gt; 0) &amp;
+                        if ( grid % CellsOnCell % array (i,cell2) &lt;= grid%nCells) &amp;
                         d2fdx2_cell2 = d2fdx2_cell2 + &amp;
                                        deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
                      end do
@@ -1305,7 +1324,7 @@
          do iEdge=1,grid%nEdges
             cell1 = cellsOnEdge(1,iEdge)
             cell2 = cellsOnEdge(2,iEdge)
-            if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+            if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
 
                do k=1,grid % nVertLevels
    
@@ -1313,12 +1332,12 @@
                      d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iScalar,k,cell1)
                      d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iScalar,k,cell2)
                      do i=1, grid % nEdgesOnCell % array (cell1)
-                        if ( grid % CellsOnCell % array (i,cell1) &gt; 0) &amp;
+                        if ( grid % CellsOnCell % array (i,cell1) &lt;= grid%nCells) &amp;
                            d2fdx2_cell1 = d2fdx2_cell1 + &amp;
                                           deriv_two(i+1,1,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell1))
                      end do
                      do i=1, grid % nEdgesOnCell % array (cell2)
-                        if ( grid % CellsOnCell % array (i,cell2) &gt; 0) &amp;
+                        if ( grid % CellsOnCell % array (i,cell2) &lt;= grid%nCells) &amp;
                         d2fdx2_cell2 = d2fdx2_cell2 + &amp;
                                        deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
                      end do
@@ -1343,7 +1362,7 @@
          do iEdge=1,grid % nEdges
             cell1 = grid % cellsOnEdge % array(1,iEdge)
             cell2 = grid % cellsOnEdge % array(2,iEdge)
-            if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+            if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
 
                do k=1,grid % nVertLevels
                   do iScalar=1,s_old % num_scalars
@@ -1363,7 +1382,7 @@
          do iEdge=1,grid % nEdges
             cell1 = grid % cellsOnEdge % array(1,iEdge)
             cell2 = grid % cellsOnEdge % array(2,iEdge)
-            if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+            if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
 
                do k=1,grid % nVertLevels
                   do iScalar=1,s_old % num_scalars
@@ -1536,7 +1555,7 @@
       scalar_new  =&gt; s_new % scalars % array
       deriv_two   =&gt; grid % deriv_two % array
 !****      uhAvg       =&gt; grid % uhAvg % array
-      uhAvg       =&gt; grid % ruAvg % array
+      uhAvg       =&gt; diag % ruAvg % array
       dvEdge      =&gt; grid % dvEdge % array
       dcEdge      =&gt; grid % dcEdge % array
       cellsOnEdge =&gt; grid % cellsOnEdge % array
@@ -1545,7 +1564,7 @@
 !****      h_new       =&gt; s_new % h % array
       h_old       =&gt; s_old % rho % array
       h_new       =&gt; s_new % rho % array
-      wwAvg       =&gt; grid % wwAvg % array
+      wwAvg       =&gt; diag % wwAvg % array
       areaCell    =&gt; grid % areaCell % array
 
 !****      fnm         =&gt; grid % fnm % array
@@ -1586,7 +1605,7 @@
 
             if (k &lt; grid % nVertLevels) then
 
-               if ((config_scalar_vadv_order == 2) .or. (k==1) .or. (k==grid % nVertLevels)) then
+               if ((config_scalar_vadv_order == 2) .or. (k==1) .or. (k==grid % nVertLevels-1)) then
                   do iScalar=1,s_old % num_scalars
                      v_flux(iScalar,iCell,km0) = dt * wwAvg(k+1,iCell) *   &amp;
                           (fnm(k+1) * scalar_new(iScalar,k+1,iCell) + fnp(k+1) * scalar_new(iScalar,k,iCell))
@@ -1634,7 +1653,7 @@
             do iEdge=1,grid%nEdges
                cell1 = cellsOnEdge(1,iEdge)
                cell2 = cellsOnEdge(2,iEdge)
-               if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+               if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
                   cell_upwind = cell2
                   if (uhAvg(k,iEdge) &gt;= 0) cell_upwind = cell1
                   do iScalar=1,s_old % num_scalars
@@ -1654,7 +1673,7 @@
             do iEdge=1,grid%nEdges
                cell1 = cellsOnEdge(1,iEdge)
                cell2 = cellsOnEdge(2,iEdge)
-               if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+               if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
                   cell_upwind = cell2
                   if (uhAvg(k,iEdge) &gt;= 0) cell_upwind = cell1
                   do iScalar=1,s_old % num_scalars
@@ -1662,12 +1681,12 @@
                      d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iScalar,k,cell1)
                      d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iScalar,k,cell2)
                      do i=1, grid % nEdgesOnCell % array (cell1)
-                        if ( grid % CellsOnCell % array (i,cell1) &gt; 0) &amp;
+                        if ( grid % CellsOnCell % array (i,cell1) &lt;= grid%nCells) &amp;
                         d2fdx2_cell1 = d2fdx2_cell1 + &amp;
                                        deriv_two(i+1,1,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell1))
                      end do
                      do i=1, grid % nEdgesOnCell % array (cell2)
-                        if ( grid % CellsOnCell % array (i,cell2) &gt; 0) &amp;
+                        if ( grid % CellsOnCell % array (i,cell2) &lt;= grid%nCells) &amp;
                         d2fdx2_cell2 = d2fdx2_cell2 + &amp;
                                        deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
                      end do
@@ -1721,7 +1740,7 @@
                end do
    
                do i = 1, grid % nEdgesOnCell % array(iCell)  ! go around the edges of each cell
-                  if (grid % cellsOnCell % array(i,iCell) &gt; 0) then
+                  if (grid % cellsOnCell % array(i,iCell) &lt;= grid%nCells) then
                      do iScalar=1,s_old % num_scalars
     
                         s_max(iScalar)  = max(scalar_old(iScalar,k,grid % cellsOnCell % array(i,iCell)), s_max(iScalar))
@@ -1777,7 +1796,7 @@
             do iEdge = 1, grid % nEdges
                cell1 = grid % cellsOnEdge % array(1,iEdge)
                cell2 = grid % cellsOnEdge % array(2,iEdge)
-               if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+               if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
                   do iScalar=1,s_old % num_scalars
                      flux = h_flux(iScalar,iEdge)
                      if (flux &gt; 0) then
@@ -1822,7 +1841,7 @@
          do iEdge=1,grid%nEdges
             cell1 = cellsOnEdge(1,iEdge)
             cell2 = cellsOnEdge(2,iEdge)
-            if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+            if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
                do iScalar=1,s_old % num_scalars
                   s_update(iScalar,cell1,km0) = s_update(iScalar,cell1,km0) - &amp;
                       h_flux(iScalar,iEdge) / grid % areaCell % array(cell1)
@@ -1939,14 +1958,14 @@
 
       rho          =&gt; s % rho % array
       rho_edge     =&gt; diag % rho_edge % array
-      rb           =&gt; grid % rho_base % array
-      rr           =&gt; s % rho_p % array
+      rb           =&gt; diag % rho_base % array
+      rr           =&gt; diag % rho_p % array
       u            =&gt; s % u % array
       v            =&gt; diag % v % array
       kdiff        =&gt; diag % kdiff % array
-      ru           =&gt; grid % ru % array
+      ru           =&gt; diag % ru % array
       w            =&gt; s % w % array
-      rw           =&gt; grid % rw % array
+      rw           =&gt; diag % rw % array
       theta        =&gt; s % theta % array
       circulation  =&gt; diag % circulation % array
       divergence   =&gt; diag % divergence % array
@@ -1954,8 +1973,8 @@
       ke           =&gt; diag % ke % array
       pv_edge      =&gt; diag % pv_edge % array
       pp           =&gt; s % pressure % array
-      pressure_b   =&gt; grid % pressure_base % array
-      h_divergence =&gt; grid % h_divergence % array
+      pressure_b   =&gt; diag % pressure_base % array
+      h_divergence =&gt; diag % h_divergence % array
 
 
       weightsOnEdge     =&gt; grid % weightsOnEdge % array
@@ -1981,7 +2000,7 @@
       tend_theta  =&gt; tend % theta % array
       tend_w      =&gt; tend % w % array
       tend_rho    =&gt; tend % rho % array
-      rt_diabatic_tend  =&gt; grid % rt_diabatic_tend % array
+      rt_diabatic_tend  =&gt; tend % rt_diabatic_tend % array
 
       t_init      =&gt; grid % t_init % array
       qv_init     =&gt; grid % qv_init % array
@@ -1991,8 +2010,8 @@
       fzm         =&gt; grid % fzm % array
       fzp         =&gt; grid % fzp % array
       zgrid       =&gt; grid % zgrid % array
-      cqw         =&gt; grid % cqw % array
-      cqu         =&gt; grid % cqu % array
+      cqw         =&gt; diag % cqw % array
+      cqu         =&gt; diag % cqu % array
 
       nCells      = grid % nCells
       nEdges      = grid % nEdges
@@ -2285,12 +2304,12 @@
 
          delsq_circulation(:,:) = 0.0
          do iEdge=1,nEdges
-            if (verticesOnEdge(1,iEdge) &gt; 0) then
+            if (verticesOnEdge(1,iEdge) &lt;= grid%nVertices) then
                do k=1,nVertLevels
                   delsq_circulation(k,verticesOnEdge(1,iEdge)) = delsq_circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * delsq_u(k,iEdge)
                end do
             end if
-            if (verticesOnEdge(2,iEdge) &gt; 0) then
+            if (verticesOnEdge(2,iEdge) &lt;= grid%nVertices) then
                do k=1,nVertLevels
                   delsq_circulation(k,verticesOnEdge(2,iEdge)) = delsq_circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * delsq_u(k,iEdge)
                end do
@@ -2447,11 +2466,11 @@
                   d2fdx2_cell1 = deriv_two(1,1,iEdge) * w(k,cell1)
                   d2fdx2_cell2 = deriv_two(1,2,iEdge) * w(k,cell2)
                   do i=1, grid % nEdgesOnCell % array (cell1)
-                     if ( grid % CellsOnCell % array (i,cell1) &gt; 0) &amp;
+                     if ( grid % CellsOnCell % array (i,cell1) &lt;= grid%nCells) &amp;
                      d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * w(k,grid % CellsOnCell % array (i,cell1))
                   end do
                   do i=1, grid % nEdgesOnCell % array (cell2)
-                     if ( grid % CellsOnCell % array (i,cell2) &gt; 0) &amp;
+                     if ( grid % CellsOnCell % array (i,cell2) &lt;= grid%nCells) &amp;
                      d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * w(k,grid % CellsOnCell % array (i,cell2))
                   end do
 
@@ -2485,11 +2504,11 @@
                   d2fdx2_cell1 = deriv_two(1,1,iEdge) * w(k,cell1)
                   d2fdx2_cell2 = deriv_two(1,2,iEdge) * w(k,cell2)
                   do i=1, grid % nEdgesOnCell % array (cell1)
-                     if ( grid % CellsOnCell % array (i,cell1) &gt; 0) &amp;
+                     if ( grid % CellsOnCell % array (i,cell1) &lt;= grid%nCells) &amp;
                      d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * w(k,grid % CellsOnCell % array (i,cell1))
                   end do
                   do i=1, grid % nEdgesOnCell % array (cell2)
-                     if ( grid % CellsOnCell % array (i,cell2) &gt; 0) &amp;
+                     if ( grid % CellsOnCell % array (i,cell2) &lt;= grid%nCells) &amp;
                      d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * w(k,grid % CellsOnCell % array (i,cell2))
                   end do
 
@@ -2720,11 +2739,11 @@
                   d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
                   d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(k,cell2)
                   do i=1, grid % nEdgesOnCell % array (cell1)
-                     if ( grid % CellsOnCell % array (i,cell1) &gt; 0) &amp;
+                     if ( grid % CellsOnCell % array (i,cell1) &lt;= grid%nCells) &amp;
                      d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta(k,grid % CellsOnCell % array (i,cell1))
                   end do
                   do i=1, grid % nEdgesOnCell % array (cell2)
-                     if ( grid % CellsOnCell % array (i,cell2) &gt; 0) &amp;
+                     if ( grid % CellsOnCell % array (i,cell2) &lt;= grid%nCells) &amp;
                      d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
                   end do
 
@@ -2761,18 +2780,18 @@
          do iEdge=1,nEdges
             cell1 = cellsOnEdge(1,iEdge)
             cell2 = cellsOnEdge(2,iEdge)
-            if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+            if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
 
                do k=1,grid % nVertLevels
 
                   d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
                   d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(k,cell2)
                   do i=1, grid % nEdgesOnCell % array (cell1)
-                     if ( grid % CellsOnCell % array (i,cell1) &gt; 0) &amp;
+                     if ( grid % CellsOnCell % array (i,cell1) &lt;= grid%nCells) &amp;
                      d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta(k,grid % CellsOnCell % array (i,cell1))
                   end do
                   do i=1, grid % nEdgesOnCell % array (cell2)
-                     if ( grid % CellsOnCell % array (i,cell2) &gt; 0) &amp;
+                     if ( grid % CellsOnCell % array (i,cell2) &lt;= grid%nCells) &amp;
                      d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
                   end do
 
@@ -3048,7 +3067,7 @@
       do iEdge=1,nEdges
          cell1 = cellsOnEdge(1,iEdge)
          cell2 = cellsOnEdge(2,iEdge)
-         if (cell1 &gt; 0 .and. cell2 &gt; 0) then
+         if (cell1 &lt;= grid%nCells .and. cell2 &lt;= grid%nCells) then
             do k=1,nVertLevels
                h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
             end do
@@ -3062,12 +3081,12 @@
       !
       circulation(:,:) = 0.0
       do iEdge=1,nEdges
-         if (verticesOnEdge(1,iEdge) &gt; 0) then
+         if (verticesOnEdge(1,iEdge) &lt;= grid%nVertices) then
             do k=1,nVertLevels
                circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
             end do
          end if
-         if (verticesOnEdge(2,iEdge) &gt; 0) then
+         if (verticesOnEdge(2,iEdge) &lt;= grid%nVertices) then
             do k=1,nVertLevels
                circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
             end do
@@ -3087,12 +3106,12 @@
       do iEdge=1,nEdges
          cell1 = cellsOnEdge(1,iEdge)
          cell2 = cellsOnEdge(2,iEdge)
-         if (cell1 &gt; 0) then
+         if (cell1 &lt;= grid%nCells) then
             do k=1,nVertLevels
               divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
             end do
          end if
-         if(cell2 &gt; 0) then
+         if(cell2 &lt;= grid%nCells) then
             do k=1,nVertLevels
               divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
             end do
@@ -3130,7 +3149,7 @@
       do iEdge = 1,nEdges
          do i=1,nEdgesOnEdge(iEdge)
             eoe = edgesOnEdge(i,iEdge)
-            if (eoe &gt; 0) then
+            if (eoe &lt;= grid%nEdges) then
                do k = 1,nVertLevels
                  v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
               end do
@@ -3182,7 +3201,7 @@
       do iVertex = 1,nVertices
         do i=1,grid % vertexDegree
           iEdge = edgesOnVertex(i,iVertex)
-          if(iEdge &gt; 0) then
+          if(iEdge &lt;= grid%nEdges) then
             do k=1,nVertLevels
               pv_edge(k,iEdge) =  pv_edge(k,iEdge)  + 0.5 * pv_vertex(k,iVertex)
             end do
@@ -3211,7 +3230,7 @@
       do iVertex = 1, nVertices
        do i=1,grid % vertexDegree
          iCell = cellsOnVertex(i,iVertex)
-         if( iCell &gt; 0) then
+         if( iCell &lt;= grid % nCells) then
            do k = 1,nVertLevels
              pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
            end do
@@ -3227,7 +3246,7 @@
       !
       gradPVn(:,:) = 0.0
       do iEdge = 1,nEdges
-        if( cellsOnEdge(1,iEdge) &gt; 0 .and. cellsOnEdge(2,iEdge) &gt; 0) then
+        if( cellsOnEdge(1,iEdge) &lt;= grid%nCells .and. cellsOnEdge(2,iEdge) &lt;= grid%nCells) then
           do k = 1,nVertLevels
             gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &amp;
                                  dcEdge(iEdge)
@@ -3249,11 +3268,12 @@
 
 !----------
 
-   subroutine init_coupled_diagnostics( state, grid )
+   subroutine init_coupled_diagnostics( state, diag, grid )
 
    implicit none
    
    type (state_type), intent(inout) :: state
+   type (diag_type), intent(inout) :: diag
    type (mesh_type), intent(inout) :: grid
 
    integer :: k,iEdge,i,iCell1,iCell2
@@ -3262,7 +3282,7 @@
         iCell1 = grid % cellsOnEdge % array(1,iEdge)
         iCell2 = grid % cellsOnEdge % array(2,iEdge)
         do k=1,grid % nVertLevels
-          grid % ru % array(k,iEdge) = 0.5 * state % u % array(k,iEdge)*(state % rho % array(k,iCell1)+state % rho % array(k,iCell2))
+          diag % ru % array(k,iEdge) = 0.5 * state % u % array(k,iEdge)*(state % rho % array(k,iCell1)+state % rho % array(k,iCell2))
         enddo
       enddo
 
@@ -3270,11 +3290,13 @@
 
 ! ------------------------
 
-   subroutine qd_kessler( state_old, state_new, grid, dt )
+   subroutine qd_kessler( state_old, state_new, diag, tend, grid, dt )
 
       implicit none
       
       type (state_type), intent(inout) :: state_old, state_new
+      type (diag_type), intent(inout) :: diag
+      type (tend_type), intent(inout) :: tend
       type (mesh_type), intent(inout) :: grid
       real (kind=RKIND), intent(in) :: dt
    
@@ -3294,11 +3316,11 @@
    
         do k = 1, grid % nVertLevels
    
-          grid % rt_diabatic_tend % array(k,iCell) = state_new % theta % array(k,iCell)
+          tend % rt_diabatic_tend % array(k,iCell) = state_new % theta % array(k,iCell)
    
           t(k) = state_new % theta % array(k,iCell)/(1. + 1.61*state_new % scalars % array(state_new % index_qv,k,iCell))
           rho(k) = grid % zz % array(k,iCell)*state_new % rho % array(k,iCell)
-          p(k) = grid % exner % array(k,iCell)
+          p(k) = diag % exner % array(k,iCell)
           qv(k) = max(0.,state_new % scalars % array(state_new % index_qv,k,iCell))
           qc(k) = max(0.,state_new % scalars % array(state_new % index_qc,k,iCell))
           qr(k) = max(0.,state_new % scalars % array(state_new % index_qr,k,iCell))
@@ -3313,24 +3335,24 @@
         do k = 1, grid % nVertLevels
    
           state_new % theta % array(k,iCell) = t(k)*(1.+1.61*qv(k))
-          grid % rt_diabatic_tend % array(k,iCell) = state_new % rho % array(k,iCell) *  &amp;
-                     (state_new % theta % array(k,iCell) - grid % rt_diabatic_tend % array(k,iCell))/dt
-          grid % rtheta_p % array(k,iCell) = state_new % rho % array(k,iCell) * state_new % theta % array(k,iCell)  &amp;
-                                         - grid % rtheta_base % array(k,iCell) 
+          tend % rt_diabatic_tend % array(k,iCell) = state_new % rho % array(k,iCell) *  &amp;
+                     (state_new % theta % array(k,iCell) - tend % rt_diabatic_tend % array(k,iCell))/dt
+          diag % rtheta_p % array(k,iCell) = state_new % rho % array(k,iCell) * state_new % theta % array(k,iCell)  &amp;
+                                         - diag % rtheta_base % array(k,iCell) 
           state_new % scalars % array(state_new % index_qv,k,iCell) = qv(k)
           state_new % scalars % array(state_new % index_qc,k,iCell) = qc(k)
           state_new % scalars % array(state_new % index_qr,k,iCell) = qr(k)
    
-          grid % exner % array(k,iCell) =                                       &amp;
+          diag % exner % array(k,iCell) =                                       &amp;
                                  ( grid % zz % array(k,iCell)*(rgas/p0) * ( &amp;
-                                     grid % rtheta_p % array(k,iCell)       &amp;
-                                   + grid % rtheta_base % array(k,iCell) ) )**rcv
+                                     diag % rtheta_p % array(k,iCell)       &amp;
+                                   + diag % rtheta_base % array(k,iCell) ) )**rcv
    
           state_new % pressure % array(k,iCell) =                                               &amp;
                grid % zz % array(k,iCell) * rgas * (                                        &amp;
-                 grid % exner % array(k,iCell)*grid % rtheta_p % array(k,iCell)             &amp;
-                                   +grid % rtheta_base % array(k,iCell) *                   &amp;
-                        (grid % exner % array(k,iCell) - grid % exner_base % array(k,iCell)) )
+                 diag % exner % array(k,iCell)*diag % rtheta_p % array(k,iCell)             &amp;
+                                   +diag % rtheta_base % array(k,iCell) *                   &amp;
+                        (diag % exner % array(k,iCell) - diag % exner_base % array(k,iCell)) )
         end do
    
       end do

Modified: branches/atmos_physics/src/core_nhyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/mpas_interface.F        2010-10-15 20:55:10 UTC (rev 561)
+++ branches/atmos_physics/src/core_nhyd_atmos/mpas_interface.F        2010-10-15 21:53:14 UTC (rev 562)
@@ -36,7 +36,7 @@
    type (mesh_type), intent(inout) :: mesh
    real (kind=RKIND), intent(in) :: dt
 
-   call init_coupled_diagnostics( block % state % time_levs(1) % state, mesh)
+   call init_coupled_diagnostics( block % state % time_levs(1) % state, block % diag, mesh)
    call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, block % diag, mesh)
 
    call rbfInterp_initialize(mesh)

Modified: branches/atmos_physics/src/core_physics/module_driver_microphysics.F
===================================================================
--- branches/atmos_physics/src/core_physics/module_driver_microphysics.F        2010-10-15 20:55:10 UTC (rev 561)
+++ branches/atmos_physics/src/core_physics/module_driver_microphysics.F        2010-10-15 21:53:14 UTC (rev 562)
@@ -116,7 +116,11 @@
  end subroutine microphysics_init
 
 !=============================================================================================
+#ifdef non_hydrostatic_core
+ subroutine microphysics_driver(state,diag,tend,grid,itimestep)
+#elif hydrostatic_core
  subroutine microphysics_driver(state,grid,itimestep)
+#endif
 !=============================================================================================
 
 !input arguments:
@@ -127,6 +131,10 @@
 !inout arguments:
 !----------------
  type(state_type),intent(inout):: state
+#ifdef non_hydrostatic_core
+ type(diag_type),intent(inout):: diag
+ type(tend_type),intent(inout):: tend
+#endif
 
 !local variables and arrays:
 !---------------------------
@@ -160,7 +168,7 @@
 
 #ifdef non_hydrostatic_core
 
- call microphysics_from_MPAS(state,grid)
+ call microphysics_from_MPAS(state,diag,grid)
 
 #elif hydrostatic_core
 
@@ -222,7 +230,7 @@
 
 #ifdef non_hydrostatic_core
 
- call microphysics_to_MPAS(state,grid,itimestep)
+ call microphysics_to_MPAS(state,diag,tend,grid,itimestep)
 
 #elif hydrostatic_core
 

Modified: branches/atmos_physics/src/core_physics/module_physics_driver.F
===================================================================
--- branches/atmos_physics/src/core_physics/module_physics_driver.F        2010-10-15 20:55:10 UTC (rev 561)
+++ branches/atmos_physics/src/core_physics/module_physics_driver.F        2010-10-15 21:53:14 UTC (rev 562)
@@ -49,7 +49,11 @@
  do while(associated(block))
 
     !compute relative humidity:
-    !call compute_relhum(block%mesh,block%state%time_levs(1)%state)
+#ifdef non_hydrostatic_core
+    !call compute_relhum(block%mesh, block%diag, block%state%time_levs(1)%state)
+#elif hydrostatic_core
+    !call compute_relhum(block%mesh, block%state%time_levs(1)%state)
+#endif
 
     !physics prep step:
 #ifdef non_hydrostatic_core

Modified: branches/atmos_physics/src/core_physics/module_physics_interface_nhyd.F
===================================================================
--- branches/atmos_physics/src/core_physics/module_physics_interface_nhyd.F        2010-10-15 20:55:10 UTC (rev 561)
+++ branches/atmos_physics/src/core_physics/module_physics_interface_nhyd.F        2010-10-15 21:53:14 UTC (rev 562)
@@ -108,10 +108,10 @@
 !initialization:
  zgrid      =&gt; grid% zgrid % array
  zz         =&gt; grid % zz % array
- exner      =&gt; grid % exner % array
- pressure_b =&gt; grid % pressure_base % array
- rtheta_p   =&gt; grid % rtheta_p % array
- rtheta_b   =&gt; grid % rtheta_base % array
+ exner      =&gt; diag % exner % array
+ pressure_b =&gt; diag % pressure_base % array
+ rtheta_p   =&gt; diag % rtheta_p % array
+ rtheta_b   =&gt; diag % rtheta_base % array
 
  rho        =&gt; vars % rho % array
  theta      =&gt; vars % theta % array
@@ -153,12 +153,13 @@
  end subroutine MPAS_to_physics
 
 !=============================================================================================
- subroutine microphysics_from_MPAS(vars,grid)
+ subroutine microphysics_from_MPAS(vars,diag,grid)
 !=============================================================================================
 
 !input variables:
+ type(state_type),intent(in):: vars
+ type(diag_type) ,intent(in):: diag
  type(mesh_type) ,intent(in):: grid
- type(state_type),intent(in):: vars
 
 !local variables:
  integer:: i,k,j
@@ -181,14 +182,14 @@
  zgrid =&gt; grid% zgrid % array
 
  zz         =&gt; grid % zz % array
- exner      =&gt; grid % exner % array
- pressure_b =&gt; grid % pressure_base % array
- rtheta_p   =&gt; grid % rtheta_p % array
- rtheta_b   =&gt; grid % rtheta_base % array
+ exner      =&gt; diag % exner % array
+ pressure_b =&gt; diag % pressure_base % array
+ rtheta_p   =&gt; diag % rtheta_p % array
+ rtheta_b   =&gt; diag % rtheta_base % array
  
  rho        =&gt; vars % rho % array
  theta      =&gt; vars % theta % array
- rh         =&gt; vars % rh % array
+ rh         =&gt; diag % rh % array
  pressure   =&gt; vars % pressure % array
 
  qv =&gt; vars % scalars % array(vars%index_qv,:,:)
@@ -253,7 +254,7 @@
  end subroutine microphysics_from_MPAS
 
 !=============================================================================================
- subroutine microphysics_to_MPAS(vars, grid, itimestep)
+ subroutine microphysics_to_MPAS(vars, diag, tend, grid, itimestep)
 !=============================================================================================
 
 !input variables:
@@ -262,6 +263,8 @@
 
 !output variables:
  type(state_type),intent(inout):: vars
+ type(diag_type),intent(inout):: diag
+ type(tend_type),intent(inout):: tend
 
  real(kind=RKIND):: min_theta,min_thp,min_tp
  real(kind=RKIND):: max_theta,max_thp,max_tp
@@ -282,17 +285,17 @@
 
 !initialization:
  zz       =&gt; grid % zz % array
- exner    =&gt; grid % exner % array
- exner_b  =&gt; grid % exner_base % array
- pressure_b =&gt; grid % pressure_base % array
- rtheta_p =&gt; grid % rtheta_p % array
- rtheta_b =&gt; grid % rtheta_base % array
+ exner    =&gt; diag % exner % array
+ exner_b  =&gt; diag % exner_base % array
+ pressure_b =&gt; diag % pressure_base % array
+ rtheta_p =&gt; diag % rtheta_p % array
+ rtheta_b =&gt; diag % rtheta_base % array
 
  rho      =&gt; vars % rho % array
  theta    =&gt; vars % theta % array
  pressure =&gt; vars % pressure % array
 
- rt_diabatic_tend =&gt; grid % rt_diabatic_tend % array
+ rt_diabatic_tend =&gt; tend % rt_diabatic_tend % array
 
 !variables common to all cloud microphysics schemes:
 
@@ -358,10 +361,11 @@
  end subroutine microphysics_to_MPAS
 
 !--------------------------------------------------------------------------------------------------
- subroutine compute_relhum(grid,vars)
+ subroutine compute_relhum(grid,diag,vars)
 !--------------------------------------------------------------------------------------------------
    
  type(state_type),intent(in) :: vars
+ type(diag_type),intent(inout) :: diag
  type(mesh_type),intent(in)  :: grid
 
  real,dimension(:,:),pointer :: theta,pressure  
@@ -371,12 +375,12 @@
  integer:: iCell,k
  real(kind=RKIND):: tt,pres,qsat
 
- exner      =&gt; grid % exner % array
- pressure_b =&gt; grid % pressure_base % array
+ exner      =&gt; diag % exner % array
+ pressure_b =&gt; diag % pressure_base % array
 
  theta    =&gt; vars % theta % array
  pressure =&gt; vars % pressure % array
- rh       =&gt; vars % rh  % array
+ rh       =&gt; diag % rh  % array
  qv       =&gt; vars % scalars % array(vars%index_qv,:,:)
 
  do iCell = 1, grid % nCells

Modified: branches/atmos_physics/src/driver/module_subdriver.F
===================================================================
--- branches/atmos_physics/src/driver/module_subdriver.F        2010-10-15 20:55:10 UTC (rev 561)
+++ branches/atmos_physics/src/driver/module_subdriver.F        2010-10-15 21:53:14 UTC (rev 562)
@@ -56,12 +56,14 @@
 
       ! Compute diagnostic fields needed in solve loop, and initialize 
       !    simulation time to 0 for all blocks
-      block_ptr =&gt; domain % blocklist
-      do while (associated(block_ptr))
-         call mpas_init(block_ptr, block_ptr % mesh, dt)
-         if (.not. config_do_restart) block_ptr % state % time_levs(1) % state % xtime % scalar = 0.0
-         block_ptr =&gt; block_ptr % next
-      end do
+      if (.not. config_do_restart) then
+         block_ptr =&gt; domain % blocklist
+         do while (associated(block_ptr))
+            call mpas_init(block_ptr, block_ptr % mesh, dt)
+            block_ptr % state % time_levs(1) % state % xtime % scalar = 0.0
+            block_ptr =&gt; block_ptr % next
+         end do
+      end if
 
       ! Before integrating, write out the initial state
       output_frame = 1

</font>
</pre>