<p><b>mpetersen@lanl.gov</b> 2010-11-15 12:50:10 -0700 (Mon, 15 Nov 2010)</p><p>BRANCH COMMIT Updating topography branch. This is the first part of the revisions after Todd's branch review. This has a bit-for-bit match with the previous commit for z-level runs.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-11-12 21:02:17 UTC (rev 613)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-11-15 19:50:10 UTC (rev 614)
@@ -168,7 +168,8 @@
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 o uReconstructZ state - -
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
-var persistent real p ( nVertLevels nCells Time ) 2 o p state - -
+var persistent real MontPot ( nVertLevels nCells Time ) 2 o MontPot state - -
+var persistent real pressure ( nVertLevels nCells Time ) 2 o pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 o wTop state - -
# Other diagnostic variables: neither read nor written to any files
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_global_diagnostics.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_global_diagnostics.F 2010-11-12 21:02:17 UTC (rev 613)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_global_diagnostics.F 2010-11-15 19:50:10 UTC (rev 614)
@@ -37,7 +37,7 @@
real (kind=RKIND) :: areaCellGlobal, areaEdgeGlobal, areaTriangleGlobal
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, areaEdge
real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, circulation, vorticity, ke, pv_edge, pv_vertex, &
- pv_cell, gradPVn, gradPVt, p, wTop, rho, tracerTemp
+ pv_cell, gradPVn, gradPVt, pressure, MontPot, wTop, rho, tracerTemp
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
@@ -80,7 +80,8 @@
pv_cell => state % pv_cell % array
gradPVn => state % gradPVn % array
gradPVt => state % gradPVt % array
- p => state % p % array
+ MontPot => state % MontPot % array
+ pressure => state % pressure % array
variableIndex = 0
! h
@@ -152,12 +153,18 @@
gradPVt(:,1:nEdgesSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! p
+ ! pressure
variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- p(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ pressure(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
+ ! MontPot
+ variableIndex = variableIndex + 1
+ call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
+ MontPot(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ verticalSumMaxes(variableIndex))
+
! wTop vertical velocity
variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels+1, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
@@ -281,10 +288,14 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeEdgeGlobal
- ! p
+ ! pressure
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
+ ! MontPot
+ variableIndex = variableIndex + 1
+ averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
+
! wTop vertical velocity
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
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-11-12 21:02:17 UTC (rev 613)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-11-15 19:50:10 UTC (rev 614)
@@ -268,16 +268,15 @@
vertex1, vertex2, eoe, i, j, kmax
integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: h_mom_eddy_visc2, h_mom_eddy_visc4
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
upstream_bias, wTopEdge, rho0Inv, r
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, p, wTop, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
tend_h, tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- divergence
+ MontPot, wTop, divergence
type (dm_info) :: dminfo
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
@@ -307,7 +306,8 @@
ke => s % ke % array
ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
- p => s % p % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -343,9 +343,6 @@
u_src => grid % u_src % array
- h_mom_eddy_visc2 = config_h_mom_eddy_visc2
- h_mom_eddy_visc4 = config_h_mom_eddy_visc4
-
!
! height tendency: start accumulating tendency terms
!
@@ -434,7 +431,7 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - (p(k,cell2) - p(k,cell1))/dcEdge(iEdge)
+ - (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
end do
enddo
elseif (config_vert_grid_type.eq.'zlevel') then
@@ -443,8 +440,8 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - rho0Inv*( p(k,cell2) &
- - p(k,cell1) )/dcEdge(iEdge)
+ - rho0Inv*( pressure(k,cell2) &
+ - pressure(k,cell1) )/dcEdge(iEdge)
end do
enddo
endif
@@ -452,9 +449,9 @@
!
! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="black">abla^2 u
! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity )
- ! strictly only valid for h_mom_eddy_visc2 == constant
+ ! strictly only valid for config_h_mom_eddy_visc2 == constant
!
- if ( h_mom_eddy_visc2 > 0.0 ) then
+ if ( config_h_mom_eddy_visc2 > 0.0 ) then
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -469,7 +466,7 @@
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = h_mom_eddy_visc2 * u_diffusion
+ u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
@@ -481,9 +478,9 @@
! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="black">abla^4 u
! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity
! applied recursively.
- ! strictly only valid for h_mom_eddy_visc4 == constant
+ ! strictly only valid for config_h_mom_eddy_visc4 == constant
!
- if ( h_mom_eddy_visc4 > 0.0 ) then
+ if ( config_h_mom_eddy_visc4 > 0.0 ) then
allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_u(nVertLevels, nEdges+1))
@@ -562,7 +559,7 @@
-( delsq_vorticity(k,vertex2) &
- delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - config_h_mom_eddy_visc4 * u_diffusion
end do
end do
@@ -662,8 +659,6 @@
fluxVertTop(k) = vertViscTop(k) &
* ( u(k-1,iEdge) - u(k,iEdge) ) &
* 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
-! mrp previous line was this, but z variables are redundant with h.
-! / (zMidEdge(k-1,iEdge) - zMidEdge(k,iEdge))
enddo
fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
@@ -671,8 +666,6 @@
tend_u(k,iEdge) = tend_u(k,iEdge) &
+ (fluxVertTop(k) - fluxVertTop(k+1)) &
/ h_edge(k,iEdge)
-! mrp previous line was this, but z variables are redundant with h.
-! /(zTopEdge(k,iEdge) - zTopEdge(k+1,iEdge))
enddo
end do
@@ -882,7 +875,6 @@
!
allocate(tracerTop(num_tracers,nVertLevels+1))
tracerTop(:,1)=0.0
- tracerTop(:,nVertLevels+1)=0.0
do iCell=1,grid % nCellsSolve
if (config_vert_tracer_adv.eq.'centered') then
@@ -1112,7 +1104,7 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, pTop, rho0Inv
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
@@ -1121,11 +1113,12 @@
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
hZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, p, wTop, &
- circulation, vorticity, ke, ke_edge, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
+ circulation, vorticity, ke, ke_edge, MontPot, wTop, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
rho, temperature, salinity
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
@@ -1133,7 +1126,8 @@
verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, maxLevelVertexBot
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot, maxLevelVertexTop
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
@@ -1157,7 +1151,8 @@
gradPVt => s % gradPVt % array
rho => s % rho % array
tracers => s % tracers % array
- p => s % p % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -1181,7 +1176,8 @@
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+ maxLevelVertexTop => grid % maxLevelVertexTop % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -1197,7 +1193,6 @@
! Namelist options control the order of accuracy of the reconstructed h_edge value
!
- ! mrp 101026 efficiency note: zlevel does not need to compute h_edge for k>1
coef_3rd_order = 0.
if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
@@ -1369,7 +1364,7 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- do k = 1,maxLevelEdgeBot(iEdge)
+ do k = 1,maxLevelEdgeTop(iEdge)
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
end do
@@ -1491,25 +1486,27 @@
! For Isopycnal model.
! Compute pressure at top of each layer, and then
! Montgomery Potential.
+ allocate(pTop(nVertLevels))
do iCell=1,nCells
! assume atmospheric pressure at the surface is zero for now.
- pTop = 0.0
+ pTop(1) = 0.0
! For isopycnal mode, p is the Montgomery Potential.
! At top layer it is g*SSH, where SSH may be off by a
! constant (ie, h_s can be relative to top or bottom)
- p(1,iCell) = gravity &
+ MontPot(1,iCell) = gravity &
* (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
do k=2,nVertLevels
- pTop = pTop + rho(k-1,iCell)*gravity* h(k-1,iCell)
+ pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
! from delta M = p delta / rho
- p(k,iCell) = p(k-1,iCell) &
- + pTop*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
+ MontPot(k,iCell) = MontPot(k-1,iCell) &
+ + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
end do
end do
+ deallocate(pTop)
elseif (config_vert_grid_type.eq.'zlevel') then
@@ -1522,11 +1519,11 @@
! compute pressure for z-level coordinates
! assume atmospheric pressure at the surface is zero for now.
- p(1,iCell) = rho(1,iCell)*gravity &
+ pressure(1,iCell) = rho(1,iCell)*gravity &
* (h(1,iCell)-0.5*hZLevel(1))
do k=2,maxLevelCell(iCell)
- p(k,iCell) = p(k-1,iCell) &
+ pressure(k,iCell) = pressure(k-1,iCell) &
+ 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
+ rho(k ,iCell)*hZLevel(k ))
end do
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-11-12 21:02:17 UTC (rev 613)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-11-15 19:50:10 UTC (rev 614)
@@ -16,43 +16,50 @@
! Initialize grid variables that are computed from the netcdf input file.
use grid_types
+ use dmpar
use configure
- use constants
- use dmpar
implicit none
type (domain_type), intent(inout) :: domain
-
- 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, 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'
+ call init_ZLevel(domain)
else
print *, ' Incorrect choice of config_vert_grid_type:',&
config_vert_grid_type
call dmpar_abort(dminfo)
endif
+
+ call compute_maxLevel(domain)
+
+end subroutine mpas_init_domain
+
+
+subroutine init_ZLevel(domain)
+! Initialize maxLevel and bouncary grid variables.
+
+ use grid_types
+ use configure
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ integer :: iTracer, cell
+ real (kind=RKIND), dimension(:), pointer :: &
+ hZLevel, zMidZLevel, zTopZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: h
+ integer :: nVertLevels
+
! Initialize z-level grid variables from h, read in from input file.
block => domain % blocklist
do while (associated(block))
@@ -61,6 +68,60 @@
hZLevel => block % mesh % hZLevel % array
zMidZLevel => block % mesh % zMidZLevel % array
zTopZLevel => block % mesh % zTopZLevel % array
+ nVertLevels = block % mesh % nVertLevels
+
+ ! These should eventually be in an input file. For now
+ ! I just read them in from h(:,1).
+ ! Upon restart, the correct hZLevel should be in restart.nc
+ if (.not. config_do_restart) hZLevel = h(:,1)
+
+ ! 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 k = 1,nVertLevels
+ zMidZLevel(k) = zTopZLevel(k)-0.5*hZLevel(k)
+ zTopZLevel(k+1) = zTopZLevel(k)- hZLevel(k)
+ end do
+
+ block => block % next
+
+ end do
+
+end subroutine init_ZLevel
+
+
+subroutine compute_maxLevel(domain)
+! Initialize maxLevel and bouncary grid variables.
+
+ use grid_types
+ use configure
+ use constants
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ 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, &
+ boundaryVertex, verticesOnEdge
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
maxLevelCell => block % mesh % maxLevelCell % array
maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
@@ -68,8 +129,10 @@
maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
cellsOnEdge => block % mesh % cellsOnEdge % array
cellsOnVertex => block % mesh % cellsOnVertex % array
+ verticesOnEdge => block % mesh % verticesOnEdge % array
boundaryEdge => block % mesh % boundaryEdge % array
boundaryCell => block % mesh % boundaryCell % array
+ boundaryVertex => block % mesh % boundaryVertex % array
nCells = block % mesh % nCells
nEdges = block % mesh % nEdges
@@ -77,24 +140,6 @@
nVertLevels = block % mesh % nVertLevels
vertexDegree = block % mesh % vertexDegree
- if (config_vert_grid_type.eq.'zlevel') then
-
- ! These should eventually be in an input file. For now
- ! I just read them in from h(:,1).
- ! Upon restart, the correct hZLevel should be in restart.nc
- if (.not. config_do_restart) hZLevel = h(:,1)
-
- ! 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
-
- endif
-
! for z-grids, maxLevelCell should be in input state
! Isopycnal grid uses all vertical cells
if (config_vert_grid_type.eq.'isopycnal') then
@@ -104,34 +149,42 @@
! 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
-
+ maxLevelEdgeTop(iEdge) = &
+ min( maxLevelCell(cellsOnEdge(1,iEdge)), &
+ maxLevelCell(cellsOnEdge(2,iEdge)) )
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(iEdge) = &
- max( maxLevelCell(cell1), &
- maxLevelCell(cell2) )
- else
- maxLevelEdgeBot(iEdge) = 0
- endif
+ maxLevelEdgeBot(iEdge) = &
+ max( maxLevelCell(cellsOnEdge(1,iEdge)), &
+ maxLevelCell(cellsOnEdge(2,iEdge)) )
end do
maxLevelEdgeBot(nEdges+1) = 0
+ ! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexBot(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
+ do i=2,vertexDegree
+ maxLevelVertexBot(iVertex) = &
+ max( maxLevelVertexBot(iVertex), &
+ maxLevelCell(cellsOnVertex(i,iVertex)))
+ end do
+ end do
+ maxLevelVertexBot(nVertices+1) = 0
+
+ ! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexTop(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
+ do i=2,vertexDegree
+ maxLevelVertexTop(iVertex) = &
+ min( maxLevelVertexTop(iVertex), &
+ maxLevelCell(cellsOnVertex(i,iVertex)))
+ end do
+ end do
+ maxLevelVertexTop(nVertices+1) = 0
+
! set boundary edge
boundaryEdge=1
do iEdge=1,nEdges
@@ -139,55 +192,27 @@
end do
!
- ! Find those cells that have an edge on the boundary
+ ! Find cells and vertices 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) )
+ do k=1,nVertLevels
+ if (boundaryEdge(k,iEdge).eq.1) then
+ boundaryCell(k,cellsOnEdge(1,iEdge)) = 1
+ boundaryCell(k,cellsOnEdge(2,iEdge)) = 1
+ boundaryVertex(k,verticesOnEdge(1,iEdge)) = 1
+ boundaryVertex(k,verticesOnEdge(2,iEdge)) = 1
endif
end do
- end do
- maxLevelVertexBot(nVertices+1) = 0
+ end do
- ! 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
-
block => block % next
end do
! Note: We do not update halos on maxLevel* variables. I want the
! outside edge of a halo to be zero on each processor.
-end subroutine mpas_init_domain
+end subroutine compute_maxLevel
subroutine mpas_init(block, mesh, dt)
</font>
</pre>