<p><b>laura@ucar.edu</b> 2010-06-28 13:27:06 -0600 (Mon, 28 Jun 2010)</p><p>Added latest changes to dynamical core<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-06-28 19:26:56 UTC (rev 365)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2010-06-28 19:27:06 UTC (rev 366)
@@ -74,6 +74,7 @@
! logical, parameter :: debug = .true.
logical, parameter :: debug_mass_conservation = .true.
logical, parameter :: do_microphysics = .false.
+ logical, parameter :: scalar_advection = .false.
real (kind=RKIND) :: domain_mass, scalar_mass, scalar_min, scalar_max
real (kind=RKIND) :: global_domain_mass, global_scalar_mass, global_scalar_min, global_scalar_max
@@ -173,6 +174,9 @@
! ************ advection of moist variables here...
+
+ if(scalar_advection) then
+
block => domain % blocklist
do while (associated(block))
!
@@ -193,6 +197,12 @@
block => block % next
end do
+ else
+
+ write(0,*) ' no scalar advection '
+
+ end if
+
block => domain % blocklist
do while (associated(block))
call compute_solve_diagnostics( dt, block % time_levs(2) % state, block % mesh )
@@ -489,7 +499,7 @@
wwAvg, rho_pp, cofwt, coftz, zx, &
a_tri, alpha_tri, gamma_tri, dss, &
tend_ru, tend_rho, tend_rt, tend_rw, &
- zgrid, cofwr, cofwz, w
+ zgrid, cofwr, cofwz, w, h_divergence
real (kind=RKIND), dimension(:), pointer :: fzm, fzp, rdzw, dcEdge, AreaCell, cofrz, dvEdge
real (kind=RKIND) :: smdiv, c2, rcv
@@ -519,6 +529,7 @@
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
@@ -608,7 +619,8 @@
pgrad = (rtheta_pp_old(k,cell2)-rtheta_pp_old(k,cell1))/dcEdge(iEdge) &
- rdzw(k)*(dpzx(k+1)-dpzx(k))
pgrad = 0.5*c2*(exner(k,cell1)+exner(k,cell2))*pgrad
- du(k) = dts*(tend_ru(k,iEdge) - cqu(k,iEdge) * pgrad)
+ du(k) = dts*(tend_ru(k,iEdge) - cqu(k,iEdge) * pgrad)
+! + (0.005/6.)*dcEdge(iEdge)*(h_divergence(k,cell2)-h_divergence(k,cell1))
ru_p(k,iEdge) = ru_p(k,iEdge) + du(k)
@@ -1546,7 +1558,8 @@
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
+ h_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
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -1568,7 +1581,10 @@
logical, parameter :: debug = .false.
logical, parameter :: mix_full = .false.
! logical, parameter :: mix_full = .true.
+ integer :: w_adv_order
+ real (kind=RKIND) :: coef_3rd_order
+
rho => s % rho % array
rho_edge => s % rho_edge % array
rb => grid % rho_base % array
@@ -1585,7 +1601,9 @@
pv_edge => s % pv_edge % array
pp => s % pressure % array
pressure_b => grid % pressure_base % array
+ h_divergence => grid % h_divergence % array
+
weightsOnEdge => grid % weightsOnEdge % array
cellsOnEdge => grid % cellsOnEdge % array
verticesOnEdge => grid % verticesOnEdge % array
@@ -1645,9 +1663,14 @@
allocate(qtot(nVertLevels, nCells))
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
@@ -1660,6 +1683,7 @@
r = 1.0 / areaCell(iCell)
do k = 1,nVertLevels
divergence_ru(k,iCell) = divergence_ru(k,iCell) * r
+ h_divergence(k,iCell) = divergence_ru(k,iCell)
tend_rho(k,iCell) = -divergence_ru(k,iCell)-rdzw(k)*(rw(k+1,iCell)-rw(k,iCell))
do iq = moist_start, moist_end
@@ -1748,7 +1772,7 @@
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)) )
+ -rdzw(k)*(dpzx(k+1)-dpzx(k)) )
end do
@@ -1967,8 +1991,10 @@
! horizontal advection for w
!
- if (config_theta_adv_order == 2) then
+ w_adv_order = 2
+ if (w_adv_order == 2) then
+
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -1982,7 +2008,7 @@
end if
end do
- else if (config_theta_adv_order == 3) then
+ else if (w_adv_order == 3) then
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -2020,7 +2046,7 @@
end if
end do
- else if (config_theta_adv_order == 4) then
+ else if (w_adv_order == 4) then
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -2227,14 +2253,26 @@
end do
! 3rd order stencil
+
+ coef_3rd_order = 0.25
+
if( u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
- else
- flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
- 0.5*(theta(k,cell1) + theta(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
+ flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ else
+ flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+! flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
+! 0.5*(theta(k,cell1) + theta(k,cell2)) &
+! -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
+! else
+! flux = dvEdge(iEdge) * ru(k,iEdge) * ( &
+! 0.5*(theta(k,cell1) + theta(k,cell2)) &
+! -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
end if
tend_theta(k,cell1) = tend_theta(k,cell1) - flux
</font>
</pre>