<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 :: &amp; 
         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) &amp; 
@@ -953,8 +953,8 @@
               cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
 
               sshEdge = 0.5 &amp;
-                 *(  block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &amp;
-                   + block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
+                 *(  block % state % time_levs(2) % state % ssh % array(cell1) &amp;
+                   + block % state % time_levs(2) % state % ssh % array(cell2) )
 
              ! This is u*
                uTemp(:) &amp;
@@ -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) &amp;
+                     /block % mesh % dvEdge % array(iEdge) &amp;
+                      - 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) &amp;
@@ -1062,8 +1069,8 @@
 
         block =&gt; 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 &lt; 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) &amp;
+                + 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) &amp;
                    * block % state % time_levs(1) % state % h % array(k,iCell) &amp;
                  + dt * block % tend % tracers % array(i,k,iCell) &amp;
-                  ) / 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)  &amp;
@@ -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) &amp;
+                + 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) &amp;
                    * block % state % time_levs(1) % state % h % array(k,iCell) &amp;
                  + dt * block % tend % tracers % array(i,k,iCell) &amp;
-                  ) / block % state % time_levs(2) % state % h % array(k,iCell)
+                  ) / hNew(k)
 
                  enddo
               end do
            end do
 
          endif ! higdon_step
+           deallocate(hNew)
 
          block =&gt; block % next
        end do

</font>
</pre>