<p><b>duda</b> 2009-08-12 12:33:24 -0600 (Wed, 12 Aug 2009)</p><p>1) Add halo update code for 3d fields (like tracers) in module_dmpar<br>
2) Make corrections to scalar transport in module_time_integration<br>
3) Initialize test tracer fields to more realistic range <br>
in module_test_cases<br>
<br>
M module_test_cases.F<br>
M module_time_integration.F<br>
M module_dmpar.F<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/module_dmpar.F
===================================================================
--- trunk/swmodel/module_dmpar.F 2009-08-11 18:57:40 UTC (rev 23)
+++ trunk/swmodel/module_dmpar.F 2009-08-12 18:33:24 UTC (rev 24)
@@ -1661,5 +1661,143 @@
end subroutine dmpar_exch_halo_field2dReal
+
+ subroutine dmpar_exch_halo_field3dReal(dminfo, array, dim1, dim2, dim3, sendList, recvList)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: dim1, dim2, dim3
+ real (kind=RKIND), dimension(dim1,dim2, dim3), intent(inout) :: array
+ type (exchange_list), pointer :: sendList, recvList
+
+ type (exchange_list), pointer :: sendListPtr, recvListPtr
+ integer :: i, j, istart, nLoop, leftNeighbor, rightNeighbor
+ integer :: nSendBuf, nRecvBuf, lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ real (kind=RKIND), allocatable, dimension(:) :: buffer
+ logical :: sendFirst
+ integer :: mpi_ierr
+ integer :: d1, d2, d3
+
+
+ nLoop = dminfo % nprocs / 2
+
+ d1 = dim1
+ d2 = dim2
+ d3 = dim3
+
+ allocate(buffer(BUFSIZE))
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID == dminfo % my_proc_id) exit
+ sendListPtr => sendListPtr % next
+ end do
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID == dminfo % my_proc_id) exit
+ recvListPtr => recvListPtr % next
+ end do
+
+ do i=1,nLoop
+ if (mod(dminfo % my_proc_id / i, 2) /= 0) then
+ sendFirst = .true.
+ else
+ sendFirst = .false.
+ end if
+
+ leftNeighbor = dminfo % my_proc_id - i + dminfo % nprocs
+ leftNeighbor = mod(leftNeighbor, dminfo % nprocs)
+ rightNeighbor = dminfo % my_proc_id + i
+ rightNeighbor = mod(rightNeighbor, dminfo % nprocs)
+
+ if (sendFirst) then
+ if (isInList(leftNeighbor, sendList, sendListPtr)) then
+ istart = 1
+ do while (istart <= sendListPtr % nlist)
+ call packSendBuf3dReal(1, d1, 1, d2, d3, array, sendListPtr, istart, BUFSIZE, buffer, nSendBuf, lastPackedIdx)
+ call MPI_Send(nSendBuf, 1, MPI_INTEGER, leftNeighbor, i, dminfo % comm, mpi_ierr)
+ call MPI_Send(buffer, nSendBuf, MPI_REALKIND, leftNeighbor, i, dminfo % comm, mpi_ierr)
+ istart = lastPackedIdx + 1
+ end do
+ end if
+ if (isInList(leftNeighbor, recvList, recvListPtr)) then
+ istart = 1
+ do while (istart <= recvListPtr % nList)
+ call MPI_Recv(nRecvBuf, 1, MPI_INTEGER, leftNeighbor, i, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call MPI_Recv(buffer, nRecvBuf, MPI_REALKIND, leftNeighbor, i, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call unpackRecvBuf3dReal(1, d1, 1, d2, d3, array, recvListPtr, istart, BUFSIZE, buffer, nUnpacked, lastUnpackedIdx)
+ istart = lastUnpackedIdx + 1
+ end do
+ end if
+ if (leftNeighbor /= rightNeighbor) then
+ if (isInList(rightNeighbor, recvList, recvListPtr)) then
+ istart = 1
+ do while (istart <= recvListPtr % nList)
+ call MPI_Recv(nRecvBuf, 1, MPI_INTEGER, rightNeighbor, i, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call MPI_Recv(buffer, nRecvBuf, MPI_REALKIND, rightNeighbor, i, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call unpackRecvBuf3dReal(1, d1, 1, d2, d3, array, recvListPtr, istart, BUFSIZE, buffer, nUnpacked, lastUnpackedIdx)
+ istart = lastUnpackedIdx + 1
+ end do
+ end if
+ if (isInList(rightNeighbor, sendList, sendListPtr)) then
+ istart = 1
+ do while (istart <= sendListPtr % nlist)
+ call packSendBuf3dReal(1, d1, 1, d2, d3, array, sendListPtr, istart, BUFSIZE, buffer, nSendBuf, lastPackedIdx)
+ call MPI_Send(nSendBuf, 1, MPI_INTEGER, rightNeighbor, i, dminfo % comm, mpi_ierr)
+ call MPI_Send(buffer, nSendBuf, MPI_REALKIND, rightNeighbor, i, dminfo % comm, mpi_ierr)
+ istart = lastPackedIdx + 1
+ end do
+ end if
+ end if
+
+ else
+ if (isInList(rightNeighbor, recvList, recvListPtr)) then
+ istart = 1
+ do while (istart <= recvListPtr % nList)
+ call MPI_Recv(nRecvBuf, 1, MPI_INTEGER, rightNeighbor, i, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call MPI_Recv(buffer, nRecvBuf, MPI_REALKIND, rightNeighbor, i, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call unpackRecvBuf3dReal(1, d1, 1, d2, d3, array, recvListPtr, istart, BUFSIZE, buffer, nUnpacked, lastUnpackedIdx)
+ istart = lastUnpackedIdx + 1
+ end do
+ end if
+ if (isInList(rightNeighbor, sendList, sendListPtr)) then
+ istart = 1
+ do while (istart <= sendListPtr % nlist)
+ call packSendBuf3dReal(1, d1, 1, d2, d3, array, sendListPtr, istart, BUFSIZE, buffer, nSendBuf, lastPackedIdx)
+ call MPI_Send(nSendBuf, 1, MPI_INTEGER, rightNeighbor, i, dminfo % comm, mpi_ierr)
+ call MPI_Send(buffer, nSendBuf, MPI_REALKIND, rightNeighbor, i, dminfo % comm, mpi_ierr)
+ istart = lastPackedIdx + 1
+ end do
+ end if
+ if (leftNeighbor /= rightNeighbor) then
+ if (isInList(leftNeighbor, sendList, sendListPtr)) then
+ istart = 1
+ do while (istart <= sendListPtr % nlist)
+ call packSendBuf3dReal(1, d1, 1, d2, d3, array, sendListPtr, istart, BUFSIZE, buffer, nSendBuf, lastPackedIdx)
+ call MPI_Send(nSendBuf, 1, MPI_INTEGER, leftNeighbor, i, dminfo % comm, mpi_ierr)
+ call MPI_Send(buffer, nSendBuf, MPI_REALKIND, leftNeighbor, i, dminfo % comm, mpi_ierr)
+ istart = lastPackedIdx + 1
+ end do
+ end if
+ if (isInList(leftNeighbor, recvList, recvListPtr)) then
+ istart = 1
+ do while (istart <= recvListPtr % nList)
+ call MPI_Recv(nRecvBuf, 1, MPI_INTEGER, leftNeighbor, i, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call MPI_Recv(buffer, nRecvBuf, MPI_REALKIND, leftNeighbor, i, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call unpackRecvBuf3dReal(1, d1, 1, d2, d3, array, recvListPtr, istart, BUFSIZE, buffer, nUnpacked, lastUnpackedIdx)
+ istart = lastUnpackedIdx + 1
+ end do
+ end if
+ end if
+
+ end if
+ end do
+
+ deallocate(buffer)
+
+ end subroutine dmpar_exch_halo_field3dReal
+
end module dmpar
Modified: trunk/swmodel/module_test_cases.F
===================================================================
--- trunk/swmodel/module_test_cases.F 2009-08-11 18:57:40 UTC (rev 23)
+++ trunk/swmodel/module_test_cases.F 2009-08-12 18:33:24 UTC (rev 24)
@@ -348,14 +348,14 @@
!
do iCell=1,grid % nCells
r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- state % tracers % array(1,1,iCell) = hs0 * (1.0 - r/rr)
+ state % tracers % array(1,1,iCell) = 1.0 - r/rr
end do
do iCell=1,grid % nCells
r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + &
(grid % latCell % array(iCell) - theta_c - pii/6.0)**2.0 &
) &
)
- state % tracers % array(2,1,iCell) = hs0 * (1.0 - r/rr)
+ state % tracers % array(2,1,iCell) = 1.0 - r/rr
end do
!
Modified: trunk/swmodel/module_time_integration.F
===================================================================
--- trunk/swmodel/module_time_integration.F 2009-08-11 18:57:40 UTC (rev 23)
+++ trunk/swmodel/module_time_integration.F 2009-08-12 18:33:24 UTC (rev 24)
@@ -59,6 +59,7 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
+ integer :: iCell, k
type (block_type), pointer :: block
integer, parameter :: TEMP = 1
@@ -94,6 +95,9 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q1) % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(Q1) % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
@@ -105,8 +109,15 @@
block % intermediate_step(Q1) % u % array(:,:)/2.0
block % intermediate_step(TEMP) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) + &
block % intermediate_step(Q1) % h % array(:,:)/2.0
- block % intermediate_step(TEMP) % tracers % array(:,:,:) = block % time_levs(1) % state % tracers % array(:,:,:) + &
- block % intermediate_step(Q1) % tracers % array(:,:,:)/2.0
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % intermediate_step(TEMP) % tracers % array(:,k,iCell) = ( &
+ block % time_levs(1) % state % h % array(k,iCell) * &
+ block % time_levs(1) % state % tracers % array(:,k,iCell) + &
+ block % intermediate_step(Q1) % tracers % array(:,k,iCell) / 2.0 &
+ ) / block % intermediate_step(TEMP) % h % array(k,iCell)
+ end do
+ end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
@@ -127,6 +138,9 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q2) % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(Q2) % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
@@ -138,8 +152,15 @@
block % intermediate_step(Q2) % u % array(:,:)/2.0
block % intermediate_step(TEMP) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) + &
block % intermediate_step(Q2) % h % array(:,:)/2.0
- block % intermediate_step(TEMP) % tracers % array(:,:,:) = block % time_levs(1) % state % tracers % array(:,:,:) + &
- block % intermediate_step(Q2) % tracers % array(:,:,:)/2.0
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % intermediate_step(TEMP) % tracers % array(:,k,iCell) = ( &
+ block % time_levs(1) % state % h % array(k,iCell) * &
+ block % time_levs(1) % state % tracers % array(:,k,iCell) + &
+ block % intermediate_step(Q2) % tracers % array(:,k,iCell) / 2.0 &
+ ) / block % intermediate_step(TEMP) % h % array(k,iCell)
+ end do
+ end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
@@ -160,6 +181,9 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q3) % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(Q3) % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
@@ -171,8 +195,15 @@
block % intermediate_step(Q3) % u % array(:,:)
block % intermediate_step(TEMP) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) + &
block % intermediate_step(Q3) % h % array(:,:)
- block % intermediate_step(TEMP) % tracers % array(:,:,:) = block % time_levs(1) % state % tracers % array(:,:,:) + &
- block % intermediate_step(Q3) % tracers % array(:,:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % intermediate_step(TEMP) % tracers % array(:,k,iCell) = ( &
+ block % time_levs(1) % state % h % array(k,iCell) * &
+ block % time_levs(1) % state % tracers % array(:,k,iCell) + &
+ block % intermediate_step(Q3) % tracers % array(:,k,iCell) / 2.0 &
+ ) / block % intermediate_step(TEMP) % h % array(k,iCell)
+ end do
+ end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % intermediate_step(TEMP) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
@@ -193,6 +224,9 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(Q4) % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(Q4) % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
@@ -212,12 +246,20 @@
2.0 * block % intermediate_step(Q3) % h % array(:,:) + &
block % intermediate_step(Q4) % h % array(:,:) &
) / 6.0
- block % time_levs(2) % state % tracers % array(:,:,:) = block % time_levs(1) % state % tracers % array(:,:,:) + &
- ( block % intermediate_step(Q1) % tracers % array(:,:,:) + &
- 2.0 * block % intermediate_step(Q2) % tracers % array(:,:,:) + &
- 2.0 * block % intermediate_step(Q3) % tracers % array(:,:,:) + &
- block % intermediate_step(Q4) % tracers % array(:,:,:) &
- ) / 6.0
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % time_levs(2) % state % tracers % array(:,k,iCell) = ( &
+ block % time_levs(1) % state % h % array(k,iCell) * &
+ block % time_levs(1) % state % tracers % array(:,k,iCell) + &
+ ( block % intermediate_step(Q1) % tracers % array(:,k,iCell) + &
+ 2.0 * block % intermediate_step(Q2) % tracers % array(:,k,iCell) + &
+ 2.0 * block % intermediate_step(Q3) % tracers % array(:,k,iCell) + &
+ block % intermediate_step(Q4) % tracers % array(:,k,iCell) &
+ ) / 6.0 &
+ ) / block % time_levs(2) % state % h % array(k,iCell)
+ end do
+ end do
+
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
end if
@@ -483,23 +525,25 @@
real (kind=RKIND) :: flux, tracer_edge
tend % tracers % array(:,:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- do k=1,grid % nVertLevelsSolve
- do iTracer=1,grid % nTracers
- tracer_edge = 0.5*(s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
- flux = s % u % array(k,iEdge) * grid % dvEdge % array(iEdge) * tend % h_edge % array(k,iEdge) * tracer_edge
- tend % tracers % array(iTracer,k,cell1) = tend % tracers % array(iTracer,k,cell1) - flux
- tend % tracers % array(iTracer,k,cell2) = tend % tracers % array(iTracer,k,cell2) + flux
- end do
- end do
- end do
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) 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))
+ flux = s % u % array(k,iEdge) * grid % dvEdge % array(iEdge) * tend % h_edge % array(k,iEdge) * tracer_edge
+ tend % tracers % array(iTracer,k,cell1) = tend % tracers % array(iTracer,k,cell1) - flux
+ tend % tracers % array(iTracer,k,cell2) = tend % tracers % array(iTracer,k,cell2) + flux
+ end do
+ end do
+ end if
+ end do
+
do iCell=1,grid % nCellsSolve
do k=1,grid % nVertLevelsSolve
do iTracer=1,grid % nTracers
- tend % tracers % array(iTracer,k,iCell) = tend % tracers % array(iTracer,k,iCell) / &
- (grid % areaCell % array(iCell) * s % h % array(k,iCell))
+ tend % tracers % array(iTracer,k,iCell) = tend % tracers % array(iTracer,k,iCell) / grid % areaCell % array(iCell)
end do
end do
end do
</font>
</pre>