<p><b>mpetersen@lanl.gov</b> 2010-06-09 11:33:50 -0600 (Wed, 09 Jun 2010)</p><p>Added variables maxLevelCell, maxLevelEdge, maxLevelVertex for z-level topography. Change from nVertLevels in compute_tend and compute_scalar_tend.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/Registry 2010-06-07 17:49:33 UTC (rev 339)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/Registry 2010-06-09 17:33:50 UTC (rev 340)
@@ -76,6 +76,10 @@
var integer nEdgesOnEdge ( nEdges ) iro nEdgesOnEdge - -
var integer edgesOnCell ( maxEdges nCells ) iro edgesOnCell - -
var integer edgesOnEdge ( maxEdges2 nEdges ) iro edgesOnEdge - -
+# var integer maxVertLevel ( nCells ) iro maxVertLevel - -
+var integer maxLevelCell ( nCells ) iro maxLevelCell - -
+var integer maxLevelEdge ( nEdges ) ro maxLevelEdge - -
+var integer maxLevelVertex ( nVertices ) ro maxLevelVertex - -
var real weightsOnEdge ( maxEdges2 nEdges ) iro weightsOnEdge - -
var real dvEdge ( nEdges ) iro dvEdge - -
Modified: branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F 2010-06-07 17:49:33 UTC (rev 339)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F 2010-06-09 17:33:50 UTC (rev 340)
@@ -26,13 +26,17 @@
type (dm_info) :: dminfo
! mrp 100507: for diagnostic output
- integer :: iTracer
+ integer :: iTracer, cell1, cell2, cell
real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
hZLevel, zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: delta_rho
- integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+ integer, dimension(:), pointer :: &
+ maxLevelCell, maxLevelEdge, maxLevelVertex
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex
! mrp 100507 end: for diagnostic output
if (config_test_case == 0) then
@@ -110,11 +114,17 @@
hZLevel => block_ptr % mesh % hZLevel % array
zMidZLevel => block_ptr % mesh % zMidZLevel % array
zTopZLevel => block_ptr % mesh % zTopZLevel % array
+ maxLevelCell => block_ptr % mesh % maxLevelCell % array
+ maxLevelEdge => block_ptr % mesh % maxLevelEdge % array
+ maxLevelVertex => block_ptr % mesh % maxLevelVertex % array
+ cellsOnEdge => block_ptr % mesh % cellsOnEdge % array
+ cellsOnVertex => block_ptr % mesh % cellsOnVertex % array
nCells = block_ptr % mesh % nCells
nEdges = block_ptr % mesh % nEdges
nVertices = block_ptr % mesh % nVertices
nVertLevels = block_ptr % mesh % nVertLevels
+ vertexDegree = block_ptr % mesh % vertexDegree
if (config_vert_grid_type.eq.'zlevel') then
! These should eventually be in an input file. For now
@@ -160,6 +170,40 @@
endif
+ ! mrp 100608 For now, use a flat bottom
+ maxLevelCell = nVertLevels
+
+ ! maxLevelEdge is the minimum of the surrounding cells
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ maxLevelEdge(iEdge) = &
+ min( maxLevelCell(cell1), &
+ maxLevelCell(cell2) )
+ else
+ maxLevelEdge(iEdge) = 0
+ endif
+
+ end do
+
+ ! maxLevelVertex is the minimum of the surrounding cells
+ VTX_LOOP: do iVtx = 1,nVertices
+ maxLevelVertex(iVtx) = nVertLevels
+ do i=1,vertexDegree
+ cell = cellsOnVertex(i,iVtx)
+ if (cell > nCells) then
+ maxLevelVertex(iVtx) = 0
+ cycle VTX_LOOP
+ else
+ maxLevelVertex(iVtx) = &
+ min( maxLevelVertex(iVtx), &
+ maxLevelCell(cell) )
+ endif
+ end do
+ end do VTX_LOOP
+
! print some diagnostics
print '(10a)', 'ilevel',&
' rho ',&
Modified: branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F 2010-06-07 17:49:33 UTC (rev 339)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F 2010-06-09 17:33:50 UTC (rev 340)
@@ -89,7 +89,7 @@
block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
block % time_levs(2) % state % h % array(:,:) = block % time_levs(1) % state % h % array(:,:)
do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
block % time_levs(2) % state % tracers % array(:,k,iCell) = block % time_levs(1) % state % tracers % array(:,k,iCell) &
* block % time_levs(1) % state % h % array(k,iCell)
end do
@@ -162,7 +162,7 @@
block % intermediate_step(PROVIS) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) &
+ rk_substep_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
block % intermediate_step(PROVIS) % tracers % array(:,k,iCell) = ( &
block % time_levs(1) % state % h % array(k,iCell) * &
block % time_levs(1) % state % tracers % array(:,k,iCell) &
@@ -191,7 +191,7 @@
+ rk_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
block % time_levs(2) % state % tracers % array(:,k,iCell) = &
block % time_levs(2) % state % tracers % array(:,k,iCell) &
+ rk_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell)
@@ -212,7 +212,7 @@
block => domain % blocklist
do while (associated(block))
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
block % time_levs(2) % state % tracers % array(:,k,iCell) = &
block % time_levs(2) % state % tracers % array(:,k,iCell) &
/ block % time_levs(2) % state % h % array(k,iCell)
@@ -263,7 +263,8 @@
divergence, MontPot, pZLevel, zMidEdge, zTopEdge
type (dm_info) :: dminfo
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdge, maxLevelVertex
integer, dimension(:,:), pointer :: &
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
edgesOnEdge, edgesOnVertex
@@ -308,6 +309,9 @@
fEdge => grid % fEdge % array
zMidZLevel => grid % zMidZLevel % array
zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdge => grid % maxLevelEdge % array
+ maxLevelVertex => grid % maxLevelVertex % array
tend_h => tend % h % array
tend_u => tend % u % array
@@ -329,20 +333,20 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelCell(cell1)
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell1) = tend_h(k,cell1) - flux
end do
end if
if (cell2 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelCell(cell2)
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell2) = tend_h(k,cell2) + flux
end do
end if
end do
do iCell=1,nCells
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
end do
end do
@@ -361,7 +365,7 @@
! by just tend_h(2:nVertLevels,:) = 0.0, and then there is
! no need to compute the horizontal tend_h term for k=2:nVertLevels
! on a zlevel grid, above.
- do k=2,nVertLevels
+ do k=2,maxLevelCell(iCell)
tend_h(k,iCell) = tend_h(k,iCell) &
- (wTop(k,iCell) - wTop(k+1,iCell))
end do
@@ -381,7 +385,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=2,nVertLevels
+ do k=2,maxLevelEdge(iEdge)
! Average w from cell center to edge
wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
@@ -391,7 +395,7 @@
end do
! Average w*du/dz from vertical interface to vertical middle of cell
- do k=1,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
enddo
enddo
@@ -406,7 +410,7 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
end do
@@ -415,7 +419,7 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- rho0Inv*( pZLevel(k,cell2) &
- pZLevel(k,cell1) )/dcEdge(iEdge)
@@ -435,7 +439,7 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
@@ -460,29 +464,40 @@
!
do iEdge=1,grid % nEdgesSolve
- ! forcing in top layer only
- tend_u(1,iEdge) = tend_u(1,iEdge) &
+ ! forcing in top layer only
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+ k = maxLevelEdge(iEdge)
+
+ ! efficiency note: it would be nice to avoid this
+ ! if within a do. This could be done with
+ ! k = max(maxLevelEdge(iEdge),1)
+ ! and then tend_u(1,iEdge) is just not used for land cells.
+
+ if (k>0) then
+
! bottom drag is the same as POP:
! -c |u| u where c is unitless and 1.0e-3.
! see POP Reference guide, section 3.4.4.
- tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- - 1.0e-3*u(nVertLevels,iEdge) &
- *sqrt(2.0*ke_edge(nVertLevels,iEdge))
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - 1.0e-3*u(k,iEdge) &
+ *sqrt(2.0*ke_edge(k,iEdge))
! mrp 100603 The following method is more straight forward,
! that the above method of computing ke_edge, but I have
! not verified that v is working correctly yet.
- !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- ! - 1.0e-3*u(nVertLevels,iEdge) &
- ! *sqrt(u(nVertLevels,iEdge)**2 + v(nVertLevels,iEdge)**2)
+ !tend_u(k,iEdge) = tend_u(k,iEdge) &
+ ! - 1.0e-3*u(k,iEdge) &
+ ! *sqrt(u(k,iEdge)**2 + v(k,iEdge)**2)
! old bottom drag, just linear friction
! du/dt = u/tau where tau=100 days.
- !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- ! - u(nVertLevels,iEdge)/(100.0*86400.0)
+ !tend_u(k,iEdge) = tend_u(k,iEdge) &
+ ! - u(k,iEdge)/(100.0*86400.0)
+ endif
+
enddo
!
@@ -513,12 +528,12 @@
fluxVertTop(1) = 0.0
fluxVertTop(nVertLevels+1) = 0.0
do iEdge=1,grid % nEdgesSolve
- do k=2,nVertLevels
+ do k=2,maxLevelEdge(iEdge)
fluxVertTop(k) = vertViscTop(k) &
* ( u(k-1,iEdge) - u(k,iEdge) ) &
/ (zMidEdge(k-1,iEdge) - zMidEdge(k,iEdge))
enddo
- do k=1,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
+ (fluxVertTop(k) - fluxVertTop(k+1)) &
/(zTopEdge(k,iEdge) - zTopEdge(k+1,iEdge))
@@ -556,7 +571,8 @@
tracers, tend_tr
type (dm_info) :: dminfo
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdge, maxLevelVertex
integer, dimension(:,:), pointer :: cellsOnEdge
real (kind=RKIND), dimension(:), pointer :: &
zTopZLevel
@@ -580,6 +596,9 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdge => grid % maxLevelEdge % array
+ maxLevelVertex => grid % maxLevelVertex % array
nEdges = grid % nEdges
nCells = grid % nCells
@@ -594,8 +613,8 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
+ !if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdge(iEdge)
do iTracer=1,num_tracers
tracer_edge = 0.5 * ( tracers(iTracer,k,cell1) &
+ tracers(iTracer,k,cell2))
@@ -605,7 +624,7 @@
tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux
end do
end do
- end if
+ !end if
end do
elseif (config_hor_tracer_adv.eq.'upwind') then
@@ -613,8 +632,8 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
+ !if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdge(iEdge)
if (u(k,iEdge)>0.0) then
upwindCell = cell1
else
@@ -627,12 +646,12 @@
tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux
end do
end do
- end if
+ !end if
end do
endif
do iCell=1,grid % nCellsSolve
- do k=1,grid % nVertLevelsSolve
+ do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) / areaCell(iCell)
end do
@@ -648,7 +667,7 @@
do iCell=1,grid % nCellsSolve
if (config_vert_tracer_adv.eq.'centered') then
- do k=2,nVertLevels
+ do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
tracerTop(iTracer,k) = ( tracers(iTracer,k-1,iCell) &
+tracers(iTracer,k ,iCell))/2.0
@@ -656,7 +675,7 @@
end do
elseif (config_vert_tracer_adv.eq.'upwind') then
- do k=2,nVertLevels
+ do k=2,maxLevelCell(iCell)
if (wTop(k,iCell)>0.0) then
upwindCell = k
else
@@ -669,7 +688,7 @@
endif
- do k=1,nVertLevels
+ 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 ) &
@@ -690,14 +709,14 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
+ !if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdge(iEdge)
do iTracer=1,num_tracers
tr_flux(iTracer,k,iEdge) = h_edge(k,iEdge)*config_hor_diffusion * &
(Tracers(iTracer,k,cell2) - Tracers(iTracer,k,cell1))/dcEdge(iEdge)
enddo
enddo
- endif
+ !endif
enddo
! Compute the divergence, div(h \kappa_h </font>
<font color="gray">abla\phi) at cell centers
@@ -707,7 +726,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelCell(cell1)
do iTracer=1,num_tracers
tr_div(iTracer,k,cell1) = tr_div(iTracer,k,cell1) &
+ tr_flux(iTracer,k,iEdge)*dvEdge(iEdge)
@@ -715,7 +734,7 @@
enddo
endif
if (cell2 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelCell(cell2)
do iTracer=1,num_tracers
tr_div(iTracer,k,cell2) = tr_div(iTracer,k,cell2) &
- tr_flux(iTracer,k,iEdge)*dvEdge(iEdge)
@@ -727,7 +746,7 @@
! add div(h \kappa_h </font>
<font color="gray">abla\phi ) to tracer tendency
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
+ do k = 1,maxLevelCell(iCell)
do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ tr_div(iTracer,k,iCell)*r
@@ -764,7 +783,7 @@
fluxVertTop(:,1) = 0.0
fluxVertTop(:,nVertLevels+1) = 0.0
do iCell=1,grid % nCellsSolve
- do k=2,nVertLevels
+ do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
! compute \kappa_v d\phi/dz
fluxVertTop(iTracer,k) = vertDiffTop(k) &
@@ -773,7 +792,7 @@
enddo
enddo
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
dist = zTop(k,iCell) - zTop(k+1,iCell)
do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
@@ -834,7 +853,8 @@
character :: c1*6
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdge, maxLevelVertex
real (kind=RKIND) :: r, h1, h2
@@ -883,6 +903,9 @@
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
hZLevel => grid % hZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdge => grid % maxLevelEdge % array
+ maxLevelVertex => grid % maxLevelVertex % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -897,19 +920,18 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- elseif(cell1 <= nCells) then
- do k=1,nVertLevels
+ ! edge with cell on both sides
+ do k=1,maxLevelEdge(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
+
+ do k=maxLevelEdge(iEdge)+1,max(maxLevelCell(cell1),maxLevelCell(cell2))
+ if (k <= maxLevelCell(cell1)) then
h_edge(k,iEdge) = h(k,cell1)
- end do
- elseif(cell2 <= nCells) then
- do k=1,nVertLevels
+ elseif (k <= maxLevelCell(cell2)) then
h_edge(k,iEdge) = h(k,cell2)
- end do
- end if
+ end if
+ end do
end do
@@ -962,7 +984,7 @@
end do
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
+ do k = 1,maxLevelCell(iCell)
divergence(k,iCell) = divergence(k,iCell) * r
enddo
enddo
@@ -975,10 +997,10 @@
do i=1,nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
do k=1,nVertLevels
- ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
end do
end do
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
end do
end do
@@ -1022,7 +1044,12 @@
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! mrp 100609 was
+! if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
+! mrp 100609 should be
+ if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end do
do k=1,nVertLevels
h_vertex = 0.0
@@ -1082,7 +1109,7 @@
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
if (iCell <= nCells) then
- do k = 1,nVertLevels
+ do k = 1,maxLevelCell(iCell)
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
enddo
endif
@@ -1125,7 +1152,7 @@
! For a isopycnal model, density should remain constant.
if (config_vert_grid_type.eq.'zlevel') then
do iCell=1,nCells
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
! Linear equation of state, for the time being
rho(k,iCell) = 1000.0*( 1.0 &
- 2.5e-4*tracers(index_temperature,k,iCell) &
@@ -1193,7 +1220,7 @@
! assume atmospheric pressure at the surface is zero for now.
pZLevel(1,iCell) = rho(1,iCell)*gravity &
* (h(1,iCell)-0.5*hZLevel(1))
- do k=2,nVertLevels
+ do k=2,maxLevelCell(iCell)
delta_p = rho(k,iCell)*gravity*hZLevel(k)
pZLevel(k,iCell) = pZLevel(k-1,iCell) + 0.5*delta_p
end do
@@ -1231,14 +1258,14 @@
end do
do iCell=1,nCells
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
end do
! Vertical velocity at bottom is zero.
! this next line can be set permanently somewhere upon startup.
wTop(nVertLevels+1,iCell) = 0.0
- do k=nVertLevels,1,-1
+ do k=maxLevelCell(iCell),1,-1
wTop(k,iCell) = wTop(k+1,iCell) - div_u(k,iCell)*h(k,iCell)
end do
Modified: branches/ocean_projects/topography_mrp/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_sw/module_time_integration.F 2010-06-07 17:49:33 UTC (rev 339)
+++ branches/ocean_projects/topography_mrp/src/core_sw/module_time_integration.F 2010-06-09 17:33:50 UTC (rev 340)
@@ -618,7 +618,12 @@
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! mrp 100609 was
+! if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
+! mrp 100609 should be
+ if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end do
do k=1,nVertLevels
h_vertex = 0.0
</font>
</pre>