<p><b>mpetersen@lanl.gov</b> 2011-03-24 15:42:29 -0600 (Thu, 24 Mar 2011)</p><p>BRANCH COMMIT Update to implicit vertical mix branch. Code works for double-gyre with inverse stratification. This commit includes print lines for debugging.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry 2011-03-24 18:27:00 UTC (rev 767)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/Registry 2011-03-24 21:42:29 UTC (rev 768)
@@ -199,9 +199,9 @@
var persistent real CFLNumberGlobal ( Time ) 2 o CFLNumberGlobal state - -
# Diagnostics fields, only one time level required
-var persistent real RiTopOfCell ( nVertLevelsP1 nCells ) 0 o RiTopOfCell diagnostics - -
-var persistent real RiTopOfEdge ( nVertLevelsP1 nEdges ) 0 o RiTopOfEdge diagnostics - -
-var persistent real vertViscTopOfEdge ( nVertLevelsP1 nEdges ) 0 o vertViscTopOfEdge diagnostics - -
-var persistent real vertDiffTopOfCell ( nVertLevelsP1 nCells ) 0 o vertDiffTopOfCell diagnostics - -
+var persistent real RiTopOfCell ( nVertLevelsP1 nCells Time ) 1 o RiTopOfCell diagnostics - -
+var persistent real RiTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 o RiTopOfEdge diagnostics - -
+var persistent real vertViscTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 o vertViscTopOfEdge diagnostics - -
+var persistent real vertDiffTopOfCell ( nVertLevelsP1 nCells Time ) 1 o vertDiffTopOfCell diagnostics - -
Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_mpas_core.F 2011-03-24 18:27:00 UTC (rev 767)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_mpas_core.F 2011-03-24 21:42:29 UTC (rev 768)
@@ -81,10 +81,38 @@
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+! mrp causes error Subscript #1 of the array INDX has value 7 which is greater than the upper bound of 6
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
call reconstruct(block % state % time_levs(1) % state, mesh)
+ ! mrp 110131 testing only
+ block % mesh % u_src % array(:,:) = &
+ block % mesh % u_src % array(:,:) *100
+ print *, 'min,maxval(u_src) ',minval( block % mesh % u_src % array(:,:)),maxval( block % mesh % u_src % array(:,:))
+
+ do iCell=1,block % mesh % nCells
+ do k=1,mesh % nVertLevels
+ !block % state % time_levs(1) % state % tracers % array( &
+ ! 1, k ,iCell) = 31-k/1000.0 !31-k
+ !block % state % time_levs(1) % state % tracers % array( &
+ ! 2, k ,iCell) = 32 !+ k*0.32
+ !block % state % time_levs(1) % state % tracers % array( &
+ ! 3, k ,iCell) = k
+ enddo
+ block % state % time_levs(1) % state % tracers % array( &
+ 3, :,iCell) = 1.
+ block % state % time_levs(1) % state % tracers % array( &
+ 4, :,iCell) = 1.
+
+ do k=1,mesh % nVertLevels,2
+ block % state % time_levs(1) % state % tracers % array( &
+ 4, k ,iCell) = 0.
+ enddo
+
+ end do
+ ! mrp 110131 testing only end
+
! initialize velocities and tracers on land to be -1e34
! The reconstructed velocity on land will have values not exactly
! -1e34 due to the interpolation of reconstruction.
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-24 18:27:00 UTC (rev 767)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F 2011-03-24 21:42:29 UTC (rev 768)
@@ -72,11 +72,20 @@
type (block_type), pointer :: block
type (state_type) :: provis
- integer :: rk_step, iEdge
+ integer :: rk_step, iEdge, cell1, cell2
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
- real (kind=RKIND), dimension(:), allocatable:: A,C,R,X
+ 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
+
+
block => domain % blocklist
call allocate_state(provis, &
block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
@@ -231,45 +240,116 @@
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(block % mesh % nVertLevels),C(block % mesh % nVertLevels),R(block % mesh % nVertLevels),X(block % mesh % nVertLevels))
+ 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)
+ call compute_vertical_mix_coefficients(provis, block % diagnostics, block % mesh)
- do iCell=1,block % mesh % nCells
+!print '(a,10es12.4)', 'u1', minval(u(:,1:nEdges)),maxval(u(:,1:nEdges))
+! should I use nEdges or nEdgesSolve here?
+ do iEdge=1,nEdges
+ if (maxLevelEdgeTop(iEdge).gt.0) then
- C=-2
- do k=1,10
- A(k) = 9.0 + k/10.0
- R(k)=k
- enddo
+ ! 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
- 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,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
- ! do k=1,maxLevelCell(iCell)
- ! Compute A(k), C(k), R(k) for momentum
- ! enddo
- !call tridiagonal_solve(u(:,iCell),A(2:N),C(1:N),A(1:N-1),R)
+ C(1) = 1 - A(1)
+ do k=2,maxLevelEdgeTop(iEdge)
+ C(k) = 1 - A(k) - A(k-1)
+ enddo
- ! do k=1,maxLevelCell(iCell)
- ! Compute A(k), C(k), R(iTracer,k) for tracers
- !enddo
- !call tridiagonal_solve_mult(tracers(:,:,iCell),A(2:N),C(1:N),A(1:N-1),R)
+ call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
+
+!print '(a,2i4)', 'iEdge, maxLevelEdgeTop(iEdge))', iEdge,maxLevelEdgeTop(iEdge)
+!print '(a,100es12.4)', 'vertViscTopOfEdge(:,iEdge)', vertViscTopOfEdge(:,iEdge)
+!print '(a,100es12.4)', 'h_edge', h_edge(:,iEdge)
+!print '(a,100es12.4)', 'dt', dt
+!print '(a,100es12.4)', 'A', A
+!print '(a,100es12.4)', 'C', C
+!print '(a,100es12.4)', 'u(:,iEdge)', u(:,iEdge)
+!print '(a,100es12.4)', 'uTemp ', uTemp
+!stop
+
+! mrp temp, turn off implicit on u:
+ u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
+ u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
+
+ end if
end do
- deallocate(A,C,R,X)
- else
- 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) &
- / block % state % time_levs(2) % state % h % array(k,iCell)
- end do
+!print '(a,10es12.4)', 'u2', minval(u(:,1:nEdges)),maxval(u(:,1:nEdges))
+
+
+ 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)
+
+!mrp temp
+ tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
+! mrp temp remove soon:
+! do k=1,maxLevelCell(iCell)
+! tracers(:,k,iCell) = tracers(:,k,iCell) / h(k,iCell)
+! end do
+! mrp temp end
+! mrp set as follows
+! tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
+ tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = 0
+
end do
+!print '(a,10es12.4)', 'tracer', minval(tracers(:,:,1:nCells)),maxval(tracers(:,:,1:nCells))
+ deallocate(A,C,uTemp,tracersTemp)
end if
if (config_test_case == 1) then ! For case 1, wind field should be fixed
@@ -278,6 +358,7 @@
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+! mrp causes error Subscript #1 of the array INDX has value 7 which is greater than the upper bound of 6
call reconstruct(block % state % time_levs(2) % state, block % mesh)
block => block % next
@@ -681,6 +762,8 @@
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
+!mrp indent later
+ if (.not.config_implicit_vertical_mix) then
allocate(fluxVertTop(nVertLevels+1))
fluxVertTop(1) = 0.0
do iEdge=1,grid % nEdgesSolve
@@ -700,7 +783,7 @@
end do
deallocate(fluxVertTop)
-
+ endif
end subroutine compute_tend
@@ -1207,6 +1290,8 @@
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
+!mrp indent later
+ if (.not.config_implicit_vertical_mix) then
allocate(fluxVertTop(num_tracers,nVertLevels+1))
fluxVertTop(:,1) = 0.0
do iCell=1,nCellsSolve
@@ -1232,6 +1317,7 @@
enddo ! iCell loop
deallocate(fluxVertTop)
+ endif
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
@@ -1638,6 +1724,11 @@
! For an isopycnal model, density should remain constant.
if (config_vert_grid_type.eq.'zlevel') then
call equation_of_state(s,grid,-1)
+
+! mrp temp, just to visualize rhoDisplaced
+ call equation_of_state(s, grid, 1)
+! mrp temp, just to visualize rhoDisplaced end
+
endif
!
@@ -1794,17 +1885,25 @@
real (kind=RKIND) :: r, h1, h2
rhoDisplaced => s % rhoDisplaced % array
+ u => s % u % array
+ h => s % h % array
+ h_edge => s % h_edge % array
vertViscTopOfEdge => d % vertViscTopOfEdge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
RiTopOfEdge => d % RiTopOfEdge % array
RiTopOfCell => d % RiTopOfCell % array
- h => s % h % array
- u => s % u % array
+
+!
+!
+! mmmmmmmmmmmmmmmmmmmmmmm mrp delete variables later
+!
+!
+
+
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
@@ -1846,6 +1945,7 @@
maxLevelVertexBot => grid % maxLevelVertexBot % array
maxLevelVertexTop => grid % maxLevelVertexTop % array
+
nCells = grid % nCells
nEdges = grid % nEdges
nVertices = grid % nVertices
@@ -1867,6 +1967,11 @@
! compute density of parcel displaced to surface, $\rho^*$,
! in state % rhoDisplaced
call equation_of_state(s, grid, 1)
+! mrp delete:
+!print '(a,10es12.4)', 'temperature', minval(tracers(1,2:10,1:nCells)),maxval(tracers(1,2:10,1:nCells))
+!print '(a,10es12.4)', 'salinity', minval(tracers(2,2:10,1:nCells)),maxval(tracers(2,2:10,1:nCells))
+!print '(a,10es12.4)', 'rhoDisplaced', minval(rhoDisplaced(2:10,1:nCells)),maxval(rhoDisplaced(2:10,1:nCells))
+!print '(a,10es12.4)', 'h_edge', minval(h_edge(2:10,1:nEdges)),maxval(h_edge(2:10,1:nEdges))
! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
drhoTopOfCell = 0.0
@@ -1915,6 +2020,8 @@
!test this interp
+! are there issues with updating proc boundaries here?
+
! compute RiTopOfEdge using drhoTopOfEdge and du2TopOfEdge
! coef = -g/rho_0/2
RiTopOfEdge = 0.0
@@ -1939,6 +2046,13 @@
end do
end do
+!print '(a,10es12.4)', 'RiTopOfCell', minval(RiTopOfCell(2:10,1:nCells)),maxval(RiTopOfCell(2:10,1:nCells))
+!print '(a,10es12.4)', 'RiTopOfEdge', minval(RiTopOfEdge(2:10,1:nEdges)),maxval(RiTopOfEdge(2:10,1:nEdges))
+!print *, 'gravity',gravity
+!print '(a,10es12.4)', 'drhoTopOfCell', minval(drhoTopOfCell(2:10,1:nCells)),maxval(drhoTopOfCell(2:10,1:nCells))
+!print '(a,10es12.4)', 'du2TopOfCell', minval(du2TopOfCell(2:10,1:nCells)),maxval(du2TopOfCell(2:10,1:nCells))
+!print '(a,10es12.4)', 'h', minval(h(2:10,1:nCells)),maxval(h(2:10,1:nCells))
+
! what about boundary update?
@@ -1947,6 +2061,9 @@
du2TopOfCell, du2TopOfEdge)
endif
+ !
+ ! Compute vertical viscosity
+ !
if (config_vert_visc_type.eq.'const') then
vertViscTopOfEdge = config_vert_viscosity
@@ -1971,8 +2088,18 @@
vertViscTopOfEdge = 0.0
do iEdge = 1,nEdges
do k = 2,maxLevelEdgeTop(iEdge)
- vertViscTopOfEdge(k,iEdge) = config_bkrd_vert_visc &
- + config_rich_mix / (1.0 + 5.0*RiTopOfEdge(k,iEdge))**2
+ ! mrp 110324 efficiency note: this if is inside iCell and k loops.
+ ! Perhaps there is a more efficient way to do this.
+ if (RiTopOfCell(k,iCell)>0.0) then
+ vertViscTopOfEdge(k,iEdge) = config_bkrd_vert_visc &
+ + config_rich_mix / (1.0 + 5.0*RiTopOfEdge(k,iEdge))**2
+ else
+ ! for Ri<0, use maximum diffusion allowed by CFL criterion
+ ! mrp 110324 efficiency note: for z-level, could use fixed
+ ! grid array hMeanTopZLevel and compute maxdiff on startup.
+ vertViscTopOfEdge(k,iEdge) = &
+ (h_edge(k-1,iEdge)+h_edge(k,iEdge))**2/config_dt/4.0
+ endif
end do
end do
@@ -1983,6 +2110,9 @@
endif
+ !
+ ! Compute vertical tracer diffusion
+ !
if (config_vert_diff_type.eq.'const') then
vertDiffTopOfCell = config_vert_diffusion
@@ -2008,10 +2138,20 @@
coef = -gravity/1000.0/2.0
do iCell = 1,nCells
do k = 2,maxLevelCell(iCell)
- vertDiffTopOfCell(k,iCell) = config_bkrd_vert_diff &
- + (config_bkrd_vert_visc &
- + config_rich_mix / (1.0 + 5.0*RiTopOfCell(k,iCell))**2) &
- / (1.0 + 5.0*RiTopOfCell(k,iCell))**2
+ ! mrp 110324 efficiency note: this if is inside iCell and k loops.
+ ! Perhaps there is a more efficient way to do this.
+ if (RiTopOfCell(k,iCell)>0.0) then
+ vertDiffTopOfCell(k,iCell) = config_bkrd_vert_diff &
+ + (config_bkrd_vert_visc &
+ + config_rich_mix / (1.0 + 5.0*RiTopOfCell(k,iCell))**2) &
+ / (1.0 + 5.0*RiTopOfCell(k,iCell))
+ else
+ ! for Ri<0, use maximum diffusion allowed by CFL criterion
+ ! mrp 110324 efficiency note: for z-level, could use fixed
+ ! grid array hMeanTopZLevel and compute maxdiff on startup.
+ vertDiffTopOfCell(k,iCell) = &
+ (h(k-1,iCell)+h(k,iCell))**2/config_dt/4.0
+ end if
end do
end do
@@ -2022,6 +2162,8 @@
endif
+!print '(a,10es12.4)', 'vertDiffTopOfCell', minval(vertDiffTopOfCell(2:10,1:nCells)),maxval(vertDiffTopOfCell(2:10,1:nCells))
+!print '(a,10es12.4)', 'vertViscTopOfEdge', minval(vertViscTopOfEdge(2:10,1:nEdges)),maxval(vertViscTopOfEdge(2:10,1:nEdges))
end subroutine compute_vertical_mix_coefficients
@@ -2349,21 +2491,36 @@
DENOMK = 1.0/(BULK_MOD - p(k))
rhoPointer(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+!mrp
+! if (k_displaced.gt.0) then
+!print '(a,2i5,10es22.14)', 'eos icell,k,T,S,rho', icell,k,TQ,SQ,rhoPointer(k,iCell)
+! endif
+!mrp
end do
end do
+!mrp
+! if (k_displaced.gt.0) then
+!print '(a,10es22.14)', 'eos temperature', minval(tracers(1,2:10,1:nCells)),maxval(tracers(1,2:10,1:nCells))
+!print '(a,10es22.14)', 'eos salinity', minval(tracers(2,2:10,1:nCells)),maxval(tracers(2,2:10,1:nCells))
+!print '(a,10es22.14)', 'eos rhoDisplaced', minval(rhoPointer(2:10,1:nCells)),maxval(rhoPointer(2:10,1:nCells))
+! stop
+! endif
+!mrp
+
deallocate(pRefEOS,p,p2)
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
+! a sub-diagonal, filled from 1:n-1 (a(1) appears on row 2)
! b diagonal, filled from 1:n
-! c sup-diagonal, filled from 1:n-1
+! c sup-diagonal, filled from 1:n-1 (c(1) apears on row 1)
!
! Input: a,b,c,r,n
!
@@ -2385,12 +2542,10 @@
! First pass: set the coefficients
do i = 2,n
- m = a(i)/bTemp(i-1)
+ m = a(i-1)/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
@@ -2401,4 +2556,56 @@
end subroutine tridiagonal_solve
+subroutine tridiagonal_solve_mult(a,b,c,r,x,n,nDim,nSystems)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Solve the matrix equation Ax=r for x, where A is tridiagonal.
+! A is an nxn matrix, with:
+! a sub-diagonal, filled from 1:n-1 (a(1) appears on row 2)
+! b diagonal, filled from 1:n
+! c sup-diagonal, filled from 1:n-1 (c(1) apears on row 1)
+!
+! Input: a,b,c,r,n
+!
+! Output: x
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer,intent(in) :: n, nDim, nSystems
+ real (KIND=RKIND), dimension(n), intent(in) :: a,b,c
+ real (KIND=RKIND), dimension(nSystems,nDim), intent(in) :: r
+ real (KIND=RKIND), dimension(nSystems,nDim), intent(out) :: x
+ real (KIND=RKIND), dimension(n) :: bTemp
+ real (KIND=RKIND), dimension(nSystems,n) :: rTemp
+ real (KIND=RKIND) :: m
+ integer i,j
+
+ ! Use work variables for b and r
+ bTemp(1) = b(1)
+ do j = 1,nSystems
+ rTemp(j,1) = r(j,1)
+ end do
+
+ ! First pass: set the coefficients
+ do i = 2,n
+ m = a(i-1)/bTemp(i-1)
+ bTemp(i) = b(i) - m*c(i-1)
+ do j = 1,nSystems
+ rTemp(j,i) = r(j,i) - m*rTemp(j,i-1)
+ end do
+ end do
+
+ do j = 1,nSystems
+ x(j,n) = rTemp(j,n)/bTemp(n)
+ end do
+ ! Second pass: back-substition
+ do i = n-1, 1, -1
+ do j = 1,nSystems
+ x(j,i) = (rTemp(j,i) - c(i)*x(j,i+1))/bTemp(i)
+ end do
+ end do
+
+end subroutine tridiagonal_solve_mult
+
+
end module time_integration
</font>
</pre>