<p><b>mpetersen@lanl.gov</b> 2011-06-06 16:40:49 -0600 (Mon, 06 Jun 2011)</p><p>Changed h divisor in tracer equation to be at time n+1. Split version is working for short times now, but still has trouble after hundreds of timesteps.<br>
</p><hr noshade><pre><font color="gray">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-06-06 22:15:39 UTC (rev 878)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F        2011-06-06 22:40:49 UTC (rev 879)
@@ -385,7 +385,7 @@
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
- real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp
+ real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, hNew
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
!print *, '1'
@@ -799,7 +799,7 @@
! mrp 110606 efficiency note: could make this a 1D integer factor instead of an if statement.
if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
- block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) = 0
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) = 0.0
else
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
@@ -953,8 +953,8 @@
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) )
+ *( block % state % time_levs(2) % state % ssh % array(cell1) &
+ + block % state % time_levs(2) % state % ssh % array(cell2) )
! This is u*
uTemp(:) &
@@ -970,12 +970,19 @@
enddo
! mrp temp 110602 testing split: remove correction term.
-! uCorr = (block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum
- uCorr = 0.0
+ uCorr = ( block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ /block % mesh % dvEdge % array(iEdge) &
+ - uhSum)/hSum
+! uCorr = 0.0
! mrp temp 110602 testing split end
! put u^{tr}, the velocity for tracer transport, in uNew
+ ! mrp 060611 not sure if boundary enforcement is needed here.
+ if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
+ block % state % time_levs(2) % state % u % array(:,iEdge) = 0.0
+ else
block % state % time_levs(2) % state % u % array(:,iEdge) = uTemp(:) + uCorr
+ endif
! 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) &
@@ -1062,8 +1069,8 @@
block => domain % blocklist
do while (associated(block))
+ allocate(hNew(nVertLevels))
-
if (trim(config_time_integration) == 'higdon_unsplit') then
! this is h_{n+1}
@@ -1090,7 +1097,10 @@
! then all the tracers needed the last time through.
if (higdon_step < config_n_ts_iter) then
+ hNew(:) = block % mesh % hZLevel % array(:)
do iCell=1,block % mesh % nCells
+ hNew(1) = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ + block % mesh % hZLevel % array(1)
do k=1,block % mesh % maxLevelCell % array(iCell)
do i=1,2
! This is Phi at n+1
@@ -1098,7 +1108,7 @@
= ( 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)
+ ) / hNew(k)
! This is Phi at n+1/2
block % state % time_levs(2) % state % tracers % array(i,k,iCell) &
@@ -1140,7 +1150,10 @@
elseif (higdon_step == config_n_ts_iter) then
+ hNew(:) = block % mesh % hZLevel % array(:)
do iCell=1,block % mesh % nCells
+ hNew(1) = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ + block % mesh % hZLevel % array(1)
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
@@ -1148,13 +1161,14 @@
= ( 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)
+ ) / hNew(k)
enddo
end do
end do
endif ! higdon_step
+ deallocate(hNew)
block => block % next
end do
</font>
</pre>