<p><b>mpetersen@lanl.gov</b> 2011-05-12 10:19:04 -0600 (Thu, 12 May 2011)</p><p>Split compute_tend into compute_tend_u and compute_tend_h. RK45 runs checked, are identical.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-11 19:25:16 UTC (rev 828)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-12 16:19:04 UTC (rev 829)
@@ -16,6 +16,12 @@
namelist real restart config_restart_time 172800.0
namelist character grid config_vert_grid_type isopycnal
namelist real grid config_rho0 1028
+namelist integer timestep_higdon config_n_ts_iter 2
+namelist integer timestep_higdon config_n_bcl_iter_beg 4
+namelist integer timestep_higdon config_n_bcl_iter_mid 4
+namelist integer timestep_higdon config_n_bcl_iter_end 4
+namelist integer timestep_higdon config_n_btr_subcycles 10
+namelist logical timestep_higdon config_compute_tr_midstage true
namelist real hmix config_h_mom_eddy_visc2 0.0
namelist real hmix config_h_mom_eddy_visc4 0.0
namelist real hmix config_h_tracer_eddy_diff2 0.0
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-11 19:25:16 UTC (rev 828)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-12 16:19:04 UTC (rev 829)
@@ -29,10 +29,14 @@
if (trim(config_time_integration) == 'RK4') then
call rk4(domain, dt)
+ elseif (trim(config_time_integration) == 'higdon_split' &
+ .or.trim(config_time_integration) == 'higdon_unsplit') then
+ call higdon_timestep(domain, dt)
else
write(0,*) 'Abort: Unknown time integration option '&
//trim(config_time_integration)
- write(0,*) 'Currently, only ''RK4'' is supported.'
+ write(0,*) 'Currently, only RK4, hidgdon_split and '// &
+ 'higdon_unsplit are supported.'
call dmpar_abort(dminfo)
end if
@@ -156,7 +160,8 @@
if (.not.config_implicit_vertical_mix) then
call compute_vertical_mix_coefficients(provis, block % diagnostics, block % mesh)
end if
- call compute_tend(block % tend, provis, block % diagnostics, block % mesh)
+ call compute_tend_h(block % tend, provis, block % diagnostics, block % mesh)
+ call compute_tend_u(block % tend, provis, block % diagnostics, block % mesh)
call compute_scalar_tend(block % tend, provis, block % diagnostics, block % mesh)
call enforce_boundaryEdge(block % tend, block % mesh)
block => block % next
@@ -348,8 +353,273 @@
end subroutine rk4
+ subroutine higdon_timestep(domain, dt)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Advance model state forward in time by the specified time step using
+ ! 4th order Runge-Kutta
+ !
+ ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
+ ! plus grid meta-data
+ ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
+ ! model state advanced forward in time by dt seconds
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine compute_tend(tend, s, d, grid)
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ real (kind=RKIND), intent(in) :: dt
+
+ integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step
+ type (block_type), pointer :: block
+
+ integer :: nCells, nEdges, nVertLevels, num_tracers
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer, dimension(:), pointer :: &
+ maxLevelCell, maxLevelEdgeTop
+ real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp
+ real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
+
+
+ !
+ ! 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 % 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
+
+!mrp 110512 need to add w, p, and eta here?
+
+ block => block % next
+ end do
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! BEGIN large iteration (predictor/corrector when config_n_ts_iter=2)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do higdon_step = 1, config_n_ts_iter
+! --- update halos for diagnostic variables
+
+ block => domain % blocklist
+ do while (associated(block))
+! mrp 110512 not sure if I need the following three. Leave be, assume I need it.
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % pv_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ end if
+
+ block => block % next
+ end do
+
+ !
+ ! Stage 1: Baroclinic velocity (3D) prediction, explicit with long timestep
+ !
+
+! --- compute tendencies
+
+ block => domain % blocklist
+ do while (associated(block))
+ if (.not.config_implicit_vertical_mix) then
+ call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+ end if
+ call compute_tend_u(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call enforce_boundaryEdge(block % tend, block % mesh)
+ block => block % next
+ end do
+
+ !
+ ! Stage 2: Barotropic velocity (2D) prediction, explicitly subcycled
+ !
+
+
+ !
+ ! Stage 3: Tracer, density, pressure, vertical velocity prediction
+ !
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call compute_tend_h(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call enforce_boundaryEdge(block % tend, block % mesh)
+ 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 % tend % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % h % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % tend % tracers % array(:,:,:), &
+ block % tend % num_tracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
+!--- accumulate update (for HIGDON_TIMESTEP)
+
+ 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 % maxLevelCell % array(iCell)
+ ! 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
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! END large iteration loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ !
+ ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
+ !
+ block => domain % blocklist
+ do while (associated(block))
+
+ u => block % state % time_levs(2) % state % u % array
+ tracers => block % state % time_levs(2) % state % tracers % array
+ h => block % state % time_levs(2) % state % h % array
+ h_edge => block % state % time_levs(2) % state % h_edge % array
+ num_tracers = block % state % time_levs(2) % state % num_tracers
+ vertViscTopOfEdge => block % diagnostics % vertViscTopOfEdge % array
+ vertDiffTopOfCell => block % diagnostics % vertDiffTopOfCell % array
+ maxLevelCell => block % mesh % maxLevelCell % array
+ maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertLevels = block % mesh % nVertLevels
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ tracers(:,k,iCell) = tracers(:,k,iCell) / h(k,iCell)
+ end do
+ end do
+
+ if (config_implicit_vertical_mix) then
+ allocate(A(nVertLevels),C(nVertLevels),uTemp(nVertLevels), &
+ tracersTemp(num_tracers,nVertLevels))
+
+ call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+
+ !
+ ! Implicit vertical solve for momentum
+ !
+ do iEdge=1,nEdges
+ if (maxLevelEdgeTop(iEdge).gt.0) then
+
+ ! Compute A(k), C(k) for momentum
+ ! mrp 110315 efficiency note: for z-level, could precompute
+ ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
+ ! h_edge is computed in compute_solve_diag, and is not available yet.
+ ! This could be removed if hZLevel used instead.
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
+
+ do k=1,maxLevelEdgeTop(iEdge)-1
+ A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &
+ / (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &
+ / h_edge(k,iEdge)
+ enddo
+ A(maxLevelEdgeTop(iEdge)) = 0.0
+
+ C(1) = 1 - A(1)
+ do k=2,maxLevelEdgeTop(iEdge)
+ C(k) = 1 - A(k) - A(k-1)
+ enddo
+
+ call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
+
+ u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
+ u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
+
+ end if
+ end do
+
+ !
+ ! Implicit vertical solve for tracers
+ !
+ do iCell=1,nCells
+ ! Compute A(k), C(k) for tracers
+ ! mrp 110315 efficiency note: for z-level, could precompute
+ ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
+ do k=1,maxLevelCell(iCell)-1
+ A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &
+ / (h(k,iCell) + h(k+1,iCell)) &
+ / h(k,iCell)
+ enddo
+
+ A(maxLevelCell(iCell)) = 0.0
+
+ C(1) = 1 - A(1)
+ do k=2,maxLevelCell(iCell)
+ C(k) = 1 - A(k) - A(k-1)
+ enddo
+
+ call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &
+ tracersTemp, maxLevelCell(iCell), &
+ nVertLevels,num_tracers)
+
+ tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
+ tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
+ end do
+ deallocate(A,C,uTemp,tracersTemp)
+ end if
+
+ if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ end if
+
+ call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+
+ call reconstruct(block % state % time_levs(2) % state, block % mesh)
+
+ block => block % next
+ end do
+
+ end subroutine higdon_timestep
+
+
+ subroutine compute_tend_h(tend, s, d, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute height and normal wind tendencies, as well as diagnostic variables
!
@@ -369,6 +639,8 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
vertex1, vertex2, eoe, i, j
+! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
+! Most of these variables can be removed, but at a later time.
integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
upstream_bias, wTopEdge, rho0Inv, r
@@ -377,7 +649,7 @@
zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
- tend_h, tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
+ tend_h, circulation, vorticity, ke, ke_edge, pv_edge, &
MontPot, wTop, divergence, vertViscTopOfEdge
type (dm_info) :: dminfo
@@ -393,10 +665,6 @@
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
- real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
-
h => s % h % array
u => s % u % array
v => s % v % array
@@ -436,7 +704,6 @@
maxLevelVertexBot => grid % maxLevelVertexBot % array
tend_h => tend % h % array
- tend_u => tend % u % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -444,8 +711,6 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- u_src => grid % u_src % array
-
!
! height tendency: start accumulating tendency terms
!
@@ -503,6 +768,106 @@
end do
endif ! coordinate type
+ end subroutine compute_tend_h
+
+ subroutine compute_tend_u(tend, s, d, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute height and normal wind tendencies, as well as diagnostic variables
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (diagnostics_type), intent(in) :: d
+ type (mesh_type), intent(in) :: grid
+
+! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
+! Some of these variables can be removed, but at a later time.
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
+ upstream_bias, wTopEdge, rho0Inv, r
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ zMidZLevel, zTopZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
+ tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
+ MontPot, wTop, divergence, vertViscTopOfEdge
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
+ edgesOnEdge, edgesOnVertex
+ real (kind=RKIND) :: u_diffusion
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+ vertViscTopOfEdge => d % vertViscTopOfEdge % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+ tend_u => tend % u % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ u_src => grid % u_src % array
+
!
! velocity tendency: start accumulating tendency terms
!
@@ -764,7 +1129,7 @@
deallocate(fluxVertTop)
endif
- end subroutine compute_tend
+ end subroutine compute_tend_u
subroutine compute_scalar_tend(tend, s, d, grid)
</font>
</pre>