<p><b>duda</b> 2012-07-12 11:03:41 -0600 (Thu, 12 Jul 2012)</p><p>BRANCH COMMIT<br>
<br>
- changes to hydrostatic balance in J&amp;W test cases<br>
- bug fix to computation of zb and zb3 in J&amp;W test cases<br>
- in test case setup, take omega_e (rotation rate) from omega <br>
  in mpas_constants.F for consistency<br>
<br>
<br>
M    src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F        2012-07-11 21:40:59 UTC (rev 2015)
+++ branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F        2012-07-12 17:03:41 UTC (rev 2016)
@@ -166,7 +166,10 @@
 
       real (kind=RKIND), parameter :: u0 = 35.0
       real (kind=RKIND), parameter :: alpha_grid = 0.  ! no grid rotation
-      real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+
+!      real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+      real (kind=RKIND) :: omega_e
+
       real (kind=RKIND), parameter :: t0b = 250., t0 = 288., delta_t = 4.8e+05, dtdz = 0.005, eta_t = 0.2
       real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
       real (kind=RKIND), parameter :: theta_c = pii/4.0
@@ -312,6 +315,7 @@
 
       etavs = (1.-0.252)*pii/2.
       r_earth = a
+      omega_e = omega
       p0 = 1.e+05
 
       write(0,*) ' point 1 in test case setup '
@@ -518,9 +522,14 @@
 
             ppi(1) = ppi(1)-ppb(1,i)
             do k=1,nz1-1
-              ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity*                        &amp;
-                            (rr(k  ,i)+(rr(k  ,i)+rb(k  ,i))*qv_2d(k  ,i)   &amp;
-                            +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*qv_2d(k+1,i))
+
+!              ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity*                        &amp;
+!                            (rr(k  ,i)+(rr(k  ,i)+rb(k  ,i))*qv_2d(k  ,i)   &amp;
+!                            +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*qv_2d(k+1,i))
+
+              ppi(k+1) = ppi(k)-dzu(k+1)*gravity*                                       &amp;
+                            ( (rr(k  ,i)+(rr(k  ,i)+rb(k  ,i))*qv_2d(k  ,i))*fzp(k+1)   &amp;
+                            + (rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*qv_2d(k+1,i))*fzm(k+1) )
             end do
 
             do k=1,nz1
@@ -651,9 +660,15 @@
 
             ppi(1) = ppi(1)-ppb(1,i)
             do k=1,nz1-1
-              ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity*                     &amp;
-                            (rr(k  ,i)+(rr(k  ,i)+rb(k  ,i))*scalars(index_qv,k  ,i)   &amp;
-                            +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))
+
+!              ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity*                     &amp;
+!                            (rr(k  ,i)+(rr(k  ,i)+rb(k  ,i))*scalars(index_qv,k  ,i)   &amp;
+!                            +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))
+
+               ppi(k+1) = ppi(k)-dzu(k+1)*gravity*                                                  &amp;
+                             ( (rr(k  ,i)+(rr(k  ,i)+rb(k  ,i))*scalars(index_qv,k  ,i))*fzp(k+1)   &amp;
+                             + (rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*scalars(index_qv,k+1,i))*fzm(k+1) )
+
             end do
 
             do k=1,nz1
@@ -778,11 +793,18 @@
 
                d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
                d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
-               do i=1, grid % nEdgesOnCell % array (cell1)
-                  d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
-                  d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
-               end do
 
+!  WCS fix 20120711
+
+                  do i=1, grid % nEdgesOnCell % array (cell1)
+                     if ( grid % CellsOnCell % array (i,cell1) &gt; 0)       &amp;
+                     d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
+                  end do
+                  do i=1, grid % nEdgesOnCell % array (cell2)
+                     if ( grid % CellsOnCell % array (i,cell2) &gt; 0)       &amp;
+                     d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
+                  end do             
+
                z_edge =  0.5*(zgrid(k,cell1) + zgrid(k,cell2))         &amp;
                              - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
 
@@ -938,11 +960,14 @@
    real (kind=RKIND), dimension(nlat-1) :: f
    real (kind=RKIND), dimension(nz1+1)  :: dpzx
 
-   real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+!   real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+   real (kind=RKIND) :: omega_e
+
    real (kind=RKIND) :: rdx, qtot, r_earth, phi
    integer :: k,i, itr
 
    r_earth  = a
+   omega_e = omega
    rdx = 1./(dlat*r_earth)
 
    do i=1,nlat-1
@@ -2179,7 +2204,10 @@
 
       real (kind=RKIND), parameter :: u0 = 35.0
       real (kind=RKIND), parameter :: alpha_grid = 0.  ! no grid rotation
-      real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+
+!      real (kind=RKIND), parameter :: omega_e = 7.29212e-05
+      real (kind=RKIND) :: omega_e
+
       real (kind=RKIND), parameter :: t0b = 250., t0 = 288., delta_t = 4.8e+05, dtdz = 0.005, eta_t = 0.2
       real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
       real (kind=RKIND), parameter :: theta_c = pii/4.0
@@ -2338,6 +2366,7 @@
       etavs = (1.-0.252)*pii/2.
       rcv = rgas/(cp-rgas)
       r_earth = a
+      omega_e = omega
       p0 = 1.e+05
 
       interp_list(1) = FOUR_POINT

</font>
</pre>