<p><b>duda</b> 2011-09-19 14:47:24 -0600 (Mon, 19 Sep 2011)</p><p>BRANCH COMMIT<br>
<br>
Bug fixes to omega computation from Bill, uncovered by Soyoung's testing.<br>
<br>
<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/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2011-09-19 19:44:04 UTC (rev 1010)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2011-09-19 20:47:24 UTC (rev 1011)
@@ -1180,9 +1180,9 @@
w(1,cell1) = w(1,cell1) + zb(1,1,iEdge)*flux/(cf1*rho_zz(1,cell1)+cf2*rho_zz(2,cell1)+cf3*rho_zz(3,cell1))
!3rd order stencil
if (config_theta_adv_order == 3) then
- w(1,cell2) = w(1,cell2) - sign(1.,ru(1,iEdge))*config_coef_3rd_order &
+ w(1,cell2) = w(1,cell2) - sign(1.,flux)*config_coef_3rd_order &
*zb3(1,2,iEdge)*flux/(cf1*rho_zz(1,cell2)+cf2*rho_zz(2,cell2)+cf3*rho_zz(3,cell2))
- w(1,cell1) = w(1,cell1) + sign(1.,ru(1,iEdge))*config_coef_3rd_order &
+ w(1,cell1) = w(1,cell1) + sign(1.,flux)*config_coef_3rd_order &
*zb3(1,1,iEdge)*flux/(cf1*rho_zz(1,cell1)+cf2*rho_zz(2,cell1)+cf3*rho_zz(3,cell1))
end if
@@ -1192,9 +1192,9 @@
w(k,cell1) = w(k,cell1) + zb(k,1,iEdge)*flux/(fzm(k)*rho_zz(k,cell1)+fzp(k)*rho_zz(k-1,cell1))
!3rd order! stencil
if (config_theta_adv_order ==3) then
- w(k,cell2) = w(k,cell2) - sign(1.,ru(k,iEdge))*config_coef_3rd_order &
+ w(k,cell2) = w(k,cell2) - sign(1.,flux)*config_coef_3rd_order &
*zb3(k,2,iEdge)*flux/(fzm(k)*rho_zz(k,cell2)+fzp(k)*rho_zz(k-1,cell2))
- w(k,cell1) = w(k,cell1) + sign(1.,ru(k,iEdge))*config_coef_3rd_order &
+ w(k,cell1) = w(k,cell1) + sign(1.,flux)*config_coef_3rd_order &
*zb3(k,1,iEdge)*flux/(fzm(k)*rho_zz(k,cell1)+fzp(k)*rho_zz(k-1,cell1))
end if
enddo
@@ -3485,8 +3485,8 @@
type (diag_type), intent(inout) :: diag
type (mesh_type), intent(inout) :: grid
- integer :: k,iCell,iEdge,i,iCell1,iCell2
- real (kind=RKIND) :: p0, rcv
+ integer :: k,iCell,iEdge,i,iCell1,iCell2, cell1, cell2
+ real (kind=RKIND) :: p0, rcv, flux
rcv = rgas / (cp-rgas)
p0 = 1.e5 ! this should come from somewhere else...
@@ -3506,16 +3506,55 @@
end do
end do
- ! Compute w from rho_zz and rw
+! ! Compute w from rho_zz and rw
+! do iCell=1,grid%nCells
+! diag % rw % array(1,iCell) = 0.
+! diag % rw % array(grid%nVertLevels+1,iCell) = 0.
+! do k=2,grid%nVertLevels
+! diag % rw % array(k,iCell) = state % w % array(k,iCell) &
+! * (grid % fzp % array(k) * state % rho_zz % array(k-1,iCell) + grid % fzm % array(k) * state % rho_zz % array(k,iCell))
+! end do
+! end do
+
+
+! WCS bug fix 20110916
+
+ ! Compute rw (i.e. rho_zz * omega) from rho_zz, w, and ru.
+ ! We are reversing the procedure we use in subroutine recover_large_step_variables.
+ ! first, the piece that depends on w.
do iCell=1,grid%nCells
diag % rw % array(1,iCell) = 0.
diag % rw % array(grid%nVertLevels+1,iCell) = 0.
do k=2,grid%nVertLevels
diag % rw % array(k,iCell) = state % w % array(k,iCell) &
- * (grid % fzp % array(k) * state % rho_zz % array(k-1,iCell) + grid % fzm % array(k) * state % rho_zz % array(k,iCell))
+ * (grid % fzp % array(k) * state % rho_zz % array(k-1,iCell) + grid % fzm % array(k) * state % rho_zz % array(k,iCell)) &
+ * (grid % fzp % array(k) * grid % zz % array(k-1,iCell) + grid % fzm % array(k) * grid % zz % array(k,iCell))
end do
end do
+ ! next, the piece that depends on ru
+ do iEdge=1,grid%nEdges
+ cell1 = grid % CellsOnEdge % array(1,iEdge)
+ cell2 = grid % CellsOnEdge % array(2,iEdge)
+ do k = 2, grid % nVertLevels
+ flux = (grid % fzm % array(k) * diag % ru % array(k,iEdge)+grid % fzp % array(k) * diag % ru % array(k-1,iEdge))
+ diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + grid % zb % array(k,2,iEdge)*flux &
+ * (grid % fzp % array(k) * grid % zz % array(k-1,cell2) + grid % fzm % array(k) * grid % zz % array(k,cell2))
+ diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - grid % zb % array(k,1,iEdge)*flux &
+ * (grid % fzp % array(k) * grid % zz % array(k-1,cell1) + grid % fzm % array(k) * grid % zz % array(k,cell1))
+!3rd order! stencil
+ if (config_theta_adv_order ==3) then
+ diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + sign(1.,flux)*config_coef_3rd_order &
+ * grid % zb3 % array(k,2,iEdge)*flux &
+ * (grid % fzp % array(k) * grid % zz % array(k-1,cell2) + grid % fzm % array(k) * grid % zz % array(k,cell2))
+ diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - sign(1.,flux)*config_coef_3rd_order &
+ * grid % zb3 % array(k,1,iEdge)*flux &
+ * (grid % fzp % array(k) * grid % zz % array(k-1,cell1) + grid % fzm % array(k) * grid % zz % array(k,cell1))
+ end if
+ enddo
+ enddo
+! end WCS bug fix
+
do iCell = 1, grid % nCells
do k=1,grid % nVertLevels
diag % rho_p % array(k,iCell) = state % rho_zz % array(k,iCell) - diag % rho_base % array(k,iCell)
</font>
</pre>