<p><b>duda</b> 2010-10-06 15:53:14 -0600 (Wed, 06 Oct 2010)</p><p>BRANCH COMMIT<br>
<br>
Move diagnostic fields from 'state' into a new structure 'diag', which<br>
has just one time level.<br>
<br>
Code was tested with TC2 (J&W baroclinic wave) and produces identical <br>
results to those of the top-of-the-branch code.<br>
<br>
Notes:<br>
<br>
- with only a single time level for diagnostic fields, there is no longer<br>
a need to copy diagnostic fields in rk_integration_setup (which was<br>
fixed in r520)<br>
<br>
- with uReconstruct* fields in diag, we need to also pass diag into<br>
the reconstruct(...) routine in module_vector_reconstruction.F; somehow<br>
we will need to reconcile the new argument list with that in the trunk<br>
code<br>
<br>
- although many differences for the qd_kessler() routine show up, these<br>
"changes" are just the addition of indentation, with no other substance;<br>
I've had the itch to indent the code in this subroutine for a while and<br>
decided to finally do it<br>
<br>
<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>
M src/operators/module_vector_reconstruction.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2010-10-06 18:22:15 UTC (rev 520)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2010-10-06 21:53:14 UTC (rev 521)
@@ -169,37 +169,33 @@
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 tracers ( nTracers nVertLevels nCells Time ) 2 iro tracers state - -
-
# state variables diagnosed from prognostic state
-# var persistent real ww ( nVertLevelsP1 nCells Time ) 2 ro ww state - -
var persistent real pressure ( nVertLevels nCells Time ) 2 ro pressure state - -
-# var persistent real pp ( nVertLevelsP1 nCells Time ) 2 - pp state - -
var persistent real u_init ( nVertLevels ) 0 iro u_init mesh - -
var persistent real t_init ( nVertLevels nCells ) 0 iro t_init mesh - -
var persistent real qv_init ( nVertLevels ) 0 iro qv_init mesh - -
# Diagnostic fields: only written to output
-var persistent real v ( nVertLevels nEdges Time ) 2 o v state - -
-var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
-var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
-var persistent real pv_edge ( nVertLevels nEdges Time ) 2 o pv_edge state - -
-var persistent real rho_edge ( nVertLevels nEdges Time ) 2 o rho_edge state - -
-var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
-var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 o pv_vertex state - -
-var persistent real pv_cell ( nVertLevels nCells Time ) 2 o pv_cell state - -
-var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
-var persistent real uReconstructY ( nVertLevels nCells Time ) 2 o uReconstructY state - -
-var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 o uReconstructZ state - -
-var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
-var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
+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 - -
# Other diagnostic variables: neither read nor written to any files
-var persistent real rv ( nVertLevels nEdges Time ) 2 - rv state - -
-var persistent real circulation ( nVertLevels nVertices Time ) 2 - circulation state - -
-var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
-var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
+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 exner ( nVertLevels nCells ) 0 - exner mesh - -
@@ -247,7 +243,7 @@
# 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 state - -
+var persistent real kdiff ( nVertLevels nCells Time ) 2 - kdiff diag - -
# Arrays required for reconstruction of velocity field
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - 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-06 18:22:15 UTC (rev 520)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_test_cases.F 2010-10-06 21:53:14 UTC (rev 521)
@@ -38,7 +38,7 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
write(0,*) ' calling test case setup '
- call nhyd_test_case_jw(block_ptr % mesh, block_ptr % state % time_levs(1) % state, config_test_case)
+ call nhyd_test_case_jw(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
write(0,*) ' returned from test case setup '
do i=2,nTimeLevs
call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
@@ -55,7 +55,7 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
write(0,*) ' calling test case setup '
- call nhyd_test_case_squall_line(domain % dminfo, block_ptr % mesh, block_ptr % state % time_levs(1) % state, config_test_case)
+ call nhyd_test_case_squall_line(domain % dminfo, block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
write(0,*) ' returned from test case setup '
do i=2,nTimeLevs
call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
@@ -70,7 +70,7 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
write(0,*) ' calling test case setup '
- call nhyd_test_case_mtn_wave(block_ptr % mesh, block_ptr % state % time_levs(1) % state, config_test_case)
+ call nhyd_test_case_mtn_wave(block_ptr % mesh, block_ptr % state % time_levs(1) % state, block_ptr % diag, config_test_case)
write(0,*) ' returned from test case setup '
do i=2,nTimeLevs
call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
@@ -90,7 +90,7 @@
!----------------------------------------------------------------------------------------------------------
- subroutine nhyd_test_case_jw(grid, state, test_case)
+ subroutine nhyd_test_case_jw(grid, state, diag, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -99,6 +99,7 @@
type (mesh_type), intent(inout) :: grid
type (state_type), intent(inout) :: state
+ type (diag_type), intent(inout) :: diag
integer, intent(in) :: test_case
real (kind=RKIND), parameter :: u0 = 35.0
@@ -729,13 +730,13 @@
!
! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
!
- state % v % array(:,:) = 0.0
+ diag % v % array(:,:) = 0.0
do iEdge = 1, grid%nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
if (eoe > 0) then
do k = 1, grid%nVertLevels
- state % v % array(k,iEdge) = state % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
+ diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
end do
end if
end do
@@ -750,7 +751,6 @@
write(0,*) ' i, psurf, lat ',i,psurf,grid%latCell%array(i)*180./3.1415828
enddo
-! stop
end subroutine nhyd_test_case_jw
@@ -806,7 +806,7 @@
!----------------------------------------------------------------------------------------------------------
- subroutine nhyd_test_case_squall_line(dminfo, grid, state, test_case)
+ subroutine nhyd_test_case_squall_line(dminfo, grid, state, diag, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup squall line and supercell test case
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -816,6 +816,7 @@
type (dm_info), intent(in) :: dminfo
type (mesh_type), intent(inout) :: grid
type (state_type), intent(inout) :: state
+ type (diag_type), intent(inout) :: diag
integer, intent(in) :: test_case
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
@@ -1310,13 +1311,13 @@
!
! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
!
- state % v % array(:,:) = 0.0
+ diag % v % array(:,:) = 0.0
do iEdge = 1, grid%nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
if (eoe > 0) then
do k = 1, grid%nVertLevels
- state % v % array(k,iEdge) = state % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
+ diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
end do
end if
end do
@@ -1333,7 +1334,7 @@
!----------------------------------------------------------------------------------------------------------
- subroutine nhyd_test_case_mtn_wave(grid, state, test_case)
+ subroutine nhyd_test_case_mtn_wave(grid, state, diag, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1342,6 +1343,7 @@
type (mesh_type), intent(inout) :: grid
type (state_type), intent(inout) :: state
+ type (diag_type), intent(inout) :: diag
integer, intent(in) :: test_case
real (kind=RKIND), parameter :: t0=288., hm=250.
@@ -1904,13 +1906,13 @@
!
! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
!
- state % v % array(:,:) = 0.0
+ diag % v % array(:,:) = 0.0
do iEdge = 1, grid%nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
if (eoe > 0) then
do k = 1, grid%nVertLevels
- state % v % array(k,iEdge) = state % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
+ diag % v % array(k,iEdge) = diag % v %array(k,iEdge) + weightsOnEdge(i,iEdge) * state % u % array(k, eoe)
end do
end if
end do
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-06 18:22:15 UTC (rev 520)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2010-10-06 21:53:14 UTC (rev 521)
@@ -113,15 +113,15 @@
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
! vorticity
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(1) % state % vorticity % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % vorticity % array(:,:), &
block % mesh % nVertLevels, block % mesh % nVertices, &
block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
! divergence
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(1) % state % divergence % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % divergence % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
! pv_edge
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(1) % state % pv_edge % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % pv_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
! rtheta_p
@@ -164,7 +164,7 @@
if (debug) write(0,*) ' compute_dyn_tend '
block => domain % blocklist
do while (associated(block))
- call compute_dyn_tend( block % tend, block % state % time_levs(2) % state, block % mesh )
+ call compute_dyn_tend( block % tend, block % state % time_levs(2) % state, block % diag, block % mesh )
block => block % next
end do
if (debug) write(0,*) ' finished compute_dyn_tend '
@@ -261,10 +261,12 @@
if (rk_step < 3 .or. (.not. config_monotonic .and. .not. config_positive_definite)) then
call advance_scalars( block % tend, &
block % state % time_levs(1) % state, block % state % time_levs(2) % state, &
+ block % diag, &
block % mesh, rk_timestep(rk_step) )
else
call advance_scalars_mono( block % tend, &
block % state % time_levs(1) % state, block % state % time_levs(2) % state, &
+ block % diag, &
block % mesh, rk_timestep(rk_step), rk_step, 3, &
domain % dminfo, block % parinfo % cellsToSend, block % parinfo % cellsToRecv )
end if
@@ -291,7 +293,7 @@
block => domain % blocklist
do while (associated(block))
- call compute_solve_diagnostics( dt, block % state % time_levs(2) % state, block % mesh )
+ call compute_solve_diagnostics( dt, block % state % time_levs(2) % state, block % diag, block % mesh )
block => block % next
end do
@@ -322,28 +324,28 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % w % array(:,:), &
block % mesh % nVertLevels+1, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % vorticity % array(:,:), &
block % mesh % nVertLevels, block % mesh % nVertices, &
block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % pv_vertex % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % pv_vertex % array(:,:), &
block % mesh % nVertLevels, block % mesh % nVertices, &
block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % pv_cell % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % pv_cell % 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 % divergence % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % divergence % 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 % ke % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % ke % 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 % v % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % v % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % pv_edge % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % pv_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % rho_edge % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % diag % rho_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
call dmpar_exch_halo_field3dReal(domain % dminfo, block % state % time_levs(2) % state % scalars % array(:,:,:), &
@@ -449,17 +451,8 @@
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
- s_old % rho_edge % array = s_new % rho_edge % array
- s_old % v % array = s_new % v % array
- s_old % circulation % array = s_new % circulation % array
- s_old % vorticity % array = s_new % vorticity % array
- s_old % divergence % array = s_new % divergence % array
- s_old % ke % array = s_new % ke % array
- s_old % pv_edge % array = s_new % pv_edge % array
-
end subroutine rk_integration_setup
!-----
@@ -624,11 +617,11 @@
end do ! loop over cells
- end subroutine compute_vert_imp_coefs
+ end subroutine compute_vert_imp_coefs
!------------------------
- subroutine set_smlstep_pert_variables( s_old, s_new, tend, grid )
+ subroutine set_smlstep_pert_variables( s_old, s_new, tend, grid )
implicit none
type (state_type) :: s_new, s_old
@@ -687,11 +680,11 @@
grid % ruAvg % array = 0.
grid % wwAvg % array = 0.
- end subroutine set_smlstep_pert_variables
+ end subroutine set_smlstep_pert_variables
!-------------------------------
- subroutine advance_acoustic_step( s, tend, grid, dts )
+ subroutine advance_acoustic_step( s, tend, grid, dts )
implicit none
@@ -1112,7 +1105,7 @@
!---------------------------------------------------------------------------------------
- subroutine advance_scalars( tend, s_old, s_new, grid, dt)
+ subroutine advance_scalars( tend, s_old, s_new, diag, grid, dt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -1126,6 +1119,7 @@
type (tend_type), intent(in) :: tend
type (state_type), intent(in) :: s_old
type (state_type), intent(inout) :: s_new
+ type (diag_type), intent(in) :: diag
type (mesh_type), intent(in) :: grid
real (kind=RKIND) :: dt
@@ -1164,7 +1158,7 @@
scalar_old => s_old % scalars % array
scalar_new => s_new % scalars % array
- kdiff => s_new % kdiff % array
+ kdiff => diag % kdiff % array
deriv_two => grid % deriv_two % array
!**** uhAvg => grid % uhAvg % array
uhAvg => grid % ruAvg % array
@@ -1190,7 +1184,7 @@
h_theta_eddy_visc2 = config_h_theta_eddy_visc2
v_theta_eddy_visc2 = config_v_theta_eddy_visc2
- rho_edge => s_new % rho_edge % array
+ rho_edge => diag % rho_edge % array
rho => s_new % rho % array
qv_init => grid % qv_init % array
zgrid => grid % zgrid % array
@@ -1441,7 +1435,7 @@
end subroutine advance_scalars
- subroutine advance_scalars_mono( tend, s_old, s_new, grid, dt, rk_step, rk_order, dminfo, cellsToSend, cellsToRecv)
+ subroutine advance_scalars_mono( tend, s_old, s_new, diag, grid, dt, rk_step, rk_order, dminfo, cellsToSend, cellsToRecv)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -1455,6 +1449,7 @@
type (tend_type), intent(in) :: tend
type (state_type), intent(in) :: s_old
type (state_type), intent(inout) :: s_new
+ type (diag_type), intent(in) :: diag
type (mesh_type), intent(in) :: grid
integer, intent(in) :: rk_step, rk_order
real (kind=RKIND), intent(in) :: dt
@@ -1825,7 +1820,7 @@
!----
- subroutine compute_dyn_tend(tend, s, grid)
+ subroutine compute_dyn_tend(tend, s, diag, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute height and normal wind tendencies, as well as diagnostic variables
!
@@ -1841,6 +1836,7 @@
type (tend_type), intent(inout) :: tend
type (state_type), intent(in) :: s
+ type (diag_type), intent(in) :: diag
type (mesh_type), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, iq
@@ -1899,21 +1895,21 @@
coef_3rd_order = config_coef_3rd_order
rho => s % rho % array
- rho_edge => s % rho_edge % array
+ rho_edge => diag % rho_edge % array
rb => grid % rho_base % array
rr => s % rho_p % array
u => s % u % array
- v => s % v % array
- kdiff => s % kdiff % array
+ v => diag % v % array
+ kdiff => diag % kdiff % array
ru => grid % ru % array
w => s % w % array
rw => grid % rw % array
theta => s % theta % array
- circulation => s % circulation % array
- divergence => s % divergence % array
- vorticity => s % vorticity % array
- ke => s % ke % array
- pv_edge => s % pv_edge % array
+ circulation => diag % circulation % array
+ divergence => diag % divergence % array
+ vorticity => diag % vorticity % array
+ 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
@@ -2936,7 +2932,7 @@
!-------
- subroutine compute_solve_diagnostics(dt, s, grid)
+ subroutine compute_solve_diagnostics(dt, s, diag, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
!
@@ -2949,6 +2945,7 @@
real (kind=RKIND), intent(in) :: dt
type (state_type), intent(inout) :: s
+ type (diag_type), intent(inout) :: diag
type (mesh_type), intent(in) :: grid
@@ -2964,23 +2961,20 @@
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
-! h => s % h % array
h => s % rho % array
u => s % u % array
- v => s % v % array
- vh => s % rv % array
- h_edge => s % rho_edge % array
-! tend_h => s % h % array
-! tend_u => s % u % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- pv_edge => s % pv_edge % array
- pv_vertex => s % pv_vertex % array
- pv_cell => s % pv_cell % array
- gradPVn => s % gradPVn % array
- gradPVt => s % gradPVt % array
+ v => diag % v % array
+ vh => diag % rv % array
+ h_edge => diag % rho_edge % array
+ circulation => diag % circulation % array
+ vorticity => diag % vorticity % array
+ divergence => diag % divergence % array
+ ke => diag % ke % array
+ pv_edge => diag % pv_edge % array
+ pv_vertex => diag % pv_vertex % array
+ pv_cell => diag % pv_cell % array
+ gradPVn => diag % gradPVn % array
+ gradPVt => diag % gradPVt % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -3235,71 +3229,71 @@
subroutine qd_kessler( state_old, state_new, grid, dt )
- implicit none
+ implicit none
+
+ type (state_type), intent(inout) :: state_old, state_new
+ type (mesh_type), intent(inout) :: grid
+ real (kind=RKIND), intent(in) :: dt
- type (state_type), intent(inout) :: state_old, state_new
- type (mesh_type), intent(inout) :: grid
- real (kind=RKIND), intent(in) :: dt
+ real (kind=RKIND), dimension( grid % nVertLevels ) :: t, rho, p, dzu, qv, qc, qr, qc1, qr1
+
+ integer :: k,iEdge,i,iCell,nz1
+ real (kind=RKIND) :: p0,rcv
+
+
+ write(0,*) ' in qd_kessler '
+
+ p0 = 1.e+05
+ rcv = rgas/(cp-rgas)
+ nz1 = grid % nVertLevels
+
+ do iCell = 1, grid % nCellsSolve
+
+ do k = 1, grid % nVertLevels
+
+ grid % 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)
+ 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))
+ qc1(k) = max(0.,state_old % scalars % array(state_old % index_qc,k,iCell))
+ qr1(k) = max(0.,state_old % scalars % array(state_old % index_qr,k,iCell))
+ dzu(k) = grid % dzu % array(k)
+
+ end do
+
+ call kessler( t,qv,qc,qc1,qr,qr1,rho,p,dt,dzu,nz1, 1)
+
+ 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)
+ 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) = &
+ ( grid % zz % array(k,iCell)*(rgas/p0) * ( &
+ grid % rtheta_p % array(k,iCell) &
+ + grid % 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)) )
+ end do
+
+ end do
+
+ write(0,*) ' exiting qd_kessler '
- real (kind=RKIND), dimension( grid % nVertLevels ) :: t, rho, p, dzu, qv, qc, qr, qc1, qr1
-
- integer :: k,iEdge,i,iCell,nz1
- real (kind=RKIND) :: p0,rcv
-
-
- write(0,*) ' in qd_kessler '
-
- p0 = 1.e+05
- rcv = rgas/(cp-rgas)
- nz1 = grid % nVertLevels
-
- do iCell = 1, grid % nCellsSolve
-
- do k = 1, grid % nVertLevels
-
- grid % 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)
- 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))
- qc1(k) = max(0.,state_old % scalars % array(state_old % index_qc,k,iCell))
- qr1(k) = max(0.,state_old % scalars % array(state_old % index_qr,k,iCell))
- dzu(k) = grid % dzu % array(k)
-
- end do
-
- call kessler( t,qv,qc,qc1,qr,qr1,rho,p,dt,dzu,nz1, 1)
-
- 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)
- 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) = &
- ( grid % zz % array(k,iCell)*(rgas/p0) * ( &
- grid % rtheta_p % array(k,iCell) &
- + grid % 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)) )
- end do
-
- end do
-
- write(0,*) ' exiting qd_kessler '
-
end subroutine qd_kessler
!-----------------------------------------------------------------------
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-10-06 18:22:15 UTC (rev 520)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-10-06 21:53:14 UTC (rev 521)
@@ -30,10 +30,9 @@
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
-! call compute_solver_constants(block % state % time_levs(1) % state, mesh)
-! call compute_state_diagnostics(block % state % time_levs(1) % state, mesh)
call init_coupled_diagnostics( block % state % time_levs(1) % state, mesh)
- call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh) ! ok for nonhydrostatic model
+ call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, block % diag, mesh)
+
!
! Note: The following initialization calls have been moved to mpas_setup_test_case()
! since values computed by these routines are needed to produce initial fields
Modified: branches/atmos_nonhydrostatic/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/atmos_nonhydrostatic/src/operators/module_vector_reconstruction.F 2010-10-06 18:22:15 UTC (rev 520)
+++ branches/atmos_nonhydrostatic/src/operators/module_vector_reconstruction.F 2010-10-06 21:53:14 UTC (rev 521)
@@ -112,7 +112,7 @@
end subroutine init_reconstruct
- subroutine reconstruct(state, grid)
+ subroutine reconstruct(state, diag, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: reconstruct vector field at cell centers based on radial basis functions
!
@@ -123,7 +123,8 @@
implicit none
- type (state_type), intent(inout) :: state
+ type (state_type), intent(in) :: state
+ type (diag_type), intent(inout) :: diag
type (mesh_type), intent(in) :: grid
! temporary arrays needed in the compute procedure
@@ -151,14 +152,14 @@
nEdgesOnCell=> grid % nEdgesOnCell % array
nCellsSolve = grid % nCellsSolve
u => state % u % array
- uReconstructX => state % uReconstructX % array
- uReconstructY => state % uReconstructY % array
- uReconstructZ => state % uReconstructZ % array
+ uReconstructX => diag % uReconstructX % array
+ uReconstructY => diag % uReconstructY % array
+ uReconstructZ => diag % uReconstructZ % array
latCell => grid % latCell % array
lonCell => grid % lonCell % array
- uReconstructZonal => state % uReconstructZonal % array
- uReconstructMeridional => state % uReconstructMeridional % array
+ uReconstructZonal => diag % uReconstructZonal % array
+ uReconstructMeridional => diag % uReconstructMeridional % array
on_a_sphere = grid % on_a_sphere
! init the intent(out)
</font>
</pre>