<p><b>duda</b> 2011-04-20 13:08:23 -0600 (Wed, 20 Apr 2011)</p><p>BRANCH COMMIT<br>
<br>
- Evaluate mixing terms on the first RK step and re-use this value for<br>
the second and third RK steps. To revert to evaluation of mixing on<br>
all RK steps, set rk_diffusion = .true. in compute_dyn_tend().<br>
<br>
- Remove NCAR_FORMULATION code.<br>
<br>
- Move computations specific to APVM into an if-test on the value of<br>
config_apvm_upwinding.<br>
<br>
<br>
M src/core_nhyd_atmos/Registry<br>
M src/core_nhyd_atmos/module_time_integration.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 2011-04-20 00:04:37 UTC (rev 795)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/Registry 2011-04-20 19:08:23 UTC (rev 796)
@@ -169,6 +169,10 @@
var persistent real tend_qr ( nVertLevels nCells Time ) 1 - qr tend scalars moist
var persistent real rt_diabatic_tend ( nVertLevels nCells Time ) 1 - rt_diabatic_tend tend - -
+var persistent real euler_tend_u ( nVertLevels nEdges Time ) 1 - u_euler tend - -
+var persistent real euler_tend_w ( nVertLevelsP1 nCells Time ) 1 - w_euler tend - -
+var persistent real euler_tend_theta ( nVertLevels nCells Time ) 1 - theta_euler tend - -
+
# state variables diagnosed from prognostic state
var persistent real pressure_p ( nVertLevels nCells Time ) 1 ro pressure_p diag - -
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2011-04-20 00:04:37 UTC (rev 795)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/module_time_integration.F 2011-04-20 19:08:23 UTC (rev 796)
@@ -158,7 +158,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 % diag, block % mesh )
+ call compute_dyn_tend( block % tend, block % state % time_levs(2) % state, block % diag, block % mesh, rk_step )
block => block % next
end do
if (debug) write(0,*) ' finished compute_dyn_tend '
@@ -1814,7 +1814,7 @@
!----
- subroutine compute_dyn_tend(tend, s, diag, grid)
+ subroutine compute_dyn_tend(tend, s, diag, grid, rk_step)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute height and normal wind tendencies, as well as diagnostic variables
!
@@ -1832,7 +1832,10 @@
type (state_type), intent(in) :: s
type (diag_type), intent(in) :: diag
type (mesh_type), intent(in) :: grid
+ integer, intent(in) :: rk_step
+ logical, parameter :: rk_diffusion = .false.
+
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, iq
real (kind=RKIND) :: flux, vorticity_abs, rho_vertex, workpv, q, upstream_bias
@@ -1845,6 +1848,9 @@
circulation, divergence, vorticity, ke, pv_edge, theta, rw, tend_rho, &
rt_diabatic_tend, tend_theta, tend_w, w, cqw, rb, rr, pp, pressure_b, zz, zx, cqu, &
h_divergence, kdiff
+
+ real (kind=RKIND), dimension(:,:), pointer :: tend_u_euler, tend_w_euler, tend_theta_euler
+
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -1934,6 +1940,10 @@
tend_rho => tend % rho % array
rt_diabatic_tend => tend % rt_diabatic_tend % array
+ tend_u_euler => tend % u_euler % array
+ tend_theta_euler => tend % theta_euler % array
+ tend_w_euler => tend % w_euler % array
+
t_init => grid % t_init % array
qv_init => grid % qv_init % array
@@ -1962,6 +1972,8 @@
! Compute u (normal) velocity tendency for each edge (cell face)
!
+ write(0,*) ' rk_step in compute_dyn_tend ',rk_step
+
tend_u(:,:) = 0.0
cf1 = grid % cf1 % scalar
@@ -2001,7 +2013,7 @@
end do
delsq_horiz_mixing = .false.
- if (config_horiz_mixing == "2d_smagorinsky") then
+ if (config_horiz_mixing == "2d_smagorinsky" .and. (rk_step == 1 .or. rk_diffusion)) then
do iCell = 1, nCells
d_diag(:) = 0.
d_off_diag(:) = 0.
@@ -2022,7 +2034,6 @@
delsq_horiz_mixing = .true.
end if
-#ifdef LANL_FORMULATION
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -2055,59 +2066,6 @@
end do
-#endif
-
-#ifdef NCAR_FORMULATION
- !
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
-
- allocate(rv(nVertLevels, nEdges+1))
- rv(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- k = 1
- dpzx(k) = .5*zx(k,iEdge)*(cf1*(pp(k ,cell2)+pp(k ,cell1)) &
- +cf2*(pp(k+1,cell2)+pp(k+1,cell1)) &
- +cf3*(pp(k+2,cell2)+pp(k+2,cell1)))
- do k = 2, nVertLevels
- dpzx(k) = .5*zx(k,iEdge)*(fzm(k)*(pp(k ,cell2)+pp(k ,cell1)) &
- +fzp(k)*(pp(k-1,cell2)+pp(k-1,cell1)))
- end do
- dpzx(nVertLevels+1) = 0.
-
- do j=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- do k=1,nVertLevels
- rv(k,iEdge) = rv(k,iEdge) + weightsOnEdge(j,iEdge) * ru(k,eoe)
- end do
- end do
- end do
-
- do iEdge=1,grid % nEdgesSolve
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,nVertLevels
- vorticity_abs = fEdge(iEdge) + (circulation(k,vertex1) + circulation(k,vertex2)) / &
- (areaTriangle(vertex1) + areaTriangle(vertex2))
-
- workpv = 2.0 * vorticity_abs / (rho(k,cell1) + rho(k,cell2))
-
- tend_u(k,iEdge) = rho_edge(k,iEdge)* (workpv * rv(k,iEdge) - (ke(k,cell2) - ke(k,cell1)) / dcEdge(iEdge)) &
- - u(k,iEdge)*0.5*(divergence_ru(k,cell1)+divergence_ru(k,cell2)) &
- - cqu(k,iEdge)*( (pp(k,Cell2)/zz(k,cell2) - pp(k,cell1)/zz(k,cell1)) / dcEdge(iEdge) &
- -rdzw(k)*(dpzx(k+1)-dpzx(k)) )
-
- end do
-
- end do
- deallocate(rv)
-#endif
deallocate(divergence_ru)
!
@@ -2157,6 +2115,10 @@
! horizontal mixing for u
!
+ if (rk_step == 1 .or. rk_diffusion) then
+
+ tend_u_euler = 0.
+
if (delsq_horiz_mixing) then
if (h_mom_eddy_visc2 > 0.0) then
@@ -2176,7 +2138,8 @@
-( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
u_diffusion = rho_edge(k,iEdge)*h_mom_eddy_visc2 * u_diffusion
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+! tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ tend_u_euler(k,iEdge) = tend_u_euler(k,iEdge) + u_diffusion
end do
end do
@@ -2197,7 +2160,8 @@
-( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
u_diffusion = rho_edge(k,iEdge)* 0.5*(kdiff(k,cell1)+kdiff(k,cell2)) * u_diffusion
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+! tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ tend_u_euler(k,iEdge) = tend_u_euler(k,iEdge) + u_diffusion
end do
end do
end if
@@ -2278,7 +2242,8 @@
-( delsq_vorticity(k,vertex2) - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge) &
)
- tend_u(k,iEdge) = tend_u(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
+! tend_u(k,iEdge) = tend_u(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
+ tend_u_euler(k,iEdge) = tend_u_euler(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
end do
end do
@@ -2312,7 +2277,10 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
- tend_u(k,iEdge) = tend_u(k,iEdge) + rho_edge(k,iEdge) * v_mom_eddy_visc2*( &
+! tend_u(k,iEdge) = tend_u(k,iEdge) + rho_edge(k,iEdge) * v_mom_eddy_visc2*( &
+! (u(k+1,iEdge)-u(k ,iEdge))/(zp-z0) &
+! -(u(k ,iEdge)-u(k-1,iEdge))/(z0-zm) )/(0.5*(zp-zm))
+ tend_u_euler(k,iEdge) = tend_u_euler(k,iEdge) + rho_edge(k,iEdge) * v_mom_eddy_visc2*( &
(u(k+1,iEdge)-u(k ,iEdge))/(zp-z0) &
-(u(k ,iEdge)-u(k-1,iEdge))/(z0-zm) )/(0.5*(zp-zm))
end do
@@ -2340,9 +2308,12 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
- tend_u(k,iEdge) = tend_u(k,iEdge) + rho_edge(k,iEdge) * v_mom_eddy_visc2*( &
+ tend_u_euler(k,iEdge) = tend_u_euler(k,iEdge) + rho_edge(k,iEdge) * v_mom_eddy_visc2*( &
(u_mix(k+1)-u_mix(k ))/(zp-z0) &
-(u_mix(k )-u_mix(k-1))/(z0-zm) )/(0.5*(zp-zm))
+! tend_u(k,iEdge) = tend_u(k,iEdge) + rho_edge(k,iEdge) * v_mom_eddy_visc2*( &
+! (u_mix(k+1)-u_mix(k ))/(zp-z0) &
+! -(u_mix(k )-u_mix(k-1))/(z0-zm) )/(0.5*(zp-zm))
end do
end do
@@ -2350,6 +2321,17 @@
end if
+ end if ! (rk_step 1 test for computing mixing terms)
+
+! add in mixing for u
+
+ do iEdge=1,grid % nEdgesSolve
+ do k=1,nVertLevels
+ tend_u(k,iEdge) = tend_u(k,iEdge) + tend_u_euler(k,iEdge)
+ end do
+ end do
+
+
!----------- rhs for w
tend_w(:,:) = 0.
@@ -2451,6 +2433,10 @@
! Note: we are using quite a bit of the theta code here - could be combined later???
+ if (rk_step == 1 .or. rk_diffusion) then
+
+ tend_w_euler = 0.
+
if (delsq_horiz_mixing) then
if (h_mom_eddy_visc2 > 0.0) then
@@ -2463,8 +2449,10 @@
do k=2,grid % nVertLevels
theta_turb_flux = h_mom_eddy_visc2*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
flux = 0.5*dvEdge (iEdge) * (rho_edge(k,iEdge)+rho_edge(k-1,iEdge)) * theta_turb_flux
- tend_w(k,cell1) = tend_w(k,cell1) + flux
- tend_w(k,cell2) = tend_w(k,cell2) - flux
+! tend_w(k,cell1) = tend_w(k,cell1) + flux
+! tend_w(k,cell2) = tend_w(k,cell2) - flux
+ tend_w_euler(k,cell1) = tend_w_euler(k,cell1) + flux/areaCell(cell1)
+ tend_w_euler(k,cell2) = tend_w_euler(k,cell2) - flux/areaCell(cell2)
end do
end if
@@ -2481,8 +2469,10 @@
theta_turb_flux = 0.25*(kdiff(k,cell1)+kdiff(k,cell2)+kdiff(k-1,cell1)+kdiff(k,cell2)) &
*(w(k,cell2) - w(k,cell1))/dcEdge(iEdge)
flux = 0.5*dvEdge (iEdge) * (rho_edge(k,iEdge)+rho_edge(k-1,iEdge)) * theta_turb_flux
- tend_w(k,cell1) = tend_w(k,cell1) + flux
- tend_w(k,cell2) = tend_w(k,cell2) - flux
+! tend_w(k,cell1) = tend_w(k,cell1) + flux
+! tend_w(k,cell2) = tend_w(k,cell2) - flux
+ tend_w_euler(k,cell1) = tend_w_euler(k,cell1) + flux/areaCell(cell1)
+ tend_w_euler(k,cell2) = tend_w_euler(k,cell2) - flux/areaCell(cell2)
end do
end if
@@ -2526,8 +2516,10 @@
theta_turb_flux = h_mom_eddy_visc4*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
- tend_w(k,cell1) = tend_w(k,cell1) - flux
- tend_w(k,cell2) = tend_w(k,cell2) + flux
+! tend_w(k,cell1) = tend_w(k,cell1) - flux
+! tend_w(k,cell2) = tend_w(k,cell2) + flux
+ tend_w_euler(k,cell1) = tend_w_euler(k,cell1) - flux/areaCell(cell1)
+ tend_w_euler(k,cell2) = tend_w_euler(k,cell2) + flux/areaCell(cell2)
end do
end if
@@ -2537,6 +2529,8 @@
end if
+ end if ! horizontal mixing for w computed in first rk_step
+
!
! vertical advection, pressure gradient and buoyancy for w
! Note: we are also dividing through by the cell area after the horizontal flux divergence
@@ -2611,17 +2605,36 @@
!
! vertical mixing for w - 2nd order
!
+
+ if (rk_step == 1 .or. rk_diffusion) then
+
if ( v_mom_eddy_visc2 > 0.0 ) then
do iCell = 1, grid % nCellsSolve
do k=2,nVertLevels
- tend_w(k,iCell) = tend_w(k,iCell) + v_mom_eddy_visc2*0.5*(rho(k,iCell)+rho(k-1,iCell))*( &
+! tend_w(k,iCell) = tend_w(k,iCell) + v_mom_eddy_visc2*0.5*(rho(k,iCell)+rho(k-1,iCell))*( &
+! (w(k+1,iCell)-w(k ,iCell))*rdzw(k) &
+! -(w(k ,iCell)-w(k-1,iCell))*rdzw(k-1) )*rdzu(k)
+
+ tend_w_euler(k,iCell) = tend_w_euler(k,iCell) + v_mom_eddy_visc2*0.5*(rho(k,iCell)+rho(k-1,iCell))*( &
(w(k+1,iCell)-w(k ,iCell))*rdzw(k) &
-(w(k ,iCell)-w(k-1,iCell))*rdzw(k-1) )*rdzu(k)
end do
end do
end if
+
+ end if ! mixing term computed first rk_step
+
+
+! add in mixing terms for w
+
+ do iCell = 1, grid % nCellsSolve
+ do k=2,nVertLevels
+ tend_w(k,iCell) = tend_w(k,iCell) + tend_w_euler(k,iCell)
+ end do
+ end do
+
deallocate(qtot)
!----------- rhs for theta
@@ -2727,12 +2740,15 @@
end do
end if
-! write(0,*) ' pt 1 tend_theta(3,1120) ',tend_theta(3,1120)/AreaCell(1120)
-
!
! horizontal mixing for theta - we could combine this with advection directly (i.e. as a turbulent flux),
! but here we can also code in hyperdiffusion if we wish (2nd order at present)
!
+
+ if (rk_step == 1 .or. rk_diffusion) then
+
+ tend_theta_euler = 0.
+
if (delsq_horiz_mixing) then
if ( h_theta_eddy_visc2 > 0.0 ) then
@@ -2744,8 +2760,10 @@
do k=1,grid % nVertLevels
theta_turb_flux = h_theta_eddy_visc2*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
- tend_theta(k,cell1) = tend_theta(k,cell1) + flux
- tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+! tend_theta(k,cell1) = tend_theta(k,cell1) + flux
+! tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+ tend_theta_euler(k,cell1) = tend_theta_euler(k,cell1) + flux/areaCell(cell1)
+ tend_theta_euler(k,cell2) = tend_theta_euler(k,cell2) - flux/areaCell(cell2)
end do
end if
@@ -2762,8 +2780,10 @@
theta_turb_flux = 0.5*(kdiff(k,cell1)+kdiff(k,cell2))*prandtl &
*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * rho_edge(k,iEdge) * theta_turb_flux
- tend_theta(k,cell1) = tend_theta(k,cell1) + flux
- tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+! tend_theta(k,cell1) = tend_theta(k,cell1) + flux
+! tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+ tend_theta_euler(k,cell1) = tend_theta_euler(k,cell1) + flux/areaCell(cell1)
+ tend_theta_euler(k,cell2) = tend_theta_euler(k,cell2) - flux/areaCell(cell2)
end do
end if
@@ -2803,8 +2823,10 @@
theta_turb_flux = h_theta_eddy_visc4*prandtl*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * theta_turb_flux
- tend_theta(k,cell1) = tend_theta(k,cell1) - flux
- tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+! tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+! tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ tend_theta_euler(k,cell1) = tend_theta_euler(k,cell1) - flux/areaCell(cell1)
+ tend_theta_euler(k,cell2) = tend_theta_euler(k,cell2) + flux/areaCell(cell2)
end do
end if
@@ -2814,6 +2836,8 @@
end if
+ end if ! theta mixing calculated first rk_step
+
!
! vertical advection plus diabatic term
! Note: we are also dividing through by the cell area after the horizontal flux divergence
@@ -2861,6 +2885,9 @@
!
! vertical mixing for theta - 2nd order
!
+
+ if (rk_step == 1 .or. rk_diffusion) then
+
if ( v_theta_eddy_visc2 > 0.0 ) then
if (config_mix_full) then
@@ -2876,7 +2903,10 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
- tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
+! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
+! (theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
+! -(theta(k ,iCell)-theta(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
+ tend_theta_euler(k,iCell) = tend_theta_euler(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
(theta(k+1,iCell)-theta(k ,iCell))/(zp-z0) &
-(theta(k ,iCell)-theta(k-1,iCell))/(z0-zm) )/(0.5*(zp-zm))
end do
@@ -2895,7 +2925,10 @@
z0 = 0.5*(z2+z3)
zp = 0.5*(z3+z4)
- tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
+! tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
+! ((theta(k+1,iCell)-t_init(k+1,iCell))-(theta(k ,iCell)-t_init(k,iCell)))/(zp-z0) &
+! -((theta(k ,iCell)-t_init(k,iCell))-(theta(k-1,iCell)-t_init(k-1,iCell)))/(z0-zm) )/(0.5*(zp-zm))
+ tend_theta_euler(k,iCell) = tend_theta_euler(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
((theta(k+1,iCell)-t_init(k+1,iCell))-(theta(k ,iCell)-t_init(k,iCell)))/(zp-z0) &
-((theta(k ,iCell)-t_init(k,iCell))-(theta(k-1,iCell)-t_init(k-1,iCell)))/(z0-zm) )/(0.5*(zp-zm))
end do
@@ -2905,6 +2938,14 @@
end if
+ end if ! compute theta mixing on first rk_step
+
+ do iCell = 1, grid % nCellsSolve
+ do k=1,nVertLevels
+ tend_theta(k,iCell) = tend_theta(k,iCell) + tend_theta_euler(k,iCell)
+ end do
+ end do
+
end subroutine compute_dyn_tend
!-------
@@ -3121,17 +3162,6 @@
!
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells )
- !
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
- dvEdge(iEdge)
- end do
- end do
-
- !
! Compute pv at the edges
! ( this computes pv_edge at all edges bounding real cells )
!
@@ -3146,16 +3176,6 @@
end do
!
- ! Modify PV edge with upstream bias.
- !
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
- end do
- end do
-
-
- !
! Compute pv at cell centers
! ( this computes pv_cell for all real cells )
!
@@ -3169,7 +3189,23 @@
end do
end do
+
+ if (config_apvm_upwinding > 0.0) then
+
!
+ ! Modify PV edge with upstream bias.
+ !
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradPVt at all edges bounding real cells )
+ !
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+ gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
+ dvEdge(iEdge)
+ end do
+ end do
+
+ !
! Compute gradient of PV in normal direction
! (tdr: 2009-10-02: this is not correct because the pv_cell in the halo is not correct)
!
@@ -3181,6 +3217,12 @@
end do
end do
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ end do
+ end do
+
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
@@ -3189,7 +3231,9 @@
end do
end do
+ end if ! apvm upwinding
+
end subroutine compute_solve_diagnostics
!----------
</font>
</pre>