<p><b>mpetersen@lanl.gov</b> 2010-11-02 06:52:48 -0600 (Tue, 02 Nov 2010)</p><p>BRANCH COMMIT: Removed more if statements from time_integration, adjusted spacing.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-10-29 20:18:00 UTC (rev 588)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-11-02 12:52:48 UTC (rev 589)
@@ -305,7 +305,7 @@
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
- ke_edge => s % ke_edge % array
+ ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
MontPot => s % MontPot % array
pZLevel => s % pZLevel % array
@@ -368,11 +368,6 @@
kmax = 1
endif
- ! mrp efficiency note: I try to remove if statements from within
- ! loops, but did not here. I could use a loop like
- ! do k=1,min(1,maxLevelEdgeTop(iedge))
- ! for z-level, but then cells in the outer halo ring do not get
- ! the proper flux, as they do in the current code.
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -382,31 +377,6 @@
tend_h(k,cell2) = tend_h(k,cell2) + flux
end do
end do
-!mrp old begin
-if (1==2) then
- ! mrp efficiency note: I try to remove if statements from within
- ! loops, but did not here. I could use a loop like
- ! do k=1,min(1,maxLevelEdgeTop(iedge))
- ! for z-level, but then cells in the outer halo ring do not get
- ! the proper flux, as they do in the current code.
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,kmax
- 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,kmax
- 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
-endif
-!mrp old end
do iCell=1,nCells
do k=1,kmax
@@ -414,9 +384,6 @@
end do
end do
- !print *, 'nEdges,nEdgesSolve,nCells,size(cellsOnEdge),maxval(cellsOnEdge(1,1:nEdges)),maxval(cellsOnEdge(1,1:nEdgesSolve))', &
- !nEdges,nEdgesSolve,nCells,size(cellsOnEdge),maxval(cellsOnEdge(1,1:nEdges)),maxval(cellsOnEdge(1,1:nEdgesSolve))
-
!
! height tendency: vertical advection term -d/dz(hw)
!
@@ -1346,10 +1313,13 @@
endif ! if(config_thickness_adv_order == 2)
!
- ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
- ! used to when reading for edges that do not exist
+ ! set the velocity and height at dummy address
+ ! used -1e34 so error clearly occurs if these values are used.
!
- u(:,nEdges+1) = 0.0
+ u(:,nEdges+1) = -1e34
+ h(:,nCells+1) = -1e34
+ tracers(s % index_temperature,:,nCells+1) = -1e34
+ tracers(s % index_salinity,:,nCells+1) = -1e34
!
! Compute circulation and relative vorticity at each vertex
@@ -1467,18 +1437,6 @@
enddo
!
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
- !
- do iEdge = 1,nEdges
- do k = 1,maxLevelEdgeBot(iEdge)
- gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
- - pv_vertex(k,verticesOnEdge(1,iEdge))) &
- /dvEdge(iEdge)
- enddo
- enddo
-
- !
! Compute pv at the edges
! ( this computes pv_edge at all edges bounding real cells )
!
@@ -1493,21 +1451,12 @@
end do
!
- ! Modify PV edge with upstream bias.
- !
- do iEdge = 1,nEdges
- do k = 1,maxLevelEdgeBot(iEdge)
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
- enddo
- enddo
-
- !
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
+ do k=1,maxLevelEdgeTop(iEdge)
gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
- pv_cell(k,cellsOnEdge(1,iEdge))) &
/ dcEdge(iEdge)
@@ -1515,11 +1464,25 @@
enddo
!
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
+ !
+ do iEdge = 1,nEdges
+ do k = 1,maxLevelEdgeBot(iEdge)
+ gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
+ - pv_vertex(k,verticesOnEdge(1,iEdge))) &
+ /dvEdge(iEdge)
+ enddo
+ enddo
+
+ !
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
- do k=1,maxLevelEdgeTop(iEdge)
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ do k = 1,maxLevelEdgeBot(iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) &
+ - 0.5 * dt* ( u(k,iEdge) * gradPVn(k,iEdge) &
+ + v(k,iEdge) * gradPVt(k,iEdge) )
enddo
enddo
@@ -1556,61 +1519,58 @@
!
if (config_vert_grid_type.eq.'isopycnal') then
- ! Compute mid- and top-depths of each layer, from bottom up.
- do iCell=1,nCells
- ! h_s is ocean topography: elevation above lowest point,
- ! and is positive. z coordinates are pos upward.
- ! h_s is only used for isopycnal coordinates.
- zTop(nVertLevels+1,iCell) = h_s(iCell)
- do k=nVertLevels,1,-1
- zMid(k,iCell) = zTop(k+1,iCell) + 0.5*h(k,iCell)
- zTop(k,iCell) = zTop(k+1,iCell) + h(k,iCell)
- end do
- end do
+ ! Compute mid- and top-depths of each layer, from bottom up.
+ do iCell=1,nCells
+ ! h_s is ocean topography: elevation above lowest point,
+ ! and is positive. z coordinates are pos upward.
+ ! h_s is only used for isopycnal coordinates.
+ zTop(nVertLevels+1,iCell) = h_s(iCell)
+ do k=nVertLevels,1,-1
+ zMid(k,iCell) = zTop(k+1,iCell) + 0.5*h(k,iCell)
+ zTop(k,iCell) = zTop(k+1,iCell) + h(k,iCell)
+ end do
+ end do
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if(cell1<=nCells .and. cell2<=nCells) then
- do k=1,nVertLevels
- zTopEdge(k,iEdge) = (zTop(k,cell1)+zTop(k,cell2))/2.0
- zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ zTopEdge(k,iEdge) = (zTop(k,cell1)+zTop(k,cell2))/2.0
+ zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
enddo
zTopEdge(nVertLevels+1,iEdge) = ( zTop(nVertLevels+1,cell1) &
+ zTop(nVertLevels+1,cell2))/2.0
- endif
- enddo
+ enddo
elseif (config_vert_grid_type.eq.'zlevel') then
- ! For z-level coordinates, z variables only change for the top
- ! layer, so z variables for k=2:nVertLevels+1 can be computed
- ! from 1D ZLevel variables
-! mrp efficiency. Move zmid ztop ztopedge zmidedge for k>1 to mpas_init, and remove from test_cases
- do iCell=1,nCells
- do k=2,nVertLevels
- zMid(k,iCell) = zMidZLevel(k)
- zTop(k,iCell) = zTopZLevel(k)
- end do
- zTop(nVertLevels+1,iCell) = zTopZLevel(nVertLevels+1)
+ ! For z-level coordinates, z variables only change for the top
+ ! layer, so z variables for k=2:nVertLevels+1 can be computed
+ ! from 1D ZLevel variables
+ do iCell=1,nCells
+ do k=2,nVertLevels
+ zMid(k,iCell) = zMidZLevel(k)
+ zTop(k,iCell) = zTopZLevel(k)
+ end do
+ zTop(nVertLevels+1,iCell) = zTopZLevel(nVertLevels+1)
- zMid(1,iCell) = zTopZLevel(2) + 0.5*h(1,iCell)
- zTop(1,iCell) = zTopZLevel(2) + h(1,iCell)
- enddo
+ zMid(1,iCell) = zTopZLevel(2) + 0.5*h(1,iCell)
+ zTop(1,iCell) = zTopZLevel(2) + h(1,iCell)
+ enddo
+ zMid(1,nCells+1)=0
+ zTop(1,nCells+1)=0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1<=nCells .and. cell2<=nCells) then
- zTopEdge(1,iEdge) = (zTop(1,cell1)+zTop(1,cell2))/2.0
- zMidEdge(1,iEdge) = (zMid(1,cell1)+zMid(1,cell2))/2.0
- endif
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ zTopEdge(1,iEdge) = (zTop(1,cell1)+zTop(1,cell2))/2.0
+ zMidEdge(1,iEdge) = (zMid(1,cell1)+zMid(1,cell2))/2.0
- do k=2,nVertLevels
- zMidEdge(k,iEdge) = zMidZLevel(k)
- zTopEdge(k,iEdge) = zTopZLevel(k)
- end do
- zTopEdge(nVertLevels+1,iEdge) = zTopZLevel(nVertLevels+1)
+ do k=2,nVertLevels
+ zMidEdge(k,iEdge) = zMidZLevel(k)
+ zTopEdge(k,iEdge) = zTopZLevel(k)
+ end do
+ zTopEdge(nVertLevels+1,iEdge) = zTopZLevel(nVertLevels+1)
enddo
endif
@@ -1680,26 +1640,16 @@
! Compute div(u) for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
!
- allocate(div_u(nVertLevels,nCells))
+ allocate(div_u(nVertLevels,nCells+1))
div_u(:,:) = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
-! mrp efficiency note: I can get rid of cell1<=nCells only if these edges
-! are not on the outside edge of a halo. I dont think they are, but
-! double check before you delete it.
- flux = u(k,iEdge) * dvEdge(iEdge)
- div_u(k,cell1) = div_u(k,cell1) + flux
- end do
- end if
- if (cell2 <= nCells) then
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
- flux = u(k,iEdge) * dvEdge(iEdge)
- div_u(k,cell2) = div_u(k,cell2) - flux
- end do
- end if
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ flux = u(k,iEdge) * dvEdge(iEdge)
+ div_u(k,cell1) = div_u(k,cell1) + flux
+ div_u(k,cell2) = div_u(k,cell2) - flux
+ end do
end do
do iCell=1,nCells
@@ -1750,6 +1700,7 @@
end do
end do
+
end subroutine compute_solve_diagnostics
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-10-29 20:18:00 UTC (rev 588)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-11-02 12:52:48 UTC (rev 589)
@@ -24,166 +24,167 @@
type (domain_type), intent(inout) :: domain
- integer :: i, iCell, iEdge, iVertex, iLevel
- type (block_type), pointer :: block
- type (dm_info) :: dminfo
+ integer :: i, iCell, iEdge, iVertex, iLevel
+ type (block_type), pointer :: block
+ type (dm_info) :: dminfo
- integer :: iTracer, cell1, cell2, cell, k
- real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
- hZLevel, zMidZLevel, zTopZLevel, latCell,LonCell
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
- real (kind=RKIND) :: centerx, centery
- integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+ integer :: iTracer, cell1, cell2, cell, k
+ real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
+ hZLevel, zMidZLevel, zTopZLevel, latCell,LonCell
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
+ real (kind=RKIND) :: centerx, centery
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
- integer, dimension(:), pointer :: &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
- maxLevelVertexTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, boundaryEdge, boundaryCell
+ integer, dimension(:), pointer :: &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, boundaryEdge, boundaryCell
- if (config_vert_grid_type.eq.'isopycnal') then
- print *, ' Using isopycnal coordinates'
- elseif (config_vert_grid_type.eq.'zlevel') then
- print *, ' Using z-level coordinates'
- else
- print *, ' Incorrect choice of config_vert_grid_type:',&
- config_vert_grid_type
- call dmpar_abort(dminfo)
- endif
+ if (config_vert_grid_type.eq.'isopycnal') then
+ print *, ' Using isopycnal coordinates'
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ print *, ' Using z-level coordinates'
+ else
+ print *, ' Incorrect choice of config_vert_grid_type:',&
+ config_vert_grid_type
+ call dmpar_abort(dminfo)
+ endif
- ! Initialize z-level grid variables from h, read in from input file.
- block => domain % blocklist
- do while (associated(block))
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
- hZLevel => block % mesh % hZLevel % array
- zMidZLevel => block % mesh % zMidZLevel % array
- zTopZLevel => block % mesh % zTopZLevel % array
- maxLevelCell => block % mesh % maxLevelCell % array
- maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
- maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
- maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
- maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
- cellsOnEdge => block % mesh % cellsOnEdge % array
- cellsOnVertex => block % mesh % cellsOnVertex % array
- boundaryEdge => block % mesh % boundaryEdge % array
- boundaryCell => block % mesh % boundaryCell % array
+ hZLevel => block % mesh % hZLevel % array
+ zMidZLevel => block % mesh % zMidZLevel % array
+ zTopZLevel => block % mesh % zTopZLevel % array
+ maxLevelCell => block % mesh % maxLevelCell % array
+ maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+ maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
+ maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
+ maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
+ cellsOnEdge => block % mesh % cellsOnEdge % array
+ cellsOnVertex => block % mesh % cellsOnVertex % array
+ boundaryEdge => block % mesh % boundaryEdge % array
+ boundaryCell => block % mesh % boundaryCell % array
- nCells = block % mesh % nCells
- nEdges = block % mesh % nEdges
- nVertices = block % mesh % nVertices
- nVertLevels = block % mesh % nVertLevels
- vertexDegree = block % mesh % vertexDegree
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertices = block % mesh % nVertices
+ nVertLevels = block % mesh % nVertLevels
+ vertexDegree = block % mesh % vertexDegree
- if (config_vert_grid_type.eq.'zlevel') then
+ if (config_vert_grid_type.eq.'zlevel') then
- ! hZLevel should be in the grid.nc and restart.nc file
+ ! hZLevel should be in the grid.nc and restart.nc file,
+ ! and h for k=1 must be specified there as well.
- zTopZLevel(1) = 0.0
- do iLevel = 1,nVertLevels
- zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
- zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
- end do
+ zTopZLevel(1) = 0.0
+ do iLevel = 1,nVertLevels
+ zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
+ zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
+ end do
- endif
+ endif
- ! for z-grids, maxLevelCell should be in input state
- ! Isopycnal grid uses all vertical cells
- if (config_vert_grid_type.eq.'isopycnal') then
- maxLevelCell(1:nCells) = nVertLevels
- endif
- maxLevelCell(nCells+1) = 0
+ ! for z-grids, maxLevelCell should be in input state
+ ! Isopycnal grid uses all vertical cells
+ if (config_vert_grid_type.eq.'isopycnal') then
+ maxLevelCell(1:nCells) = nVertLevels
+ endif
+ maxLevelCell(nCells+1) = 0
- ! mrp insert topography mesa for testing
- ! do iCell = 1,nCells
- ! if (latCell(iCell)> 0.0/180.*3.14.and.&
- ! latCell(iCell)< 30.0/180.*3.14.and. &
- ! lonCell(iCell)>180.0/180.*3.14.and. &
- ! lonCell(iCell)<210.0/180.*3.14) then
- ! maxLevelCell(iCell) = 10
- ! endif
- ! end do
+ ! mrp insert topography mesa for testing
+ ! do iCell = 1,nCells
+ ! if (latCell(iCell)> 0.0/180.*3.14.and.&
+ ! latCell(iCell)< 30.0/180.*3.14.and. &
+ ! lonCell(iCell)>180.0/180.*3.14.and. &
+ ! lonCell(iCell)<210.0/180.*3.14) then
+ ! maxLevelCell(iCell) = 10
+ ! endif
+ ! end do
- ! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ ! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- maxLevelEdgeTop(iEdge) = &
- min( maxLevelCell(cell1), &
- maxLevelCell(cell2) )
- else
- maxLevelEdgeTop(iEdge) = 0
- endif
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ maxLevelEdgeTop(iEdge) = &
+ min( maxLevelCell(cell1), &
+ maxLevelCell(cell2) )
+ else
+ maxLevelEdgeTop(iEdge) = 0
+ endif
- end do
- maxLevelEdgeTop(nEdges+1) = 0
+ end do
+ maxLevelEdgeTop(nEdges+1) = 0
- ! maxLevelEdgeBot is the maximum (deepest) of the surrounding cells
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
+ ! maxLevelEdgeBot is the maximum (deepest) of the surrounding cells
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= nCells .and. cell2 <= nCells) then
maxLevelEdgeBot(iEdge) = &
- max( maxLevelCell(cell1), &
- maxLevelCell(cell2) )
- else
+ max( maxLevelCell(cell1), &
+ maxLevelCell(cell2) )
+ else
maxLevelEdgeBot(iEdge) = 0
- endif
- end do
- maxLevelEdgeBot(nEdges+1) = 0
+ endif
+ end do
+ maxLevelEdgeBot(nEdges+1) = 0
- ! set boundary edge
- boundaryEdge=1
- do iEdge=1,nEdges
- boundaryEdge(1:maxLevelEdgeTop(iEdge),iEdge)=0
- end do
+ ! set boundary edge
+ boundaryEdge=1
+ do iEdge=1,nEdges
+ boundaryEdge(1:maxLevelEdgeTop(iEdge),iEdge)=0
+ end do
- !
- ! Find those cells that have an edge on the boundary
- !
- boundaryCell(:,:) = 0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if (boundaryEdge(k,iEdge).eq.1) then
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- boundaryCell(k,cell1) = 1
- boundaryCell(k,cell2) = 1
- endif
- end do
+ !
+ ! Find those cells that have an edge on the boundary
+ !
+ boundaryCell(:,:) = 0
+ do iEdge=1,nEdges
+ do k=1,nVertLevels
+ if (boundaryEdge(k,iEdge).eq.1) then
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ boundaryCell(k,cell1) = 1
+ boundaryCell(k,cell2) = 1
+ endif
end do
+ end do
- ! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
- do iVertex = 1,nVertices
- maxLevelVertexBot(iVertex) = 0
- do i=1,vertexDegree
- cell = cellsOnVertex(i,iVertex)
- if (cell <= nCells) then
- maxLevelVertexBot(iVertex) = &
- max( maxLevelVertexBot(iVertex), &
- maxLevelCell(cell) )
- endif
- end do
- end do
- maxLevelVertexBot(nVertices+1) = 0
+ ! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexBot(iVertex) = 0
+ do i=1,vertexDegree
+ cell = cellsOnVertex(i,iVertex)
+ if (cell <= nCells) then
+ maxLevelVertexBot(iVertex) = &
+ max( maxLevelVertexBot(iVertex), &
+ maxLevelCell(cell) )
+ endif
+ end do
+ end do
+ maxLevelVertexBot(nVertices+1) = 0
- ! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
- do iVertex = 1,nVertices
- maxLevelVertexTop(iVertex) = 0
- do i=1,vertexDegree
- cell = cellsOnVertex(i,iVertex)
- if (cell <= nCells) then
- maxLevelVertexTop(iVertex) = &
- min( maxLevelVertexTop(iVertex), &
- maxLevelCell(cell) )
- endif
- end do
- end do
- maxLevelVertexTop(nVertices+1) = 0
+ ! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexTop(iVertex) = 0
+ do i=1,vertexDegree
+ cell = cellsOnVertex(i,iVertex)
+ if (cell <= nCells) then
+ maxLevelVertexTop(iVertex) = &
+ min( maxLevelVertexTop(iVertex), &
+ maxLevelCell(cell) )
+ endif
+ end do
+ end do
+ maxLevelVertexTop(nVertices+1) = 0
! mrp temp printing, can remove later
print *, 'min/max MaxLevelCell, min/max MaxLevelEdge, ',&
</font>
</pre>