[mpas-developers] /home/subversion/mpas/model revision 557
Laura Fowler
laura at ucar.edu
Fri Oct 15 10:43:09 MDT 2010
Hi Michael:
If you are happy with your changes below, can you add them to the
atmos_physics branch? Thanks.
Laura
On 10/14/10 3:48 PM, mpas-developers at ucar.edu wrote:
>
> *duda* 2010-10-14 15:48:16 -0600 (Thu, 14 Oct 2010)
>
> BRANCH COMMIT
>
> Add changes that, among other things, enable bit-for-bit restarts
> in the non-hydrostatic model.
>
> Specifically:
>
> - move time-varying diagnostic fields from mesh into diag so that
> they can be properly written as time-varying fields
>
> - output every field to the restart file (not optimal, but we can
> whittle this set down)
>
> - move rho_p from state to diag
>
> - add halo exchanges at the beginning of each time step for
> cqu, cqw, and rw
>
> - add loop checks to prevent "dangerous" accesses of the garbage
> cell/edge/vertex
>
>
> M src/driver/module_subdriver.F
> M src/core_nhyd_atmos/mpas_interface.F
> M src/core_nhyd_atmos/module_test_cases.F
> M src/core_nhyd_atmos/Registry
> M src/core_nhyd_atmos/module_time_integration.F
>
> ------------------------------------------------------------------------
> 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
>
>
>
> _______________________________________________
> mpas-developers mailing list
> mpas-developers at mailman.ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/mpas-developers
--
!----------------------------------------------------
Laura D. Fowler
Mesoscale and Microscale Meteorology Division (MMM)
National Center for Atmospheric Research
P.O. Box 3000, Boulder CO 80307-3000
e-mail: laura at ucar.edu
phone : 303-497-1628
!----------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/mpas-developers/attachments/20101015/3984eb04/attachment-0001.html
More information about the mpas-developers
mailing list