<p><b>mpetersen@lanl.gov</b> 2010-05-06 15:22:50 -0600 (Thu, 06 May 2010)</p><p>Merging lateral_boundary_conditions branch onto trunk.<br>
<br>
I tested the compile of all three cores, and validated bit-for-bit<br>
matching on the ocean core before and after the trunk-> branch and <br>
branch->trunk merges. Todd validated bit-for-bit runs on the<br>
sw and hyd_atmos cores on the trunk->branch merge.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/namelist.input.ocean 2010-05-06 21:22:50 UTC (rev 258)
@@ -1,21 +1,21 @@
&sw_model
- config_test_case = 5
+ config_test_case = 0
config_time_integration = 'RK4'
- config_dt = 300.0
- config_ntimesteps = 3000
- config_output_interval = 300
- config_stats_interval = 10
- config_visc = 4.0
+ config_dt = 120.0
+ config_ntimesteps = 30
+ config_output_interval = 30
+ config_stats_interval = 10
+ config_visc = 100.0
/
&io
- config_input_name = grid.nc
- config_output_name = output.nc
- config_restart_name = restart.nc
+ config_input_name = 'grid.quad.20km.nc'
+ config_output_name = 'output.quad.20km.nc'
+ config_restart_name = 'restart.nc'
/
&restart
- config_restart_interval = 864000
+ config_restart_interval = 3000
config_do_restart = .false.
config_restart_time = 1036800.0
/
Modified: trunk/mpas/src/core_hyd_atmos/module_advection.F
===================================================================
--- trunk/mpas/src/core_hyd_atmos/module_advection.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/core_hyd_atmos/module_advection.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -96,7 +96,7 @@
do_the_cell = .true.
do i=1,n
- if (cell_list(i) <= 0) do_the_cell = .false.
+ if (cell_list(i) > grid % nCells) do_the_cell = .false.
end do
Modified: trunk/mpas/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_hyd_atmos/module_time_integration.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/core_hyd_atmos/module_time_integration.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -671,10 +671,10 @@
if ( h_mom_eddy_visc4 > 0.0 ) then
- allocate(delsq_divergence(nVertLevels, nCells))
- allocate(delsq_u(nVertLevels, nEdges))
- allocate(delsq_circulation(nVertLevels, nVertices))
- allocate(delsq_vorticity(nVertLevels, nVertices))
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
delsq_u(:,:) = 0.0
@@ -849,7 +849,7 @@
if ( h_theta_eddy_visc4 > 0.0 ) then
- allocate(delsq_theta(nVertLevels, nCells))
+ allocate(delsq_theta(nVertLevels, nCells+1))
delsq_theta(:,:) = 0.
@@ -1519,9 +1519,9 @@
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension( num_scalars, grid % nEdges) :: h_flux
- real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: v_flux, v_flux_upwind, s_update
- real (kind=RKIND), dimension( num_scalars, grid % nCells, 2 ) :: scale_out, scale_in
+ real (kind=RKIND), dimension( num_scalars, grid % nEdges+1) :: h_flux
+ real (kind=RKIND), dimension( num_scalars, grid % nCells+1, 2 ) :: v_flux, v_flux_upwind, s_update
+ real (kind=RKIND), dimension( num_scalars, grid % nCells+1, 2 ) :: scale_out, scale_in
real (kind=RKIND), dimension( num_scalars ) :: s_max, s_min, s_max_update, s_min_update
integer :: nVertLevels, km0, km1, ktmp, kcp1, kcm1
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/core_ocean/Registry 2010-05-06 21:22:50 UTC (rev 258)
@@ -6,9 +6,7 @@
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
-# mrp 100120:
namelist integer sw_model config_stats_interval 100
-# mrp 100120 end
namelist real sw_model config_visc 0.0
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -85,7 +83,8 @@
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
# Boundary conditions: read from input, saved in restart and written to output
-var integer uBC ( nVertLevels nEdges ) iro uBC - -
+var integer boundaryEdge ( nVertLevels nEdges ) iro boundaryEdge - -
+var integer boundaryVertex ( nVertLevels nVertices ) iro boundaryVertex - -
var real u_src ( nVertLevels nEdges ) iro u_src - -
# Prognostic variables: read from input, saved in restart, and written to output
@@ -105,14 +104,15 @@
var real uReconstructX ( nVertLevels nCells Time ) o uReconstructX - -
var real uReconstructY ( nVertLevels nCells Time ) o uReconstructY - -
var real uReconstructZ ( nVertLevels nCells Time ) o uReconstructZ - -
-# mrp 100112:
-var real zmid ( nVertLevels nCells Time ) o zmid - -
-var real zbot ( nVertLevels nCells Time ) o zbot - -
+var real zMid ( nVertLevels nCells Time ) o zMid - -
+var real zBot ( nVertLevels nCells Time ) o zBot - -
var real zSurface ( nCells Time ) o zSurface - -
+var real zMidEdge ( nVertLevels nEdges Time ) o zMidEdge - -
+var real zBotEdge ( nVertLevels nEdges Time ) o zBotEdge - -
+var real zSurfaceEdge ( nEdges Time ) o zSurfaceEdge - -
var real pmid ( nVertLevels nCells Time ) o pmid - -
var real pbot ( nVertLevels nCells Time ) o pbot - -
var real MontPot ( nVertLevels nCells Time ) o MontPot - -
-# mrp 100112 end
# Other diagnostic variables: neither read nor written to any files
Modified: trunk/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- trunk/mpas/src/core_ocean/module_test_cases.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/core_ocean/module_test_cases.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -23,9 +23,6 @@
integer :: i, iCell, iEdge, iVtx, iLevel
type (block_type), pointer :: block_ptr
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
- real (kind=RKIND) :: delta_rho
- integer :: nCells, nEdges, nVertices, nVertLevels
if (config_test_case == 0) then
write(0,*) 'Using initial conditions supplied in input file'
@@ -76,84 +73,16 @@
stop
end if
- ! mrp 100121:
- !
- ! Initialize u_src, rho, alpha
- ! This is a temporary fix until everything is specified correctly
- ! in an initial conditions file.
- !
block_ptr => domain % blocklist
do while (associated(block_ptr))
- h => block_ptr % time_levs(1) % state % h % array
- u => block_ptr % time_levs(1) % state % u % array
- rho => block_ptr % time_levs(1) % state % rho % array
- u_src => block_ptr % mesh % u_src % array
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % time_levs(1) % state, &
+ block_ptr % time_levs(i) % state)
+ end do
- nCells = block_ptr % mesh % nCells
- nEdges = block_ptr % mesh % nEdges
- nVertices = block_ptr % mesh % nVertices
- nVertLevels = block_ptr % mesh % nVertLevels
-
- ! Momentum forcing u_src
- if (config_test_case > 0) then
- ! for shallow water test cases:
- u_src = 0.0
- elseif (config_test_case == 0 ) then
- ! for rectangular basin:
- do iEdge=1,nEdges
- u_src(1:nVertLevels, iEdge) = u_src(1, iEdge) / nVertLevels
- end do
- endif
-
- ! define the density for multiple layers
- delta_rho=0.0
- do iLevel = 1,nVertLevels
- rho(iLevel,1) = delta_rho*(iLevel-1)
- enddo
- rho(:,1) = rho(:,1) - sum(rho(:,1))/nVertLevels + 1000.0
- do iLevel = 1,nVertLevels
- rho(iLevel,:) = rho(iLevel,1)
- enddo
-
-#ifdef EXPAND_LEVELS
- print '(10a)', 'EXPAND_LEVELS compiler flag is on.', &
- ' Copying h and u from k=1 to other levels.'
- print '(a,i5)', 'EXPAND_LEVELS =', EXPAND_LEVELS
- print '(a,i5)', 'nVertLevels =', nVertLevels
- do iCell=1,nCells
- ! make the total thickness equal to the single-layer thickness:
- h(1:nVertLevels, iCell) = h(1,iCell) / nVertLevels
- end do
-
- do iEdge=1,nEdges
- u(2:nVertLevels, iEdge) = u(1,iEdge)
- end do
-#else
- print '(10a)', 'EXPAND_LEVELS compiler flag is off.'
-#endif
-
- do i=2,nTimeLevs
- call copy_state(block_ptr % time_levs(1) % state, &
- block_ptr % time_levs(i) % state)
- end do
-
- ! print some diagnostics
- print '(10a)', 'ilevel',&
- ' rho ',&
- ' min u max u ',&
- ' min h max h ',&
- ' min u_src max u_src '
- do iLevel = 1,nVertLevels
- print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
- minval(u(iLevel,:)), maxval(u(iLevel,:)), &
- minval(h(iLevel,:)), maxval(h(iLevel,:)), &
- minval(u_src(iLevel,:)), maxval(u_src(iLevel,:))
- enddo
-
- block_ptr => block_ptr % next
+ block_ptr => block_ptr % next
end do
- ! mrp 100121 end
end subroutine setup_sw_test_case
Modified: trunk/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_ocean/module_time_integration.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/core_ocean/module_time_integration.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -38,19 +38,18 @@
stop
end if
- block => domain % blocklist
- do while (associated(block))
- block % time_levs(2) % state % xtime % scalar = block % time_levs(1) % state % xtime % scalar + dt
- ! mrp 100310 I added this to avoid running with NaNs
- if (isNaN(sum(block % time_levs(2) % state % u % array))) then
- print *, 'Stopping: NaN detected'
- call MPI_abort(MPI_COMM_WORLD,errorcode,ierr)
- endif
- ! mrp 100310 end
+ block => domain % blocklist
+ do while (associated(block))
+ block % time_levs(2) % state % xtime % scalar = block % time_levs(1) % state % xtime % scalar + dt
+ ! ! mrp 100310 I added this to avoid running with NaNs
+ ! if (isNaN(sum(block % time_levs(2) % state % u % array))) then
+ ! print *, 'Stopping: NaN detected'
+ ! call MPI_abort(MPI_COMM_WORLD,errorcode,ierr)
+ ! endif
+ ! ! mrp 100310 end
+ block => block % next
+ end do
- block => block % next
- end do
-
end subroutine timestep
@@ -136,7 +135,7 @@
call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call enforce_uBC(block % intermediate_step(TEND), block % mesh)
+ call enforce_boundaryEdge(block % intermediate_step(TEND), block % mesh)
block => block % next
end do
@@ -253,12 +252,13 @@
type (grid_meta), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias, dist
+ real (kind=RKIND), allocatable, dimension(:) :: fluxVert
integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zSurfaceEdge
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, divergence
+ circulation, vorticity, ke, pv_edge, divergence, zBotEdge, zMidEdge
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: u_diffusion, visc
@@ -267,10 +267,8 @@
real (kind=RKIND), dimension(:,:), pointer :: MontPot
!mrp 100112 end
-!ocean
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
-!ocean
visc = config_visc
@@ -284,9 +282,10 @@
ke => s % ke % array
pv_edge => s % pv_edge % array
vh => s % vh % array
- !mrp 100112:
MontPot => s % MontPot % array
- !mrp 100112 end
+ zBotEdge => s % zBotEdge % array
+ zSurfaceEdge=> s % zSurfaceEdge % array
+ zMidEdge => s % zMidEdge % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -314,9 +313,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
-!ocean
u_src => grid % u_src % array
-!ocean
!
@@ -326,13 +323,13 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
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 > 0) then
+ if (cell2 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell2) = tend_h(k,cell2) + flux
@@ -372,28 +369,45 @@
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
end do
- ! mrp 100112, this is the original shallow water formulation with grad H:
- !tend_u(k,iEdge) = &
- ! q &
- ! + u_diffusion &
- ! - ( ke(k,cell2) - ke(k,cell1) &
- ! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ! ) / dcEdge(iEdge)
- ! mrp 100112, changed to grad Montgomery potential:
+
tend_u(k,iEdge) = &
q &
+ u_diffusion &
- ( ke(k,cell2) - ke(k,cell1) &
+ MontPot(k,cell2) - MontPot(k,cell1) &
) / dcEdge(iEdge)
- ! mrp 100112 end
+ end do
+ end do
-!ocean
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_src(k,iEdge)/rho_ref/h_edge(k,iEdge)
-!ocean
+ do iEdge=1,grid % nEdgesSolve
+ ! surface forcing
+ tend_u(1,iEdge) = tend_u(1,iEdge) + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
- end do
- end do
+ ! bottom drag
+ tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) - u(nVertLevels,iEdge)/(100.0*86400.0)
+
+ enddo
+
+! vertical mixing
+ allocate(fluxVert(0:nVertLevels))
+ do iEdge=1,grid % nEdgesSolve
+ fluxVert(0) = 0.0
+ fluxVert(nVertLevels) = 0.0
+ do k=nVertLevels-1,1,-1
+ fluxVert(k) = ( u(k,iEdge) - u(k+1,iEdge) ) / (zMidEdge(k,iEdge) - zMidEdge(k+1,iEdge))
+ enddo
+ fluxVert = 1.0e-4 * fluxVert
+ do k=1,nVertLevels
+ if(k.eq.1) then
+ dist = zSurfaceEdge(iEdge) - zBotEdge(k,iEdge)
+ else
+ dist = zBotEdge(k-1,iEdge) - zBotEdge(k,iEdge)
+ endif
+ tend_u(k,iEdge) = tend_u(k,iEdge) + (fluxVert(k-1) - fluxVert(k))/dist
+ enddo
+ enddo
+ deallocate(fluxVert)
+
#endif
#ifdef NCAR_FORMULATION
@@ -447,7 +461,7 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
tracer_edge = 0.5 * (s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
@@ -493,15 +507,13 @@
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
- !mrp 100112:
real (kind=RKIND), dimension(:,:), pointer :: &
- zmid, zbot, pmid, pbot, MontPot, rho
- real (kind=RKIND), dimension(:), pointer :: zSurface
+ zMid, zBot, pmid, pbot, MontPot, rho, zBotEdge, zMidEdge
+ real (kind=RKIND), dimension(:), pointer :: zSurface, zSurfaceEdge
real (kind=RKIND) :: delta_p
character :: c1*6
- !mrp 100112 end
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -524,13 +536,15 @@
gradPVt => s % gradPVt % array
!mrp 100112:
rho => s % rho % array
- zmid => s % zmid % array
- zbot => s % zbot % array
+ zMid => s % zMid % array
+ zBot => s % zBot % array
+ zMidEdge => s % zMidEdge % array
+ zBotEdge => s % zBotEdge % array
+ zSurfaceEdge=> s % zSurfaceEdge % array
pmid => s % pmid % array
pbot => s % pbot % array
MontPot => s % MontPot % array
zSurface => s % zSurface % array
- !mrp 100112 end
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -555,7 +569,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
!
! Compute height on cell edges at velocity locations
@@ -563,24 +577,39 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ 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
+ h_edge(k,iEdge) = h(k,cell1)
+ end do
+ elseif(cell2 <= nCells) then
+ do k=1,nVertLevels
+ h_edge(k,iEdge) = h(k,cell2)
+ end do
end if
end do
+
!
+ ! 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
+ !
+ u(:,nEdges+1) = 0.0
+
+ !
! Compute circulation and relative vorticity at each vertex
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
+ if (verticesOnEdge(1,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
end do
end if
- if (verticesOnEdge(2,iEdge) > 0) then
+ if (verticesOnEdge(2,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
end do
@@ -592,6 +621,7 @@
end do
end do
+
!
! Compute the divergence at each cell center
!
@@ -599,12 +629,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
enddo
endif
- if(cell2 > 0) then
+ if(cell2 <= nCells) then
do k=1,nVertLevels
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
@@ -640,7 +670,7 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
+ if (eoe <= nEdges) then
do k = 1,nVertLevels
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
@@ -664,14 +694,13 @@
#endif
- ! tdr
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) <= 0) cycle VTX_LOOP
+ if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
end do
do k=1,nVertLevels
h_vertex = 0.0
@@ -679,14 +708,12 @@
h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
h_vertex = h_vertex / areaTriangle(iVertex)
-
+
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
end do VTX_LOOP
- ! tdr
- ! tdr
!
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
@@ -698,7 +725,6 @@
enddo
enddo
- ! tdr
!
! Compute pv at the edges
! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
@@ -707,16 +733,14 @@
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
+ if(iEdge <= nEdges) then
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
endif
end do
end do
- ! tdr
- ! tdr
!
! Modify PV edge with upstream bias.
!
@@ -727,7 +751,6 @@
enddo
- ! tdr
!
! Compute pv at cell centers
! ( this computes pv_cell for all real cells and distance-1 ghost cells )
@@ -736,31 +759,29 @@
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
- if( iCell > 0) then
+ if (iCell <= nCells) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
enddo
endif
enddo
enddo
- ! tdr
- ! tdr
!
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
+ if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
enddo
endif
enddo
- ! tdr
+
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
@@ -770,25 +791,6 @@
enddo
!
- ! set pv_edge = fEdge / h_edge at boundary points
- !
- if (maxval(uBC).gt.0) then
- do iEdge = 1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
- if(cell1.gt.0) h1 = h(k,cell1)
- if(cell2.gt.0) h2 = h(k,cell2)
- pv_edge(k,iEdge) = fEdge(iEdge) / ( max(h1,h2) )
- v(k,iEdge) = 0.0
- endif
- enddo
- enddo
- endif
-
- !mrp 100112:
- !
! Compute mid- and bottom-depth of each layer, from bottom up.
! Then compute mid- and bottom-pressure of each layer, and
! Montgomery Potential, from top down
@@ -798,14 +800,14 @@
! h_s is ocean topography: elevation above lowest point,
! and is positive. z coordinates are pos upward, and z=0
! is at lowest ocean point.
- zbot(nVertLevels,iCell) = h_s(iCell)
- zmid(nVertLevels,iCell) = zbot(nVertLevels,iCell) + 0.5*h(nVertLevels,iCell)
+ zBot(nVertLevels,iCell) = h_s(iCell)
+ zMid(nVertLevels,iCell) = zBot(nVertLevels,iCell) + 0.5*h(nVertLevels,iCell)
do k=nVertLevels-1,1,-1
- zbot(k,iCell) = zbot(k+1,iCell) + h(k+1,iCell)
- zmid(k,iCell) = zbot(k,iCell) + 0.5*h(k,iCell)
+ zBot(k,iCell) = zBot(k+1,iCell) + h(k+1,iCell)
+ zMid(k,iCell) = zBot(k,iCell) + 0.5*h(k,iCell)
end do
- ! rather than using zbot(0,iCell), I am adding another 2D variable.
- zSurface(iCell) = zbot(1,iCell) + h(1,iCell)
+ ! rather than using zBot(0,iCell), I am adding another 2D variable.
+ zSurface(iCell) = zBot(1,iCell) + h(1,iCell)
! assume pressure at the surface is zero for now.
pmid(1,iCell) = 0.5*rho(1,iCell)*gravity* h(1,iCell) ! + psurf(iCell)
@@ -821,18 +823,30 @@
+ pbot(k-1,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
end do
end do
- !mrp 100112 end
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if(cell1<=nCells .and. cell2<=nCells) then
+ do k=1,nVertLevels
+ zBotEdge(k,iEdge) = (zBot(k,cell1)+zBot(k,cell2))/2.0
+ zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
+ enddo
+ zSurfaceEdge(iEdge) = (zSurface(cell1) + zSurface(cell2))/2.0
+ endif
+ enddo
+
+
end subroutine compute_solve_diagnostics
- subroutine enforce_uBC(tend, grid)
+ subroutine enforce_boundaryEdge(tend, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
!
! Input: grid - grid metadata
!
- ! Output: tend_u set to zero at uBC == 1 locations
+ ! Output: tend_u set to zero at boundaryEdge == 1 locations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -841,7 +855,7 @@
type (grid_state), intent(inout) :: tend
type (grid_meta), intent(in) :: grid
- integer, dimension(:,:), pointer :: uBC
+ integer, dimension(:,:), pointer :: boundaryEdge
real (kind=RKIND), dimension(:,:), pointer :: tend_u
integer :: nCells, nEdges, nVertices, nVertLevels
integer :: iEdge, k
@@ -851,21 +865,21 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
tend_u => tend % u % array
- if(maxval(uBC).le.0) return
+ if(maxval(boundaryEdge).le.0) return
do iEdge = 1,nEdges
do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
+ if(boundaryEdge(k,iEdge).eq.1) then
tend_u(k,iEdge) = 0.0
endif
enddo
enddo
- end subroutine enforce_uBC
+ end subroutine enforce_boundaryEdge
end module time_integration
Modified: trunk/mpas/src/core_sw/Registry
===================================================================
--- trunk/mpas/src/core_sw/Registry 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/core_sw/Registry 2010-05-06 21:22:50 UTC (rev 258)
@@ -81,7 +81,8 @@
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
# Boundary conditions: read from input, saved in restart and written to output
-var integer uBC ( nVertLevels nEdges ) iro uBC - -
+var integer boundaryEdge ( nVertLevels nEdges ) iro boundaryEdge - -
+var integer boundaryVertex ( nVertLevels nVertices ) iro boundaryVertex - -
# Prognostic variables: read from input, saved in restart, and written to output
var real u ( nVertLevels nEdges Time ) iro u - -
Modified: trunk/mpas/src/core_sw/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_sw/module_time_integration.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/core_sw/module_time_integration.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -125,7 +125,7 @@
do while (associated(block))
call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call enforce_uBC(block % intermediate_step(TEND), block % mesh)
+ call enforce_boundaryEdge(block % intermediate_step(TEND), block % mesh)
block => block % next
end do
@@ -297,13 +297,13 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
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 > 0) then
+ if (cell2 <= nCells) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell2) = tend_h(k,cell2) + flux
@@ -343,6 +343,7 @@
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
end do
+
tend_u(k,iEdge) = &
q &
+ u_diffusion &
@@ -404,7 +405,7 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
tracer_edge = 0.5 * (s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
@@ -450,7 +451,7 @@
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -495,7 +496,7 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
+ boundaryEdge => grid % boundaryEdge % array
!
! Compute height on cell edges at velocity locations
@@ -503,7 +504,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
+ 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
@@ -511,16 +512,22 @@
end do
!
+ ! 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
+ !
+ u(:,nEdges+1) = 0.0
+
+ !
! Compute circulation and relative vorticity at each vertex
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) > 0) then
+ if (verticesOnEdge(1,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
end do
end if
- if (verticesOnEdge(2,iEdge) > 0) then
+ if (verticesOnEdge(2,iEdge) <= nVertices) then
do k=1,nVertLevels
circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
end do
@@ -532,6 +539,7 @@
end do
end do
+
!
! Compute the divergence at each cell center
!
@@ -539,12 +547,12 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0) then
+ if (cell1 <= nCells) then
do k=1,nVertLevels
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
enddo
endif
- if(cell2 > 0) then
+ if(cell2 <= nCells) then
do k=1,nVertLevels
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
@@ -580,7 +588,7 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe > 0) then
+ if (eoe <= nEdges) then
do k = 1,nVertLevels
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
@@ -604,14 +612,13 @@
#endif
- ! tdr
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) <= 0) cycle VTX_LOOP
+ if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
end do
do k=1,nVertLevels
h_vertex = 0.0
@@ -623,10 +630,8 @@
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
end do VTX_LOOP
- ! tdr
- ! tdr
!
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
@@ -638,7 +643,6 @@
enddo
enddo
- ! tdr
!
! Compute pv at the edges
! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
@@ -647,16 +651,14 @@
do iVertex = 1,nVertices
do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
- if(iEdge > 0) then
+ if(iEdge <= nEdges) then
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
endif
end do
end do
- ! tdr
- ! tdr
!
! Modify PV edge with upstream bias.
!
@@ -667,7 +669,6 @@
enddo
- ! tdr
!
! Compute pv at cell centers
! ( this computes pv_cell for all real cells and distance-1 ghost cells )
@@ -676,30 +677,28 @@
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
- if( iCell > 0) then
+ if (iCell <= nCells) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
enddo
endif
enddo
enddo
- ! tdr
- ! tdr
+
!
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) > 0 .and. cellsOnEdge(2,iEdge) > 0) then
+ if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
enddo
endif
enddo
- ! tdr
! Modify PV edge with upstream bias.
!
@@ -712,32 +711,38 @@
!
! set pv_edge = fEdge / h_edge at boundary points
!
- if (maxval(uBC).gt.0) then
- do iEdge = 1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
- if(cell1.gt.0) h1 = h(k,cell1)
- if(cell2.gt.0) h2 = h(k,cell2)
- pv_edge(k,iEdge) = fEdge(iEdge) / ( max(h1,h2) )
- v(k,iEdge) = 0.0
- endif
- enddo
- enddo
- endif
+ ! if (maxval(boundaryEdge).ge.0) then
+ ! do iEdge = 1,nEdges
+ ! cell1 = cellsOnEdge(1,iEdge)
+ ! cell2 = cellsOnEdge(2,iEdge)
+ ! do k = 1,nVertLevels
+ ! if(boundaryEdge(k,iEdge).eq.1) then
+ ! v(k,iEdge) = 0.0
+ ! if(cell1.gt.0) then
+ ! h1 = h(k,cell1)
+ ! pv_edge(k,iEdge) = fEdge(iEdge) / h1
+ ! h_edge(k,iEdge) = h1
+ ! else
+ ! h2 = h(k,cell2)
+ ! pv_edge(k,iEdge) = fEdge(iEdge) / h2
+ ! h_edge(k,iEdge) = h2
+ ! endif
+ ! endif
+ ! enddo
+ ! enddo
+ ! endif
end subroutine compute_solve_diagnostics
- subroutine enforce_uBC(tend, grid)
+ subroutine enforce_boundaryEdge(tend, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Enforce any boundary conditions on the normal velocity at each edge
!
! Input: grid - grid metadata
!
- ! Output: tend_u set to zero at uBC == 1 locations
+ ! Output: tend_u set to zero at boundaryEdge == 1 locations
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -746,7 +751,7 @@
type (grid_state), intent(inout) :: tend
type (grid_meta), intent(in) :: grid
- integer, dimension(:,:), pointer :: uBC
+ integer, dimension(:,:), pointer :: boundaryEdge
real (kind=RKIND), dimension(:,:), pointer :: tend_u
integer :: nCells, nEdges, nVertices, nVertLevels
integer :: iEdge, k
@@ -756,22 +761,22 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- uBC => grid % uBC % array
- tend_u => tend % u % array
+ boundaryEdge => grid % boundaryEdge % array
+ tend_u => tend % u % array
- if(maxval(uBC).le.0) return
+ if(maxval(boundaryEdge).le.0) return
do iEdge = 1,nEdges
do k = 1,nVertLevels
- if(uBC(k,iEdge).eq.1) then
+ if(boundaryEdge(k,iEdge).eq.1) then
tend_u(k,iEdge) = 0.0
endif
enddo
enddo
- end subroutine enforce_uBC
+ end subroutine enforce_boundaryEdge
end module time_integration
Modified: trunk/mpas/src/framework/module_block_decomp.F
===================================================================
--- trunk/mpas/src/framework/module_block_decomp.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/framework/module_block_decomp.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -142,7 +142,11 @@
edgeIDListLocal(:) = edgeIDList(:)
do i=1,nEdges
- if (hash_search(h, cellsOnEdge(1,i))) then
+ do j=1,maxCells
+ if (cellsOnEdge(j,i) /= 0) exit
+ end do
+if (j > maxCells) write(0,*) 'Error in block_decomp_partitioned_edge_list: edge/vertex is not adjacent to any valid cells'
+ if (hash_search(h, cellsOnEdge(j,i))) then
lastEdge = lastEdge + 1
edgeIDList(lastEdge) = edgeIDListLocal(i)
else
@@ -231,8 +235,10 @@
do i=1,local_graph_info % nVertices
do j=1,local_graph_info % nAdjacent(i)
- if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
- call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ if (local_graph_info % adjacencyList(j,i) /= 0) then
+ if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
+ call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ end if
end if
end do
end do
@@ -264,10 +270,12 @@
write(0,*) 'block_decomp_add_halo: Somehow we don''t have the right number of non-ghost cells.'
do i=1,local_graph_info % nVertices
do j=1,local_graph_info % nAdjacent(i)
- if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
- call hash_insert(h, local_graph_info % adjacencyList(j,i))
- local_graph_with_halo % vertexID(k) = local_graph_info % adjacencyList(j,i)
- k = k + 1
+ if (local_graph_info % adjacencyList(j,i) /= 0) then
+ if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
+ call hash_insert(h, local_graph_info % adjacencyList(j,i))
+ local_graph_with_halo % vertexID(k) = local_graph_info % adjacencyList(j,i)
+ k = k + 1
+ end if
end if
end do
end do
Modified: trunk/mpas/src/framework/module_dmpar.F
===================================================================
--- trunk/mpas/src/framework/module_dmpar.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/framework/module_dmpar.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -707,8 +707,8 @@
implicit none
type (dm_info), intent(in) :: dminfo
- integer, dimension(nOwnedList), intent(in) :: arrayIn
- integer, dimension(nNeededList), intent(inout) :: arrayOut
+ integer, dimension(*), intent(in) :: arrayIn
+ integer, dimension(*), intent(inout) :: arrayOut
integer, intent(in) :: nOwnedList, nNeededList
type (exchange_list), pointer :: sendList, recvList
@@ -784,7 +784,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:) = arrayIn(:)
+ arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
end if
#endif
@@ -797,8 +797,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, nOwnedList, nNeededList
- integer, dimension(dim1,nOwnedList), intent(in) :: arrayIn
- integer, dimension(dim1,nNeededList), intent(inout) :: arrayOut
+ integer, dimension(dim1,*), intent(in) :: arrayIn
+ integer, dimension(dim1,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -876,7 +876,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:) = arrayIn(:,:)
+ arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
end if
#endif
@@ -888,8 +888,8 @@
implicit none
type (dm_info), intent(in) :: dminfo
- real (kind=RKIND), dimension(nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(*), intent(inout) :: arrayOut
integer, intent(in) :: nOwnedList, nNeededList
type (exchange_list), pointer :: sendList, recvList
@@ -965,7 +965,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:) = arrayIn(:)
+ arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
end if
#endif
@@ -978,8 +978,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, nOwnedList, nNeededList
- real (kind=RKIND), dimension(dim1,nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(dim1,nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(dim1,*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(dim1,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1057,7 +1057,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:) = arrayIn(:,:)
+ arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
end if
#endif
@@ -1070,8 +1070,8 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2, nOwnedList, nNeededList
- real (kind=RKIND), dimension(dim1,dim2,nOwnedList), intent(in) :: arrayIn
- real (kind=RKIND), dimension(dim1,dim2,nNeededList), intent(inout) :: arrayOut
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(in) :: arrayIn
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(inout) :: arrayOut
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1149,7 +1149,7 @@
write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
- arrayOut(:,:,:) = arrayIn(:,:,:)
+ arrayOut(:,:,1:nNeededList) = arrayIn(:,:,1:nOwnedList)
end if
#endif
@@ -1161,7 +1161,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startPackIdx
- integer, dimension(nField), intent(in) :: field
+ integer, dimension(*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
integer, dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1188,7 +1188,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startPackIdx
- integer, dimension(ds:de,1:nField), intent(in) :: field
+ integer, dimension(ds:de,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
integer, dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1222,7 +1222,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(nField), intent(in) :: field
+ real (kind=RKIND), dimension(*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1249,7 +1249,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(ds:de,1:nField), intent(in) :: field
+ real (kind=RKIND), dimension(ds:de,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1283,7 +1283,7 @@
implicit none
integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startPackIdx
- real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,1:nField), intent(in) :: field
+ real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,*), intent(in) :: field
type (exchange_list), intent(in) :: sendList
real (kind=RKIND), dimension(nBuffer), intent(out) :: buffer
integer, intent(inout) :: nPacked, lastPackedIdx
@@ -1321,7 +1321,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startUnpackIdx
- integer, dimension(nField), intent(inout) :: field
+ integer, dimension(*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
integer, dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1348,7 +1348,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startUnpackIdx
- integer, dimension(ds:de,1:nField), intent(inout) :: field
+ integer, dimension(ds:de,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
integer, dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1377,7 +1377,7 @@
implicit none
integer, intent(in) :: nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(nField), intent(inout) :: field
+ real (kind=RKIND), dimension(*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1404,7 +1404,7 @@
implicit none
integer, intent(in) :: ds, de, nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(ds:de,1:nField), intent(inout) :: field
+ real (kind=RKIND), dimension(ds:de,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1434,7 +1434,7 @@
implicit none
integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startUnpackIdx
- real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,1:nField), intent(inout) :: field
+ real (kind=RKIND), dimension(d1s:d1e,d2s:d2e,*), intent(inout) :: field
type (exchange_list), intent(in) :: recvList
real (kind=RKIND), dimension(nBuffer), intent(in) :: buffer
integer, intent(inout) :: nUnpacked, lastUnpackedIdx
@@ -1468,7 +1468,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1
- real (kind=RKIND), dimension(dim1), intent(inout) :: array
+ real (kind=RKIND), dimension(*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1528,7 +1528,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2
- real (kind=RKIND), dimension(dim1,dim2), intent(inout) :: array
+ real (kind=RKIND), dimension(dim1,*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
@@ -1592,7 +1592,7 @@
type (dm_info), intent(in) :: dminfo
integer, intent(in) :: dim1, dim2, dim3
- real (kind=RKIND), dimension(dim1,dim2,dim3), intent(inout) :: array
+ real (kind=RKIND), dimension(dim1,dim2,*), intent(inout) :: array
type (exchange_list), pointer :: sendList, recvList
type (exchange_list), pointer :: sendListPtr, recvListPtr
Modified: trunk/mpas/src/framework/module_io_input.F
===================================================================
--- trunk/mpas/src/framework/module_io_input.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/framework/module_io_input.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -805,7 +805,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnCell % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnCell % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnCell % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnCell % array(j,i) = 0
end if
k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
@@ -813,7 +814,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnCell % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnCell % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnCell % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnCell % array(j,i) = 0
end if
k = binary_search(vertexIDSorted, 2, 1, domain % blocklist % mesh % nVertices, &
@@ -821,7 +823,8 @@
if (k <= domain % blocklist % mesh % nVertices) then
domain % blocklist % mesh % verticesOnCell % array(j,i) = vertexIDSorted(2,k)
else
- domain % blocklist % mesh % verticesOnCell % array(j,i) = 0
+ domain % blocklist % mesh % verticesOnCell % array(j,i) = domain % blocklist % mesh % nVertices + 1
+! domain % blocklist % mesh % verticesOnCell % array(j,i) = 0
end if
end do
@@ -835,7 +838,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnEdge % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnEdge % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnEdge % array(j,i) = 0
end if
k = binary_search(vertexIDSorted, 2, 1, domain % blocklist % mesh % nVertices, &
@@ -843,7 +847,8 @@
if (k <= domain % blocklist % mesh % nVertices) then
domain % blocklist % mesh % verticesOnEdge % array(j,i) = vertexIDSorted(2,k)
else
- domain % blocklist % mesh % verticesOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % verticesOnEdge % array(j,i) = domain % blocklist % mesh % nVertices + 1
+! domain % blocklist % mesh % verticesOnEdge % array(j,i) = 0
end if
end do
@@ -855,7 +860,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnEdge % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnEdge % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnEdge % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnEdge % array(j,i) = 0
end if
end do
@@ -869,7 +875,8 @@
if (k <= domain % blocklist % mesh % nCells) then
domain % blocklist % mesh % cellsOnVertex % array(j,i) = cellIDSorted(2,k)
else
- domain % blocklist % mesh % cellsOnVertex % array(j,i) = 0
+ domain % blocklist % mesh % cellsOnVertex % array(j,i) = domain % blocklist % mesh % nCells + 1
+! domain % blocklist % mesh % cellsOnVertex % array(j,i) = 0
end if
k = binary_search(edgeIDSorted, 2, 1, domain % blocklist % mesh % nEdges, &
@@ -877,7 +884,8 @@
if (k <= domain % blocklist % mesh % nEdges) then
domain % blocklist % mesh % edgesOnVertex % array(j,i) = edgeIDSorted(2,k)
else
- domain % blocklist % mesh % edgesOnVertex % array(j,i) = 0
+ domain % blocklist % mesh % edgesOnVertex % array(j,i) = domain % blocklist % mesh % nEdges + 1
+! domain % blocklist % mesh % edgesOnVertex % array(j,i) = 0
end if
end do
Modified: trunk/mpas/src/framework/module_io_output.F
===================================================================
--- trunk/mpas/src/framework/module_io_output.F 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/framework/module_io_output.F 2010-05-06 21:22:50 UTC (rev 258)
@@ -188,8 +188,12 @@
domain % blocklist % mesh % edgesOnEdge % array(j,i))
end do
do j=domain % blocklist % mesh % nEdgesOnEdge % array(i)+1,2*domain % blocklist % mesh % maxEdges
- edgesOnEdge(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
+ if(domain % blocklist % mesh % nEdgesOnEdge % array(i) .eq. 0) then
+ edgesOnEdge(j,i) = domain % blocklist % mesh % nEdgesSolve + 1
+ else
+ edgesOnEdge(j,i) = domain % blocklist % mesh % indexToEdgeID % array( &
domain % blocklist % mesh % nEdgesOnEdge % array(i))
+ endif
end do
end do
do i=1,domain % blocklist % mesh % nVerticesSolve
Modified: trunk/mpas/src/registry/gen_inc.c
===================================================================
--- trunk/mpas/src/registry/gen_inc.c 2010-05-06 20:48:16 UTC (rev 257)
+++ trunk/mpas/src/registry/gen_inc.c 2010-05-06 21:22:50 UTC (rev 258)
@@ -387,10 +387,20 @@
fortprintf(fd, " allocate(g %% %s %% ioinfo)</font>
<font color="red">", var_ptr2->super_array);
fortprintf(fd, " allocate(g %% %s %% array(%i, ", var_ptr2->super_array, i);
dimlist_ptr = var_ptr2->dimlist;
- fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -417,10 +427,20 @@
if (var_ptr->ndims > 0) {
fortprintf(fd, " allocate(g %% %s %% array(", var_ptr->name_in_code);
dimlist_ptr = var_ptr->dimlist;
- fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -516,10 +536,20 @@
fortprintf(fd, " allocate(s %% %s %% ioinfo)</font>
<font color="red">", var_ptr2->super_array);
fortprintf(fd, " allocate(s %% %s %% array(%i, ", var_ptr2->super_array, i);
dimlist_ptr = var_ptr2->dimlist;
- fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, "b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, ", b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="gray">");
@@ -546,10 +576,20 @@
fortprintf(fd, " allocate(s %% %s %% ioinfo)</font>
<font color="red">", var_ptr->name_in_code);
fortprintf(fd, " allocate(s %% %s %% array(", var_ptr->name_in_code);
dimlist_ptr = var_ptr->dimlist;
- fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, "b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, "b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
while (dimlist_ptr) {
- fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ fortprintf(fd, ", b %% mesh %% %s + 1", dimlist_ptr->dim->name_in_code);
+ else
+ fortprintf(fd, ", b %% mesh %% %s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
}
fortprintf(fd, "))</font>
<font color="black">");
</font>
</pre>