<p><b>laura@ucar.edu</b> 2010-07-23 15:57:43 -0600 (Fri, 23 Jul 2010)</p><p>updated from Michael Duda; added call to cloud microphysics<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2010-07-23 21:55:43 UTC (rev 428)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2010-07-23 21:57:43 UTC (rev 429)
@@ -5,11 +5,14 @@
use constants
use dmpar
+#ifdef DO_PHYSICS
+ use module_microphysics
+#endif
contains
- subroutine timestep(domain, dt)
+ subroutine timestep(domain, dt, itimestep)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step
!
@@ -23,11 +26,12 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
+ integer, intent(in) :: itimestep
type (block_type), pointer :: block
if (trim(config_time_integration) == 'SRK3') then
- call srk3(domain, dt)
+ call srk3(domain, dt, itimestep)
else
write(0,*) 'Unknown time integration option '//trim(config_time_integration)
write(0,*) 'Currently, only ''SRK3'' is supported.'
@@ -43,7 +47,7 @@
end subroutine timestep
- subroutine srk3(domain, dt)
+ subroutine srk3(domain, dt, itimestep)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step using
! time-split RK3 scheme
@@ -60,7 +64,9 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
+ integer, intent(in) :: itimestep
+ integer :: iq
integer :: iCell, k, iEdge
type (block_type), pointer :: block
@@ -100,6 +106,43 @@
block => domain % blocklist
do while (associated(block))
+! theta
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! scalars
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(1) % state % scalars % array(:,:,:), &
+ num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! pressure
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % pressure % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! vorticity
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+! divergence
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! pv_edge
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(1) % state % pv_edge % array(:,:), &
+ 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(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! ru
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % ru % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ block => block % next
+ end do
+
+ block => domain % blocklist
+ 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 % time_levs(2) % state, block % time_levs(1) % state, block % mesh )
@@ -136,6 +179,22 @@
! we will need to communicate the momentum tendencies here - we want tendencies for all edges of owned cells
! because we are solving for all edges of owned cells
!***********************************
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % rho % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % w % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
block => domain % blocklist
do while (associated(block))
@@ -159,10 +218,30 @@
if(debug) write(0,*) ' acoustic step complete '
! 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(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % 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(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ block => block % next
+ end do
end do ! end of small stimestep loop
! 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(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
block => domain % blocklist
do while (associated(block))
@@ -197,6 +276,18 @@
block => block % next
end do
+ block => domain % blocklist
+ do while (associated(block))
+! For now, we do scalar halo updates later on...
+! call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(TEND) % scalars % array(:,:,:), &
+! num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % scalars % array(:,:,:), &
+! num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
else
write(0,*) ' no scalar advection '
@@ -214,11 +305,74 @@
! need communications here to fill out u, w, theta, p, and pp, scalars, etc
! so that they are available for next RK step or the first rk substep of the next timestep
+ block => domain % blocklist
+ do while (associated(block))
+! NB: if any of these cause differences, better to uncomment the version after qd_kessler?
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % 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 % 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 % time_levs(2) % state % rho_p % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % 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 % time_levs(2) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_vertex % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_cell % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % ke % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % v % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % rho_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % scalars % array(:,:,:), &
+ num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
end do ! rk_step loop
! microphysics here...
+#ifdef DO_PHYSICS
+ if(config_microp_scheme .ne. 'off') then
+ block => domain % blocklist
+ do while(associated(block))
+
+ call microphysics_driver ( block % intermediate_step(TEND), block % time_levs(2) % state, &
+ block % mesh, itimestep )
+
+ block => block % next
+ end do
+ endif
+#endif
if(do_microphysics) then
block => domain % blocklist
do while (associated(block))
@@ -227,6 +381,29 @@
end do
end if
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % 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(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % 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(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % scalars % array(:,:,:), &
+ num_scalars, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
! if(debug) then
block => domain % blocklist
do while (associated(block))
@@ -320,7 +497,7 @@
grid % cqw % array(k,iCell) = 1./(1.+qtot)
end do
end do
-
+
do iEdge = 1, nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
@@ -371,7 +548,8 @@
nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
! epssm = grid % epssm ! this should come in through the namelist ******************
- epssm = 0.1
+! epssm = 0.1
+ epssm = config_epssm
rdzu => grid % rdzu % array
rdzw => grid % rdzw % array
@@ -505,8 +683,8 @@
real (kind=RKIND) :: smdiv, c2, rcv
real (kind=RKIND), dimension( grid % nVertLevels ) :: du
real (kind=RKIND), dimension( grid % nVertLevels + 1 ) :: dpzx
- real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells ) :: ts, rs
- real (kind=RKIND), dimension( grid % nVertLevels + 1 , grid % nCells ) :: ws
+ real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells+1 ) :: ts, rs
+ real (kind=RKIND), dimension( grid % nVertLevels + 1 , grid % nCells+1 ) :: ws
integer :: cell1, cell2, iEdge, iCell, k
real (kind=RKIND) :: pgrad, flux1, flux2, flux, resm, epssm
@@ -569,17 +747,13 @@
nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
-! cf1, cf2 and cf3 should come from the initialization *************
+ cf1 = grid % cf1 % scalar
+ cf2 = grid % cf2 % scalar
+ cf3 = grid % cf3 % scalar
- cf1 = 1.5
- cf2 = -0.5
- cf3 = 0.
+ epssm = config_epssm
+ smdiv = config_smdiv
-! these values should come from the namelist *****************
-
- epssm = 0.1
- smdiv = 0.1
-
rcv = rgas/(cp-rgas)
c2 = cp*rcv
resm = (1.-epssm)/(1.+epssm)
@@ -723,7 +897,7 @@
real (kind=RKIND), intent(in) :: dt
real (kind=RKIND), dimension(:,:), pointer :: wwAvg, rw_save, w, rw, rw_p, rtheta_p, rtheta_pp, &
- rtheta_p_save, rt_diabatic_tend, rho_p, rho_p_save, &
+ rtheta_p_save, rt_diabatic, rho_p, rho_p_save, &
rho_pp, rho, rho_base, ruAvg, ru_save, ru_p, u, ru, &
exner, exner_base, rtheta_base, pressure_p, &
zz, theta, zgrid
@@ -749,7 +923,7 @@
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
+ rt_diabatic => grid % rt_diabatic_tend % array
theta => s % theta % array
rho => s % rho % array
@@ -785,10 +959,11 @@
rcv = rgas/(cp-rgas)
p0 = 1.e+05 ! this should come from somewhere else...
- cf1 = 1.5
- cf2 = -0.5
- cf3 = 0.
+ cf1 = grid % cf1 % scalar
+ cf2 = grid % cf2 % scalar
+ cf3 = grid % cf3 % scalar
+
! compute new density everywhere so we can compute u from ru.
! we will also need it to compute theta below
@@ -812,7 +987,7 @@
! recover owned-cell values in block
- if( iCell <= nCellsSolve ) then
+! if( iCell <= nCellsSolve ) then ! If using this test, more halo exchanges will be needed
if(debug) then
if( iCell == 479 ) then
@@ -847,9 +1022,8 @@
do k = 1, nVertLevels
- rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) ! - dt * rt_diabatic_tend(k,iCell)
+ rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) - dt*rho(k,iCell)*rt_diabatic(k,iCell)
-
theta(k,iCell) = (rtheta_p(k,iCell) + rtheta_base(k,iCell))/rho(k,iCell)
exner(k,iCell) = (zz(k,iCell)*(rgas/p0)*(rtheta_p(k,iCell)+rtheta_base(k,iCell)))**rcv
! pressure below is perturbation pressure - perhaps we should rename it in the Registry????
@@ -857,7 +1031,7 @@
* (exner(k,iCell)-exner_base(k,iCell)))
end do
- end if
+! end if
end do
@@ -870,7 +1044,7 @@
cell1 = CellsOnEdge(1,iEdge)
cell2 = CellsOnEdge(2,iEdge)
- if( cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+! if( cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then ! If using this test, more halo exchanges will be needed
do k = 1, nVertLevels
ruAvg(k,iEdge) = ru(k,iEdge) + (ruAvg(k,iEdge) / float(ns))
@@ -890,7 +1064,7 @@
w(k,cell1) = w(k,cell1)+flux/AreaCell(cell1)
enddo
- end if
+! end if
enddo
@@ -1210,9 +1384,9 @@
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension( num_scalars, grid % nEdges) :: h_flux
- real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: v_flux, v_flux_upwind, s_update
- real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: scale_out, scale_in
+ real (kind=RKIND), dimension( num_scalars, grid % nEdges+1) :: h_flux
+ real (kind=RKIND), dimension( num_scalars, grid % nCells+1, 2 ) :: v_flux, v_flux_upwind, s_update
+ real (kind=RKIND), dimension( num_scalars, grid % nCells+1, 2 ) :: scale_out, scale_in
real (kind=RKIND), dimension( num_scalars ) :: s_max, s_min, s_max_update, s_min_update
integer :: nVertLevels, km0, km1, ktmp, kcp1, kcm1
@@ -1558,7 +1732,7 @@
real (kind=RKIND), dimension(:), pointer :: fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, zgrid, rho_edge, rho, ru, u, v, tend_u, &
circulation, divergence, vorticity, ke, pv_edge, theta, rw, tend_rho, &
- h_diabatic, tend_theta, tend_w, w, cqw, rb, rr, pp, pressure_b, zz, zx, cqu, &
+ rt_diabatic, tend_theta, tend_w, w, cqw, rb, rr, pp, pressure_b, zz, zx, cqu, &
h_divergence
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
@@ -1622,7 +1796,7 @@
tend_theta => tend % theta % array
tend_w => tend % w % array
tend_rho => tend % rho % array
- h_diabatic => grid % rt_diabatic_tend % array
+ rt_diabatic => grid % rt_diabatic_tend % array
t_init => grid % t_init % array
@@ -1653,24 +1827,20 @@
tend_u(:,:) = 0.0
- cf1 = 1.5
- cf2 = -.5
- cf3 = 0.
+ cf1 = grid % cf1 % scalar
+ cf2 = grid % cf2 % scalar
+ cf3 = grid % cf3 % scalar
! tendency for density
! divergence_ru may calculated in the diagnostic subroutine - it is temporary
- allocate(divergence_ru(nVertLevels, nCells))
- allocate(qtot(nVertLevels, nCells))
+ allocate(divergence_ru(nVertLevels, nCells+1))
+ allocate(qtot(nVertLevels, nCells+1))
divergence_ru(:,:) = 0.0
h_divergence(:,:) = 0.
do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- write(6,*) '--- nCells=', nCells
- write(6,*) '--- cell1 =', cell1
- write(6,*) '--- cell2 =', cell2
-
do k=1,nVertLevels
flux = ru(k,iEdge)*dvEdge(iEdge)
divergence_ru(k,cell1) = divergence_ru(k,cell1) + flux
@@ -1733,7 +1903,7 @@
! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
!
- allocate(rv(nVertLevels, nEdges))
+ allocate(rv(nVertLevels, nEdges+1))
rv(:,:) = 0.0
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
@@ -1826,10 +1996,10 @@
if ( h_mom_eddy_visc4 > 0.0 ) then
- allocate(delsq_divergence(nVertLevels, nCells))
- allocate(delsq_u(nVertLevels, nEdges))
- allocate(delsq_circulation(nVertLevels, nVertices))
- allocate(delsq_vorticity(nVertLevels, nVertices))
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
delsq_u(:,:) = 0.0
@@ -2107,7 +2277,7 @@
if ( h_mom_eddy_visc4 > 0.0 ) then
- allocate(delsq_theta(nVertLevels, nCells))
+ allocate(delsq_theta(nVertLevels, nCells+1))
delsq_theta(:,:) = 0.
@@ -2342,7 +2512,7 @@
if ( h_theta_eddy_visc4 > 0.0 ) then
- allocate(delsq_theta(nVertLevels, nCells))
+ allocate(delsq_theta(nVertLevels, nCells+1))
delsq_theta(:,:) = 0.
@@ -2398,7 +2568,7 @@
wdtz(nVertLevels+1) = 0.
do k=1,nVertLevels
tend_theta(k,iCell) = tend_theta(k,iCell)/areaCell(iCell) -rdzw(k)*(wdtz(k+1)-wdtz(k))
-!! tend_theta(k,iCell) = tend_theta(k) + rho(k,iCell)*h_diabatic(k,iCell)
+ tend_theta(k,iCell) = tend_theta(k,iCell) + rho(k,iCell)*rt_diabatic(k,iCell)
end do
end do
@@ -2718,11 +2888,11 @@
! Modify PV edge with upstream bias.
!
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) *dt * gradPVn(k,iEdge)
- end do
- end do
+! do iEdge = 1,nEdges
+! do k = 1,nVertLevels
+! pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) *dt * gradPVn(k,iEdge)
+! end do
+! end do
end subroutine compute_solve_diagnostics
@@ -2801,7 +2971,9 @@
grid % rt_diabatic_tend % array(k,iCell) = state_new % theta % array(k,iCell)
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) * &
+! 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 % rt_diabatic_tend % 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)
</font>
</pre>