<p><b>mpetersen@lanl.gov</b> 2011-06-02 07:06:27 -0600 (Thu, 02 Jun 2011)</p><p>All code in place for Higdon Split. It does not work yet, I am still debugging.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/Makefile
===================================================================
--- branches/ocean_projects/timesplitting_mrp/Makefile 2011-05-31 22:39:08 UTC (rev 865)
+++ branches/ocean_projects/timesplitting_mrp/Makefile 2011-06-02 13:06:27 UTC (rev 866)
@@ -86,7 +86,7 @@
"CC = gcc" \
"SFC = ifort" \
"SCC = gcc" \
- "FFLAGS = -real-size 64 -O3 -convert big_endian" \
+ "FFLAGS = -real-size 64 -convert big_endian -check all" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-31 22:39:08 UTC (rev 865)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-06-02 13:06:27 UTC (rev 866)
@@ -188,7 +188,7 @@
var persistent real sshSubcycle ( nCells Time ) 2 o sshSubcycle state - -
var persistent real FBtr ( nEdges Time ) 1 o FBtr state - -
var persistent real GBtrForcing ( nEdges Time ) 1 o GBtrForcing state - -
-var persistent real uBcl ( nVertLevels nEdges Time ) 3 o uBcl state - -
+var persistent real uBcl ( nVertLevels nEdges Time ) 2 o uBcl state - -
# Diagnostic fields: only written to output
var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-05-31 22:39:08 UTC (rev 865)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-06-02 13:06:27 UTC (rev 866)
@@ -282,7 +282,8 @@
integer :: i, iCell, iEdge, iVertex, k
type (block_type), pointer :: block
- integer :: iTracer, cell
+ integer :: iTracer, cell, cell1, cell2
+ real (kind=RKIND) :: uhSum, hSum, sshEdge
real (kind=RKIND), dimension(:), pointer :: &
hZLevel, zMidZLevel, zTopZLevel, &
hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
@@ -325,6 +326,53 @@
hRatioZLevelKm1(k) = 0.5*hZLevel(k-1)/hMeanTopZLevel(k)
end do
+! mrp delete later - probably dont need HTotZLevelEdge variable
+! do iEdge=1,block % mesh % nEdges
+! block % mesh % HTotZLevelEdge % array(iEdge) = 0.0
+! do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+! block % mesh % HTotZLevelEdge % array(iEdge) &
+! = block % mesh % HTotZLevelEdge % array(iEdge) &
+! + block % mesh % hZLevel % array(k)
+! enddo
+! enddo
+
+ ! mrp 110601 For now, h is the variable saved in the restart file
+ ! I am computing SSH here. In the future, could make smaller
+ ! restart files for z-Level runs by saving SSH only.
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(1) % state % ssh % array(iCell) &
+ = block % state % time_levs(1) % state % h % array(1,iCell) &
+ - block % mesh % hZLevel % array(1)
+ enddo
+
+ ! Compute barotropic velocity at first timestep
+ ! This is only done upon start-up.
+ if (trim(config_time_integration) == 'higdon_unsplit') then
+ block % state % time_levs(1) % state % uBtr % array(:) = 0.0
+ elseif (trim(config_time_integration) == 'higdon_split') then
+ do iEdge=1,block % mesh % nEdges
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ sshEdge = 0.5*( &
+ block % state % time_levs(1) % state % ssh % array(cell1) &
+ + block % state % time_levs(1) % state % ssh % array(cell2) )
+
+ uhSum = (sshEdge + block % mesh % hZLevel % array(1)) &
+ * block % state % time_levs(1) % state % u % array(1,iEdge)
+ hSum = sshEdge + block % mesh % hZLevel % array(1)
+
+ do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
+ uhSum = uhSum &
+ + block % mesh % hZLevel % array(k) &
+ *block % state % time_levs(1) % state % u % array(k,iEdge)
+ hSum = hSum &
+ + block % mesh % hZLevel % array(k)
+ enddo
+ block % state % time_levs(1) % state % uBtr % array(iEdge) = uhSum/hsum
+ enddo
+ endif
+
block => block % next
end do
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-31 22:39:08 UTC (rev 865)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-06-02 13:06:27 UTC (rev 866)
@@ -372,10 +372,12 @@
real (kind=RKIND), intent(in) :: dt
integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step, split, &
- eoe, oldBtrSubcycleTime, newBtrSubcycleTime
+ eoe, oldBtrSubcycleTime, newBtrSubcycleTime, uPerpTime, BtrCorIter, &
+ n_bcl_iter(config_n_ts_iter)
type (block_type), pointer :: block
- real (kind=RKIND) :: uSum, hSum, GSum, grav_rho0Inv, sshEdge, flux, uPerp
+ real (kind=RKIND) :: uhSum, hSum, grav_rho0Inv, sshEdge, flux, &
+ uPerp, uCorr, tracerTemp
integer :: nCells, nEdges, nVertLevels, num_tracers
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -396,24 +398,6 @@
block => domain % blocklist
do while (associated(block))
- ! Compute barotropic velocity at the old timestep
- if (trim(config_time_integration) == 'higdon_unsplit') then
- block % state % time_levs(1) % state % uBtr % array(:) = 0.0
- elseif (trim(config_time_integration) == 'higdon_split') then
- do iEdge=1,block % mesh % nEdges
- uSum = 0.0
- hSum = 0.0
- ! 110518 mrp efficiency note:
- ! in a z-level model, could change hSum to use sum of hZLevel, or have another
- ! 1-D z array which is total thickness between the surface and that level.
- do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
- uSum = uSum + block % state % time_levs(1) % state % u % array(k,iEdge)
- hSum = hSum + block % state % time_levs(1) % state % h_edge % array(k,iEdge)
- enddo
- block % state % time_levs(1) % state % uBtr % array(iEdge) = uSum/hSum
- enddo
- endif
-
!mrp 110512 need to add w, p, and ssh here?
do iEdge=1,block % mesh % nEdges
@@ -474,6 +458,10 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ n_bcl_iter = config_n_bcl_iter_mid
+ n_bcl_iter(1) = config_n_bcl_iter_beg
+ n_bcl_iter(config_n_ts_iter) = config_n_bcl_iter_end
+
do higdon_step = 1, config_n_ts_iter
! --- update halos for diagnostic variables
@@ -514,9 +502,10 @@
block => block % next
end do
- ! baroclinic iterations on linear Coriolis term
-! change this from beg later
- do j=1,config_n_bcl_iter_beg
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! BEGIN baroclinic iterations on linear Coriolis term
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do j=1,n_bcl_iter(higdon_step)
! mrp soon: add linear Coriolis to tend_u maybe at end of this iteration?
!compute {\bf u}_{k,n+1}^{'\;\perp} from {\bf u}_{k,n+1}'
@@ -562,17 +551,21 @@
enddo
- GSum = 0.0
- hSum = 0.0
- ! 110518 mrp efficiency note:
- ! in a z-level model, could change hSum to use sum of hZLevel, or have another
- ! 1-D z array which is total thickness between the surface and that level.
- do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
- GSum = GSum + block % state % time_levs(2) % state % h_edge % array(k,iEdge)*uTemp(k)
- hSum = hSum + block % state % time_levs(2) % state % h_edge % array(k,iEdge)
+ ! Compute GBtrForcing, the vertically averaged forcing
+ sshEdge = 0.5*( &
+ block % state % time_levs(1) % state % ssh % array(cell1) &
+ + block % state % time_levs(1) % state % ssh % array(cell2) )
+
+ uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * uTemp(1)
+ hSum = sshEdge + block % mesh % hZLevel % array(1)
+
+ do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
+ uhSum = uhSum + block % mesh % hZLevel % array(k) * uTemp(k)
+ hSum = hSum + block % mesh % hZLevel % array(k)
enddo
- block % state % time_levs(1) % state % GBtrForcing % array(iEdge) = split*GSum/hSum/dt
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge) = split*uhSum/hSum/dt
+
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
! These two steps are together here:
!{\bf u}'_{k,n+1} = {\bf u}'_{k,n} - \Delta t {\overline {\bf G}}
@@ -608,7 +601,10 @@
end do
enddo ! do j=1,config_n_bcl_iter
-!print *, '5'
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! END baroclinic iterations on linear Coriolis term
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+print *, '5'
@@ -625,9 +621,12 @@
do while (associated(block))
block % state % time_levs(2) % state % uBtr % array(:) = 0.0
- ! mrp 110531 change following lines once I have pointers set up for subcycling
+
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:) = 0.0
+ block % state % time_levs(2) % state % u % array(:,:) &
+ = block % state % time_levs(2) % state % uBcl % array(:,:)
+
block => block % next
end do ! block
@@ -646,6 +645,7 @@
block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(1) % state % ssh % array(iCell)
enddo
+print *, '6'
do iEdge=1,block % mesh % nEdges
@@ -663,9 +663,11 @@
block => block % next
end do ! block
+print *, '7'
-
- ! Barotropic subcycle loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! BEGIN Barotropic subcycle loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do j=1,config_n_btr_subcycles
! Compute thickness flux and new SSH
@@ -674,6 +676,7 @@
allocate(tend_ssh(block % mesh % nCells))
+print *, '7.1'
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -693,8 +696,7 @@
! \left(\zeta_{n+j/J}^{edge} + \sum_{k=1}^{N^{edge}} \Delta z_k \right)
flux = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
* block % mesh % dvEdge % array(iEdge) &
- * h_edge(k,iEdge) &
- * (hSum + sshEdge)
+ * (sshEdge + hSum)
tend_ssh(cell1) = tend_ssh(cell1) - flux
tend_ssh(cell2) = tend_ssh(cell2) + flux
@@ -702,6 +704,7 @@
= block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ flux
end do
+print *, '7.2'
do iCell=1,block % mesh % nCells
@@ -720,24 +723,28 @@
end do
+print *, '7.3'
deallocate(tend_ssh)
+print *, '7.4'
block => block % next
end do ! block
+print *, '8'
-
-! |boundary update on \zeta_{n+(j+1)/J}
+! boundary update on \zeta_{n+(j+1)/J}
block => domain % blocklist
do while (associated(block))
- call dmpar_exch_halo_field1dReal(domain % dminfo, &
- block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
- block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
+! call dmpar_exch_halo_field1dReal(domain % dminfo, &
+! block % state % time_levs(1) % state % sshSubcycle % array(:), &
+! block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
block => block % next
end do ! block
+print *, '9'
-
! compute new barotropic velocity
! do iEdge=1,nEdges
@@ -752,6 +759,10 @@
do while (associated(block))
+ ! Iterate on u two times.
+ do BtrCorIter=1,2
+ if (BtrCorIter==1) uPerpTime = oldBtrSubcycleTime
+ if (BtrCorIter==2) uPerpTime = newBtrSubcycleTime
do iEdge=1,block % mesh % nEdges
@@ -763,7 +774,7 @@
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
eoe = block % mesh % edgesOnEdge % array(i,iEdge)
uPerp = uPerp + block % mesh % weightsOnEdge % array(i,iEdge) &
- * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(eoe) &
+ * block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
* block % mesh % fEdge % array(eoe)
end do
@@ -790,11 +801,13 @@
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
end do
+ enddo !do BtrCorIter=1,2
+
block => block % next
end do ! block
+print *, '10'
-
-! |boundary update on \ubtr_{n+(j+1)/J}
+! boundary update on \ubtr_{n+(j+1)/J}
block => domain % blocklist
do while (associated(block))
@@ -809,27 +822,18 @@
! advance time pointers
oldBtrSubcycleTime = mod(oldBtrSubcycleTime,2)+1
newBtrSubcycleTime = mod(newBtrSubcycleTime,2)+1
+print *, '11 oldBtrSubcycleTime, newBtrSubcycleTime',oldBtrSubcycleTime, newBtrSubcycleTime
- end do ! barotropic subcycling on j
+ end do ! j=1,config_n_btr_subcycles
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! END Barotropic subcycle loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Note: Normalize so that \sum_{j=0}^{2J} a_j=1 and \sum_{j=1}^{2J} b_j=1. Then the following
-! three lines are already done, so that\\
-! \zeta^* = \left. \sum_{j=0}^{2J} a_j \zeta_{n+j/J} \right/ \sum_{j=0}^{2J} a_j is ! sshNew
-! \ubtr^* = \left. \sum_{j=0}^{2J} a_j \ubtr_{n+j/J} \right/ \sum_{j=0}^{2J} a_j is ! uBtrNew
-! {\overline {\bf F}} = \left. \sum_{j=1}^{2J} b_j {\bf F}_j \right/ \sum_{j=1}^{2J} b_j is ! FBtr
-! {\bf u}^*_k = \ubtr^* + {\bf u}^{'*}_k
-! uNew(iEdge) = uBtrNew(iEdge) + 0.5*(uBclOld(iEdge) + uBclNew(iEdge))
-! h_1^* = \Delta z_1 + \zeta^*
-! hNew(iCell) = hZlevel(1) + sshNew(iCell)
+
+ ! Normalize Barotropic subcycle sums: ssh, uBtr, and F
block => domain % blocklist
do while (associated(block))
-
-
-
-
-
-
do iEdge=1,block % mesh % nEdges
block % state % time_levs(1) % state % FBtr % array(iEdge) &
= block % state % time_levs(1) % state % FBtr % array(iEdge) &
@@ -846,23 +850,112 @@
/ (config_n_btr_subcycles + 1)
end do
+ block => block % next
+ end do ! block
+ ! boundary update on F
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field1dReal(domain % dminfo, &
+ block % state % time_levs(1) % state % FBtr % array(:), &
+ block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
block => block % next
end do ! block
-! boundary update on {\overline {\bf F}}
+ ! Check that \zeta_{n} + \Delta t \left( - </font>
<font color="gray">abla \cdot {\overline {\bf F}} \right)
block => domain % blocklist
do while (associated(block))
+ allocate(tend_ssh(block % mesh % nCells))
+
+ do iEdge=1,block % mesh % nEdges
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ tend_ssh(cell1) = tend_ssh(cell1) - block % state % time_levs(1) % state % FBtr % array(iEdge)
+ tend_ssh(cell2) = tend_ssh(cell2) + block % state % time_levs(1) % state % FBtr % array(iEdge)
+
+ end do
+
+ do iCell=1,block % mesh % nCells
+
+ block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ = block % state % time_levs(1) % state % ssh % array(iCell) &
+ + dt &
+ * tend_ssh(iCell) / block % mesh % areaCell % array (iCell)
+ end do
+ ! Now can compare sshSubcycleNEW (big step using summed fluxes) with
+ ! sshSubcycleOLD (individual steps to get there)
+
+ ! Correction velocity, u^{corr}:
+!u^{corr} = \left( {\overline {\bf F}}
+! - \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right) u_k^* \right)
+!\left/ \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right) \right.
+
+ do iEdge=1,block % mesh % nEdges
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ sshEdge = 0.5 &
+ *( block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
+
+ ! This is u*
+ uTemp(:) &
+ = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
+
+ uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * uTemp(1)
+ hSum = sshEdge + block % mesh % hZLevel % array(1)
+
+ do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
+ uhSum = uhSum + block % mesh % hZLevel % array(k) * uTemp(k)
+ hSum = hSum + block % mesh % hZLevel % array(k)
+ enddo
+ uCorr = (block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum
+
+ ! put u^{tr}, the velocity for tracer transport, in uNew
+ block % state % time_levs(2) % state % u % array(:,iEdge) = uTemp(:) + uCorr
+
+ ! Put new SSH values in h array, for the compute_scalar_tend call below.
+ block % state % time_levs(2) % state % h_edge % array(1,iEdge) &
+ = sshEdge + block % mesh % hZLevel % array(1)
+
+ do k=2,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % h_edge % array(k,iEdge) &
+ = block % mesh % hZLevel % array(k)
+ enddo
+
+
+ end do ! iEdge
+
+ ! Put new SSH values in h array, for the compute_scalar_tend call below.
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(2) % state % h % array(1,iCell) &
+ = block % state % time_levs(2) % state % ssh % array(iCell) &
+ + block % mesh % hZLevel % array(1)
+
+ ! mrp 110601 efficiency note: Since h just moves back and forth between pointers,
+ ! this is not necessary once initialized.
+ do k=2,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % h % array(k,iCell) &
+ = block % mesh % hZLevel % array(k)
+ enddo
+ enddo ! iCell
+
+
+
+ deallocate(tend_ssh)
block => block % next
end do ! block
- endif
+ endif ! higdon_split
@@ -875,15 +968,6 @@
block => domain % blocklist
do while (associated(block))
- do iEdge=1,block % mesh % nEdges
- ! compute u*, the velocity for tendency terms. Put in uNew.
- ! uBclNew is at time n+1/2 here.
- block % state % time_levs(2) % state % u % array(:,iEdge) &
- = block % state % time_levs(2) % state % uBtr % array(iEdge) &
- + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
-
- end do ! iEdge
-
call compute_wTop(block % state % time_levs(2) % state, block % mesh)
if (trim(config_time_integration) == 'higdon_unsplit') then
@@ -922,13 +1006,15 @@
block => domain % blocklist
do while (associated(block))
+
if (trim(config_time_integration) == 'higdon_unsplit') then
+ ! this is h_{n+1}
block % state % time_levs(2) % state % h % array(:,:) &
= block % state % time_levs(1) % state % h % array(:,:) &
+ dt* block % tend % h % array(:,:)
- endif
+ endif ! higdon_unsplit
! printing:
nCells = block % mesh % nCells
@@ -943,17 +1029,73 @@
! enddo
! printing end
- ! mrp 110512 at a later time, only need T & S for iterations
+ ! Only need T & S for earlier iterations,
+ ! then all the tracers needed the last time through.
+ if (higdon_step < config_n_ts_iter) then
+
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(1) % state % tracers % array(:,k,iCell) &
+ do i=1,2
+ ! This is Phi at n+1
+ tracerTemp &
+ = ( block % state % time_levs(1) % state % tracers % array(i,k,iCell) &
* block % state % time_levs(1) % state % h % array(k,iCell) &
- + dt * block % tend % tracers % array(:,k,iCell) &
+ + dt * block % tend % tracers % array(i,k,iCell) &
) / block % state % time_levs(2) % state % h % array(k,iCell)
+
+ ! This is Phi at n+1/2
+ block % state % time_levs(2) % state % tracers % array(i,k,iCell) &
+ = 0.5*( &
+ block % state % time_levs(1) % state % tracers % array(i,k,iCell) &
+ + tracerTemp )
+ enddo
end do
+ end do ! iCell
+
+
+ if (trim(config_time_integration) == 'higdon_unsplit') then
+
+ ! compute h*, which is h at n+1/2 and put into array hNew
+ ! on last iteration, hNew remains at n+1
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(2) % state % h % array(1,iCell) &
+ = 0.5*( &
+ block % state % time_levs(2) % state % h % array(1,iCell) &
+ + block % state % time_levs(1) % state % h % array(1,iCell) )
+
+ end do ! iCell
+ endif ! higdon_unsplit
+
+
+ elseif (higdon_step == config_n_ts_iter) then
+
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % maxLevelCell % array(iCell)
+ do i=1,block % state % time_levs(1) % state % num_tracers
+ ! This is Phi at n+1
+ block % state % time_levs(2) % state % tracers % array(i,k,iCell) &
+ = ( block % state % time_levs(1) % state % tracers % array(i,k,iCell) &
+ * block % state % time_levs(1) % state % h % array(k,iCell) &
+ + dt * block % tend % tracers % array(i,k,iCell) &
+ ) / block % state % time_levs(2) % state % h % array(k,iCell)
+
+ enddo
+ end do
end do
+ endif ! higdon_step
+
+
+ ! compute u*, the velocity for tendency terms. Put in uNew.
+ ! uBclNew is at time n+1/2 here.
+ ! This overwrites u^{tr}, the tracer transport velocity, which was in uNew.
+ ! The following must occur after call compute_scalar_tend
+ do iEdge=1,block % mesh % nEdges
+ block % state % time_levs(2) % state % u % array(:,iEdge) &
+ = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
+ end do ! iEdge
+
! mrp 110512 I really only need this to compute h_edge, density, pressure.
! I can par this down later.
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
@@ -1015,13 +1157,38 @@
= block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+2*block % state % time_levs(2) % state % uBcl % array(:,iEdge) &
- block % state % time_levs(1) % state % uBcl % array(:,iEdge)
- enddo
+ end do
+ if (trim(config_time_integration) == 'higdon_split') then
do iCell=1,block % mesh % nCells
+ ! SSH for the next step is from the end of the barotropic subcycle.
block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell)
- enddo
+ ! Put new SSH values in h array, for the compute_scalar_tend call below.
+ block % state % time_levs(2) % state % h % array(1,iCell) &
+ = block % state % time_levs(2) % state % ssh % array(iCell) &
+ + block % mesh % hZLevel % array(1)
+
+ ! mrp 110601 efficiency note: Since h just moves back and forth between pointers,
+ ! this is not necessary once initialized.
+ do k=2,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % h % array(k,iCell) &
+ = block % mesh % hZLevel % array(k)
+ end do
+ end do ! iCell
+ end if ! higdon_split
+
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % maxLevelCell % array(iCell)
+ block % state % time_levs(2) % state % u % array(k,iCell) &
+ = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iCell) &
+ +2*block % state % time_levs(2) % state % uBcl % array(k,iCell) &
+ - block % state % time_levs(1) % state % uBcl % array(k,iCell)
+ enddo
+
+ enddo ! iCell
+
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
</font>
</pre>