<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&W test cases<br>
- bug fix to computation of zb and zb3 in J&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* &
- (rr(k ,i)+(rr(k ,i)+rb(k ,i))*qv_2d(k ,i) &
- +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* &
+! (rr(k ,i)+(rr(k ,i)+rb(k ,i))*qv_2d(k ,i) &
+! +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* &
+ ( (rr(k ,i)+(rr(k ,i)+rb(k ,i))*qv_2d(k ,i))*fzp(k+1) &
+ + (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* &
- (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i) &
- +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* &
+! (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i) &
+! +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* &
+ ( (rr(k ,i)+(rr(k ,i)+rb(k ,i))*scalars(index_qv,k ,i))*fzp(k+1) &
+ + (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) > 0) &
+ 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) > 0) &
+ 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)) &
- (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>