<p><b>laura@ucar.edu</b> 2010-10-07 16:05:31 -0600 (Thu, 07 Oct 2010)</p><p>time_integration is now completely consistent with trunk, except for stuff related to reorganization of Registry<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F 2010-10-07 22:01:57 UTC (rev 530)
+++ branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F 2010-10-07 22:05:31 UTC (rev 531)
@@ -7,19 +7,14 @@
use vector_reconstruction
#ifdef DO_PHYSICS
- use module_microphysics_driver
+ use module_microphysics
use module_physics_todynamics
#endif
contains
-#ifdef DO_PHYSICS
- subroutine timestep(domain, dt, itimestep, config_ntimesteps)
-#else
- subroutine timestep(domain, dt)
-#endif
-
+ subroutine timestep(domain, dt, itimestep)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step
!
@@ -32,20 +27,13 @@
implicit none
type (domain_type), intent(inout) :: domain
+ integer, intent(in):: itimestep
real (kind=RKIND), intent(in) :: dt
type (block_type), pointer :: block
-#ifdef DO_PHYSICS
- integer, intent(in):: itimestep, config_ntimesteps
-#endif
-
if (trim(config_time_integration) == 'SRK3') then
-#ifdef DO_PHYSICS
- call srk3(domain, dt, itimestep, config_ntimesteps)
-#else
- call srk3(domain, dt)
-#endif
+ call srk3(domain, dt, itimestep)
else
write(0,*) 'Unknown time integration option '//trim(config_time_integration)
write(0,*) 'Currently, only ''SRK3'' is supported.'
@@ -61,11 +49,7 @@
end subroutine timestep
-#ifdef DO_PHYSICS
- subroutine srk3(domain, dt, itimestep , config_ntimesteps)
-#else
- subroutine srk3(domain, dt)
-#endif
+ subroutine srk3(domain, dt, itimestep)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step using
! time-split RK3 scheme
@@ -82,12 +66,14 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
+ integer, intent(in):: itimestep
integer :: iCell, k
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
@@ -98,10 +84,6 @@
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
-#ifdef DO_PHYSICS
- integer, intent(in):: itimestep, config_ntimesteps
-#endif
-
!
! Initialize time_levs(2) with state at current time
! Initialize RK weights
@@ -176,15 +158,26 @@
block => domain % blocklist
do while (associated(block))
-#ifdef DO_PHYSICS
-! call physics_addtend( block % intermediate_step(TEND), block % time_levs(2) % state, block % mesh )
-#endif
call compute_dyn_tend( block % intermediate_step(TEND), block % time_levs(2) % state, block % mesh )
block => block % next
end do
if(debug) write(0,*) ' returned from 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 % h % array, block % mesh )
+! 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)
+! block => block % next
+! end do
+! if (debug) write(0,*) ' finished add physics tendencies '
+!#endif
+
!
! --- update halos for tendencies
!
@@ -195,7 +188,7 @@
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)
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
@@ -224,8 +217,8 @@
block => domain % blocklist
do while (associated(block))
call advance_dynamics( block % intermediate_step(TEND), block % time_levs(2) % state, &
- block % mesh, &
- small_step, number_sub_steps(rk_step), rk_sub_timestep(rk_step) )
+ block % mesh,small_step, number_sub_steps(rk_step), &
+ rk_sub_timestep(rk_step), rk_step )
block => block % next
end do
@@ -286,7 +279,6 @@
if(debug) write(0,*) ' advance scalars '
-
! --- advance scalars with time integrated mass fluxes
block => domain % blocklist
@@ -297,15 +289,18 @@
! so we keep the advance_scalars routine as well
!
if (rk_step < 3 .or. (.not. config_monotonic .and. .not. config_positive_definite)) then
+
call advance_scalars( block % intermediate_step(TEND), &
block % time_levs(1) % state, block % time_levs(2) % state, &
block % mesh, rk_timestep(rk_step) )
else
+
call advance_scalars_mono( block % intermediate_step(TEND), &
block % time_levs(1) % state, block % time_levs(2) % state, &
block % mesh, rk_timestep(rk_step), rk_step, 3, &
domain % dminfo, block % parinfo % cellsToSend, block % parinfo % cellsToRecv )
end if
+
block => block % next
end do
@@ -353,7 +348,6 @@
block => block % next
end do
-
if(debug) write(0,*) ' rk step complete - mass diagnostics '
if(debug .or. debug_mass_conservation) then
@@ -368,6 +362,7 @@
block % mesh % areaCell % array (iCell) &
- block % time_levs(2) % state % pressure % array (block % mesh % nVertLevels + 1, 1) * &
block % mesh % areaCell % array (iCell)
+
do k=1, block % mesh % nVertLevelsSolve
scalar_mass = scalar_mass - block % time_levs(2) % state % scalars % array (2,k,iCell) * &
block % time_levs(2) % state % h % array (k,iCell) * &
@@ -387,82 +382,32 @@
write(0,*) ' scalar mass in the domain = ',global_scalar_mass
write(0,*) ' scalar_min, scalar_max ',global_scalar_min, global_scalar_max
end if
-
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! LDF begin (04-26-2010):
- ! CALL TO CLOUD MICROPHYISCS SCHEMES STARTS HERE:
- ! LDF end.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+
#ifdef DO_PHYSICS
- if(config_mp_physics /=0) then
-
+ !.. CALL TO CLOUD MICROPHYSICS DRIVER STARTS HERE:
block => domain % blocklist
do while(associated(block))
- call microphysics_driver ( block % intermediate_step(TEND), block % time_levs(2) % state, &
- block % mesh, itimestep , config_ntimesteps )
- !to be moved into a module:
- do iCell = 1, block % mesh % nCellsSolve
- block % time_levs(2) % state % qv_col % array(iCell) = 0.
- block % time_levs(2) % state % qc_col % array(iCell) = 0.
- block % time_levs(2) % state % qr_col % array(iCell) = 0.
- block % time_levs(2) % state % qi_col % array(iCell) = 0.
- block % time_levs(2) % state % qs_col % array(iCell) = 0.
+ !call microphysics schemes:
+ if(config_microp_scheme .ne. 'off') &
+ call microphysics_driver ( block % time_levs(2) % state, block % mesh, itimestep )
- do k=1, block % mesh % nVertLevelsSolve
- block % time_levs(2) % state % qv_col % array(iCell) = &
- block % time_levs(2) % state % qv_col % array(iCell) &
- - block % time_levs(2) % state % scalars % array (index_qv,k,iCell) &
- * block % time_levs(2) % state % h % array (k,iCell) &
- * block % mesh % dnw % array (k)
-
- block % time_levs(2) % state % qc_col % array(iCell) = &
- block % time_levs(2) % state % qc_col % array(iCell) &
- - block % time_levs(2) % state % scalars % array (index_qc,k,iCell) &
- * block % time_levs(2) % state % h % array (k,iCell) &
- * block % mesh % dnw % array (k)
-
- block % time_levs(2) % state % qr_col % array(iCell) = &
- block % time_levs(2) % state % qr_col % array(iCell) &
- - block % time_levs(2) % state % scalars % array (index_qr,k,iCell) &
- * block % time_levs(2) % state % h % array (k,iCell) &
- * block % mesh % dnw % array (k)
-
- block % time_levs(2) % state % qi_col % array(iCell) = &
- block % time_levs(2) % state % qi_col % array(iCell) &
- - block % time_levs(2) % state % scalars % array (index_qi,k,iCell) &
- * block % time_levs(2) % state % h % array (k,iCell) &
- * block % mesh % dnw % array (k)
-
- block % time_levs(2) % state % qs_col % array(iCell) = &
- block % time_levs(2) % state % qs_col % array(iCell) &
- - block % time_levs(2) % state % scalars % array (index_qs,k,iCell) &
- * block % time_levs(2) % state % h % array (k,iCell) &
- * block % mesh % dnw % array (k)
-
- block % time_levs(2) % state % qg_col % array(iCell) = &
- block % time_levs(2) % state % qg_col % array(iCell) &
- - block % time_levs(2) % state % scalars % array (index_qg,k,iCell) &
- * block % time_levs(2) % state % h % array (k,iCell) &
- * block % mesh % dnw % array (k)
+ 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
block => block % next
end do
-
- endif
#endif
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! LDF begin (04-26-2010):
- ! CALL TO CLOUD MICROPHYISCS SCHEMES ENDS HERE:
- ! LDF end.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
end subroutine srk3
!------------------------------------------------------------------------------------------------------------------
@@ -680,7 +625,6 @@
vh => tend % vh % array
tend_u => tend % u % array
tend_theta => tend % theta % array
-! h_diabatic => grid % h_diabatic % array
ww => s % ww % array
rdnu => grid % rdnu % array
@@ -920,10 +864,7 @@
!----------- rhs for theta
-!#ifndef DO_PHYSICS
tend_theta(:,:) = 0.
-!#endif
-
!
! 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)
@@ -1080,7 +1021,9 @@
wdtn(nVertLevels+1) = 0.
do k=1,nVertLevels
tend_theta(k,iCell) = tend_theta(k,iCell)/areaCell(iCell) -rdnw(k)*(wdtn(k+1)-wdtn(k))
-!! tend_theta(k,iCell) = tend_theta(k) + h(k,iCell)*h_diabatic(k,iCell)
+
+ !ldf added theta tendency due to cloud microphysics only (09-03-2010)::
+ tend_theta(k,iCell) = tend_theta(k,iCell) + h(k,iCell)*h_diabatic(k,iCell)
end do
end do
@@ -1113,7 +1056,7 @@
!---------------------------------------------------------------------------------------------------------
- subroutine advance_dynamics(tend, s, grid, small_step, number_small_steps, dt)
+ subroutine advance_dynamics(tend, s, grid, small_step, number_small_steps, dt, rk_step)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance the dry dynamics a small timestep (forward-backward integration)
!
@@ -1131,7 +1074,7 @@
type (grid_state), intent(inout) :: s
type (grid_meta), intent(in) :: grid
real (kind=RKIND), intent(in) :: dt
- integer, intent(in) :: small_step, number_small_steps
+ integer, intent(in) :: small_step, number_small_steps, rk_step
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
@@ -1142,7 +1085,7 @@
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, geopotential, alpha, theta, &
pressure, pressure_old, tend_theta, uhAvg, wwAvg, ww, u_old, &
- theta_old, h_edge_old, qtot, ww_old, cqu, h_old
+ theta_old, h_edge_old, qtot, ww_old, cqu, h_old, h_diabatic
real (kind=RKIND), dimension(:,:,:), pointer :: scalar
! real (kind=RKIND), pointer :: smext, p0, ptop
@@ -1187,7 +1130,7 @@
tend_h => tend % h % array
tend_u => tend % u % array
tend_theta => tend % theta % array
-
+ h_diabatic => s % h_diabatic % array
uhAvg => grid % uhAvg % array
wwAvg => grid % wwAvg % array
@@ -1355,15 +1298,27 @@
end do
! if last small step of a set, decouple theta
+!ldf (09-20-2010): modified to include the diabatic heating due to cloud mcirophysics processes:
if(small_step == number_small_steps) then
- do iCell=1,grid % nCells
- do k=1,nVertLevels
- theta(k,iCell) = theta(k,iCell)/h(k,iCell)
- end do
- end do
- uhAvg = uhAvg/real(number_small_steps)
- wwAvg = wwAvg/real(number_small_steps)
+
+ if(rk_step < 3) then
+ do iCell=1,grid % nCells
+ do k=1,nVertLevels
+ theta(k,iCell) = theta(k,iCell)/h(k,iCell)
+ end do
+ end do
+ else
+ do iCell=1,grid % nCells
+ do k=1,nVertLevels
+ theta(k,iCell) = (theta(k,iCell)-dt*number_small_steps*h_diabatic(k,iCell)) &
+ / h(k,iCell)
+ end do
+ end do
+ endif
+
+ uhAvg = uhAvg/real(number_small_steps)
+ wwAvg = wwAvg/real(number_small_steps)
end if
end subroutine advance_dynamics
@@ -1424,9 +1379,7 @@
nVertLevels = grid % nVertLevels
-#ifndef DO_PHYSICS
scalar_tend = 0. ! testing purposes - we have no sources or sinks
-#endif
!
! Runge Kutta integration, so we compute fluxes from scalar_new values, update starts form scalar_old
!
</font>
</pre>