<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta content="text/html; charset=UTF-8" http-equiv="Content-Type">
</head>
<body bgcolor="#ffffff" text="#000000">
Hi Michael:<br>
If you are happy with your changes below, can you add them to the
atmos_physics branch? Thanks.<br>
Laura<br>
<br>
On 10/14/10 3:48 PM, <a class="moz-txt-link-abbreviated" href="mailto:mpas-developers@ucar.edu">mpas-developers@ucar.edu</a> wrote:
<blockquote cite="mid:20101014214817.0240D50043@svn02.cgd.ucar.edu"
type="cite">
<p><b>duda</b> 2010-10-14 15:48:16 -0600 (Thu, 14 Oct 2010)</p>
<p>BRANCH COMMIT<br>
<br>
Add changes that, among other things, enable bit-for-bit
restarts <br>
in the non-hydrostatic model.<br>
<br>
Specifically:<br>
<br>
- move time-varying diagnostic fields from mesh into diag so
that<br>
they can be properly written as time-varying fields<br>
<br>
- output every field to the restart file (not optimal, but we
can<br>
whittle this set down)<br>
<br>
- move rho_p from state to diag<br>
<br>
- add halo exchanges at the beginning of each time step for<br>
cqu, cqw, and rw<br>
<br>
- add loop checks to prevent "dangerous" accesses of the garbage<br>
cell/edge/vertex<br>
<br>
<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="noshade">
<pre><font color="gray">Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2010-10-14 21:39:11 UTC (rev 556)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2010-10-14 21:48:16 UTC (rev 557)
@@ -97,9 +97,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 - -
@@ -137,14 +137,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
@@ -154,20 +154,20 @@
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 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
# 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 - -
@@ -177,74 +177,74 @@
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 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 - -
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2010-10-14 21:39:11 UTC (rev 556)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2010-10-14 21:48:16 UTC (rev 557)
@@ -204,19 +204,19 @@
hx => grid % hx % array
dss => grid % dss % array
- pb => grid % exner_base % array
- rb => grid % rho_base % array
- tb => grid % theta_base % array
- rtb => grid % rtheta_base % array
- p => grid % exner % array
+ pb => diag % exner_base % array
+ rb => diag % rho_base % array
+ tb => diag % theta_base % array
+ rtb => diag % rtheta_base % array
+ p => diag % exner % array
- ppb => grid % pressure_base % array
+ ppb => diag % pressure_base % array
pp => state % pressure % array
rho => state % rho % array
- rr => state % rho_p % array
+ rr => diag % rho_p % array
t => state % theta % array
- rt => grid % rtheta_p % array
+ rt => diag % rtheta_p % array
scalars => state % scalars % array
@@ -620,7 +620,7 @@
cell2 = grid % CellsOnEdge % array(2,i)
if(cell1 <= nCellsSolve .or. cell2 <= 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
@@ -701,7 +701,7 @@
end do
! for including terrain
- grid % rw % array = 0.
+ diag % rw % array = 0.
state % w % array = 0.
do iEdge = 1,grid % nEdges
@@ -710,15 +710,15 @@
if (cell1 <= nCellsSolve .or. cell2 <= 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) &
- - 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) &
- + 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) &
+ - 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) &
+ + sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
end if
end do
@@ -892,22 +892,22 @@
hx => grid % hx % array
dss => grid % dss % array
- ppb => grid % pressure_base % array
- pb => grid % exner_base % array
- rb => grid % rho_base % array
- tb => grid % theta_base % array
- rtb => grid % rtheta_base % array
- p => grid % exner % array
- cqw => grid % cqw % array
+ ppb => diag % pressure_base % array
+ pb => diag % exner_base % array
+ rb => diag % rho_base % array
+ tb => diag % theta_base % array
+ rtb => diag % rtheta_base % array
+ p => diag % exner % array
+ cqw => diag % cqw % array
rho => state % rho % array
pp => state % pressure % array
- rr => state % rho_p % array
+ rr => diag % rho_p % array
t => state % theta % array
- rt => grid % rtheta_p % array
+ rt => diag % rtheta_p % array
u => state % u % array
- ru => grid % ru % array
+ ru => diag % ru % array
scalars => state % scalars % array
@@ -1289,7 +1289,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.
@@ -1435,22 +1435,22 @@
xCell => grid % xCell % array
yCell => grid % yCell % array
- ppb => grid % pressure_base % array
- pb => grid % exner_base % array
- rb => grid % rho_base % array
- tb => grid % theta_base % array
- rtb => grid % rtheta_base % array
- p => grid % exner % array
- cqw => grid % cqw % array
+ ppb => diag % pressure_base % array
+ pb => diag % exner_base % array
+ rb => diag % rho_base % array
+ tb => diag % theta_base % array
+ rtb => diag % rtheta_base % array
+ p => diag % exner % array
+ cqw => diag % cqw % array
rho => state % rho % array
pp => state % pressure % array
- rr => state % rho_p % array
+ rr => diag % rho_p % array
t => state % theta % array
- rt => grid % rtheta_p % array
+ rt => diag % rtheta_p % array
u => state % u % array
- ru => grid % ru % array
+ ru => diag % ru % array
scalars => state % scalars % array
@@ -1865,7 +1865,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
@@ -1879,13 +1879,13 @@
if (cell1 <= nCellsSolve .or. cell2 <= 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) &
+ diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
- sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
- grid % rw % array(k,cell1) = grid % rw % array(k,cell1) &
+ diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
+ sign(1.,ru(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
end if
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-10-14 21:39:11 UTC (rev 556)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-10-14 21:48:16 UTC (rev 557)
@@ -125,13 +125,27 @@
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
! rtheta_p
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_p % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rtheta_p % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
! ru
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % ru % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % ru % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ call dmpar_exch_halo_field2dReal(domain % dminfo, &
+ block % diag % cqu % array, &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % EdgesToSend, block % parinfo % EdgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, &
+ block % diag % cqw % array, &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % CellsToSend, block % parinfo % CellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, &
+ block % diag % rw % array, &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % CellsToSend, block % parinfo % CellsToRecv)
+
block => block % next
end do
@@ -139,7 +153,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 => block % next
end do
@@ -156,7 +170,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 => block % next
end do
@@ -193,8 +207,8 @@
block => domain % blocklist
do while (associated(block))
call set_smlstep_pert_variables( block % state % time_levs(1) % state, block % state % time_levs(2) % state, &
- 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 => block % next
end do
@@ -204,8 +218,8 @@
block => domain % blocklist
do while (associated(block))
- call advance_acoustic_step( block % state % time_levs(2) % state, block % tend, &
- block % mesh, rk_sub_timestep(rk_step) )
+ call advance_acoustic_step( block % state % time_levs(2) % state, block % diag, block % tend, &
+ block % mesh, rk_sub_timestep(rk_step) )
block => block % next
end do
@@ -214,13 +228,13 @@
! will need communications here for rtheta_pp
block => domain % blocklist
do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_pp % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rtheta_pp % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rw_p % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rw_p % array(:,:), &
block % mesh % nVertLevels+1, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % ru_p % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % ru_p % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
block => block % next
@@ -231,7 +245,7 @@
! will need communications here for rho_pp
block => domain % blocklist
do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rho_pp % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rho_pp % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
@@ -239,9 +253,9 @@
block => domain % blocklist
do while (associated(block))
- call recover_large_step_variables( block % state % time_levs(2) % state, &
- block % mesh, rk_timestep(rk_step), &
- number_sub_steps(rk_step), rk_step )
+ call recover_large_step_variables( block % state % time_levs(2) % state, &
+ block % diag, block % tend, block % mesh, &
+ rk_timestep(rk_step), number_sub_steps(rk_step), rk_step )
block => block % next
end do
@@ -318,7 +332,7 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % rho % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % rho_p % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rho_p % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % w % array(:,:), &
@@ -361,17 +375,17 @@
if(do_microphysics) then
block => domain % blocklist
do while (associated(block))
- call qd_kessler( block % state % time_levs(1) % state, block % state % time_levs(2) % state, block % mesh, dt )
+ call qd_kessler( block % state % time_levs(1) % state, block % state % time_levs(2) % state, block % diag, block % tend, block % mesh, dt )
block => block % next
end do
end if
block => domain % blocklist
do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rtheta_p % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rtheta_p % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % exner % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % exner % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % theta % array(:,:), &
@@ -380,7 +394,7 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % pressure % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % rt_diabatic_tend % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % rt_diabatic_tend % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
call dmpar_exch_halo_field3dReal(domain % dminfo, block % state % time_levs(2) % state % scalars % array(:,:,:), &
@@ -433,22 +447,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
@@ -457,11 +470,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
@@ -478,7 +492,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
@@ -491,7 +505,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
@@ -500,7 +514,7 @@
!---
- 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
!
@@ -513,7 +527,8 @@
implicit none
type (state_type), intent(in) :: s
- type (mesh_type), intent(inout) :: grid
+ type (mesh_type), intent(in) :: grid
+ type (diag_type), intent(in) :: diag
real (kind=RKIND), intent(in) :: dts
integer :: i, k, iq
@@ -540,22 +555,22 @@
fzm => grid % fzm % array
fzp => grid % fzp % array
zz => grid % zz % array
- cqw => grid % cqw % array
+ cqw => diag % cqw % array
- p => grid % exner % array
- pb => grid % exner_base % array
- rt => grid % rtheta_p % array
- rtb => grid % rtheta_base % array
- rb => grid % rho_base % array
+ p => diag % exner % array
+ pb => diag % exner_base % array
+ rt => diag % rtheta_p % array
+ rtb => diag % rtheta_base % array
+ rb => diag % rho_base % array
- alpha_tri => grid % alpha_tri % array
- gamma_tri => grid % gamma_tri % array
- a_tri => grid % a_tri % array
- cofwr => grid % cofwr % array
- cofwz => grid % cofwz % array
- coftz => grid % coftz % array
- cofwt => grid % cofwt % array
- cofrz => grid % cofrz % array
+ alpha_tri => diag % alpha_tri % array
+ gamma_tri => diag % gamma_tri % array
+ a_tri => diag % a_tri % array
+ cofwr => diag % cofwr % array
+ cofwz => diag % cofwz % array
+ coftz => diag % coftz % array
+ cofwt => diag % cofwt % array
+ cofrz => diag % cofrz % array
t => s % theta % array
@@ -621,11 +636,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
@@ -641,12 +657,12 @@
areaCell => grid % areaCell % array
cellsOnEdge => 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
@@ -677,18 +693,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
@@ -725,24 +742,24 @@
theta => s % theta % array
w => s % w % array
- rtheta_pp => grid % rtheta_pp % array
- rtheta_pp_old => grid % rtheta_pp_old % array
- h_divergence => grid % h_divergence % array
- ru_p => grid % ru_p % array
- rw_p => grid % rw_p % array
- exner => grid % exner % array
- cqu => grid % cqu % array
- ruAvg => grid % ruAvg % array
- wwAvg => grid % wwAvg % array
- rho_pp => grid % rho_pp % array
- cofwt => grid % cofwt % array
- coftz => grid % coftz % array
- cofrz => grid % cofrz % array
- cofwr => grid % cofwr % array
- cofwz => grid % cofwz % array
- a_tri => grid % a_tri % array
- alpha_tri => grid % alpha_tri % array
- gamma_tri => grid % gamma_tri % array
+ rtheta_pp => diag % rtheta_pp % array
+ rtheta_pp_old => diag % rtheta_pp_old % array
+ h_divergence => diag % h_divergence % array
+ ru_p => diag % ru_p % array
+ rw_p => diag % rw_p % array
+ exner => diag % exner % array
+ cqu => diag % cqu % array
+ ruAvg => diag % ruAvg % array
+ wwAvg => diag % wwAvg % array
+ rho_pp => diag % rho_pp % array
+ cofwt => diag % cofwt % array
+ coftz => diag % coftz % array
+ cofrz => diag % cofrz % array
+ cofwr => diag % cofwr % array
+ cofwz => diag % cofwz % array
+ a_tri => diag % a_tri % array
+ alpha_tri => diag % alpha_tri % array
+ gamma_tri => diag % gamma_tri % array
dss => grid % dss % array
tend_ru => tend % u % array
@@ -899,10 +916,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
@@ -925,33 +944,33 @@
!---
- wwAvg => grid % wwAvg % array
- rw_save => grid % rw_save % array
- rw => grid % rw % array
- rw_p => grid % rw_p % array
+ wwAvg => diag % wwAvg % array
+ rw_save => diag % rw_save % array
+ rw => diag % rw % array
+ rw_p => diag % rw_p % array
w => s % w % array
- rtheta_p => grid % rtheta_p % array
- rtheta_p_save => grid % rtheta_p_save % array
- rtheta_pp => grid % rtheta_pp % array
- rtheta_base => grid % rtheta_base % array
- rt_diabatic_tend => grid % rt_diabatic_tend % array
+ rtheta_p => diag % rtheta_p % array
+ rtheta_p_save => diag % rtheta_p_save % array
+ rtheta_pp => diag % rtheta_pp % array
+ rtheta_base => diag % rtheta_base % array
+ rt_diabatic_tend => tend % rt_diabatic_tend % array
theta => s % theta % array
rho => s % rho % array
- rho_p => s % rho_p % array
- rho_p_save => grid % rho_p_save % array
- rho_pp => grid % rho_pp % array
- rho_base => grid % rho_base % array
+ rho_p => diag % rho_p % array
+ rho_p_save => diag % rho_p_save % array
+ rho_pp => diag % rho_pp % array
+ rho_base => diag % rho_base % array
- ruAvg => grid % ruAvg % array
- ru_save => grid % ru_save % array
- ru_p => grid % ru_p % array
- ru => grid % ru % array
+ ruAvg => diag % ruAvg % array
+ ru_save => diag % ru_save % array
+ ru_p => diag % ru_p % array
+ ru => diag % ru % array
u => s % u % array
- exner => grid % exner % array
- exner_base => grid % exner_base % array
+ exner => diag % exner % array
+ exner_base => diag % exner_base % array
pressure_p => s % pressure % array
@@ -1062,8 +1081,8 @@
cell1 = CellsOnEdge(1,iEdge)
cell2 = CellsOnEdge(2,iEdge)
-! if( cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then ! If using this test, more halo exchanges will be needed
-
+ if( cell1 <= nCellsSolve .or. cell2 <= 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))
@@ -1097,7 +1116,7 @@
end if
enddo
-! end if
+ end if
enddo
@@ -1161,7 +1180,7 @@
kdiff => diag % kdiff % array
deriv_two => grid % deriv_two % array
!**** uhAvg => grid % uhAvg % array
- uhAvg => grid % ruAvg % array
+ uhAvg => diag % ruAvg % array
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
cellsOnEdge => grid % cellsOnEdge % array
@@ -1170,7 +1189,7 @@
!**** h_new => s_new % h % array
h_old => s_old % rho % array
h_new => s_new % rho % array
- wwAvg => grid % wwAvg % array
+ wwAvg => diag % wwAvg % array
areaCell => grid % areaCell % array
!**** fnm => grid % fnm % array
@@ -1202,7 +1221,7 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= 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))
@@ -1219,7 +1238,7 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
@@ -1227,12 +1246,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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
d2fdx2_cell1 = d2fdx2_cell1 + &
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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + &
deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
end do
@@ -1262,7 +1281,7 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= 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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
d2fdx2_cell1 = d2fdx2_cell1 + &
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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + &
deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
end do
@@ -1300,7 +1319,7 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
do iScalar=1,s_old % num_scalars
@@ -1320,7 +1339,7 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
do iScalar=1,s_old % num_scalars
@@ -1493,7 +1512,7 @@
scalar_new => s_new % scalars % array
deriv_two => grid % deriv_two % array
!**** uhAvg => grid % uhAvg % array
- uhAvg => grid % ruAvg % array
+ uhAvg => diag % ruAvg % array
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
cellsOnEdge => grid % cellsOnEdge % array
@@ -1502,7 +1521,7 @@
!**** h_new => s_new % h % array
h_old => s_old % rho % array
h_new => s_new % rho % array
- wwAvg => grid % wwAvg % array
+ wwAvg => diag % wwAvg % array
areaCell => grid % areaCell % array
!**** fnm => grid % fnm % array
@@ -1591,7 +1610,7 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
cell_upwind = cell2
if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
do iScalar=1,s_old % num_scalars
@@ -1611,7 +1630,7 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
cell_upwind = cell2
if (uhAvg(k,iEdge) >= 0) cell_upwind = cell1
do iScalar=1,s_old % num_scalars
@@ -1619,12 +1638,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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
d2fdx2_cell1 = d2fdx2_cell1 + &
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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + &
deriv_two(i+1,2,iEdge) * scalar_new(iScalar,k,grid % CellsOnCell % array (i,cell2))
end do
@@ -1678,7 +1697,7 @@
end do
do i = 1, grid % nEdgesOnCell % array(iCell) ! go around the edges of each cell
- if (grid % cellsOnCell % array(i,iCell) > 0) then
+ if (grid % cellsOnCell % array(i,iCell) <= 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))
@@ -1734,7 +1753,7 @@
do iEdge = 1, grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do iScalar=1,s_old % num_scalars
flux = h_flux(iScalar,iEdge)
if (flux > 0) then
@@ -1779,7 +1798,7 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do iScalar=1,s_old % num_scalars
s_update(iScalar,cell1,km0) = s_update(iScalar,cell1,km0) - &
h_flux(iScalar,iEdge) / grid % areaCell % array(cell1)
@@ -1896,14 +1915,14 @@
rho => s % rho % array
rho_edge => diag % rho_edge % array
- rb => grid % rho_base % array
- rr => s % rho_p % array
+ rb => diag % rho_base % array
+ rr => diag % rho_p % array
u => s % u % array
v => diag % v % array
kdiff => diag % kdiff % array
- ru => grid % ru % array
+ ru => diag % ru % array
w => s % w % array
- rw => grid % rw % array
+ rw => diag % rw % array
theta => s % theta % array
circulation => diag % circulation % array
divergence => diag % divergence % array
@@ -1911,8 +1930,8 @@
ke => diag % ke % array
pv_edge => diag % pv_edge % array
pp => s % pressure % array
- pressure_b => grid % pressure_base % array
- h_divergence => grid % h_divergence % array
+ pressure_b => diag % pressure_base % array
+ h_divergence => diag % h_divergence % array
weightsOnEdge => grid % weightsOnEdge % array
@@ -1938,7 +1957,7 @@
tend_theta => tend % theta % array
tend_w => tend % w % array
tend_rho => tend % rho % array
- h_diabatic => grid % rt_diabatic_tend % array
+ h_diabatic => tend % rt_diabatic_tend % array
t_init => grid % t_init % array
qv_init => grid % qv_init % array
@@ -1948,8 +1967,8 @@
fzm => grid % fzm % array
fzp => grid % fzp % array
zgrid => grid % zgrid % array
- cqw => grid % cqw % array
- cqu => grid % cqu % array
+ cqw => diag % cqw % array
+ cqu => diag % cqu % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -2242,12 +2261,12 @@
delsq_circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
+ if (verticesOnEdge(1,iEdge) <= 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) > 0) then
+ if (verticesOnEdge(2,iEdge) <= 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
@@ -2404,11 +2423,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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * w(k,grid % CellsOnCell % array (i,cell2))
end do
@@ -2442,11 +2461,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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * w(k,grid % CellsOnCell % array (i,cell2))
end do
@@ -2677,11 +2696,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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
end do
@@ -2718,18 +2737,18 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= 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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell1) <= grid%nCells) &
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) > 0) &
+ if ( grid % CellsOnCell % array (i,cell2) <= grid%nCells) &
d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
end do
@@ -3005,7 +3024,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,nVertLevels
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
end do
@@ -3019,12 +3038,12 @@
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
+ if (verticesOnEdge(1,iEdge) <= 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) > 0) then
+ if (verticesOnEdge(2,iEdge) <= grid%nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
end do
@@ -3044,12 +3063,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= grid%nCells) then
do k=1,nVertLevels
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
end do
end if
- if(cell2 > 0) then
+ if(cell2 <= grid%nCells) then
do k=1,nVertLevels
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
end do
@@ -3087,7 +3106,7 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
+ if (eoe <= grid%nEdges) then
do k = 1,nVertLevels
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
@@ -3139,7 +3158,7 @@
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
+ if(iEdge <= grid%nEdges) then
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
end do
@@ -3168,7 +3187,7 @@
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
- if( iCell > 0) then
+ if( iCell <= 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
@@ -3184,7 +3203,7 @@
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
+ if( cellsOnEdge(1,iEdge) <= grid%nCells .and. cellsOnEdge(2,iEdge) <= grid%nCells) then
do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
@@ -3206,11 +3225,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
@@ -3219,7 +3239,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
@@ -3227,11 +3247,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
@@ -3251,11 +3273,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))
@@ -3270,24 +3292,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) * &
- (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) &
- - grid % rtheta_base % array(k,iCell)
+ tend % rt_diabatic_tend % array(k,iCell) = state_new % rho % array(k,iCell) * &
+ (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) &
+ - 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) = &
+ diag % exner % array(k,iCell) = &
( grid % zz % array(k,iCell)*(rgas/p0) * ( &
- grid % rtheta_p % array(k,iCell) &
- + grid % rtheta_base % array(k,iCell) ) )**rcv
+ diag % rtheta_p % array(k,iCell) &
+ + diag % rtheta_base % array(k,iCell) ) )**rcv
state_new % pressure % array(k,iCell) = &
grid % zz % array(k,iCell) * rgas * ( &
- grid % exner % array(k,iCell)*grid % rtheta_p % array(k,iCell) &
- +grid % rtheta_base % array(k,iCell) * &
- (grid % exner % array(k,iCell) - grid % exner_base % array(k,iCell)) )
+ diag % exner % array(k,iCell)*diag % rtheta_p % array(k,iCell) &
+ +diag % rtheta_base % array(k,iCell) * &
+ (diag % exner % array(k,iCell) - diag % exner_base % array(k,iCell)) )
end do
end do
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-10-14 21:39:11 UTC (rev 556)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-10-14 21:48:16 UTC (rev 557)
@@ -30,7 +30,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)
!
Modified: branches/atmos_nonhydrostatic/src/driver/module_subdriver.F
===================================================================
--- branches/atmos_nonhydrostatic/src/driver/module_subdriver.F 2010-10-14 21:39:11 UTC (rev 556)
+++ branches/atmos_nonhydrostatic/src/driver/module_subdriver.F 2010-10-14 21:48:16 UTC (rev 557)
@@ -56,12 +56,14 @@
! Compute diagnostic fields needed in solve loop, and initialize
! simulation time to 0 for all blocks
- block_ptr => 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 => block_ptr % next
- end do
+ if (.not. config_do_restart) then
+ block_ptr => 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 => block_ptr % next
+ end do
+ end if
! Before integrating, write out the initial state
output_frame = 1
</font>
</pre>
<pre wrap="">
<fieldset class="mimeAttachmentHeader"></fieldset>
_______________________________________________
mpas-developers mailing list
<a class="moz-txt-link-abbreviated" href="mailto:mpas-developers@mailman.ucar.edu">mpas-developers@mailman.ucar.edu</a>
<a class="moz-txt-link-freetext" href="http://mailman.ucar.edu/mailman/listinfo/mpas-developers">http://mailman.ucar.edu/mailman/listinfo/mpas-developers</a>
</pre>
</blockquote>
<br>
<pre class="moz-signature" cols="72">--
!----------------------------------------------------
Laura D. Fowler
Mesoscale and Microscale Meteorology Division (MMM)
National Center for Atmospheric Research
P.O. Box 3000, Boulder CO 80307-3000
e-mail: <a class="moz-txt-link-abbreviated" href="mailto:laura@ucar.edu">laura@ucar.edu</a>
phone : 303-497-1628
!----------------------------------------------------
</pre>
</body>
</html>