<p><b>duda</b> 2009-09-04 16:10:17 -0600 (Fri, 04 Sep 2009)</p><p>Switch to a loop-based RK4 implementation with a storage factor of <br>
two; from Bill S.<br>
<br>
M src/module_time_integration.F<br>
M src/module_grid_types.F<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/src/module_grid_types.F
===================================================================
--- trunk/swmodel/src/module_grid_types.F 2009-09-04 17:07:13 UTC (rev 45)
+++ trunk/swmodel/src/module_grid_types.F 2009-09-04 22:10:17 UTC (rev 46)
@@ -3,7 +3,7 @@
use dmpar
integer, parameter :: nTimeLevs = 2
- integer, parameter :: storageFactor = 5 ! Additional storage used by time integration scheme
+ integer, parameter :: storageFactor = 2 ! Additional storage used by time integration scheme
! Derived type describing info for doing I/O specific to a field
Modified: trunk/swmodel/src/module_time_integration.F
===================================================================
--- trunk/swmodel/src/module_time_integration.F 2009-09-04 17:07:13 UTC (rev 45)
+++ trunk/swmodel/src/module_time_integration.F 2009-09-04 22:10:17 UTC (rev 46)
@@ -62,198 +62,137 @@
integer :: iCell, k
type (block_type), pointer :: block
- integer, parameter :: TEMP = 1
- integer, parameter :: Q1 = 2
- integer, parameter :: Q2 = 3
- integer, parameter :: Q3 = 4
- integer, parameter :: Q4 = 5
+ integer, parameter :: PROVIS = 1
+ integer, parameter :: TEND = 2
+ integer :: rk_step
+ real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
- ! Take first RK step
- block => domain % blocklist
- do while (associated(block))
- call compute_tend(block % intermediate_step(Q1), block % time_levs(1) % state, block % mesh)
- call compute_scalar_tend(block % intermediate_step(Q1), block % time_levs(1) % state, block % mesh)
- block % intermediate_step(Q1) % u % array(:,:) = dt * block % intermediate_step(Q1) % u % array(:,:)
- block % intermediate_step(Q1) % h % array(:,:) = dt * block % intermediate_step(Q1) % h % array(:,:)
- block % intermediate_step(Q1) % tracers % array(:,:,:) = dt * block % intermediate_step(Q1) % tracers % array(:,:,:)
- block => block % next
- end do
-
+ !
+ ! Initialize time_levs(2) with state at current time
+ ! Initialize first RK state
+ ! Couple tracers time_levs(2) with h in time-levels
+ ! Initialize RK weights
+ !
block => domain % blocklist
do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q1) % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q1) % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(Q1) % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
-
- ! Take second RK step
- block => domain % blocklist
- do while (associated(block))
- block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:) + &
- block % intermediate_step(Q1) % u % array(:,:)/2.0
- block % intermediate_step(TEMP) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) + &
- block % intermediate_step(Q1) % h % array(:,:)/2.0
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % intermediate_step(TEMP) % tracers % array(:,k,iCell) = ( &
- block % time_levs(1) % state % h % array(k,iCell) * &
- block % time_levs(1) % state % tracers % array(:,k,iCell) + &
- block % intermediate_step(Q1) % tracers % array(:,k,iCell) / 2.0 &
- ) / block % intermediate_step(TEMP) % h % array(k,iCell)
+ block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
+ block % time_levs(2) % state % h % array(:,:) = block % time_levs(1) % state % h % array(:,:)
+ do iCell=1,block % mesh % nCells ! couple tracers to h
+ do k=1,block % mesh % nVertLevels
+ block % time_levs(2) % state % tracers % array(:,k,iCell) = block % time_levs(1) % state % tracers % array(:,k,iCell) &
+ * block % time_levs(1) % state % h % array(k,iCell)
end do
end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
- end if
- call compute_solve_diagnostics(block % intermediate_step(TEMP), block % mesh)
- call compute_tend(block % intermediate_step(Q2), block % intermediate_step(TEMP), block % mesh)
- call compute_scalar_tend(block % intermediate_step(Q2), block % intermediate_step(TEMP), block % mesh)
- block % intermediate_step(Q2) % u % array(:,:) = dt * block % intermediate_step(Q2) % u % array(:,:)
- block % intermediate_step(Q2) % h % array(:,:) = dt * block % intermediate_step(Q2) % h % array(:,:)
- block % intermediate_step(Q2) % tracers % array(:,:,:) = dt * block % intermediate_step(Q2) % tracers % array(:,:,:)
- block => block % next
- end do
+ call copy_state(block % time_levs(1) % state, block % intermediate_step(PROVIS))
- block => domain % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q2) % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q2) % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(Q2) % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
+ rk_weights(1) = dt/6.
+ rk_weights(2) = dt/3.
+ rk_weights(3) = dt/3.
+ rk_weights(4) = dt/6.
- ! Take third RK step
- block => domain % blocklist
- do while (associated(block))
- block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:) + &
- block % intermediate_step(Q2) % u % array(:,:)/2.0
- block % intermediate_step(TEMP) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) + &
- block % intermediate_step(Q2) % h % array(:,:)/2.0
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % intermediate_step(TEMP) % tracers % array(:,k,iCell) = ( &
- block % time_levs(1) % state % h % array(k,iCell) * &
- block % time_levs(1) % state % tracers % array(:,k,iCell) + &
- block % intermediate_step(Q2) % tracers % array(:,k,iCell) / 2.0 &
- ) / block % intermediate_step(TEMP) % h % array(k,iCell)
- end do
- end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
- end if
- call compute_solve_diagnostics(block % intermediate_step(TEMP), block % mesh)
- call compute_tend(block % intermediate_step(Q3), block % intermediate_step(TEMP), block % mesh)
- call compute_scalar_tend(block % intermediate_step(Q3), block % intermediate_step(TEMP), block % mesh)
- block % intermediate_step(Q3) % u % array(:,:) = dt * block % intermediate_step(Q3) % u % array(:,:)
- block % intermediate_step(Q3) % h % array(:,:) = dt * block % intermediate_step(Q3) % h % array(:,:)
- block % intermediate_step(Q3) % tracers % array(:,:,:) = dt * block % intermediate_step(Q3) % tracers % array(:,:,:)
+ rk_substep_weights(1) = dt/2.
+ rk_substep_weights(2) = dt/2.
+ rk_substep_weights(3) = dt
+ rk_substep_weights(4) = 0.
- block => block % next
- end do
- block => domain % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q3) % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q3) % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(Q3) % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! BEGIN RK loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do rk_step = 1, 4
+! --- compute tendencies
- ! Take fourth RK step
- block => domain % blocklist
- do while (associated(block))
- block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:) + &
- block % intermediate_step(Q3) % u % array(:,:)
- block % intermediate_step(TEMP) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) + &
- block % intermediate_step(Q3) % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % intermediate_step(TEMP) % tracers % array(:,k,iCell) = ( &
- block % time_levs(1) % state % h % array(k,iCell) * &
- block % time_levs(1) % state % tracers % array(:,k,iCell) + &
- block % intermediate_step(Q3) % tracers % array(:,k,iCell) &
- ) / block % intermediate_step(TEMP) % h % array(k,iCell)
- end do
- end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
- end if
- call compute_solve_diagnostics(block % intermediate_step(TEMP), block % mesh)
- call compute_tend(block % intermediate_step(Q4), block % intermediate_step(TEMP), block % mesh)
- call compute_scalar_tend(block % intermediate_step(Q4), block % intermediate_step(TEMP), block % mesh)
- block % intermediate_step(Q4) % u % array(:,:) = dt * block % intermediate_step(Q4) % u % array(:,:)
- block % intermediate_step(Q4) % h % array(:,:) = dt * block % intermediate_step(Q4) % h % array(:,:)
- block % intermediate_step(Q4) % tracers % array(:,:,:) = dt * block % intermediate_step(Q4) % tracers % array(:,:,:)
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
+ call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
+ block => block % next
+ end do
- block => block % next
- end do
+! --- update halos for prognostic variables
- block => domain % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q4) % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q4) % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(Q4) % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
+ 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) % h % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(TEND) % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
+! --- compute next substep state
+
+ if (rk_step < 4) then
+ block => domain % blocklist
+ do while (associated(block))
+ block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:) &
+ + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
+ block % intermediate_step(PROVIS) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) &
+ + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % intermediate_step(PROVIS) % tracers % array(:,k,iCell) = ( &
+ block % time_levs(1) % state % h % array(k,iCell) * &
+ block % time_levs(1) % state % tracers % array(:,k,iCell) &
+ + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell) &
+ ) / block % intermediate_step(PROVIS) % h % array(k,iCell)
+ end do
+ end do
+ if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
+ end if
+ call compute_solve_diagnostics(block % intermediate_step(PROVIS), block % mesh)
+ block => block % next
+ end do
+ end if
+
+!--- accumulate update (for RK4)
+
+ block => domain % blocklist
+ do while (associated(block))
+ block % time_levs(2) % state % u % array(:,:) = block % time_levs(2) % state % u % array(:,:) &
+ + rk_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
+ block % time_levs(2) % state % h % array(:,:) = block % time_levs(2) % state % h % array(:,:) &
+ + rk_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ block % time_levs(2) % state % tracers % array(:,k,iCell) &
+ + rk_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell)
+ end do
+ end do
+ block => block % next
+ end do
+
end do
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! END RK loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ! Compute final step and add to starting model state
+
+ !
+ ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
+ !
block => domain % blocklist
do while (associated(block))
- block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:) + &
- ( block % intermediate_step(Q1) % u % array(:,:) + &
- 2.0 * block % intermediate_step(Q2) % u % array(:,:) + &
- 2.0 * block % intermediate_step(Q3) % u % array(:,:) + &
- block % intermediate_step(Q4) % u % array(:,:) &
- ) / 6.0
- block % time_levs(2) % state % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) + &
- ( block % intermediate_step(Q1) % h % array(:,:) + &
- 2.0 * block % intermediate_step(Q2) % h % array(:,:) + &
- 2.0 * block % intermediate_step(Q3) % h % array(:,:) + &
- block % intermediate_step(Q4) % h % array(:,:) &
- ) / 6.0
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
- block % time_levs(2) % state % tracers % array(:,k,iCell) = ( &
- block % time_levs(1) % state % h % array(k,iCell) * &
- block % time_levs(1) % state % tracers % array(:,k,iCell) + &
- ( block % intermediate_step(Q1) % tracers % array(:,k,iCell) + &
- 2.0 * block % intermediate_step(Q2) % tracers % array(:,k,iCell) + &
- 2.0 * block % intermediate_step(Q3) % tracers % array(:,k,iCell) + &
- block % intermediate_step(Q4) % tracers % array(:,k,iCell) &
- ) / 6.0 &
- ) / block % time_levs(2) % state % h % array(k,iCell)
+ block % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ block % time_levs(2) % state % tracers % array(:,k,iCell) &
+ / block % time_levs(2) % state % h % array(k,iCell)
end do
end do
</font>
</pre>