<p><b>mpetersen@lanl.gov</b> 2010-06-11 11:52:52 -0600 (Fri, 11 Jun 2010)</p><p>Adding topography. All loops in vertical now end at maxLevelCell, maxLevelEdge, or maxLevelVertex, rather than nVertLevels.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/Registry 2010-06-10 20:40:08 UTC (rev 342)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/Registry 2010-06-11 17:52:52 UTC (rev 343)
@@ -76,7 +76,6 @@
var integer nEdgesOnEdge ( nEdges ) iro nEdgesOnEdge - -
var integer edgesOnCell ( maxEdges nCells ) iro edgesOnCell - -
var integer edgesOnEdge ( maxEdges2 nEdges ) iro edgesOnEdge - -
-# var integer maxVertLevel ( nCells ) iro maxVertLevel - -
var integer maxLevelCell ( nCells ) iro maxLevelCell - -
var integer maxLevelEdge ( nEdges ) ro maxLevelEdge - -
var integer maxLevelVertex ( nVertices ) ro maxLevelVertex - -
Modified: branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F 2010-06-10 20:40:08 UTC (rev 342)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F 2010-06-11 17:52:52 UTC (rev 343)
@@ -3,8 +3,8 @@
use grid_types
use configure
use constants
+ use dmpar
-
contains
@@ -21,22 +21,23 @@
type (domain_type), intent(inout) :: domain
- integer :: i, iCell, iEdge, iVtx, iLevel
- type (block_type), pointer :: block_ptr
+ integer :: i, iCell, iEdge, iVertex, iLevel
+ type (block_type), pointer :: block
type (dm_info) :: dminfo
! mrp 100507: for diagnostic output
integer :: iTracer, cell1, cell2, cell
real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
- hZLevel, zMidZLevel, zTopZLevel
+ 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
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
integer, dimension(:), pointer :: &
- maxLevelCell, maxLevelEdge, maxLevelVertex
+ maxLevelCell, maxLevelEdge, maxLevelVertex
integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex
+ cellsOnEdge, cellsOnVertex, boundaryEdge
+
! mrp 100507 end: for diagnostic output
if (config_test_case == 0) then
@@ -46,40 +47,40 @@
write(0,*) ' Setting up shallow water test case 1:'
write(0,*) ' Advection of Cosine Bell over the Pole'
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_1(block_ptr % mesh, block_ptr % time_levs(1) % state)
- block_ptr => block_ptr % next
+ block => domain % blocklist
+ do while (associated(block))
+ call sw_test_case_1(block % mesh, block % time_levs(1) % state)
+ block => block % next
end do
else if (config_test_case == 2) then
write(0,*) ' Setup shallow water test case 2: '// &
'Global Steady State Nonlinear Zonal Geostrophic Flow'
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_2(block_ptr % mesh, block_ptr % time_levs(1) % state)
- block_ptr => block_ptr % next
+ block => domain % blocklist
+ do while (associated(block))
+ call sw_test_case_2(block % mesh, block % time_levs(1) % state)
+ block => block % next
end do
else if (config_test_case == 5) then
write(0,*) ' Setup shallow water test case 5:'// &
' Zonal Flow over an Isolated Mountain'
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_5(block_ptr % mesh, block_ptr % time_levs(1) % state)
- block_ptr => block_ptr % next
+ block => domain % blocklist
+ do while (associated(block))
+ call sw_test_case_5(block % mesh, block % time_levs(1) % state)
+ block => block % next
end do
else if (config_test_case == 6) then
write(0,*) ' Set up shallow water test case 6:'
write(0,*) ' Rossby-Haurwitz Wave'
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_6(block_ptr % mesh, block_ptr % time_levs(1) % state)
- block_ptr => block_ptr % next
+ block => domain % blocklist
+ do while (associated(block))
+ call sw_test_case_6(block % mesh, block % time_levs(1) % state)
+ block => block % next
end do
else
@@ -89,43 +90,92 @@
call dmpar_abort(dminfo)
end if
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
+ block => domain % blocklist
+ do while (associated(block))
do i=2,nTimeLevs
- call copy_state(block_ptr % time_levs(1) % state, &
- block_ptr % time_levs(i) % state)
+ call copy_state(block % time_levs(1) % state, &
+ block % time_levs(i) % state)
end do
- block_ptr => block_ptr % next
+ block => block % next
end do
! Initialize z-level grid variables from h, read in from input 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
- tracers => block_ptr % time_levs(1) % state % tracers % array
- u_src => block_ptr % mesh % u_src % array
- xCell => block_ptr % mesh % xCell % array
- yCell => block_ptr % mesh % yCell % array
+ block => domain % blocklist
+ do while (associated(block))
+ h => block % time_levs(1) % state % h % array
+ u => block % time_levs(1) % state % u % array
+ rho => block % time_levs(1) % state % rho % array
+ tracers => block % time_levs(1) % state % tracers % array
+ u_src => block % mesh % u_src % array
+ xCell => block % mesh % xCell % array
+ yCell => block % mesh % yCell % array
+ latCell => block % mesh % latCell % array
+ lonCell => block % mesh % lonCell % array
- hZLevel => block_ptr % mesh % hZLevel % array
- zMidZLevel => block_ptr % mesh % zMidZLevel % array
- zTopZLevel => block_ptr % mesh % zTopZLevel % array
- maxLevelCell => block_ptr % mesh % maxLevelCell % array
- maxLevelEdge => block_ptr % mesh % maxLevelEdge % array
- maxLevelVertex => block_ptr % mesh % maxLevelVertex % array
- cellsOnEdge => block_ptr % mesh % cellsOnEdge % array
- cellsOnVertex => block_ptr % mesh % cellsOnVertex % array
+ hZLevel => block % mesh % hZLevel % array
+ zMidZLevel => block % mesh % zMidZLevel % array
+ zTopZLevel => block % mesh % zTopZLevel % array
+ maxLevelCell => block % mesh % maxLevelCell % array
+ maxLevelEdge => block % mesh % maxLevelEdge % array
+ maxLevelVertex => block % mesh % maxLevelVertex % array
+ cellsOnEdge => block % mesh % cellsOnEdge % array
+ cellsOnVertex => block % mesh % cellsOnVertex % array
+ boundaryEdge => block % mesh % boundaryEdge % array
- nCells = block_ptr % mesh % nCells
- nEdges = block_ptr % mesh % nEdges
- nVertices = block_ptr % mesh % nVertices
- nVertLevels = block_ptr % mesh % nVertLevels
- vertexDegree = block_ptr % mesh % vertexDegree
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertices = block % mesh % nVertices
+ nVertLevels = block % mesh % nVertLevels
+ vertexDegree = block % mesh % vertexDegree
+ ! mrp 100608 For now, use a flat bottom
+! maxLevelCell = nVertLevels
+
+ ! 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
+ ! enddo
+
+ ! maxLevelEdge is the minimum of the surrounding cells
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ maxLevelEdge(iEdge) = &
+ min( maxLevelCell(cell1), &
+ maxLevelCell(cell2) )
+ else
+ maxLevelEdge(iEdge) = 0
+ endif
+
+ end do
+
+ ! set boundary edge
+ boundaryEdge=1
+ do iEdge=1,nEdges
+ boundaryEdge(1:maxLevelEdge(iEdge),iEdge)=0
+ end do
+
+ ! maxLevelVertex is the maximum of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertex(iVertex) = 0
+ do i=1,vertexDegree
+ cell = cellsOnVertex(i,iVertex)
+ if (cell <= nCells) then
+ maxLevelVertex(iVertex) = &
+ max( maxLevelVertex(iVertex), &
+ maxLevelCell(cell) )
+ endif
+ end do
+ end do
+
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).
@@ -146,9 +196,12 @@
endif
! Set tracers, if not done in grid.nc file
- !tracers = 0.0
+ tracers = -2.0
do iCell = 1,nCells
- do iLevel = 1,nVertLevels
+ tracers(index_tracer2,1,iCell) = maxLevelCell(iCell)
+! tracers(index_tracer1,1,iCell) = latCell(iCell)
+! tracers(index_tracer2,1,iCell) = lonCell(iCell)
+ do iLevel = 1,maxLevelCell(iCell)
! for 20 layer test
! tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
! tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
@@ -158,8 +211,9 @@
tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
tracers(index_tracer1,iLevel,iCell) = 1.0
- tracers(index_tracer2,iLevel,iCell) = &
- (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
+! tracers(index_tracer1,iLevel,iCell) = maxLevelCell(iCell)
+! tracers(index_tracer2,iLevel,iCell) = &
+! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
rho(iLevel,iCell) = 1000.0*( 1.0 &
- 2.5e-4*tracers(index_temperature,iLevel,iCell) &
@@ -170,40 +224,7 @@
endif
- ! mrp 100608 For now, use a flat bottom
- maxLevelCell = nVertLevels
- ! maxLevelEdge is the minimum of the surrounding cells
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- if (cell1 <= nCells .and. cell2 <= nCells) then
- maxLevelEdge(iEdge) = &
- min( maxLevelCell(cell1), &
- maxLevelCell(cell2) )
- else
- maxLevelEdge(iEdge) = 0
- endif
-
- end do
-
- ! maxLevelVertex is the minimum of the surrounding cells
- VTX_LOOP: do iVtx = 1,nVertices
- maxLevelVertex(iVtx) = nVertLevels
- do i=1,vertexDegree
- cell = cellsOnVertex(i,iVtx)
- if (cell > nCells) then
- maxLevelVertex(iVtx) = 0
- cycle VTX_LOOP
- else
- maxLevelVertex(iVtx) = &
- min( maxLevelVertex(iVtx), &
- maxLevelCell(cell) )
- endif
- end do
- end do VTX_LOOP
-
! print some diagnostics
print '(10a)', 'ilevel',&
' rho ',&
@@ -228,10 +249,30 @@
enddo
enddo
- block_ptr => block_ptr % next
+ block => block % next
end do
+! --- update halos for maxLevel variables
+! mrp 100610 Integer halo updates currently cause a seg fault.
+ block => domain % blocklist
+ do while (associated(block))
+ maxLevelCell => block % mesh % maxLevelCell % array
+ maxLevelEdge => block % mesh % maxLevelEdge % array
+ maxLevelVertex => block % mesh % maxLevelVertex % array
+! call dmpar_exch_halo_field2dInteger(domain % dminfo, maxLevelCell, &
+! block % mesh % nVertLevels, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! call dmpar_exch_halo_field2dInteger(domain % dminfo, maxLevelEdge, &
+! block % mesh % nVertLevels, block % mesh % nEdges, &
+! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+! call dmpar_exch_halo_field2dInteger(domain % dminfo, maxLevelVertex, &
+! block % mesh % nVertLevels, block % mesh % nVertices, &
+! block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+
+ block => block % next
+ end do
+
end subroutine setup_sw_test_case
Modified: branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F 2010-06-10 20:40:08 UTC (rev 342)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F 2010-06-11 17:52:52 UTC (rev 343)
@@ -115,6 +115,7 @@
! BEGIN RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do rk_step = 1, 4
+
! --- update halos for diagnostic variables
block => domain % blocklist
@@ -323,12 +324,18 @@
u_src => grid % u_src % array
+
!
+ ! Initialize tendencies to zero
+ !
+ tend_h(:,:) = 0.0
+ tend_u(:,:) = 0.0
+
+ !
! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
!
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
! for explanation of divergence operator.
- tend_h(:,:) = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -376,32 +383,37 @@
!
! velocity tendency: vertical advection term -w du/dz
!
- allocate(w_dudzTopEdge(nVertLevels+1))
- w_dudzTopEdge(1) = 0.0
- w_dudzTopEdge(nVertLevels+1) = 0.0
- tend_u(:,:) = 0.0
if (config_vert_grid_type.eq.'zlevel') then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ allocate(w_dudzTopEdge(nVertLevels+1))
+ w_dudzTopEdge(1) = 0.0
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdge(iEdge)
- ! Average w from cell center to edge
- wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
+ do k=2,maxLevelEdge(iEdge)
+ ! Average w from cell center to edge
+ wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
- ! compute dudz at vertical interface with first order derivative.
- w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
- / (zMidZLevel(k-1) - zMidZLevel(k))
- end do
+ ! compute dudz at vertical interface with first order derivative.
+ w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
+ / (zMidZLevel(k-1) - zMidZLevel(k))
+ end do
+ w_dudzTopEdge(maxLevelEdge(iEdge)+1) = 0.0
- ! Average w*du/dz from vertical interface to vertical middle of cell
- do k=1,maxLevelEdge(iEdge)
- tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
- enddo
- enddo
+ ! Average w*du/dz from vertical interface to vertical middle of cell
+ do k=1,maxLevelEdge(iEdge)
+ tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
+ enddo
+ enddo
+ deallocate(w_dudzTopEdge)
endif
- deallocate(w_dudzTopEdge)
+ if (1==2) then !if (isNaN(sum(tend_u))) then
+ write(0,*) '1 Abort: NaN detected in vertical advection'
+ call dmpar_abort(dminfo)
+ endif
+
+
!
! velocity tendency: pressure gradient
!
@@ -427,6 +439,10 @@
enddo
endif
+ if (1==2) then !if (isNaN(sum(tend_u))) then
+ write(0,*) '2 Abort: NaN detected in p grad'
+ call dmpar_abort(dminfo)
+ endif
!
! velocity tendency: -q(h u^\perp) - </font>
<font color="black">abla K
! +</font>
<font color="black">u_h(</font>
<font color="black">abla \delta + {\bf k}\times </font>
<font color="gray">abla \xi)
@@ -459,6 +475,10 @@
end do
end do
+ if (1==2) then !if (isNaN(sum(tend_u))) then
+ write(0,*) '3 Abort: NaN detected in hor dif'
+ call dmpar_abort(dminfo)
+ endif
!
! velocity tendency: forcing and bottom drag
!
@@ -500,6 +520,11 @@
enddo
+ if (1==2) then !if (isNaN(sum(tend_u))) then
+ write(0,*) '4 Abort: NaN detected in forcing'
+ call dmpar_abort(dminfo)
+ endif
+
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
@@ -524,23 +549,31 @@
call dmpar_abort(dminfo)
endif
- allocate(fluxVertTop(1:nVertLevels+1))
+ allocate(fluxVertTop(nVertLevels+1))
fluxVertTop(1) = 0.0
- fluxVertTop(nVertLevels+1) = 0.0
do iEdge=1,grid % nEdgesSolve
+
do k=2,maxLevelEdge(iEdge)
fluxVertTop(k) = vertViscTop(k) &
* ( u(k-1,iEdge) - u(k,iEdge) ) &
/ (zMidEdge(k-1,iEdge) - zMidEdge(k,iEdge))
enddo
+ fluxVertTop(maxLevelEdge(iEdge)+1) = 0.0
+
do k=1,maxLevelEdge(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
+ (fluxVertTop(k) - fluxVertTop(k+1)) &
/(zTopEdge(k,iEdge) - zTopEdge(k+1,iEdge))
enddo
+
end do
deallocate(fluxVertTop, vertViscTop)
+ if (1==2) then !if (isNaN(sum(tend_u))) then
+ write(0,*) '3 Abort: NaN detected in vertical mix'
+ call dmpar_abort(dminfo)
+ endif
+
end subroutine compute_tend
@@ -613,7 +646,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !if (cell1 <= nCells .and. cell2 <= nCells) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,maxLevelEdge(iEdge)
do iTracer=1,num_tracers
tracer_edge = 0.5 * ( tracers(iTracer,k,cell1) &
@@ -624,7 +657,7 @@
tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux
end do
end do
- !end if
+ end if
end do
elseif (config_hor_tracer_adv.eq.'upwind') then
@@ -632,7 +665,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !if (cell1 <= nCells .and. cell2 <= nCells) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,maxLevelEdge(iEdge)
if (u(k,iEdge)>0.0) then
upwindCell = cell1
@@ -646,7 +679,7 @@
tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux
end do
end do
- !end if
+ end if
end do
endif
@@ -658,12 +691,16 @@
end do
end do
+ if (1==2) then !if (isNaN(sum(tend_tr))) then
+ write(0,*) '1 Abort: NaN detected tend_tr'
+ call dmpar_abort(dminfo)
+ endif
+
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
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
@@ -687,6 +724,7 @@
end do
endif
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
@@ -699,6 +737,11 @@
enddo ! iCell
deallocate(tracerTop)
+ if (1==2) then !if (isNaN(sum(tend_tr))) then
+ write(0,*) ' 2 Abort: NaN detected tend_tr'
+ call dmpar_abort(dminfo)
+ endif
+
!
! tracer tendency: horizontal tracer diffusion
! div(h \kappa_h </font>
<font color="gray">abla\phi )
@@ -709,14 +752,14 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !if (cell1 <= nCells .and. cell2 <= nCells) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
do k=1,maxLevelEdge(iEdge)
do iTracer=1,num_tracers
tr_flux(iTracer,k,iEdge) = h_edge(k,iEdge)*config_hor_diffusion * &
(Tracers(iTracer,k,cell2) - Tracers(iTracer,k,cell1))/dcEdge(iEdge)
enddo
enddo
- !endif
+ endif
enddo
! Compute the divergence, div(h \kappa_h </font>
<font color="gray">abla\phi) at cell centers
@@ -755,6 +798,11 @@
enddo
deallocate(tr_flux, tr_div)
+ if (1==2) then !if (isNaN(sum(tend_tr))) then
+ write(0,*) '3 Abort: NaN detected tend_tr'
+ call dmpar_abort(dminfo)
+ endif
+
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
@@ -781,8 +829,8 @@
allocate(fluxVertTop(num_tracers,nVertLevels+1))
fluxVertTop(:,1) = 0.0
- fluxVertTop(:,nVertLevels+1) = 0.0
do iCell=1,grid % nCellsSolve
+
do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
! compute \kappa_v d\phi/dz
@@ -791,6 +839,7 @@
/ (zMid(k-1,iCell) -zMid(k,iCell))
enddo
enddo
+ fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
do k=1,maxLevelCell(iCell)
dist = zTop(k,iCell) - zTop(k+1,iCell)
@@ -803,6 +852,10 @@
enddo ! iCell loop
deallocate(fluxVertTop, vertDiffTop)
+ if (1==2) then !if (isNaN(sum(tend_tr))) then
+ write(0,*) '4 Abort: NaN detected in tend tr vert diff'
+ call dmpar_abort(dminfo)
+ endif
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
@@ -834,6 +887,7 @@
type (grid_meta), intent(in) :: grid
+ type (dm_info) :: dminfo
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
@@ -921,20 +975,32 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
! edge with cell on both sides
- do k=1,maxLevelEdge(iEdge)
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdge(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
- do k=maxLevelEdge(iEdge)+1,max(maxLevelCell(cell1),maxLevelCell(cell2))
- if (k <= maxLevelCell(cell1)) then
+ do k=maxLevelEdge(iEdge)+1,max(maxLevelCell(cell1),maxLevelCell(cell2))
+ if (k <= maxLevelCell(cell1)) then
+ h_edge(k,iEdge) = h(k,cell1)
+ elseif (k <= maxLevelCell(cell2)) then
+ h_edge(k,iEdge) = h(k,cell2)
+ end if
+ end do
+
+ ! edge with cell on one side
+ elseif(cell1 <= nCells) then
+ do k=1,maxLevelCell(cell1)
h_edge(k,iEdge) = h(k,cell1)
- elseif (k <= maxLevelCell(cell2)) then
+ end do
+ elseif(cell2 <= nCells) then
+ do k=1,maxLevelCell(cell2)
h_edge(k,iEdge) = h(k,cell2)
- end if
- end do
+ 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
@@ -946,19 +1012,23 @@
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- 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)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ if (vertex1 <= nVertices) then
+ do k=1,maxLevelVertex(vertex1)
+ circulation(k,vertex1) = circulation(k,vertex1) &
+ - dcEdge(iEdge) * u(k,iEdge)
end do
end if
- 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)
+ if (vertex2 <= nVertices) then
+ do k=1,maxLevelVertex(vertex2)
+ circulation(k,vertex2) = circulation(k,vertex2) &
+ + dcEdge(iEdge) * u(k,iEdge)
end do
end if
end do
do iVertex=1,nVertices
- do k=1,nVertLevels
+ do k=1,maxLevelVertex(iVertex)
vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
end do
end do
@@ -972,12 +1042,12 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelCell(cell1)
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
enddo
endif
if(cell2 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelCell(cell2)
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
end if
@@ -996,7 +1066,7 @@
do iCell=1,nCells
do i=1,nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
end do
end do
@@ -1014,7 +1084,7 @@
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
if (eoe <= nEdges) then
- do k = 1,nVertLevels
+ do k = 1,maxLevelEdge(iEdge)
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
end if
@@ -1024,17 +1094,14 @@
!
! Compute ke on cell edges at velocity locations for quadratic bottom drag.
!
+ ke_edge = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
end do
- else
- do k=1,nVertLevels
- ke_edge(k,iEdge) = 0.0
- end do
end if
end do
@@ -1044,14 +1111,9 @@
!
VTX_LOOP: do iVertex = 1,nVertices
do i=1,grid % vertexDegree
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! mrp 100609 was
-! if (cellsOnVertex(i,iVertex) > nVertices) cycle VTX_LOOP
-! mrp 100609 should be
if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end do
- do k=1,nVertLevels
+ do k=1,maxLevelVertex(iVertex)
h_vertex = 0.0
do i=1,grid % vertexDegree
h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
@@ -1068,7 +1130,7 @@
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
do iEdge = 1,nEdges
- do k = 1,nVertLevels
+ do k = 1,maxLevelEdge(iEdge)
gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
dvEdge(iEdge)
enddo
@@ -1083,7 +1145,7 @@
do i=1,grid % vertexDegree
iEdge = edgesOnVertex(i,iVertex)
if(iEdge <= nEdges) then
- do k=1,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
endif
@@ -1094,7 +1156,7 @@
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
- do k = 1,nVertLevels
+ do k = 1,maxLevelEdge(iEdge)
pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
enddo
enddo
@@ -1123,7 +1185,7 @@
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
- do k = 1,nVertLevels
+ do k = 1,maxLevelEdge(iEdge)
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
enddo
@@ -1134,7 +1196,7 @@
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
- do k = 1,nVertLevels
+ do k = 1,maxLevelEdge(iEdge)
pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
enddo
enddo
@@ -1157,7 +1219,9 @@
rho(k,iCell) = 1000.0*( 1.0 &
- 2.5e-4*tracers(index_temperature,k,iCell) &
+ 7.6e-4*tracers(index_salinity,k,iCell))
+
end do
+
end do
endif
@@ -1227,6 +1291,7 @@
end do
+
! compute vertical velocity
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
@@ -1243,28 +1308,31 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+
if (cell1 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
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,nVertLevels
+ do k=1,maxLevelEdge(iEdge)
flux = u(k,iEdge) * dvEdge(iEdge)
div_u(k,cell2) = div_u(k,cell2) - flux
end do
end if
+
end do
+
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
+
end do
! Vertical velocity at bottom is zero.
- ! this next line can be set permanently somewhere upon startup.
- wTop(nVertLevels+1,iCell) = 0.0
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0
do k=maxLevelCell(iCell),1,-1
wTop(k,iCell) = wTop(k+1,iCell) - div_u(k,iCell)*h(k,iCell)
end do
@@ -1274,7 +1342,6 @@
endif
-
end subroutine compute_solve_diagnostics
Modified: branches/ocean_projects/topography_mrp/src/framework/module_dmpar.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/framework/module_dmpar.F 2010-06-10 20:40:08 UTC (rev 342)
+++ branches/ocean_projects/topography_mrp/src/framework/module_dmpar.F 2010-06-11 17:52:52 UTC (rev 343)
@@ -4,6 +4,7 @@
#ifdef _MPI
include 'mpif.h'
+ integer, parameter :: MPI_INTEGERKIND = MPI_INTEGER
#if (RKIND == 8)
integer, parameter :: MPI_REALKIND = MPI_DOUBLE_PRECISION
@@ -64,7 +65,8 @@
dminfo % nprocs = mpi_size
dminfo % my_proc_id = mpi_rank
- write(0,'(a,i5,a,i5,a)') 'task ', mpi_rank, ' of ', mpi_size, ' is running'
+ write(0,'(a,i5,a,i5,a)') 'task ', mpi_rank, ' of ', mpi_size, &
+ ' is running'
call open_streams(dminfo % my_proc_id)
@@ -781,7 +783,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
@@ -873,7 +876,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
@@ -962,7 +966,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
@@ -1054,7 +1059,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
@@ -1146,7 +1152,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(:,:,1:nNeededList) = arrayIn(:,:,1:nOwnedList)
@@ -1198,7 +1205,8 @@
n = de-ds+1
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dInteger: Not enough space in buffer to fit a single slice.'
+ write(0,*) 'packSendBuf2dInteger: Not enough space in buffer', &
+ ' to fit a single slice.'
return
end if
@@ -1217,6 +1225,45 @@
end subroutine packSendBuf2dInteger
+ subroutine packSendBuf3dInteger(d1s, d1e, d2s, d2e, nField, field, sendList, startPackIdx, nBuffer, buffer, nPacked, lastPackedIdx)
+
+ implicit none
+
+ integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startPackIdx
+ integer, dimension(d1s:d1e,d2s:d2e,*), intent(in) :: field
+ type (exchange_list), intent(in) :: sendList
+ integer, dimension(nBuffer), intent(out) :: buffer
+ integer, intent(inout) :: nPacked, lastPackedIdx
+
+ integer :: i, j, k, n
+
+ n = (d1e-d1s+1) * (d2e-d2s+1)
+
+ if (n > nBuffer) then
+ write(0,*) 'packSendBuf2dInteger: Not enough space in buffer', &
+ ' to fit a single slice.'
+ return
+ end if
+
+ nPacked = 0
+ do i=startPackIdx, sendList % nlist
+ nPacked = nPacked + n
+ if (nPacked > nBuffer) then
+ nPacked = nPacked - n
+ lastPackedIdx = i - 1
+ return
+ end if
+ k = nPacked-n+1
+ do j=d2s,d2e
+ buffer(k:k+d1e-d1s) = field(d1s:d1e,j,sendList % list(i))
+ k = k + d1e-d1s+1
+ end do
+ end do
+ lastPackedIdx = sendList % nlist
+
+ end subroutine packSendBuf3dInteger
+
+
subroutine packSendBuf1dReal(nField, field, sendList, startPackIdx, nBuffer, buffer, nPacked, lastPackedIdx)
implicit none
@@ -1259,7 +1306,8 @@
n = de-ds+1
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dReal: Not enough space in buffer to fit a single slice.'
+ write(0,*) 'packSendBuf2dInteger: Not enough space in buffer', &
+ ' to fit a single slice.'
return
end if
@@ -1293,7 +1341,8 @@
n = (d1e-d1s+1) * (d2e-d2s+1)
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dReal: Not enough space in buffer to fit a single slice.'
+ write(0,*) 'packSendBuf2dInteger: Not enough space in buffer', &
+ ' to fit a single slice.'
return
end if
@@ -1372,6 +1421,230 @@
end subroutine unpackRecvBuf2dInteger
+ subroutine unpackRecvBuf3dInteger(d1s, d1e, d2s, d2e, nField, field, recvList, startUnpackIdx, nBuffer, buffer, &
+ nUnpacked, lastUnpackedIdx)
+
+ implicit none
+
+ integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startUnpackIdx
+ integer, dimension(d1s:d1e,d2s:d2e,*), intent(inout) :: field
+ type (exchange_list), intent(in) :: recvList
+ integer, dimension(nBuffer), intent(in) :: buffer
+ integer, intent(inout) :: nUnpacked, lastUnpackedIdx
+
+ integer :: i, j, k, n
+
+ n = (d1e-d1s+1) * (d2e-d2s+1)
+
+ nUnpacked = 0
+ do i=startUnpackIdx, recvList % nlist
+ nUnpacked = nUnpacked + n
+ if (nUnpacked > nBuffer) then
+ nUnpacked = nUnpacked - n
+ lastUnpackedIdx = i - 1
+ return
+ end if
+ k = nUnpacked-n+1
+ do j=d2s,d2e
+ field(d1s:d1e,j,recvList % list(i)) = buffer(k:k+d1e-d1s)
+ k = k + d1e-d1s+1
+ end do
+ end do
+ lastUnpackedIdx = recvList % nlist
+
+ end subroutine unpackRecvBuf3dInteger
+
+
+ subroutine dmpar_exch_halo_field1dInteger(dminfo, array, dim1, sendList, recvList)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: dim1
+ integer, dimension(*), intent(inout) :: array
+ type (exchange_list), pointer :: sendList, recvList
+
+ type (exchange_list), pointer :: sendListPtr, recvListPtr
+ integer :: lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ integer :: mpi_ierr
+
+#ifdef _MPI
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ allocate(recvListPtr % ibuffer(recvListPtr % nlist))
+ call MPI_Irecv(recvListPtr % ibuffer, recvListPtr % nlist, MPI_INTEGERKIND, &
+ recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ allocate(sendListPtr % ibuffer(sendListPtr % nlist))
+ call packSendBuf1dInteger(dim1, array, sendListPtr, 1, sendListPtr % nlist, sendListPtr % ibuffer, nPacked, lastPackedIdx)
+ call MPI_Isend(sendListPtr % ibuffer, sendListPtr % nlist, MPI_INTEGERKIND, &
+ sendListPtr % procID, dminfo % my_proc_id, dminfo % comm, sendListPtr % reqID, mpi_ierr)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(recvListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ call unpackRecvBuf1dInteger(dim1, array, recvListPtr, 1, recvListPtr % nlist, recvListPtr % ibuffer, nUnpacked, lastUnpackedIdx)
+ deallocate(recvListPtr % ibuffer)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(sendListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ deallocate(sendListPtr % ibuffer)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+#endif
+
+ end subroutine dmpar_exch_halo_field1dInteger
+
+
+ subroutine dmpar_exch_halo_field2dInteger(dminfo, array, dim1, dim2, sendList, recvList)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: dim1, dim2
+ integer, dimension(dim1,*), intent(inout) :: array
+ type (exchange_list), pointer :: sendList, recvList
+
+ type (exchange_list), pointer :: sendListPtr, recvListPtr
+ integer :: lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ integer :: mpi_ierr
+ integer :: d2
+
+#ifdef _MPI
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ d2 = dim1 * recvListPtr % nlist
+ allocate(recvListPtr % ibuffer(d2))
+ call MPI_Irecv(recvListPtr % ibuffer, d2, MPI_INTEGERKIND, &
+ recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ d2 = dim1 * sendListPtr % nlist
+ allocate(sendListPtr % ibuffer(d2))
+ call packSendBuf2dInteger(1, dim1, dim2, array, sendListPtr, 1, d2, sendListPtr % ibuffer, nPacked, lastPackedIdx)
+ call MPI_Isend(sendListPtr % ibuffer, d2, MPI_INTEGERKIND, &
+ sendListPtr % procID, dminfo % my_proc_id, dminfo % comm, sendListPtr % reqID, mpi_ierr)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(recvListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ d2 = dim1 * recvListPtr % nlist
+ call unpackRecvBuf2dInteger(1, dim1, dim2, array, recvListPtr, 1, d2, recvListPtr % ibuffer, nUnpacked, lastUnpackedIdx)
+ deallocate(recvListPtr % ibuffer)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(sendListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ deallocate(sendListPtr % ibuffer)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+#endif
+
+ end subroutine dmpar_exch_halo_field2dInteger
+
+
+ subroutine dmpar_exch_halo_field3dInteger(dminfo, array, dim1, dim2, dim3, sendList, recvList)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: dim1, dim2, dim3
+ integer, dimension(dim1,dim2,*), intent(inout) :: array
+ type (exchange_list), pointer :: sendList, recvList
+
+ type (exchange_list), pointer :: sendListPtr, recvListPtr
+ integer :: lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ integer :: mpi_ierr
+ integer :: d3
+
+#ifdef _MPI
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ d3 = dim1 * dim2 * recvListPtr % nlist
+ allocate(recvListPtr % ibuffer(d3))
+ call MPI_Irecv(recvListPtr % ibuffer, d3, MPI_INTEGERKIND, &
+ recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ d3 = dim1 * dim2 * sendListPtr % nlist
+ allocate(sendListPtr % ibuffer(d3))
+ call packSendBuf3dInteger(1, dim1, 1, dim2, dim3, array, sendListPtr, 1, d3, &
+ sendListPtr % ibuffer, nPacked, lastPackedIdx)
+ call MPI_Isend(sendListPtr % ibuffer, d3, MPI_INTEGERKIND, &
+ sendListPtr % procID, dminfo % my_proc_id, dminfo % comm, sendListPtr % reqID, mpi_ierr)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(recvListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ d3 = dim1 * dim2 * recvListPtr % nlist
+ call unpackRecvBuf3dInteger(1, dim1, 1, dim2, dim3, array, recvListPtr, 1, d3, &
+ recvListPtr % ibuffer, nUnpacked, lastUnpackedIdx)
+ deallocate(recvListPtr % ibuffer)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(sendListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ deallocate(sendListPtr % ibuffer)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+#endif
+
+ end subroutine dmpar_exch_halo_field3dInteger
+
+
subroutine unpackRecvBuf1dReal(nField, field, recvList, startUnpackIdx, nBuffer, buffer, nUnpacked, lastUnpackedIdx)
implicit none
@@ -1651,5 +1924,5 @@
end subroutine dmpar_exch_halo_field3dReal
-
+
end module dmpar
</font>
</pre>