<p><b>dwj07@fsu.edu</b> 2012-06-05 14:23:24 -0600 (Tue, 05 Jun 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Updating ocean (RK4) and sw (RK4) cores to work with multiple blocks.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_mpas_core.F 2012-06-05 19:31:39 UTC (rev 1962)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_mpas_core.F 2012-06-05 20:23:24 UTC (rev 1963)
@@ -360,7 +360,11 @@
call mpas_timer_stop("time integration", timeIntTimer)
! Move time level 2 fields back into time level 1 for next time step
- call mpas_shift_time_levels_state(domain % blocklist % state)
+ block_ptr => domain % blocklist
+ do while(associated(block_ptr))
+ call mpas_shift_time_levels_state(block_ptr % state)
+ block_ptr => block_ptr % next
+ end do
if (mpas_is_alarm_ringing(clock, outputAlarmID, ierr=ierr)) then
call mpas_reset_clock_alarm(clock, outputAlarmID, ierr=ierr)
Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-06-05 19:31:39 UTC (rev 1962)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-06-05 20:23:24 UTC (rev 1963)
@@ -78,8 +78,6 @@
integer :: iCell, k, i, err
type (block_type), pointer :: block
- type (state_type), target :: provis
- type (state_type), pointer :: provis_ptr
integer :: rk_step, iEdge, cell1, cell2
@@ -97,13 +95,29 @@
block => domain % blocklist
- call mpas_allocate_state(block, provis, &
- block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
- block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels )
+ do while(associated(block))
+ allocate(block % provis)
+ call mpas_allocate_state(block, block % provis, &
+ block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
+ block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels )
- provis_ptr => provis
- call mpas_create_state_links(provis_ptr)
+ block => block % next
+ end do
+ block => domain % blocklist
+ do while(associated(block))
+ if(associated(block % prev) .and. associated(block % next)) then
+ call mpas_create_state_links(block % provis, prev = block % prev % provis, next = block % next % provis)
+ else if(associated(block % prev)) then
+ call mpas_create_state_links(block % provis, prev = block % prev % provis)
+ else if(associated(block % next)) then
+ call mpas_create_state_links(block % provis, next = block % next % provis)
+ else
+ call mpas_create_state_links(block % provis)
+ end if
+ block => block % next
+ end do
+
!
! Initialize time_levs(2) with state at current time
! Initialize first RK state
@@ -112,19 +126,18 @@
!
block => domain % blocklist
do while (associated(block))
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ do iCell=1,block % mesh % nCells ! couple tracers to h
+ do k=1,block % mesh % maxLevelCell % array(iCell)
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ * block % state % time_levs(1) % state % h % array(k,iCell)
+ end do
+ end do
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
- do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % maxLevelCell % array(iCell)
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % h % array(k,iCell)
- end do
- end do
+ call mpas_copy_state(block % provis, block % state % time_levs(1) % state)
- call mpas_copy_state(provis, block % state % time_levs(1) % state)
-
- block => block % next
+ block => block % next
end do
rk_weights(1) = dt/6.
@@ -146,10 +159,10 @@
! --- update halos for diagnostic variables
call mpas_timer_start("RK4-diagnostic halo update")
- call mpas_dmpar_exch_halo_field(provis % Vor_edge)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % provis % Vor_edge)
if (config_h_mom_eddy_visc4 > 0.0) then
- call mpas_dmpar_exch_halo_field(provis % divergence)
- call mpas_dmpar_exch_halo_field(provis % vorticity)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % provis % divergence)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % provis % vorticity)
end if
call mpas_timer_stop("RK4-diagnostic halo update")
@@ -160,21 +173,21 @@
do while (associated(block))
! mrp 111206 put ocn_wtop call at top for ALE
- call ocn_wtop(provis, provis, block % mesh)
+ call ocn_wtop(block % provis, block % provis, block % mesh)
if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, provis, block % diagnostics, err)
+ call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, err)
end if
- call ocn_tend_h(block % tend, provis, block % mesh)
- call ocn_tend_u(block % tend, provis, block % diagnostics, block % mesh)
+ call ocn_tend_h(block % tend, block % provis, block % mesh)
+ call ocn_tend_u(block % tend, block % provis, block % diagnostics, block % mesh)
! mrp 110718 filter btr mode out of u_tend
! still got h perturbations with just this alone. Try to set uBtr=0 after full u computation
if (config_rk_filter_btr_mode) then
- call filter_btr_mode_tend_u(block % tend, provis, block % diagnostics, block % mesh)
+ call filter_btr_mode_tend_u(block % tend, block % provis, block % diagnostics, block % mesh)
endif
- call ocn_tend_scalar(block % tend, provis, block % diagnostics, block % mesh, dt)
+ call ocn_tend_scalar(block % tend, block % provis, block % diagnostics, block % mesh, dt)
block => block % next
end do
call mpas_timer_stop("RK4-tendency computations")
@@ -194,39 +207,38 @@
block => domain % blocklist
do while (associated(block))
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
+ block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
- provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % h % array(:,:)
+ block % provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % h % array(:,:)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
- provis % tracers % array(:,k,iCell) = ( &
- block % state % time_levs(1) % state % h % array(k,iCell) * &
- block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
- ) / provis % h % array(k,iCell)
+ block % provis % tracers % array(:,k,iCell) = ( block % state % time_levs(1) % state % h % array(k,iCell) * &
+ block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
+ ) / block % provis % h % array(k,iCell)
end do
end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
if (config_prescribe_velocity) then
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
if (config_prescribe_thickness) then
- provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ block % provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
end if
- call ocn_diagnostic_solve(dt, provis, block % mesh)
+ call ocn_diagnostic_solve(dt, block % provis, block % mesh)
! Compute velocity transport, used in advection terms of h and tracer tendancy
- provis % uTransport % array(:,:) &
- = provis % u % array(:,:) &
- + provis % uBolusGM % array(:,:)
+ block % provis % uTransport % array(:,:) &
+ = block % provis % u % array(:,:) &
+ + block % provis % uBolusGM % array(:,:)
block => block % next
end do
@@ -349,7 +361,12 @@
end do
call mpas_timer_stop("RK4-cleaup phase")
- call mpas_deallocate_state(provis)
+ block => domain % blocklist
+ do while(associated(block))
+ call mpas_deallocate_state(block % provis)
+ deallocate(block % provis)
+ block => block % next
+ end do
end subroutine ocn_time_integrator_rk4!}}}
Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tracer_advection.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tracer_advection.F 2012-06-05 19:31:39 UTC (rev 1962)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tracer_advection.F 2012-06-05 20:23:24 UTC (rev 1963)
@@ -150,6 +150,7 @@
end do
nAdvCellsForEdge(iEdge) = n
+
do iCell = 1, nAdvCellsForEdge(iEdge)
advCellsForEdge(iCell,iEdge) = ordered_cell_list(iCell)
end do
Modified: branches/omp_blocks/multiple_blocks/src/core_sw/mpas_sw_mpas_core.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_sw/mpas_sw_mpas_core.F 2012-06-05 19:31:39 UTC (rev 1962)
+++ branches/omp_blocks/multiple_blocks/src/core_sw/mpas_sw_mpas_core.F 2012-06-05 20:23:24 UTC (rev 1963)
@@ -186,7 +186,11 @@
call mpas_timer_stop("time integration")
! Move time level 2 fields back into time level 1 for next time step
- call mpas_shift_time_levels_state(domain % blocklist % state)
+ block_ptr => domain % blocklist
+ do while(associated(block_ptr))
+ call mpas_shift_time_levels_state(block_ptr % state)
+ block_ptr => block_ptr % next
+ end do
!TODO: mpas_get_clock_ringing_alarms is probably faster than multiple mpas_is_alarm_ringing...
Modified: branches/omp_blocks/multiple_blocks/src/core_sw/mpas_sw_time_integration.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_sw/mpas_sw_time_integration.F 2012-06-05 19:31:39 UTC (rev 1962)
+++ branches/omp_blocks/multiple_blocks/src/core_sw/mpas_sw_time_integration.F 2012-06-05 20:23:24 UTC (rev 1963)
@@ -71,121 +71,136 @@
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
block => domain % blocklist
- call mpas_allocate_state(block, provis, &
- block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
- block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels, &
- block % mesh % nTracers)
+ do while(associated(block))
+ allocate(block % provis)
+ call mpas_allocate_state(block, block % provis, &
+ block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
+ block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels, &
+ block % mesh % nTracers)
+
+ block => block % next
+ end do
- provis_ptr => provis
- call mpas_create_state_links(provis_ptr)
+ block => domain % blocklist
+ do while(associated(block))
+ if(associated(block % prev) .and. associated(block % next)) then
+ call mpas_create_state_links(block % provis, prev = block % prev % provis, next = block % next % provis)
+ else if(associated(block % prev)) then
+ call mpas_create_state_links(block % provis, prev = block % prev % provis)
+ else if(associated(block % next)) then
+ call mpas_create_state_links(block % provis, next = block % next % provis)
+ else
+ call mpas_create_state_links(block % provis)
+ end if
+ 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))
+ !
+ ! 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))
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
- do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % h % array(k,iCell)
- end do
- end do
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ do iCell=1,block % mesh % nCells ! couple tracers to h
+ do k=1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ * block % state % time_levs(1) % state % h % array(k,iCell)
+ end do
+ end do
- call mpas_copy_state(provis, block % state % time_levs(1) % state)
+ call mpas_copy_state(block % provis, block % state % time_levs(1) % state)
- block => block % next
- end do
+ 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.
+ rk_weights(1) = dt/6.
+ rk_weights(2) = dt/3.
+ rk_weights(3) = dt/3.
+ rk_weights(4) = dt/6.
- rk_substep_weights(1) = dt/2.
- rk_substep_weights(2) = dt/2.
- rk_substep_weights(3) = dt
- rk_substep_weights(4) = 0.
+ rk_substep_weights(1) = dt/2.
+ rk_substep_weights(2) = dt/2.
+ rk_substep_weights(3) = dt
+ rk_substep_weights(4) = 0.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! BEGIN RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do rk_step = 1, 4
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! BEGIN RK loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do rk_step = 1, 4
-! --- update halos for diagnostic variables
+! --- update halos for diagnostic variables
- call mpas_dmpar_exch_halo_field(provis % pv_edge)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % provis % pv_edge)
if (config_h_mom_eddy_visc4 > 0.0) then
- call mpas_dmpar_exch_halo_field(provis % divergence)
- call mpas_dmpar_exch_halo_field(provis % vorticity)
- end if
+ call mpas_dmpar_exch_halo_field(domain % blocklist % provis % divergence)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % provis % vorticity)
+ end if
-! --- compute tendencies
+! --- compute tendencies
- block => domain % blocklist
- do while (associated(block))
- call sw_compute_tend(block % tend, provis, block % mesh)
- call sw_compute_scalar_tend(block % tend, provis, block % mesh)
- call sw_enforce_boundary_edge(block % tend, block % mesh)
- block => block % next
- end do
+ block => domain % blocklist
+ do while (associated(block))
+ call sw_compute_tend(block % tend, block % provis, block % mesh)
+ call sw_compute_scalar_tend(block % tend, block % provis, block % mesh)
+ call sw_enforce_boundary_edge(block % tend, block % mesh)
+ block => block % next
+ end do
-! --- update halos for prognostic variables
+! --- update halos for prognostic variables
- call mpas_dmpar_exch_halo_field(domain % blocklist % tend % u)
- call mpas_dmpar_exch_halo_field(domain % blocklist % tend % h)
- call mpas_dmpar_exch_halo_field(domain % blocklist % tend % tracers)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % tend % u)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % tend % h)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % tend % tracers)
-! --- compute next substep state
+! --- compute next substep state
- if (rk_step < 4) then
- block => domain % blocklist
- do while (associated(block))
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
- provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- provis % tracers % array(:,k,iCell) = ( &
- block % state % time_levs(1) % state % h % array(k,iCell) * &
- block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
- ) / provis % h % array(k,iCell)
- end do
- end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- end if
- call sw_compute_solve_diagnostics(dt, provis, block % mesh)
- block => block % next
- end do
- end if
+ if (rk_step < 4) then
+ block => domain % blocklist
+ do while (associated(block))
+ block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
+ block % provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % provis % tracers % array(:,k,iCell) = ( block % state % time_levs(1) % state % h % array(k,iCell) * &
+ block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
+ ) / block % provis % h % array(k,iCell)
+ end do
+ end do
+ if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ end if
+ call sw_compute_solve_diagnostics(dt, block % provis, block % mesh)
+ block => block % next
+ end do
+ end if
!--- accumulate update (for RK4)
- block => domain % blocklist
- do while (associated(block))
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
- + rk_weights(rk_step) * block % tend % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
- + rk_weights(rk_step) * block % tend % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
- end do
- end do
- block => block % next
- end do
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
+ + rk_weights(rk_step) * block % tend % u % array(:,:)
+ block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
+ + rk_weights(rk_step) * block % tend % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
+ end do
+ end do
+ block => block % next
+ end do
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -223,7 +238,12 @@
block => block % next
end do
- call mpas_deallocate_state(provis)
+ block => domain % blocklist
+ do while(associated(block))
+ call mpas_deallocate_state(block % provis)
+ deallocate(block % provis)
+ block => block % next
+ end do
end subroutine sw_rk4
</font>
</pre>