<p><b>laura@ucar.edu</b> 2010-10-07 12:52:52 -0600 (Thu, 07 Oct 2010)</p><p>reconciled with module_time_integration.F from /branches/atmos_nonhydrostatic version 520, except for stuff related to the reorganized registry, and stuff related to physics<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-10-07 18:50:32 UTC (rev 524)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F 2010-10-07 18:52:52 UTC (rev 525)
@@ -4,9 +4,11 @@
use configure
use constants
use dmpar
+ use vector_reconstruction
#ifdef DO_PHYSICS
use module_microphysics
+ use module_physics_todynamics
#endif
contains
@@ -25,8 +27,8 @@
implicit none
type (domain_type), intent(inout) :: domain
+ integer, intent(in):: itimestep
real (kind=RKIND), intent(in) :: dt
- integer, intent(in) :: itimestep
type (block_type), pointer :: block
@@ -64,14 +66,14 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
- integer, intent(in) :: itimestep
+ integer, intent(in):: itimestep
- integer :: iq
integer :: iCell, k, iEdge
type (block_type), pointer :: block
integer, parameter :: TEND = 1
integer :: rk_step, number_of_sub_steps
+ integer :: iScalar
real (kind=RKIND), dimension(3) :: rk_timestep, rk_sub_timestep
integer, dimension(3) :: number_sub_steps
@@ -79,8 +81,9 @@
logical, parameter :: debug = .false.
! logical, parameter :: debug = .true.
logical, parameter :: debug_mass_conservation = .true.
- logical, parameter :: do_microphysics = .false.
- logical, parameter :: scalar_advection = .true.
+! logical, parameter :: do_microphysics = .true.
+ logical, parameter :: do_microphysics = .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
@@ -174,6 +177,18 @@
end do
if (debug) write(0,*) ' finished compute_dyn_tend '
+
+!#ifdef DO_PHYSICS
+! if (debug) write(0,*) ' add physics tendencies '
+! block => domain % blocklist
+! do while (associated(block))
+! call physics_addtend( block % intermediate_step(TEND), block % time_levs(2) % state, &
+! block % time_levs(2) % state % rho % array(:,:), block % mesh )
+! block => block % next
+! end do
+! if (debug) write(0,*) ' finished add physics tendencies '
+!#endif
+
!***********************************
! 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
@@ -245,8 +260,9 @@
block => domain % blocklist
do while (associated(block))
call recover_large_step_variables( block % time_levs(2) % state, &
- block % mesh, rk_sub_timestep(rk_step), &
- number_sub_steps(rk_step) )
+ block % mesh, rk_timestep(rk_step), &
+ number_sub_steps(rk_step), rk_step )
+
block => block % next
end do
@@ -358,28 +374,38 @@
end do ! rk_step loop
-! microphysics here...
+!... compute full velocity vectors at cell centers:
+ block => domain % blocklist
+ do while (associated(block))
+ call reconstruct(block % time_levs(2) % state, block % mesh)
+ block => block % next
+ end do
+!... call to parameterizations of cloud microphysics:
#ifdef DO_PHYSICS
- if(config_microp_scheme .ne. 'off') then
- block => domain % blocklist
- do while(associated(block))
+ block => domain % blocklist
+ do while(associated(block))
- call microphysics_driver ( block % intermediate_step(TEND), block % time_levs(2) % state, &
- block % mesh, itimestep )
+ !simply set to zero negative mixing ratios of different water species (for now):
+ where ( block % time_levs(2) % state % scalars % array(:,:,:) .lt. 0.) &
+ block % time_levs(2) % state % scalars % array(:,:,:) = 0.
- block => block % next
- end do
- endif
+ !call microphysics schemes:
+ if(config_microp_scheme .ne. 'off') &
+ call microphysics_driver ( block % time_levs(2) % state, block % mesh, itimestep )
+
+ block => block % next
+ end do
#endif
- if(do_microphysics) then
- block => domain % blocklist
- do while (associated(block))
- call qd_kessler( block % time_levs(1) % state, block % time_levs(2) % state, block % mesh, dt )
- block => block % next
- end do
- end if
+! if(do_microphysics) then
+! block => domain % blocklist
+! do while (associated(block))
+! call qd_kessler( block % time_levs(1) % state, block % time_levs(2) % state, 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(:,:), &
@@ -426,16 +452,17 @@
enddo
write(0,*) ' min, max u ',scalar_min, scalar_max
- scalar_min = 0.
- scalar_max = 0.
- do iCell = 1, block % mesh % nCellsSolve
- do k = 1, block % mesh % nVertLevels
- scalar_min = min(scalar_min, block % time_levs(2) % state % scalars % array(index_qc,k,iCell))
- scalar_max = max(scalar_max, block % time_levs(2) % state % scalars % array(index_qc,k,iCell))
+ do iScalar = 1, num_scalars
+ scalar_min = 0.
+ scalar_max = 0.
+ do iCell = 1, block % mesh % nCellsSolve
+ do k = 1, block % mesh % nVertLevels
+ scalar_min = min(scalar_min, block % time_levs(2) % state % scalars % array(iScalar,k,iCell))
+ scalar_max = max(scalar_max, block % time_levs(2) % state % scalars % array(iScalar,k,iCell))
+ enddo
+ enddo
+ write(0,*) ' min, max scalar ',iScalar, scalar_min, scalar_max
enddo
- enddo
- write(0,*) ' min, max qc ',scalar_min, scalar_max
-
block => block % next
end do
@@ -465,9 +492,16 @@
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
!-----
@@ -496,7 +530,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)
@@ -587,11 +621,13 @@
do k=2,nVertLevels
cofwr(k,i) =.5*dtseps*gravity*(fzm(k)*zz(k,i)+fzp(k)*zz(k-1,i))
end do
+ coftz(1,i) = 0.0
do k=2,nVertLevels
cofwz(k,i) = dtseps*c2*(fzm(k)*zz(k,i)+fzp(k)*zz(k-1,i)) &
*rdzu(k)*cqw(k,i)*(fzm(k)*p (k,i)+fzp(k)*p (k-1,i))
coftz(k,i) = dtseps* (fzm(k)*t (k,i)+fzp(k)*t (k-1,i))
end do
+ coftz(nVertLevels+1,i) = 0.0
do k=1,nVertLevels
qtot = 0.
@@ -639,8 +675,20 @@
implicit none
type (grid_state) :: s_new, s_old, tend
type (grid_meta) :: grid
- integer :: iCell, k
+ integer :: iCell, iEdge, k, cell1, cell2
+ integer, dimension(:,:), pointer :: cellsOnEdge
+ real, dimension(:,:,:), pointer :: zf, zf3
+ real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, areaCell
+ real (kind=RKIND) :: flux
+ zf => grid % zf % array
+ zf3 => grid % zf3 % array
+ fzm => grid % fzm % array
+ fzp => grid % fzp % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+
grid % rho_pp % array = grid % rho_p_save % array - s_new % rho_p % array
grid % ru_p % array = grid % ru_save % array - grid % ru % array
@@ -650,12 +698,33 @@
do iCell = 1, grid % nCellsSolve
do k = 2, grid % nVertLevels
- tend % w % array(k,iCell) = ( grid % fzm % array (k) * grid % zz % array(k ,iCell) + &
- grid % fzp % array (k) * grid % zz % array(k-1,iCell) ) &
+ tend % w % array(k,iCell) = ( fzm(k) * grid % zz % array(k ,iCell) + &
+ fzp(k) * grid % zz % array(k-1,iCell) ) &
* tend % w % array(k,iCell)
end do
end do
+ do iEdge = 1,grid % nEdges
+
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k = 2, grid%nVertLevels
+ flux = fzm(k) * tend % u % array(k,iEdge) + fzp(k) * tend % u % array(k-1,iEdge)
+ tend % w % array(k,cell2) = tend % w % array(k,cell2) + zf(k,2,iEdge)*flux
+ tend % w % array(k,cell1) = tend % w % array(k,cell1) - zf(k,1,iEdge)*flux
+!3rd order stencil
+ if (config_theta_adv_order == 3) then
+ tend % w % array(k,cell2) = tend % w % array(k,cell2) + sign(1.,tend % u % array(k,iEdge)) &
+ *config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ tend % w % array(k,cell1) = tend % w % array(k,cell1) - sign(1.,tend % u % array(k,iEdge)) &
+ *config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ end if
+
+ end do
+
+ end do
+
grid % ruAvg % array = 0.
grid % wwAvg % array = 0.
@@ -683,7 +752,6 @@
real (kind=RKIND), dimension( grid % nVertLevels ) :: du
real (kind=RKIND), dimension( grid % nVertLevels + 1 ) :: dpzx
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
@@ -759,7 +827,6 @@
ts = 0.
rs = 0.
- ws = 0.
! acoustic step divergence damping - forward weight rtheta_pp
rtheta_pp_old = rtheta_pp + smdiv*(rtheta_pp - rtheta_pp_old)
@@ -817,14 +884,6 @@
end do
- do k=2,nVertLevels
- flux = dts*0.5*dvEdge(iEdge)*((zgrid(k,cell2)-zgrid(k,cell1))*(fzm(k)*du(k)+fzp(k)*du(k-1)) )
- flux2 = - (fzm(k)*zz(k ,cell2) +fzp(k)*zz(k-1,cell2))*flux/AreaCell(cell2)
- flux1 = - (fzm(k)*zz(k ,cell1) +fzp(k)*zz(k-1,cell1))*flux/AreaCell(cell1)
- ws(k,cell2) = ws(k,cell2) + flux2
- ws(k,cell1) = ws(k,cell1) + flux1
- enddo
-
end if ! end test for block-owned cells
end do ! end loop over edges
@@ -846,7 +905,7 @@
wwavg(k,iCell) = wwavg(k,iCell) + 0.5*(1.-epssm)*rw_p(k,iCell)
- rw_p(k,iCell) = rw_p(k,iCell) + ws(k,iCell) + dts*tend_rw(k,iCell) &
+ rw_p(k,iCell) = rw_p(k,iCell) + dts*tend_rw(k,iCell) &
- cofwz(k,iCell)*((zz(k ,iCell)*ts (k ,iCell) &
-zz(k-1,iCell)*ts (k-1,iCell)) &
+resm*(zz(k ,iCell)*rtheta_pp(k ,iCell) &
@@ -887,21 +946,22 @@
!------------------------
- subroutine recover_large_step_variables( s, grid, dt, ns )
+ subroutine recover_large_step_variables( s, grid, dt, ns, rk_step )
implicit none
type (grid_state) :: s
type (grid_meta) :: grid
- integer, intent(in) :: ns
+ integer, intent(in) :: ns, rk_step
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, rho_p, rho_p_save, &
+ rtheta_p_save, rt_diabatic_tend, 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
+ zz, theta
real (kind=RKIND), dimension(:), pointer :: fzm, fzp, dvEdge, AreaCell
- integer, dimension(:,:), pointer :: CellsOnEdge
+ real (kind=RKIND), dimension(:,:,:), pointer :: zb, zb3
+ integer, dimension(:,:), pointer :: cellsOnEdge
integer :: iCell, iEdge, k, cell1, cell2
integer :: nVertLevels, nCells, nCellsSolve, nEdges, nEdgesSolve
@@ -922,7 +982,7 @@
rtheta_p_save => grid % rtheta_p_save % array
rtheta_pp => grid % rtheta_pp % array
rtheta_base => grid % rtheta_base % array
- rt_diabatic => grid % rt_diabatic_tend % array
+ rt_diabatic_tend => grid % rt_diabatic_tend % array
theta => s % theta % array
rho => s % rho % array
@@ -943,7 +1003,8 @@
pressure_p => s % pressure % array
zz => grid % zz % array
- zgrid => grid % zgrid % array
+ zb => grid % zb % array
+ zb3 => grid % zb3 % array
fzm => grid % fzm % array
fzp => grid % fzp % array
dvEdge => grid % dvEdge % array
@@ -1021,8 +1082,14 @@
do k = 1, nVertLevels
- rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) - dt*rho(k,iCell)*rt_diabatic(k,iCell)
+ if (rk_step==3) then
+ rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell) &
+ - dt * rho(k,iCell) * rt_diabatic_tend(k,iCell)
+ else
+ rtheta_p(k,iCell) = rtheta_p(k,iCell) + rtheta_pp(k,iCell)
+ end if
+
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????
@@ -1053,14 +1120,29 @@
u(k,iEdge) = 2.*ru(k,iEdge)/(rho(k,cell1)+rho(k,cell2))
enddo
- flux = dvEdge(iEdge)*0.5*(cf1*u(1,iEdge)+cf2*u(2,iEdge)+cf3*u(3,iEdge))*(zgrid(1,cell2)-zgrid(1,cell1))
- w(1,cell2) = w(1,cell2)+flux/AreaCell(cell2)
- w(1,cell1) = w(1,cell1)+flux/AreaCell(cell1)
+ !SHP-mtn
+ flux = cf1*ru(1,iEdge) + cf2*ru(2,iEdge) + cf3*ru(3,iEdge)
+ w(1,cell2) = w(1,cell2) - zb(1,2,iEdge)*flux/(cf1*rho(1,cell2)+cf2*rho(2,cell2)+cf3*rho(3,cell2))
+ w(1,cell1) = w(1,cell1) + zb(1,1,iEdge)*flux/(cf1*rho(1,cell1)+cf2*rho(2,cell1)+cf3*rho(3,cell1))
+!3rd order stencil
+ if (config_theta_adv_order == 3) then
+ w(1,cell2) = w(1,cell2) - sign(1.,ru(1,iEdge))*config_coef_3rd_order &
+ *zb3(1,2,iEdge)*flux/(cf1*rho(1,cell2)+cf2*rho(2,cell2)+cf3*rho(3,cell2))
+ w(2,cell1) = w(2,cell1) + sign(1.,ru(1,iEdge))*config_coef_3rd_order &
+ *zb3(1,1,iEdge)*flux/(cf1*rho(1,cell1)+cf2*rho(2,cell1)+cf3*rho(3,cell1))
+ end if
do k = 2, nVertLevels
- flux = dvEdge(iEdge)*0.5*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))*(zgrid(k,cell2)-zgrid(k,cell1))
- w(k,cell2) = w(k,cell2)+flux/AreaCell(cell2)
- w(k,cell1) = w(k,cell1)+flux/AreaCell(cell1)
+ flux = (fzm(k)*ru(k,iEdge)+fzp(k)*ru(k-1,iEdge))
+ w(k,cell2) = w(k,cell2) - zb(k,2,iEdge)*flux/(fzm(k)*rho(k,cell2)+fzp(k)*rho(k-1,cell2))
+ w(k,cell1) = w(k,cell1) + zb(k,1,iEdge)*flux/(fzm(k)*rho(k,cell1)+fzp(k)*rho(k-1,cell1))
+!3rd order! stencil
+ if (config_theta_adv_order ==3) then
+ w(k,cell2) = w(k,cell2) - sign(1.,ru(k,iEdge))*config_coef_3rd_order &
+ *zb3(k,2,iEdge)*flux/(fzm(k)*rho(k,cell2)+fzp(k)*rho(k-1,cell2))
+ w(k,cell1) = w(k,cell1) + sign(1.,ru(k,iEdge))*config_coef_3rd_order &
+ *zb3(k,1,iEdge)*flux/(fzm(k)*rho(k,cell1)+fzp(k)*rho(k-1,cell1))
+ end if
enddo
! end if
@@ -1113,7 +1195,7 @@
flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
- coef3*ua*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+ coef3*abs(ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
! coef_3rd_order = 0.
! if (config_scalar_adv_order == 3) coef_3rd_order = 1.0
@@ -1449,7 +1531,7 @@
flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
- coef3*ua*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+ coef3*abs(ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
!------
@@ -1812,7 +1894,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, &
- rt_diabatic, tend_theta, tend_w, w, cqw, rb, rr, pp, pressure_b, zz, zx, cqu, &
+ 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 :: deriv_two
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
@@ -1823,7 +1905,8 @@
real (kind=RKIND) :: theta_edge, theta_turb_flux, z1, z2, z3, z4, zm, z0, zp, r
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2, pgrad
- real (kind=RKIND), dimension(:), pointer :: rdzu, rdzw, fzm, fzp, t_init
+ real (kind=RKIND), dimension(:), pointer :: rdzu, rdzw, fzm, fzp, qv_init
+ real (kind=RKIND), dimension(:,:), pointer :: t_init
real (kind=RKIND), allocatable, dimension(:,:) :: rv, divergence_ru, qtot
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_theta, delsq_divergence
@@ -1850,7 +1933,7 @@
flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
- coef3*ua*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
+ coef3*abs(ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
!-----------
@@ -1900,9 +1983,10 @@
tend_theta => tend % theta % array
tend_w => tend % w % array
tend_rho => tend % rho % array
- rt_diabatic => grid % rt_diabatic_tend % array
+ rt_diabatic_tend => grid % rt_diabatic_tend % array
t_init => grid % t_init % array
+ qv_init => grid % qv_init % array
rdzu => grid % rdzu % array
rdzw => grid % rdzw % array
@@ -2115,12 +2199,6 @@
wduz(nVertLevels+1) = 0.
- wduz(1) = 0.
- do k=2,nVertLevels
- wduz(k) = 0.5*( rw(k,cell1)+rw(k,cell2) )*(fzm(k)*u(k,iEdge)+fzp(k)*u(k-1,iEdge))
- end do
- wduz(nVertLevels+1) = 0.
-
do k=1,nVertLevels
tend_u(k,iEdge) = tend_u(k,iEdge) - rdzw(k)*(wduz(k+1)-wduz(k))
end do
@@ -2564,23 +2642,25 @@
do k=2,nVertLevels
tend_w(k,iCell) = tend_w(k,iCell)/areaCell(iCell) -rdzu(k)*(wdwz(k+1)-wdwz(k)) &
+!SHP-w
- cqw(k,iCell)*( rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
+ gravity* &
-!shpark
- ( fzm(k)*rr(k,iCell) + fzm(k)*(rb(k,iCell)+rr(k,iCell))*qtot(k,iCell) &
- +fzp(k)*rr(k-1,iCell) + fzp(k)*(rb(k-1,iCell)+rr(k-1,iCell))*qtot(k-1,iCell) ))
-
+ ( fzm(k)*(rb(k,iCell)*(qtot(k,iCell)) + &
+ rr(k,iCell)*(1.+qtot(k,iCell))) &
+ +fzp(k)*(rb(k-1,iCell)*(qtot(k-1,iCell)) + &
+ rr(k-1,iCell)*(1.+qtot(k-1,iCell))) ))
+!SHP-old version
+! ( fzm(k)*rr(k,iCell) + fzm(k)*(rb(k,iCell)+rr(k,iCell))*qtot(k,iCell) &
+! +fzp(k)*rr(k-1,iCell) + fzp(k)*(rb(k-1,iCell)+rr(k-1,iCell))*qtot(k-1,iCell) ))
+
+
! - gravity*(fzm(k)*rb(k,iCell)+fzp(k)*rb(k-1,iCell)) ) &
! - gravity*( fzm(k)*(rr(k,iCell)+rb(k,iCell)) + fzp(k)*(rr(k-1,iCell)+rb(k-1,iCell)) )
-
-
! - cqw(k,iCell)*rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
! - gravity*( fzm(k)*rr(k,iCell)+fzp(k)*rr(k-1,iCell) &
! +(1.-cqw(k,iCell))*(fzm(k)*rb(k,iCell)+fzp(k)*rb(k-1,iCell)))
-
-
! WCS version - cqw(k,iCell)*rdzu(k)*(pp(k,iCell)-pp(k-1,iCell)) &
! - gravity*0.5*(rr(k,iCell)+rr(k-1,iCell)+(1.-cqw(k,iCell))*(rb(k,iCell)+rb(k-1,iCell)))
@@ -2598,7 +2678,7 @@
if ( v_mom_eddy_visc2 > 0.0 ) then
do iCell = 1, grid % nCellsSolve
- do k=2,nVertLevels-1
+ 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))*( &
(w(k+1,iCell)-w(k ,iCell))*rdzw(k) &
-(w(k ,iCell)-w(k-1,iCell))*rdzw(k-1) )*rdzu(k)
@@ -2842,7 +2922,7 @@
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,iCell) + rho(k,iCell)*rt_diabatic(k,iCell)
+ tend_theta(k,iCell) = tend_theta(k,iCell) + rho(k,iCell)*rt_diabatic_tend(k,iCell)
end do
end do
@@ -2884,8 +2964,8 @@
zp = 0.5*(z3+z4)
tend_theta(k,iCell) = tend_theta(k,iCell) + v_theta_eddy_visc2*prandtl*rho(k,iCell)*(&
- ((theta(k+1,iCell)-t_init(k+1))-(theta(k ,iCell)-t_init(k)))/(zp-z0) &
- -((theta(k ,iCell)-t_init(k))-(theta(k-1,iCell)-t_init(k-1)))/(z0-zm) )/(0.5*(zp-zm))
+ ((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
end do
@@ -3190,12 +3270,6 @@
enddo
enddo
- do i=1,grid%nCellsSolve
- do k=1,grid % nVertLevels + 1
- grid % rw % array (k,i) = 0.
- enddo
- enddo
-
end subroutine init_coupled_diagnostics
! ------------------------
@@ -3242,11 +3316,7 @@
do k = 1, grid % nVertLevels
- 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) * &
-! (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) &
</font>
</pre>