<p><b>duda</b> 2011-03-24 11:43:19 -0600 (Thu, 24 Mar 2011)</p><p>BRANCH COMMIT<br>
<br>
Fixes for the J&W test-case initialization:<br>
<br>
- Fix bug in array indices ('i' should be 'iEdge') that was leading to <br>
the use of the wrong densities when computing ru.<br>
<br>
- Remove if-tests that were causing spurious differences in initial<br>
conditions when running on different numbers of processors.<br>
<br>
- Remove if-tests for cell or edge indices greater than zero, which<br>
have no effect since the introduction of the garbage cell/edge.<br>
<br>
<br>
M src/core_nhyd_atmos/module_test_cases.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2011-03-23 22:27:25 UTC (rev 763)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2011-03-24 17:43:19 UTC (rev 764)
@@ -696,13 +696,11 @@
state % u % array(k,iEdge) = fluxk + u_pert
end do
- cell1 = grid % CellsOnEdge % array(1,i)
- cell2 = grid % CellsOnEdge % array(2,i)
- if(cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then
- do k=1,nz1
- diag % ru % array (k,iEdge) = 0.5*(rho(k,cell1)+rho(k,cell2))*state % u % array (k,iEdge)
- end do
- end if
+ cell1 = grid % CellsOnEdge % array(1,iEdge)
+ cell2 = grid % CellsOnEdge % array(2,iEdge)
+ do k=1,nz1
+ diag % ru % array (k,iEdge) = 0.5*(rho(k,cell1)+rho(k,cell2))*state % u % array (k,iEdge)
+ end do
!
! Generate rotated Coriolis field
@@ -731,55 +729,51 @@
do iEdge = 1,grid % nEdges
cell1 = CellsOnEdge(1,iEdge)
cell2 = CellsOnEdge(2,iEdge)
- if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
- do k = 1, grid%nVertLevels
+ do k = 1, grid%nVertLevels
- if (config_theta_adv_order == 2) then
+ if (config_theta_adv_order == 2) then
- z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
+ z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
- else if (config_theta_adv_order == 3 .or. config_theta_adv_order ==4) then !theta_adv_order == 3 or 4
+ else if (config_theta_adv_order == 3 .or. config_theta_adv_order ==4) then !theta_adv_order == 3 or 4
- 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)
- 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
+ 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))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ 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.
+ z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
+ - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
- if (config_theta_adv_order == 3) then
- z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
- else
- z_edge3 = 0.
- end if
-
+ if (config_theta_adv_order == 3) then
+ z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ else
+ z_edge3 = 0.
end if
- zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
- zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
- zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
- zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
+ end if
- if (k /= 1) then
- zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb(k,1,iEdge)
- zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb(k,2,iEdge)
- zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb3(k,1,iEdge)
- zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb3(k,2,iEdge)
- end if
+ zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
+ zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
+ zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
+ zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
- end do
+ if (k /= 1) then
+ zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb(k,1,iEdge)
+ zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb(k,2,iEdge)
+ zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1))*zb3(k,1,iEdge)
+ zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2))*zb3(k,2,iEdge)
+ end if
- end if
- end do
+ end do
+ end do
+
! for including terrain
diag % rw % array = 0.
state % w % array = 0.
@@ -788,7 +782,6 @@
cell1 = CellsOnEdge(1,iEdge)
cell2 = CellsOnEdge(2,iEdge)
- if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
do k = 2, grid%nVertLevels
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
@@ -802,7 +795,6 @@
end if
end do
- end if
end do
@@ -822,11 +814,9 @@
do iEdge = 1, grid%nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
- do k = 1, grid%nVertLevels
- diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
- end do
- end if
+ do k = 1, grid%nVertLevels
+ diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
+ end do
end do
end do
</font>
</pre>