<p><b>mpetersen@lanl.gov</b> 2011-02-08 09:40:00 -0700 (Tue, 08 Feb 2011)</p><p>BRANCH COMMIT: consolidate do loops for vertical advection.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-07 23:04:29 UTC (rev 716)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-08 16:40:00 UTC (rev 717)
@@ -700,7 +700,7 @@
type (mesh_type), intent(in) :: grid
integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&
- nEdges, nCells, nVertLevels, num_tracers
+ nEdges, nCells, nCellsSolve, nVertLevels, num_tracers
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND), dimension(:), pointer :: &
@@ -719,8 +719,8 @@
hRatioZLevelK, hRatioZLevelKm1
real (kind=RKIND), dimension(:), allocatable:: vertDiffTop, tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
posZTopZLevel
- real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop, boundaryMask
- real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer
+ real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer, tracerTop
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
@@ -751,6 +751,7 @@
nEdges = grid % nEdges
nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
num_tracers = s % num_tracers
@@ -894,30 +895,21 @@
!
if (config_vert_grid_type.eq.'zlevel') then
- allocate(tracerTop(num_tracers,nVertLevels+1))
- tracerTop(:,1)=0.0
+ allocate(tracerTop(num_tracers,nVertLevels+1,nCells))
if (config_vert_tracer_adv.eq.'stencil'.and. &
config_vert_tracer_adv_order.eq.2) then
! Compute tracerTop using centered stencil, a simple average.
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
- tracerTop(iTracer,k) = &
+ tracerTop(iTracer,k,iCell) = &
( tracers(iTracer,k-1,iCell) &
+tracers(iTracer,k ,iCell))/2.0
end do
end do
- tracerTop(:,maxLevelCell(iCell)+1)=0.0
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
- end do
- end do
end do
elseif (config_vert_tracer_adv.eq.'stencil'.and. &
@@ -929,17 +921,17 @@
! Hardwire flux3Coeff at 1.0 for now. Could add this to the
! namelist, if desired.
flux3Coef = 1.0
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
k=2
do iTracer=1,num_tracers
- tracerTop(iTracer,k) = &
+ tracerTop(iTracer,k,iCell) = &
hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
end do
do k=3,maxLevelCell(iCell)-1
cSignWTop = sign(flux3Coef,wTop(k,iCell))
do iTracer=1,num_tracers
- tracerTop(iTracer,k) = &
+ tracerTop(iTracer,k,iCell) = &
( (-1.+ cSignWTop)*tracers(iTracer,k-2,iCell) &
+( 7.-3.*cSignWTop)*tracers(iTracer,k-1,iCell) &
+( 7.+3.*cSignWTop)*tracers(iTracer,k ,iCell) &
@@ -949,18 +941,10 @@
end do
k=maxLevelCell(iCell)
do iTracer=1,num_tracers
- tracerTop(iTracer,k) = &
+ tracerTop(iTracer,k,iCell) = &
hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
end do
- tracerTop(:,maxLevelCell(iCell)+1)=0.0
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
- end do
- end do
end do
elseif (config_vert_tracer_adv.eq.'stencil'.and. &
@@ -968,16 +952,16 @@
! Compute tracerTop using 4rd order stencil [-1 7 7 -1]
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
k=2
do iTracer=1,num_tracers
- tracerTop(iTracer,k) = &
+ tracerTop(iTracer,k,iCell) = &
hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
end do
do k=3,maxLevelCell(iCell)-1
do iTracer=1,num_tracers
- tracerTop(iTracer,k) = &
+ tracerTop(iTracer,k,iCell) = &
(- tracers(iTracer,k-2,iCell) &
+7.*tracers(iTracer,k-1,iCell) &
+7.*tracers(iTracer,k ,iCell) &
@@ -987,18 +971,10 @@
end do
k=maxLevelCell(iCell)
do iTracer=1,num_tracers
- tracerTop(iTracer,k) = &
+ tracerTop(iTracer,k,iCell) = &
hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
end do
- tracerTop(:,maxLevelCell(iCell)+1)=0.0
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
- end do
- end do
end do
elseif (config_vert_tracer_adv.eq.'spline'.and. &
@@ -1006,24 +982,16 @@
! Compute tracerTop using linear interpolation.
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
! Note hRatio on the k side is multiplied by tracer at k-1
! and hRatio on the Km1 (k-1) side is mult. by tracer at k.
- tracerTop(iTracer,k) = &
+ tracerTop(iTracer,k,iCell) = &
hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
end do
end do
- tracerTop(:,maxLevelCell(iCell)+1)=0.0
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
- end do
- end do
end do
elseif (config_vert_tracer_adv.eq.'spline'.and. &
@@ -1041,7 +1009,7 @@
posZMidZLevel = -zMidZLevel(1:nVertLevels)
posZTopZLevel = -zTopZLevel(2:nVertLevels)
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
! mrp 110201 efficiency note: push tracer loop down
! into spline subroutines to improve efficiency
do iTracer=1,num_tracers
@@ -1059,17 +1027,9 @@
posZMidZLevel, posZTopZLevel, &
maxLevelCell(iCell), maxLevelCell(iCell)-1 )
- tracerTop(iTracer,2:maxLevelCell(iCell)) = tracersOut(1:maxLevelCell(iCell)-1)
+ tracerTop(iTracer,2:maxLevelCell(iCell),iCell) = tracersOut(1:maxLevelCell(iCell)-1)
end do
- tracerTop(:,maxLevelCell(iCell)+1)=0.0
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
- end do
- end do
end do
deallocate(tracer2ndDer)
@@ -1086,6 +1046,18 @@
endif ! vertical tracer advection method
+ do iCell=1,nCellsSolve
+ tracerTop(:,1,iCell)=0.0
+ tracerTop(:,maxLevelCell(iCell)+1,iCell)=0.0
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ - ( wTop(k ,iCell)*tracerTop(iTracer,k ,iCell) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1,iCell))
+ end do
+ end do
+ end do
+
deallocate(tracerTop)
endif ! ZLevel
@@ -1228,7 +1200,7 @@
allocate(fluxVertTop(num_tracers,nVertLevels+1))
fluxVertTop(:,1) = 0.0
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
</font>
</pre>