<p><b>duda</b> 2009-11-18 18:36:54 -0700 (Wed, 18 Nov 2009)</p><p>Parallelize hydrostatic model and bring other infrastructure code<br>
in line with repository trunk.<br>
<br>
M src/module_io_input.F<br>
M src/module_io_output.F<br>
M src/module_time_integration.F<br>
M src/module_dmpar.F<br>
M src/module_sw_solver.F<br>
M Registry/Registry<br>
M Makefile<br>
</p><hr noshade><pre><font color="gray">Modified: branches/hyd_model/Makefile
===================================================================
--- branches/hyd_model/Makefile 2009-11-19 01:23:33 UTC (rev 76)
+++ branches/hyd_model/Makefile 2009-11-19 01:36:54 UTC (rev 77)
@@ -28,14 +28,23 @@
#CFLAGS = -O3
#LDFLAGS = -O3
-FC = g95
-CC = gcc
-SFC = g95
-SCC = gcc
-FFLAGS = -r8 -O2 -ffree-line-length-huge
-CFLAGS = -O3
-LDFLAGS = -O3
+#FC = pgf90
+#CC = pgcc
+#SFC = pgf90
+#SCC = pgcc
+#FFLAGS = -r8 -O0 -g -Mbounds -Mchkptr -Ktrap=divz,inv,fp
+#CFLAGS = -O0 -g
+#LDFLAGS = -O0 -g -Mbounds -Mchkptr -Ktrap=divz,inv,fp
+#CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE
+FC = mpif90
+CC = mpicc
+SFC = pgf90
+SCC = pgcc
+FFLAGS = -r8 -O0 -g -Mbounds -Mchkptr
+CFLAGS = -O0 -g
+LDFLAGS = -O0 -g -Mbounds -Mchkptr
+CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE -D_MPI # -DOFFSET64BIT
#ifort (for Intel Macs and Quad-Core AMD Opteron Clusters)
#FC = mpif90
@@ -48,7 +57,6 @@
#CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE # -DOFFSET64BIT
#CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE # -DOFFSET64BIT
-CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE # -DOFFSET64BIT
CPPINCLUDES = -I../inc -I$(NETCDF)/include
FCINCLUDES = -I../inc -I$(NETCDF)/include
Modified: branches/hyd_model/Registry/Registry
===================================================================
--- branches/hyd_model/Registry/Registry 2009-11-19 01:23:33 UTC (rev 76)
+++ branches/hyd_model/Registry/Registry 2009-11-19 01:36:54 UTC (rev 77)
@@ -123,4 +123,4 @@
var real theta_old ( nVertLevels nCells ) - theta_old
var real h_edge_old ( nVertLevels nEdges ) - h_edge_old
var real h_old ( nVertLevels nCells ) - h_old
-var real pressure_old ( nVertLevelsP1 nCells ) - pressure_old
\ No newline at end of file
+var real pressure_old ( nVertLevelsP1 nCells ) - pressure_old
Modified: branches/hyd_model/src/module_dmpar.F
===================================================================
--- branches/hyd_model/src/module_dmpar.F 2009-11-19 01:23:33 UTC (rev 76)
+++ branches/hyd_model/src/module_dmpar.F 2009-11-19 01:36:54 UTC (rev 77)
@@ -1228,6 +1228,66 @@
end subroutine unpackRecvBuf3dReal
+ subroutine dmpar_exch_halo_field1dReal(dminfo, array, dim1, sendList, recvList)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: dim1
+ real (kind=RKIND), 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
+
+#ifdef _MPI
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ allocate(recvListPtr % rbuffer(recvListPtr % nlist))
+ call MPI_Irecv(recvListPtr % rbuffer, recvListPtr % nlist, MPI_REALKIND, &
+ 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 % rbuffer(sendListPtr % nlist))
+ call packSendBuf1dReal(dim1, array, sendListPtr, 1, sendListPtr % nlist, sendListPtr % rbuffer, nPacked, lastPackedIdx)
+ call MPI_Isend(sendListPtr % rbuffer, sendListPtr % nlist, MPI_REALKIND, &
+ 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 unpackRecvBuf1dReal(dim1, array, recvListPtr, 1, recvListPtr % nlist, recvListPtr % rbuffer, nUnpacked, lastUnpackedIdx)
+ deallocate(recvListPtr % rbuffer)
+ 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 % rbuffer)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+#endif
+
+ end subroutine dmpar_exch_halo_field1dReal
+
+
subroutine dmpar_exch_halo_field2dReal(dminfo, array, dim1, dim2, sendList, recvList)
implicit none
Modified: branches/hyd_model/src/module_io_input.F
===================================================================
--- branches/hyd_model/src/module_io_input.F 2009-11-19 01:23:33 UTC (rev 76)
+++ branches/hyd_model/src/module_io_input.F 2009-11-19 01:36:54 UTC (rev 77)
@@ -75,6 +75,7 @@
integer, dimension(:,:), pointer :: vertexIDSorted
integer, dimension(:), pointer :: local_cell_list, local_edge_list, local_vertex_list
+ integer, dimension(:), pointer :: local_vertlevel_list, needed_vertlevel_list
integer :: nlocal_edges, nlocal_vertices
type (exchange_list), pointer :: sendCellList, recvCellList
type (exchange_list), pointer :: sendEdgeList, recvEdgeList
@@ -411,7 +412,28 @@
indexToVertexIDField % array, local_vertex_list, &
sendVertexList, recvVertexList)
+ if (domain % dminfo % my_proc_id == 0) then
+ allocate(local_vertlevel_list(nVertLevels))
+ do i=1,nVertLevels
+ local_vertlevel_list(i) = i
+ end do
+ else
+ allocate(local_vertlevel_list(0))
+ end if
+ allocate(needed_vertlevel_list(nVertLevels))
+ do i=1,nVertLevels
+ needed_vertlevel_list(i) = i
+ end do
+ call dmpar_get_owner_list(domain % dminfo, &
+ size(local_vertlevel_list), size(needed_vertlevel_list), &
+ local_vertlevel_list, needed_vertlevel_list, &
+ sendVertLevelList, recvVertLevelList)
+
+ deallocate(local_vertlevel_list)
+ deallocate(needed_vertlevel_list)
+
+
!
! Read and distribute all fields given ownership lists and exchange lists (maybe already in block?)
!
@@ -464,9 +486,6 @@
end if
- nullify(sendVertLevelList)
- nullify(recvVertLevelList)
-
call read_and_distribute_fields(domain % dminfo, input_obj, domain % blocklist, &
readCellStart, nReadCells, readEdgeStart, nReadEdges, readVertexStart, nReadVertices, &
readVertLevelStart, nReadVertLevels, &
@@ -620,8 +639,33 @@
!
! Deallocate fields, graphs, and other memory
!
+ deallocate(indexToCellIDField % ioinfo)
+ deallocate(indexToCellIDField % array)
+ deallocate(indexToEdgeIDField % ioinfo)
+ deallocate(indexToEdgeIDField % array)
+ deallocate(indexToVertexIDField % ioinfo)
+ deallocate(indexToVertexIDField % array)
+ deallocate(cellsOnCellField % ioinfo)
+ deallocate(cellsOnCellField % array)
+ deallocate(edgesOnCellField % ioinfo)
+ deallocate(edgesOnCellField % array)
+ deallocate(verticesOnCellField % ioinfo)
+ deallocate(verticesOnCellField % array)
+ deallocate(cellsOnEdgeField % ioinfo)
+ deallocate(cellsOnEdgeField % array)
+ deallocate(cellsOnVertexField % ioinfo)
+ deallocate(cellsOnVertexField % array)
+ deallocate(cellsOnCell_0Halo)
+ deallocate(nEdgesOnCell_0Halo)
+ deallocate(indexToCellID_0Halo)
+ deallocate(cellsOnEdge_2Halo)
+ deallocate(cellsOnVertex_2Halo)
+ deallocate(edgesOnCell_2Halo)
+ deallocate(verticesOnCell_2Halo)
+ deallocate(block_graph_0Halo % vertexID)
+ deallocate(block_graph_0Halo % nAdjacent)
+ deallocate(block_graph_0Halo % adjacencyList)
-
end subroutine input_state_for_domain
Modified: branches/hyd_model/src/module_io_output.F
===================================================================
--- branches/hyd_model/src/module_io_output.F 2009-11-19 01:23:33 UTC (rev 76)
+++ branches/hyd_model/src/module_io_output.F 2009-11-19 01:36:54 UTC (rev 77)
@@ -107,6 +107,7 @@
integer, dimension(:), pointer :: neededCellList
integer, dimension(:), pointer :: neededEdgeList
integer, dimension(:), pointer :: neededVertexList
+ integer, dimension(:), pointer :: neededVertLevelList
integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, verticesOnCell, &
cellsOnEdge, verticesOnEdge, edgesOnEdge, cellsOnVertex, edgesOnVertex
integer, dimension(:,:), pointer :: cellsOnCell_save, edgesOnCell_save, verticesOnCell_save, &
@@ -193,6 +194,7 @@
allocate(neededCellList(nCellsGlobal))
allocate(neededEdgeList(nEdgesGlobal))
allocate(neededVertexList(nVerticesGlobal))
+ allocate(neededVertLevelList(nVertLevelsGlobal))
do i=1,nCellsGlobal
neededCellList(i) = i
end do
@@ -202,10 +204,14 @@
do i=1,nVerticesGlobal
neededVertexList(i) = i
end do
+ do i=1,nVertLevelsGlobal
+ neededVertLevelList(i) = i
+ end do
else
allocate(neededCellList(0))
allocate(neededEdgeList(0))
allocate(neededVertexList(0))
+ allocate(neededVertLevelList(0))
end if
if (.not. output_obj % validExchangeLists) then
@@ -224,6 +230,11 @@
domain % blocklist % mesh % indexToVertexID % array, neededVertexList, &
output_obj % sendVerticesList, output_obj % recvVerticesList)
+ call dmpar_get_owner_list(domain % dminfo, &
+ size(neededVertLevelList), size(neededVertLevelList), &
+ neededVertLevelList, neededVertLevelList, &
+ output_obj % sendVertLevelsList, output_obj % recvVertLevelsList)
+
output_obj % validExchangeLists = .true.
end if
Modified: branches/hyd_model/src/module_sw_solver.F
===================================================================
--- branches/hyd_model/src/module_sw_solver.F 2009-11-19 01:23:33 UTC (rev 76)
+++ branches/hyd_model/src/module_sw_solver.F 2009-11-19 01:36:54 UTC (rev 77)
@@ -72,9 +72,7 @@
call write_output_frame(output_obj, domain)
end if
if (mod(itimestep, config_restart_interval) == 0 .and. config_restart_interval > 0) then
- if (restart_frame == 1) then
- call output_state_init(restart_obj, domain, "RESTART")
- end if
+ if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
call output_state_for_domain(restart_obj, domain, restart_frame)
restart_frame = restart_frame + 1
end if
Modified: branches/hyd_model/src/module_time_integration.F
===================================================================
--- branches/hyd_model/src/module_time_integration.F 2009-11-19 01:23:33 UTC (rev 76)
+++ branches/hyd_model/src/module_time_integration.F 2009-11-19 01:36:54 UTC (rev 77)
@@ -109,8 +109,34 @@
do rk_step = 1, 3
- ! --- compute tendencies for dynamics
+ if(debug) write(0,*) ' rk substep ', rk_step
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % qtot % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % cqu % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % h % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pressure % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % geopotential % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % alpha % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ block => block % next
+ end do
+
if(debug) write(0,*) ' rk substep ', rk_step
block => domain % blocklist
@@ -121,6 +147,21 @@
if(debug) write(0,*) ' returned from dyn_tend '
+ !
+ ! --- update halos for tendencies
+ !
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
+
! --- advance over sub_steps
block => domain % blocklist
@@ -133,12 +174,12 @@
do small_step = 1, number_sub_steps(rk_step)
- if(debug) write(0,*) ' small step ',small_step
+ if(debug) write(0,*) ' small step ',small_step
block => domain % blocklist
do while (associated(block))
call advance_dynamics( block % intermediate_step(TEND), block % time_levs(2) % state, &
- block % time_levs(1) % state, block % mesh, &
+ block % mesh, &
small_step, number_sub_steps(rk_step), rk_sub_timestep(rk_step) )
block => block % next
end do
@@ -146,6 +187,55 @@
if(debug) write(0,*) ' dynamics advance complete '
! will need communications here?
+ !
+ ! --- update halos for prognostic variables
+ !
+ block => domain % blocklist
+ do while (associated(block))
+!! call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % u % array(:,:), &
+!! block % mesh % nVertLevels, block % mesh % nEdges, &
+!! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % h_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+!! call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % uhAvg % array(:,:), &
+!! block % mesh % nVertLevels, block % mesh % nEdges, &
+!! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+!! call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % wwAvg % array(:,:), &
+!! block % mesh % nVertLevels+1, block % mesh % nCells, &
+!! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % theta % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+!! call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % h % array(:,:), &
+!! block % mesh % nVertLevels, block % mesh % nCells, &
+!! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+!! call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(TEND) % h % array(:,:), &
+!! block % mesh % nVertLevels, block % mesh % nCells, &
+!! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field1dReal(domain % dminfo, block % mesh % dpsdt % array(:), &
+ block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field1dReal(domain % dminfo, block % time_levs(2) % state % surface_pressure % array(:), &
+ block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % alpha % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % ww % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pressure % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+!! call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % pressure_old % array(:,:), &
+!! block % mesh % nVertLevels+1, block % mesh % nCells, &
+!! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % geopotential % array(:,:), &
+ block % mesh % nVertLevels+1, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
end do
@@ -161,6 +251,18 @@
block % mesh, rk_timestep(rk_step) )
block => block % next
end do
+
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % intermediate_step(TEND) % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % time_levs(2) % state % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
if(debug) write(0,*) ' advance scalars complete '
@@ -184,35 +286,34 @@
! END RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- if(debug) write(0,*) ' rk step complete - mass diagnostics '
+ if(debug) write(0,*) ' rk step complete - mass diagnostics '
- if(debug) then
-
- domain_mass = 0.
- tracer_mass = 0.
- block => domain % blocklist
- tracer_min = block % time_levs(2) % state % tracers % array (2,1,1)
- tracer_max = block % time_levs(2) % state % tracers % array (2,1,1)
- do while(associated(block))
- do iCell = 1, block % mesh % nCells
- domain_mass = domain_mass + block % time_levs(2) % state % surface_pressure % array (iCell) * &
- block % mesh % areaCell % array (iCell) &
- - block % time_levs(2) % state % pressure % array (block % mesh % nVertLevels + 1, 1) * &
- block % mesh % areaCell % array (iCell)
- do k=1, block % mesh % nVertLevels
- tracer_mass = tracer_mass - block % time_levs(2) % state % tracers % array (2,k,iCell) * &
- block % time_levs(2) % state % h % array (k,iCell) * &
- block % mesh % dnw % array (k) * &
- block % mesh % areaCell % array (iCell)
- tracer_min = min(tracer_min,block % time_levs(2) % state % tracers % array (2,k,iCell))
- tracer_max = max(tracer_max,block % time_levs(2) % state % tracers % array (2,k,iCell))
- enddo
- enddo
- block => block % next
- enddo
- write(0,*) ' mass in the domain = ',domain_mass
- write(0,*) ' tracer mass in the domain = ',tracer_mass
- write(0,*) ' tracer_min, tracer_max ',tracer_min, tracer_max
+ if(debug) then
+ domain_mass = 0.
+ tracer_mass = 0.
+ block => domain % blocklist
+ tracer_min = block % time_levs(2) % state % tracers % array (2,1,1)
+ tracer_max = block % time_levs(2) % state % tracers % array (2,1,1)
+ do while(associated(block))
+ do iCell = 1, block % mesh % nCells
+ domain_mass = domain_mass + block % time_levs(2) % state % surface_pressure % array (iCell) * &
+ block % mesh % areaCell % array (iCell) &
+ - block % time_levs(2) % state % pressure % array (block % mesh % nVertLevels + 1, 1) * &
+ block % mesh % areaCell % array (iCell)
+ do k=1, block % mesh % nVertLevels
+ tracer_mass = tracer_mass - block % time_levs(2) % state % tracers % array (2,k,iCell) * &
+ block % time_levs(2) % state % h % array (k,iCell) * &
+ block % mesh % dnw % array (k) * &
+ block % mesh % areaCell % array (iCell)
+ tracer_min = min(tracer_min,block % time_levs(2) % state % tracers % array (2,k,iCell))
+ tracer_max = max(tracer_max,block % time_levs(2) % state % tracers % array (2,k,iCell))
+ enddo
+ enddo
+ block => block % next
+ enddo
+ write(0,*) ' mass in the domain = ',domain_mass
+ write(0,*) ' tracer mass in the domain = ',tracer_mass
+ write(0,*) ' tracer_min, tracer_max ',tracer_min, tracer_max
endif
@@ -264,7 +365,9 @@
do k = 1, nVertLevels
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- grid % cqu % array(k,iEdge) = 1./( 1. + 0.5*(grid % qtot % array(k,cell1)+grid % qtot % array(k,cell2)) )
+ if (cell1 > 0 .and. cell2 > 0) then
+ grid % cqu % array(k,iEdge) = 1./( 1. + 0.5*(grid % qtot % array(k,cell1)+grid % qtot % array(k,cell2)) )
+ end if
enddo
enddo
@@ -501,16 +604,18 @@
!
! vertical advection for u
!
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
wdun(1) = 0.
+ if (cell1 > 0 .and. cell2 > 0) then
do k=2,nVertLevels
wdun(k) = &
(ww(k,cell1)/(h(k,cell1)+h(k-1,cell1)) +ww(k,cell2)/(h(k,cell2)+h(k-1,cell2)))* &
rdnu(k)*(u(k,iEdge)-u(k-1,iEdge))
end do
+ end if
wdun(nVertLevels+1) = 0.
do k=1,nVertLevels
@@ -524,11 +629,12 @@
! ---- vertical mixing - 2nd order
!
- do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
do k=2,nVertLevels-1
z1 = 0.5*(geopotential(k-1,cell1)+geopotential(k-1,cell2))/gravity
@@ -544,6 +650,7 @@
(u(k+1,iEdge)-u(k ,iEdge))/(zp-z0) &
-(u(k ,iEdge)-u(k-1,iEdge))/(z0-zm) )/(0.5*(zp-zm))
end do
+ end if
end do
!----------- rhs for theta
@@ -560,6 +667,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
do k=1,grid % nVertLevels
theta_turb_flux = hnu*prandtl*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
flux = dvEdge (iEdge) * h_edge(k,iEdge) * theta_turb_flux
@@ -570,6 +678,7 @@
! delsq_theta(k,cell2) = delsq_theta(k,iCell) - dvEdge(iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
end do
+ end if
end do
!
@@ -579,7 +688,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
-! if (cell1 > 0 .and. cell2 > 0) then
+ if (cell1 > 0 .and. cell2 > 0) then
do k=1,grid % nVertLevels
theta_edge = 0.5 * (theta(k,cell1) + theta(k,cell2))
@@ -593,7 +702,7 @@
tend_theta(k,cell1) = tend_theta(k,cell1) - flux
tend_theta(k,cell2) = tend_theta(k,cell2) + flux
end do
-! end if
+ end if
end do
!
@@ -638,7 +747,7 @@
!---------------------------------------------------------------------------------------------------------
- subroutine advance_dynamics(tend, s, s_old, grid, small_step, number_small_steps, dt)
+ subroutine advance_dynamics(tend, s, grid, small_step, number_small_steps, dt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance the dry dynamics a small timestep (forward-backward integration)
!
@@ -653,7 +762,6 @@
implicit none
type (grid_state), intent(in) :: tend
- type (grid_state), intent(in) :: s_old
type (grid_state), intent(inout) :: s
type (grid_meta), intent(in) :: grid
real (kind=RKIND), intent(in) :: dt
@@ -750,7 +858,7 @@
if(small_step == 1) then
- do iCell=1,grid % nCellsSolve
+ do iCell=1,grid % nCells
do k=1,nVertLevels
theta(k,iCell) = theta(k,iCell)*h(k,iCell)
end do
@@ -770,16 +878,18 @@
do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
- u(k,iEdge) = u(k,iEdge) + dt*tend_u(k,iEdge) &
- -(0.5*dt/dcEdge(iEdge))*( &
- (geopotential(k+1,cell2)-geopotential(k+1,cell1)) &
- +(geopotential(k ,cell2)-geopotential(k ,cell1)) &
- +cqu(k,iEdge)*(alpha(k,cell2)+alpha(k,cell1))* &
- 0.5*(pressure(k+1,cell2)-pressure(k+1,cell1) &
- +pressure(k ,cell2)-pressure(k ,cell1))) &
- -smext*dcEdge(iEdge)*(dpsdt(cell2)-dpsdt(cell1))/h_edge(k,iEdge)
- end do
+ if (cell1 > 0 .and. cell2 > 0) then
+ do k=1,nVertLevels
+ u(k,iEdge) = u(k,iEdge) + dt*tend_u(k,iEdge) &
+ -(0.5*dt/dcEdge(iEdge))*( &
+ (geopotential(k+1,cell2)-geopotential(k+1,cell1)) &
+ +(geopotential(k ,cell2)-geopotential(k ,cell1)) &
+ +cqu(k,iEdge)*(alpha(k,cell2)+alpha(k,cell1))* &
+ 0.5*(pressure(k+1,cell2)-pressure(k+1,cell1) &
+ +pressure(k ,cell2)-pressure(k ,cell1))) &
+ -smext*dcEdge(iEdge)*(dpsdt(cell2)-dpsdt(cell1))/h_edge(k,iEdge)
+ end do
+ end if
end do
!
@@ -790,15 +900,19 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
do k=1,nVertLevels
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell1) = tend_h(k,cell1) + flux
tend_h(k,cell2) = tend_h(k,cell2) - flux
- uhAvg(k,iEdge) = uhAvg(k,iEdge) + u(k,iEdge) * h_edge(k,iEdge)
end do
+ end if
+ do k=1,nVertLevels
+ uhAvg(k,iEdge) = uhAvg(k,iEdge) + u(k,iEdge) * h_edge(k,iEdge)
+ end do
end do
- do iCell=1, grid % nCellsSolve
+ do iCell=1, grid % nCells
dpsdt(iCell) = 0.
do k=1,nVertLevels
dpsdt(iCell) = dpsdt(iCell) + dnw(k)*tend_h(k,iCell)
@@ -841,6 +955,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
do k=1,nVertLevels
h_edge(k,iEdge) = 0.5*(h(k,cell1)+h(k,cell2)) ! here is update of h_edge
he_old = 0.5*(h_old(k,cell1)+h_old(k,cell2))
@@ -849,11 +964,13 @@
theta(k,cell1) = theta(k,cell1) - dt*flux/areaCell(cell1)
theta(k,cell2) = theta(k,cell2) + dt*flux/areaCell(cell2)
end do
+ end if
end do
+
! compute some diagnostics using the new state
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, grid % nCells
do k = nVertLevels,1,-1
pressure(k,iCell) = pressure(k+1,iCell) - dnw(k)*(1.+qtot(k,iCell))*h(k,iCell)
@@ -879,7 +996,7 @@
! if last small step of a set, decouple theta
if(small_step == number_small_steps) then
- do iCell=1,grid % nCellsSolve
+ do iCell=1,grid % nCells
do k=1,nVertLevels
theta(k,iCell) = theta(k,iCell)/h(k,iCell)
end do
@@ -963,7 +1080,7 @@
! vertical flux divergence
!
- do iCell=1,grid % nCellsSolve
+ do iCell=1,grid % nCells
wdtn(:,1) = 0.
do k = 2, nVertLevels
</font>
</pre>