<p><b>mpetersen@lanl.gov</b> 2011-03-15 09:04:22 -0600 (Tue, 15 Mar 2011)</p><p>BRANCH COMMIT: Added and tested tridiagonal solver.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F        2011-03-14 21:57:46 UTC (rev 758)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F        2011-03-15 15:04:22 UTC (rev 759)
@@ -75,6 +75,7 @@
       integer :: rk_step, iEdge
 
       real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
+      real (kind=RKIND), dimension(:), allocatable:: A,C,R,X
 
       block =&gt; domain % blocklist
       call allocate_state(provis, &amp;
@@ -231,12 +232,25 @@
       do while (associated(block))
 
          if (config_implicit_vertical_mix) then
+            allocate(A(block % mesh % nVertLevels),C(block % mesh % nVertLevels),R(block % mesh % nVertLevels),X(block % mesh % nVertLevels))
 
             call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
 
             do iCell=1,block % mesh % nCells
-               !N=maxLevelCell(iCell)
 
+ C=-2
+ do k=1,10
+    A(k) = 9.0 + k/10.0
+    R(k)=k
+ enddo
+
+ call tridiagonal_solve(A,C,A,R,X,10)
+print '(a,10f10.4)', 'A', A(1:10)
+print '(a,10f10.4)', 'C', C(1:10)
+print '(a,10f10.4)', 'R', R(1:10)
+print '(a,10f10.4)', 'X', X(1:10)
+stop
+
             !   do k=1,maxLevelCell(iCell)
                ! Compute A(k), C(k), R(k) for momentum
              !  enddo
@@ -247,6 +261,7 @@
                !enddo
                !call tridiagonal_solve_mult(tracers(:,:,iCell),A(2:N),C(1:N),A(1:N-1),R)
             end do
+            deallocate(A,C,R,X)
          else
             do iCell=1,block % mesh % nCells
                do k=1,block % mesh % maxLevelCell % array(iCell)
@@ -2342,4 +2357,48 @@
 
    end subroutine equation_of_state_jm
 
+subroutine tridiagonal_solve(a,b,c,r,x,n)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Solve the matrix equation Ax=r for x, where A is tridiagonal.
+! A is an nxn matrix, with:
+!   a sub-diagonal, filled from 2:n
+!   b diagonal, filled from 1:n
+!   c sup-diagonal, filled from 1:n-1
+!
+! Input: a,b,c,r,n
+!
+! Output: x
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ 
+   implicit none
+
+   integer,intent(in) :: n
+   real (KIND=RKIND), dimension(n), intent(in) :: a,b,c,r
+   real (KIND=RKIND), dimension(n), intent(out) :: x
+   real (KIND=RKIND), dimension(n) :: bTemp,rTemp
+   real (KIND=RKIND) :: m
+   integer i
+ 
+   ! Use work variables for b and r
+   bTemp(1) = b(1)
+   rTemp(1) = r(1)
+ 
+   ! First pass: set the coefficients
+   do i = 2,n
+      m = a(i)/bTemp(i-1)
+      bTemp(i) = b(i) - m*c(i-1)
+      rTemp(i) = r(i) - m*rTemp(i-1)
+   end do 
+print '(a,10f8.2)', 'bTemp', bTemp(1:n)
+print '(a,10f8.2)', 'rTemp', rTemp(1:n)
+ 
+   x(n) = rTemp(n)/bTemp(n)
+   ! Second pass: back-substition
+   do i = n-1, 1, -1
+      x(i) = (rTemp(i) - c(i)*x(i+1))/bTemp(i)
+   end do
+ 
+end subroutine tridiagonal_solve
+
+
 end module time_integration

</font>
</pre>