<p><b>mpetersen@lanl.gov</b> 2010-10-27 08:59:50 -0600 (Wed, 27 Oct 2010)</p><p>BRANCH COMMIT<br>
<br>
Change loops like<br>
do k = 1,nVertLevels<br>
to<br>
do k=1,maxLevelEdgeTop(iEdge)<br>
<br>
where I have added new variables:<br>
maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells<br>
maxLevelEdgeBot is the maximum (deepest) of the surrounding cells<br>
maxLevelVertexTop is the minimum (shallowest) of the surrounding cells<br>
maxLevelVertexBot is the maximum (deepest) of the surrounding cells<br>
<br>
Also added Jacket and McDougal EOS and merged module_dmpar.F from the trunk.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography2_mrp/Makefile
===================================================================
--- branches/ocean_projects/topography2_mrp/Makefile 2010-10-25 18:34:39 UTC (rev 581)
+++ branches/ocean_projects/topography2_mrp/Makefile 2010-10-27 14:59:50 UTC (rev 582)
@@ -80,6 +80,7 @@
"CORE = $(CORE)" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+# mrp change from -O3 to -check -g
ifort:
( make all \
"FC = mpif90" \
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-10-25 18:34:39 UTC (rev 581)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-10-27 14:59:50 UTC (rev 582)
@@ -120,8 +120,10 @@
# Arrays for z-level version of mpas-ocean
var persistent integer maxLevelCell ( nCells ) 0 iro maxLevelCell mesh - -
-var persistent integer maxLevelEdge ( nEdges ) 0 - maxLevelEdge mesh - -
-var persistent integer maxLevelVertex ( nVertices ) 0 - maxLevelVertex mesh - -
+var persistent integer maxLevelEdgeTop ( nEdges ) 0 - maxLevelEdgeTop mesh - -
+var persistent integer maxLevelEdgeBot ( nEdges ) 0 - maxLevelEdgeBot mesh - -
+var persistent integer maxLevelVertexTop ( nVertices ) 0 - maxLevelVertexTop mesh - -
+var persistent integer maxLevelVertexBot ( nVertices ) 0 - maxLevelVertexBot mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
var persistent real zMidZLevel ( nVertLevels ) 0 - zMidZLevel mesh - -
var persistent real zTopZLevel ( nVertLevelsP1 ) 0 - zTopZLevel mesh - -
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F 2010-10-25 18:34:39 UTC (rev 581)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F 2010-10-27 14:59:50 UTC (rev 582)
@@ -120,6 +120,11 @@
nVertices = block_ptr % mesh % nVertices
nVertLevels = block_ptr % mesh % nVertLevels
+!mrp temp, to agree with previous test
+h(1,:) = 500.0
+h(2,:) = 1250.0
+h(3,:) = 3250.0
+
pi=3.1415
! Tracer blob in Central Pacific, away from boundaries:
!latCenter=pi/16; lonCenter=9./8.*pi
@@ -185,11 +190,6 @@
endif
-!mrp temp, to agree with previous test
-!h(1,:) = 500.0
-!h(2,:) = 1250.0
-!h(3,:) = 3250.0
-
! print some diagnostics
print '(10a)', 'ilevel',&
' rho ',&
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-10-25 18:34:39 UTC (rev 581)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-10-27 14:59:50 UTC (rev 582)
@@ -264,9 +264,10 @@
type (state_type), intent(in) :: s
type (mesh_type), intent(in) :: grid
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j, kmax
- integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
real (kind=RKIND) :: h_mom_eddy_visc2, h_mom_eddy_visc4
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
upstream_bias, wTopEdge, rho0Inv, r
@@ -280,7 +281,7 @@
type (dm_info) :: dminfo
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdge, maxLevelVertex
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: &
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
edgesOnEdge, edgesOnVertex
@@ -331,14 +332,15 @@
zMidZLevel => grid % zMidZLevel % array
zTopZLevel => grid % zTopZLevel % array
maxLevelCell => grid % maxLevelCell % array
- maxLevelEdge => grid % maxLevelEdge % array
- maxLevelVertex => grid % maxLevelVertex % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
tend_h => tend % h % array
tend_u => tend % u % array
nCells = grid % nCells
nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
@@ -348,29 +350,40 @@
h_mom_eddy_visc4 = config_h_mom_eddy_visc4
!
+ ! height tendency: start accumulating tendency terms
+ !
+ tend_h = 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
+ !
+ ! for z-level, only compute height tendency for top layer.
+ if (config_vert_grid_type.eq.'isopycnal') then
+ kmax = nVertLevels
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ kmax = 1
+ endif
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells) then
- do k=1,nVertLevels
+ do k=1,kmax
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 <= nCells) then
- do k=1,nVertLevels
+ do k=1,kmax
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
tend_h(k,cell2) = tend_h(k,cell2) + flux
end do
end if
end do
do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
+ do k=1,kmax
tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
end do
end do
@@ -378,22 +391,10 @@
!
! height tendency: vertical advection term -d/dz(hw)
!
+ ! Vertical advection computed for top layer of a z grid only.
if (config_vert_grid_type.eq.'zlevel') then
-
do iCell=1,nCells
-
tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
-
- ! This next loop is to verify that h for levels below the first
- ! remain constant. At a later time this could be replaced
- ! by just tend_h(2:nVertLevels,:) = 0.0, and then there is
- ! no need to compute the horizontal tend_h term for k=2:nVertLevels
- ! on a zlevel grid, above.
- do k=2,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) &
- - (wTop(k,iCell) - wTop(k+1,iCell))
- end do
-
end do
endif ! coordinate type
@@ -405,30 +406,30 @@
!
! velocity tendency: vertical advection term -w du/dz
!
- allocate(w_dudzTopEdge(nVertLevels+1))
- w_dudzTopEdge(1) = 0.0
- w_dudzTopEdge(nVertLevels+1) = 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,nVertLevels
- ! Average w from cell center to edge
- wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
+ do k=2,maxLevelEdgeTop(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(maxLevelEdgeTop(iEdge)+1) = 0.0
- ! Average w*du/dz from vertical interface to vertical middle of cell
- do k=1,nVertLevels
- 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,maxLevelEdgeTop(iEdge)
+ tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
+ enddo
+ enddo
+ deallocate(w_dudzTopEdge)
endif
- deallocate(w_dudzTopEdge)
!
! velocity tendency: pressure gradient
@@ -438,7 +439,7 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
end do
@@ -447,7 +448,7 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- rho0Inv*( pZLevel(k,cell2) &
- pZLevel(k,cell1) )/dcEdge(iEdge)
@@ -467,7 +468,7 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
! Here -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
! is - </font>
<font color="gray">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
@@ -505,7 +506,7 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
@@ -520,7 +521,7 @@
do iEdge=1,nEdges
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
- dcEdge(iEdge) * delsq_u(k,iEdge)
delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
@@ -529,7 +530,7 @@
end do
do iVertex=1,nVertices
r = 1.0 / areaTriangle(iVertex)
- do k=1,nVertLevels
+ do k=1,maxLevelVertexBot(iVertex)
delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
end do
end do
@@ -539,7 +540,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+ delsq_u(k,iEdge)*dvEdge(iEdge)
delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
@@ -561,7 +562,7 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
u_diffusion = ( delsq_divergence(k,cell2) &
- delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
@@ -586,7 +587,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
@@ -606,22 +607,34 @@
!
do iEdge=1,grid % nEdgesSolve
- ! forcing in top layer only
- tend_u(1,iEdge) = tend_u(1,iEdge) &
+ ! forcing in top layer only
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+ k = maxLevelEdgeTop(iEdge)
+
+ ! efficiency note: it would be nice to avoid this
+ ! if within a do. This could be done with
+ ! k = max(maxLevelEdgeTop(iEdge),1)
+ ! and then tend_u(1,iEdge) is just not used for land cells.
+
+ if (k>0) then
+
! bottom drag is the same as POP:
! -c |u| u where c is unitless and 1.0e-3.
! see POP Reference guide, section 3.4.4.
- tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- - 1.0e-3*u(nVertLevels,iEdge) &
- *sqrt(2.0*ke_edge(nVertLevels,iEdge))
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - 1.0e-3*u(k,iEdge) &
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+
! old bottom drag, just linear friction
! du/dt = u/tau where tau=100 days.
- !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- ! - u(nVertLevels,iEdge)/(100.0*86400.0)
+ !tend_u(k,iEdge) = tend_u(k,iEdge) &
+ ! - u(k,iEdge)/(100.0*86400.0)
+ endif
+
enddo
!
@@ -648,20 +661,23 @@
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,nVertLevels
+
+ do k=2,maxLevelEdgeTop(iEdge)
fluxVertTop(k) = vertViscTop(k) &
* ( u(k-1,iEdge) - u(k,iEdge) ) &
/ (zMidEdge(k-1,iEdge) - zMidEdge(k,iEdge))
enddo
- do k=1,nVertLevels
+ fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
+
+ do k=1,maxLevelEdgeTop(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)
@@ -697,7 +713,7 @@
type (dm_info) :: dminfo
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdge, maxLevelVertex
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
real (kind=RKIND), dimension(:), pointer :: zTopZLevel
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop, boundaryMask
@@ -728,8 +744,8 @@
zTopZLevel => grid % zTopZLevel % array
boundaryEdge => grid % boundaryEdge % array
maxLevelCell => grid % maxLevelCell % array
- maxLevelEdge => grid % maxLevelEdge % array
- maxLevelVertex => grid % maxLevelVertex % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
nEdges = grid % nEdges
nCells = grid % nCells
@@ -755,8 +771,8 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,nVertLevels
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
do iTracer=1,num_tracers
tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
@@ -774,9 +790,9 @@
cell2 = cellsOnEdge(2,iEdge)
!-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
d2fdx2_cell1 = 0.0
d2fdx2_cell2 = 0.0
@@ -832,9 +848,9 @@
cell2 = cellsOnEdge(2,iEdge)
!-- if an edge is not on the outer-most ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
d2fdx2_cell1 = 0.0
d2fdx2_cell2 = 0.0
@@ -936,7 +952,7 @@
invAreaCell1 = 1.0/areaCell(cell1)
invAreaCell2 = 1.0/areaCell(cell2)
- do k=1,grid % nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
! \kappa_2 </font>
<font color="gray">abla \phi on edge
tracer_turb_flux = config_h_tracer_eddy_diff2 &
@@ -979,7 +995,7 @@
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- do k=1,grid % nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ dvEdge(iEdge)*h_edge(k,iEdge) &
@@ -994,7 +1010,7 @@
end do
- do iCell = 1, nCells
+ do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
@@ -1010,21 +1026,20 @@
invAreaCell1 = 1.0 / areaCell(cell1)
invAreaCell2 = 1.0 / areaCell(cell2)
- do k=1,grid % nVertLevels
- do iTracer=1,num_tracers
- tracer_turb_flux = config_h_tracer_eddy_diff4 &
- *( delsq_tracer(iTracer,k,cell2) &
- - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
- flux = dvEdge (iEdge) * tracer_turb_flux
+ do k=1,maxLevelEdgeTop(iEdge)
+ do iTracer=1,num_tracers
+ tracer_turb_flux = config_h_tracer_eddy_diff4 &
+ *( delsq_tracer(iTracer,k,cell2) &
+ - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * tracer_turb_flux
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) &
- - flux * invAreaCell1 * boundaryMask(k,iEdge)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) &
- + flux * invAreaCell2 * boundaryMask(k,iEdge)
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) &
+ - flux * invAreaCell1 * boundaryMask(k,iEdge)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) &
+ + flux * invAreaCell2 * boundaryMask(k,iEdge)
- end do
+ enddo
enddo
-
end do
deallocate(delsq_tracer)
@@ -1057,8 +1072,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
@@ -1080,7 +1095,6 @@
enddo ! iCell loop
deallocate(fluxVertTop, vertDiffTop)
-
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
! 'iTracer,k, minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))'
@@ -1114,12 +1128,12 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
- integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel, ssh
+ hZLevel, ssh, zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wTop, &
circulation, vorticity, ke, ke_edge, &
@@ -1129,9 +1143,11 @@
real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
+ verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
+ boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdge, maxLevelVertex
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, maxLevelVertexBot
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
@@ -1185,49 +1201,38 @@
hZLevel => grid % hZLevel % array
deriv_two => grid % deriv_two % array
maxLevelCell => grid % maxLevelCell % array
- maxLevelEdge => grid % maxLevelEdge % array
- maxLevelVertex => grid % maxLevelVertex % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
nCells = grid % nCells
nEdges = grid % nEdges
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
boundaryEdge => grid % boundaryEdge % array
boundaryCell => grid % boundaryCell % array
!
- ! Find those cells that have an edge on the boundary
- !
- boundaryCell(:,:) = 0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if(boundaryEdge(k,iEdge).eq.1) then
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- boundaryCell(k,cell1) = 1
- boundaryCell(k,cell2) = 1
- endif
- enddo
- enddo
-
-
- !
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
!
+ ! mrp 101026 efficiency note: zlevel does not need to compute h_edge for k>1
coef_3rd_order = 0.
if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
if (config_thickness_adv_order == 2) then
- do iEdge=1,grid % nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
end do
end if
@@ -1235,14 +1240,14 @@
else if (config_thickness_adv_order == 3) then
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
!-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,grid % nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
d2fdx2_cell1 = 0.0
d2fdx2_cell2 = 0.0
@@ -1287,14 +1292,14 @@
else if (config_thickness_adv_order == 4) then
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
!-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,grid % nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
d2fdx2_cell1 = 0.0
d2fdx2_cell2 = 0.0
@@ -1340,24 +1345,27 @@
!
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,maxLevelVertexBot(vertex1) !mrp remove if for efficiency
+ 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,maxLevelVertexBot(vertex2) !mrp remove if for efficiency
+ 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,maxLevelVertexBot(iVertex)
vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
end do
end do
-
!
! Compute the divergence at each cell center
!
@@ -1365,20 +1373,14 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- enddo
- endif
- if(cell2 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- enddo
- end if
+ do k=1,maxLevelEdgeTop(iEdge)
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ enddo
end do
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
+ do k = 1,maxLevelCell(iCell)
divergence(k,iCell) = divergence(k,iCell) * r
enddo
enddo
@@ -1400,7 +1402,6 @@
end do
!
- !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
@@ -1408,7 +1409,7 @@
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
if (eoe <= nEdges) then
- do k = 1,nVertLevels
+ do k = 1,maxLevelEdgeBot(iEdge)
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
end if
@@ -1418,17 +1419,16 @@
!
! Compute ke on cell edges at velocity locations for quadratic bottom drag.
!
+ ! mrp 101025 efficiency note: we could get rid of ke_edge completely by
+ ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
+ ke_edge = 0.0 !mrp remove 0 for efficiency
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,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
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
@@ -1437,12 +1437,12 @@
! ( 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
+ do i=1,vertexDegree
if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
end do
- do k=1,nVertLevels
+ do k=1,maxLevelVertexBot(iVertex)
h_vertex = 0.0
- do i=1,grid % vertexDegree
+ do i=1,vertexDegree
h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
h_vertex = h_vertex / areaTriangle(iVertex)
@@ -1451,15 +1451,33 @@
end do
end do VTX_LOOP
+ !
+ ! Compute pv at cell centers
+ ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
+ !
+ pv_cell(:,:) = 0.0
+ do iVertex = 1,nVertices
+ do i=1,vertexDegree
+ iCell = cellsOnVertex(i,iVertex)
+ if (iCell <= nCells) then
+ do k = 1,maxLevelCell(iCell) !mrp remove if for efficiency
+ pv_cell(k,iCell) = pv_cell(k,iCell) &
+ + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
+ / areaCell(iCell)
+ enddo
+ endif
+ enddo
+ enddo
!
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
do iEdge = 1,nEdges
- do k = 1,nVertLevels
- gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
- dvEdge(iEdge)
+ do k = 1,maxLevelEdgeBot(iEdge)
+ gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
+ - pv_vertex(k,verticesOnEdge(1,iEdge))) &
+ /dvEdge(iEdge)
enddo
enddo
@@ -1469,10 +1487,10 @@
!
pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
+ do i=1,vertexDegree
iEdge = edgesOnVertex(i,iVertex)
if(iEdge <= nEdges) then
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeBot(iEdge)
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
endif
@@ -1483,47 +1501,31 @@
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
- do k = 1,nVertLevels
+ do k = 1,maxLevelEdgeBot(iEdge)
pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
enddo
enddo
-
!
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
- !
- pv_cell(:,:) = 0.0
- do iVertex = 1, nVertices
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- 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
-
- !
! 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) <= 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)
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
+ gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
+ - pv_cell(k,cellsOnEdge(1,iEdge))) &
+ / dcEdge(iEdge)
enddo
endif
enddo
-
+ !
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
- do k = 1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
enddo
enddo
@@ -1540,6 +1542,12 @@
!
! For a isopycnal model, density should remain constant.
if (config_vert_grid_type.eq.'zlevel') then
+
+! Jacket and McDougal: need to add flag and organize in subroutine.
+! dont forget that eos is called after initial conditions in test_cases right now.
+ !call equation_of_state(s,grid)
+
+! linear is current default
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
! Linear equation of state, for the time being
@@ -1550,73 +1558,124 @@
end do
endif
-
- ! For Isopycnal model.
- ! 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
!
- do iCell=1,nCells
+ ! Z locations
+ !
+ if (config_vert_grid_type.eq.'isopycnal') then
- ! 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.
- zTop(nVertLevels+1,iCell) = h_s(iCell)
- do k=nVertLevels,1,-1
+ ! Compute mid- and top-depths of each layer, from bottom up.
+ do iCell=1,nCells
+ ! h_s is ocean topography: elevation above lowest point,
+ ! and is positive. z coordinates are pos upward.
+ ! h_s is only used for isopycnal coordinates.
+ zTop(nVertLevels+1,iCell) = h_s(iCell)
+ do k=nVertLevels,1,-1
zMid(k,iCell) = zTop(k+1,iCell) + 0.5*h(k,iCell)
zTop(k,iCell) = zTop(k+1,iCell) + h(k,iCell)
- end do
+ end do
+ end do
- ! assume atmospheric pressure at the surface is zero for now.
- pTop(1,iCell) = 0.0
- do k=1,nVertLevels
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if(cell1<=nCells .and. cell2<=nCells) then
+ do k=1,nVertLevels
+ zTopEdge(k,iEdge) = (zTop(k,cell1)+zTop(k,cell2))/2.0
+ zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
+ enddo
+ zTopEdge(nVertLevels+1,iEdge) = ( zTop(nVertLevels+1,cell1) &
+ + zTop(nVertLevels+1,cell2))/2.0
+ endif
+ enddo
+
+ elseif (config_vert_grid_type.eq.'zlevel') then
+
+ ! For z-level coordinates, z variables only change for the top
+ ! layer, so z variables for k=2:nVertLevels+1 can be computed
+ ! from 1D ZLevel variables
+! mrp efficiency. Move zmid ztop ztopedge zmidedge for k>1 to mpas_init, and remove from test_cases
+ do iCell=1,nCells
+ do k=2,nVertLevels
+ zMid(k,iCell) = zMidZLevel(k)
+ zTop(k,iCell) = zTopZLevel(k)
+ end do
+ zTop(nVertLevels+1,iCell) = zTopZLevel(nVertLevels+1)
+
+ zMid(1,iCell) = zTopZLevel(2) + 0.5*h(1,iCell)
+ zTop(1,iCell) = zTopZLevel(2) + h(1,iCell)
+ enddo
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1<=nCells .and. cell2<=nCells) then
+ zTopEdge(1,iEdge) = (zTop(1,cell1)+zTop(1,cell2))/2.0
+ zMidEdge(1,iEdge) = (zMid(1,cell1)+zMid(1,cell2))/2.0
+ endif
+
+ do k=2,nVertLevels
+ zMidEdge(k,iEdge) = zMidZLevel(k)
+ zTopEdge(k,iEdge) = zTopZLevel(k)
+ end do
+ zTopEdge(nVertLevels+1,iEdge) = zTopZLevel(nVertLevels+1)
+ enddo
+
+ endif
+
+ !
+ ! Pressure.
+ ! This section must be after computing rho and zTop.
+ !
+ if (config_vert_grid_type.eq.'isopycnal') then
+
+ ! For Isopycnal model.
+ ! 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
+ do iCell=1,nCells
+
+ ! assume atmospheric pressure at the surface is zero for now.
+ pTop(1,iCell) = 0.0
+ do k=1,nVertLevels
delta_p = rho(k,iCell)*gravity* h(k,iCell)
p(k ,iCell) = pTop(k,iCell) + 0.5*delta_p
pTop(k+1,iCell) = pTop(k,iCell) + delta_p
- end do
+ end do
- MontPot(1,iCell) = gravity * zTop(1,iCell)
- do k=2,nVertLevels
+ MontPot(1,iCell) = gravity * zTop(1,iCell)
+ do k=2,nVertLevels
! from delta M = p delta / rho
MontPot(k,iCell) = MontPot(k-1,iCell) &
+ pTop(k,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
- end do
+ end do
- end do
+ end do
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if(cell1<=nCells .and. cell2<=nCells) then
- do k=1,nVertLevels
- zTopEdge(k,iEdge) = (zTop(k,cell1)+zTop(k,cell2))/2.0
- zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
- enddo
- zTopEdge(nVertLevels+1,iEdge) = ( zTop(nVertLevels+1,cell1) &
- + zTop(nVertLevels+1,cell2))/2.0
- endif
- enddo
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ ! For z-level model.
+ ! Compute pressure at middle of each level.
+ ! pZLevel and p should only differ at k=1, where p is
+ ! pressure at middle of layer including SSH, and pZLevel is
+ ! pressure at a depth of hZLevel(1)/2.
- ! For z-level model.
- ! Compute pressure at middle of each level.
- ! pZLevel and p should only differ at k=1, where p is
- ! pressure at middle of layer including SSH, and pZLevel is
- ! pressure at a depth of hZLevel(1)/2.
- !
- do iCell=1,nCells
- ! compute pressure for z-level coordinates
- ! assume atmospheric pressure at the surface is zero for now.
- pZLevel(1,iCell) = rho(1,iCell)*gravity &
- * (h(1,iCell)-0.5*hZLevel(1))
- do k=2,nVertLevels
- pZLevel(k,iCell) = pZLevel(k-1,iCell) &
- + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
- + rho(k ,iCell)*hZLevel(k ))
- end do
+ do iCell=1,nCells
+ ! compute pressure for z-level coordinates
+ ! assume atmospheric pressure at the surface is zero for now.
- end do
+ pZLevel(1,iCell) = rho(1,iCell)*gravity &
+ * (h(1,iCell)-0.5*hZLevel(1))
+ do k=2,maxLevelCell(iCell)
+ pZLevel(k,iCell) = pZLevel(k-1,iCell) &
+ + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
+ + rho(k ,iCell)*hZLevel(k ))
+ end do
+
+ end do
+
+ endif
+
! compute vertical velocity
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
@@ -1634,13 +1693,16 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells) then
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
+! mrp efficiency note: I can get rid of cell1<=nCells only if these edges
+! are not on the outside edge of a halo. I dont think they are, but
+! double check before you delete it.
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,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
flux = u(k,iEdge) * dvEdge(iEdge)
div_u(k,cell2) = div_u(k,cell2) - flux
end do
@@ -1709,4 +1771,202 @@
end subroutine enforce_boundaryEdge
+ subroutine equation_of_state(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields used in the tendency computations
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+
+
+ type (dm_info) :: dminfo
+ integer :: iEdge, iCell, iVertex, k
+
+ integer :: nCells, nEdges, nVertices, nVertLevels
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ zMidZLevel, pRefEOS
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: &
+ TQ,SQ, &! adjusted T,S
+ BULK_MOD, &! Bulk modulus
+ RHO_S, &! density at the surface
+ DRDT0, &! d(density)/d(temperature), for surface
+ DRDS0, &! d(density)/d(salinity ), for surface
+ DKDT, &! d(bulk modulus)/d(pot. temp.)
+ DKDS, &! d(bulk modulus)/d(salinity )
+ SQR,DENOMK, &! work arrays
+ WORK1, WORK2, WORK3, WORK4, T2, depth
+
+ real (kind=RKIND) :: &
+ tmin, tmax, &! valid temperature range for level k
+ smin, smax ! valid salinity range for level k
+
+ real (kind=RKIND) :: p, p2 ! temporary pressure scalars
+
+!-----------------------------------------------------------------------
+!
+! UNESCO EOS constants and JMcD bulk modulus constants
+!
+!-----------------------------------------------------------------------
+
+ !*** for density of fresh water (standard UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ unt0 = 999.842594, &
+ unt1 = 6.793952e-2, &
+ unt2 = -9.095290e-3, &
+ unt3 = 1.001685e-4, &
+ unt4 = -1.120083e-6, &
+ unt5 = 6.536332e-9
+
+ !*** for dependence of surface density on salinity (UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ uns1t0 = 0.824493 , &
+ uns1t1 = -4.0899e-3, &
+ uns1t2 = 7.6438e-5, &
+ uns1t3 = -8.2467e-7, &
+ uns1t4 = 5.3875e-9, &
+ unsqt0 = -5.72466e-3, &
+ unsqt1 = 1.0227e-4, &
+ unsqt2 = -1.6546e-6, &
+ uns2t0 = 4.8314e-4
+
+ !*** from Table A1 of Jackett and McDougall
+
+ real (kind=RKIND), parameter :: &
+ bup0s0t0 = 1.965933e+4, &
+ bup0s0t1 = 1.444304e+2, &
+ bup0s0t2 = -1.706103 , &
+ bup0s0t3 = 9.648704e-3, &
+ bup0s0t4 = -4.190253e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0s1t0 = 5.284855e+1, &
+ bup0s1t1 = -3.101089e-1, &
+ bup0s1t2 = 6.283263e-3, &
+ bup0s1t3 = -5.084188e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0sqt0 = 3.886640e-1, &
+ bup0sqt1 = 9.085835e-3, &
+ bup0sqt2 = -4.619924e-4
+
+ real (kind=RKIND), parameter :: &
+ bup1s0t0 = 3.186519 , &
+ bup1s0t1 = 2.212276e-2, &
+ bup1s0t2 = -2.984642e-4, &
+ bup1s0t3 = 1.956415e-6
+
+ real (kind=RKIND), parameter :: &
+ bup1s1t0 = 6.704388e-3, &
+ bup1s1t1 = -1.847318e-4, &
+ bup1s1t2 = 2.059331e-7, &
+ bup1sqt0 = 1.480266e-4
+
+ real (kind=RKIND), parameter :: &
+ bup2s0t0 = 2.102898e-4, &
+ bup2s0t1 = -1.202016e-5, &
+ bup2s0t2 = 1.394680e-7, &
+ bup2s1t0 = -2.040237e-6, &
+ bup2s1t1 = 6.128773e-8, &
+ bup2s1t2 = 6.207323e-10
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+
+ nCells = grid % nCells
+ maxLevelCell => grid % maxLevelCell % array
+ nVertLevels = grid % nVertLevels
+ zMidZLevel => grid % zMidZLevel % array
+
+
+! Jackett and McDougall
+ tmin = -2.0 ! valid pot. temp. range
+ tmax = 40.0
+ smin = 0.0 ! valid salinity, in psu
+ smax = 42.0
+
+ ! This could be put in a startup routine.
+ ! Note I am using zMidZlevel, so pressure on top level does
+ ! not include SSH contribution. I am not sure if that matters.
+
+! This function computes pressure in bars from depth in meters
+! using a mean density derived from depth-dependent global
+! average temperatures and salinities from Levitus 1994, and
+! integrating using hydrostatic balance.
+
+ allocate(pRefEOS(nVertLevels))
+ do k = 1,nVertLevels
+ depth = -zMidZLevel(k)
+ pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ + 0.100766*depth + 2.28405e-7*depth**2
+ enddo
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+
+ p = pRefEOS(k)
+ p2 = pRefEOS(k)*pRefEOS(k)
+
+ SQ = max(min(tracers(s%index_salinity,k,iCell),smax),smin)
+ TQ = max(min(tracers(s%index_temperature,k,iCell),tmax),tmin)
+
+ SQR = sqrt(SQ)
+ T2 = TQ*TQ
+
+ !***
+ !*** first calculate surface (p=0) values from UNESCO eqns.
+ !***
+
+ WORK1 = uns1t0 + uns1t1*TQ + &
+ (uns1t2 + uns1t3*TQ + uns1t4*T2)*T2
+ WORK2 = SQR*(unsqt0 + unsqt1*TQ + unsqt2*T2)
+
+ RHO_S = unt1*TQ + (unt2 + unt3*TQ + (unt4 + unt5*TQ)*T2)*T2 &
+ + (uns2t0*SQ + WORK1 + WORK2)*SQ
+
+ !***
+ !*** now calculate bulk modulus at pressure p from
+ !*** Jackett and McDougall formula
+ !***
+
+ WORK3 = bup0s1t0 + bup0s1t1*TQ + &
+ (bup0s1t2 + bup0s1t3*TQ)*T2 + &
+ p *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
+ p2*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
+ WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
+ bup1sqt0*p)
+
+ BULK_MOD = bup0s0t0 + bup0s0t1*TQ + &
+ (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
+ p *(bup1s0t0 + bup1s0t1*TQ + &
+ (bup1s0t2 + bup1s0t3*TQ)*T2) + &
+ p2*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
+ SQ*(WORK3 + WORK4)
+
+ DENOMK = 1.0/(BULK_MOD - p)
+
+ rho(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+
+ end do
+ end do
+
+ deallocate(pRefEOS)
+ end subroutine equation_of_state
+
end module time_integration
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-10-25 18:34:39 UTC (rev 581)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-10-27 14:59:50 UTC (rev 582)
@@ -38,9 +38,10 @@
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
integer, dimension(:), pointer :: &
- maxLevelCell, maxLevelEdge, maxLevelVertex
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, boundaryEdge
+ cellsOnEdge, cellsOnVertex, boundaryEdge, boundaryCell
if (config_vert_grid_type.eq.'isopycnal') then
print *, ' Using isopycnal coordinates'
@@ -61,11 +62,15 @@
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
+ maxLevelCell => block % mesh % maxLevelCell % array
+ maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+ maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
+ maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
+ maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
+ cellsOnEdge => block % mesh % cellsOnEdge % array
cellsOnVertex => block % mesh % cellsOnVertex % array
- boundaryEdge => block % mesh % boundaryEdge % array
+ boundaryEdge => block % mesh % boundaryEdge % array
+ boundaryCell => block % mesh % boundaryCell % array
nCells = block % mesh % nCells
nEdges = block % mesh % nEdges
@@ -85,10 +90,10 @@
endif
- ! Isopycnal grid uses all vertical cells
- if (config_vert_grid_type.eq.'isopycnal') then
+ ! Isopycnal grid uses all vertical cells
+ if (config_vert_grid_type.eq.'isopycnal') then
maxLevelCell = nVertLevels
- endif
+ endif
! mrp insert topography mesa for testing
! do iCell = 1,nCells
@@ -100,47 +105,92 @@
! endif
! enddo
- ! maxLevelEdge is the minimum of the surrounding cells
+ ! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells .and. cell2 <= nCells) then
- maxLevelEdge(iEdge) = &
+ maxLevelEdgeTop(iEdge) = &
min( maxLevelCell(cell1), &
maxLevelCell(cell2) )
else
- maxLevelEdge(iEdge) = 0
+ maxLevelEdgeTop(iEdge) = 0
endif
end do
+ ! maxLevelEdgeBot is the maximum (deepest) of the surrounding cells
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ maxLevelEdgeBot(iEdge) = &
+ max( maxLevelCell(cell1), &
+ maxLevelCell(cell2) )
+ else
+ maxLevelEdgeBot(iEdge) = 0
+ endif
+
+ end do
+
! set boundary edge
boundaryEdge=1
do iEdge=1,nEdges
- boundaryEdge(1:maxLevelEdge(iEdge),iEdge)=0
+ boundaryEdge(1:maxLevelEdgeTop(iEdge),iEdge)=0
end do
- ! maxLevelVertex is the maximum of the surrounding cells
+ !
+ ! Find those cells that have an edge on the boundary
+ !
+ boundaryCell(:,:) = 0
+ do iEdge=1,nEdges
+ do k=1,nVertLevels
+ if (boundaryEdge(k,iEdge).eq.1) then
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ boundaryCell(k,cell1) = 1
+ boundaryCell(k,cell2) = 1
+ endif
+ enddo
+ enddo
+
+ ! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
do iVertex = 1,nVertices
- maxLevelVertex(iVertex) = 0
+ maxLevelVertexBot(iVertex) = 0
do i=1,vertexDegree
cell = cellsOnVertex(i,iVertex)
if (cell <= nCells) then
- maxLevelVertex(iVertex) = &
- max( maxLevelVertex(iVertex), &
+ maxLevelVertexBot(iVertex) = &
+ max( maxLevelVertexBot(iVertex), &
maxLevelCell(cell) )
endif
end do
end do
+ maxLevelVertexBot(nVertices+1) = 0
+ ! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexTop(iVertex) = 0
+ do i=1,vertexDegree
+ cell = cellsOnVertex(i,iVertex)
+ if (cell <= nCells) then
+ maxLevelVertexTop(iVertex) = &
+ min( maxLevelVertexTop(iVertex), &
+ maxLevelCell(cell) )
+ endif
+ end do
+ end do
+ maxLevelVertexTop(nVertices+1) = 0
+
! mrp temp printing, can remove later
print *, 'min/max MaxLevelCell, min/max MaxLevelEdge, ',&
'min/max MaxLevelVertex'
print '(20i5)', &
minval(maxLevelCell(1:nCells)), maxval(maxLevelCell(1:nCells)), &
- minval(maxLevelEdge(1:nEdges)), maxval(maxLevelEdge(1:nEdges)), &
- minval(maxLevelVertex(1:nVertices)), maxval(maxLevelVertex(1:nVertices))
+ minval(maxLevelEdgeTop(1:nEdges)), maxval(maxLevelEdgeTop(1:nEdges)), &
+ minval(maxLevelVertexBot(1:nVertices)), maxval(maxLevelVertexBot(1:nVertices))
block => block % next
end do
@@ -148,19 +198,22 @@
! update halos for maxLevel variables
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_field1dInteger(domain % dminfo, &
- maxLevelCell, block % mesh % nCells, &
+ block % mesh % maxLevelCell % array, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
call dmpar_exch_halo_field1dInteger(domain % dminfo, &
- maxLevelEdge, block % mesh % nEdges, &
+ block % mesh % maxLevelEdgeTop % array, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
call dmpar_exch_halo_field1dInteger(domain % dminfo, &
- maxLevelVertex, block % mesh % nVertices, &
+ block % mesh % maxLevelEdgeBot % array, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+ block % mesh % maxLevelVertexTop % array, block % mesh % nVertices, &
block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+ block % mesh % maxLevelVertexBot % array, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
block => block % next
end do
Modified: branches/ocean_projects/topography2_mrp/src/framework/module_dmpar.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/framework/module_dmpar.F 2010-10-25 18:34:39 UTC (rev 581)
+++ branches/ocean_projects/topography2_mrp/src/framework/module_dmpar.F 2010-10-27 14:59:50 UTC (rev 582)
@@ -141,7 +141,7 @@
#ifdef _MPI
integer :: mpi_ierr
- call MPI_Bcast(i, 1, MPI_INTEGER, IO_NODE, dminfo % comm, mpi_ierr)
+ call MPI_Bcast(i, 1, MPI_INTEGERKIND, IO_NODE, dminfo % comm, mpi_ierr)
#endif
end subroutine dmpar_bcast_int
@@ -158,7 +158,7 @@
#ifdef _MPI
integer :: mpi_ierr
- call MPI_Bcast(iarray, n, MPI_INTEGER, IO_NODE, dminfo % comm, mpi_ierr)
+ call MPI_Bcast(iarray, n, MPI_INTEGERKIND, IO_NODE, dminfo % comm, mpi_ierr)
#endif
end subroutine dmpar_bcast_ints
@@ -216,7 +216,7 @@
end if
end if
- call MPI_Bcast(itemp, 1, MPI_INTEGER, IO_NODE, dminfo % comm, mpi_ierr)
+ call MPI_Bcast(itemp, 1, MPI_INTEGERKIND, IO_NODE, dminfo % comm, mpi_ierr)
if (itemp == 1) then
l = .true.
@@ -255,7 +255,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(i, isum, 1, MPI_INTEGER, MPI_SUM, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(i, isum, 1, MPI_INTEGERKIND, MPI_SUM, dminfo % comm, mpi_ierr)
#else
isum = i
#endif
@@ -293,7 +293,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(i, imin, 1, MPI_INTEGER, MPI_MIN, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(i, imin, 1, MPI_INTEGERKIND, MPI_MIN, dminfo % comm, mpi_ierr)
#else
imin = i
#endif
@@ -331,7 +331,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(i, imax, 1, MPI_INTEGER, MPI_MAX, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(i, imax, 1, MPI_INTEGERKIND, MPI_MAX, dminfo % comm, mpi_ierr)
#else
imax = i
#endif
@@ -370,7 +370,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGER, MPI_SUM, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGERKIND, MPI_SUM, dminfo % comm, mpi_ierr)
#else
outArray = inArray
#endif
@@ -390,7 +390,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGER, MPI_MIN, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGERKIND, MPI_MIN, dminfo % comm, mpi_ierr)
#else
outArray = inArray
#endif
@@ -410,7 +410,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGER, MPI_MAX, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGERKIND, MPI_MAX, dminfo % comm, mpi_ierr)
#else
outArray = inArray
#endif
@@ -491,7 +491,7 @@
#ifdef _MPI
integer :: mpi_ierr
- call MPI_Scatterv(inlist, counts, displs, MPI_INTEGER, outlist, noutlist, MPI_INTEGER, IO_NODE, dminfo % comm, mpi_ierr)
+ call MPI_Scatterv(inlist, counts, displs, MPI_INTEGERKIND, outlist, noutlist, MPI_INTEGERKIND, IO_NODE, dminfo % comm, mpi_ierr)
#endif
end subroutine dmpar_scatter_ints
@@ -534,13 +534,13 @@
#ifdef _MPI
else if (dminfo % my_proc_id == dminfo % nprocs - 1) then
- call MPI_Recv(global_start, 1, MPI_INTEGER, dminfo % my_proc_id - 1, 0, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call MPI_Recv(global_start, 1, MPI_INTEGERKIND, dminfo % my_proc_id - 1, 0, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
global_end = global_start + n - 1
else
- call MPI_Recv(global_start, 1, MPI_INTEGER, dminfo % my_proc_id - 1, 0, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call MPI_Recv(global_start, 1, MPI_INTEGERKIND, dminfo % my_proc_id - 1, 0, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
global_end = global_start + n
- call MPI_Send(global_end, 1, MPI_INTEGER, dminfo % my_proc_id + 1, 0, dminfo % comm, mpi_ierr)
+ call MPI_Send(global_end, 1, MPI_INTEGERKIND, dminfo % my_proc_id + 1, 0, dminfo % comm, mpi_ierr)
global_end = global_end - 1
#endif
@@ -587,7 +587,7 @@
end do
call quicksort(nOwnedList, ownedListSorted)
- call MPI_Allreduce(nNeededList, totalSize, 1, MPI_INTEGER, MPI_MAX, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(nNeededList, totalSize, 1, MPI_INTEGERKIND, MPI_MAX, dminfo % comm, mpi_ierr)
allocate(ownerListIn(totalSize))
allocate(ownerListOut(totalSize))
@@ -636,12 +636,12 @@
end if
nMesgSend = nMesgRecv
- call MPI_Irecv(nMesgRecv, 1, MPI_INTEGER, recvNeighbor, i, dminfo % comm, mpi_rreq, mpi_ierr)
- call MPI_Isend(nMesgSend, 1, MPI_INTEGER, sendNeighbor, i, dminfo % comm, mpi_sreq, mpi_ierr)
+ call MPI_Irecv(nMesgRecv, 1, MPI_INTEGERKIND, recvNeighbor, i, dminfo % comm, mpi_rreq, mpi_ierr)
+ call MPI_Isend(nMesgSend, 1, MPI_INTEGERKIND, sendNeighbor, i, dminfo % comm, mpi_sreq, mpi_ierr)
call MPI_Wait(mpi_rreq, MPI_STATUS_IGNORE, mpi_ierr)
call MPI_Wait(mpi_sreq, MPI_STATUS_IGNORE, mpi_ierr)
- call MPI_Irecv(ownerListIn, nMesgRecv, MPI_INTEGER, recvNeighbor, i, dminfo % comm, mpi_rreq, mpi_ierr)
- call MPI_Isend(ownerListOut, nMesgSend, MPI_INTEGER, sendNeighbor, i, dminfo % comm, mpi_sreq, mpi_ierr)
+ call MPI_Irecv(ownerListIn, nMesgRecv, MPI_INTEGERKIND, recvNeighbor, i, dminfo % comm, mpi_rreq, mpi_ierr)
+ call MPI_Isend(ownerListOut, nMesgSend, MPI_INTEGERKIND, sendNeighbor, i, dminfo % comm, mpi_sreq, mpi_ierr)
call MPI_Wait(mpi_rreq, MPI_STATUS_IGNORE, mpi_ierr)
call MPI_Wait(mpi_sreq, MPI_STATUS_IGNORE, mpi_ierr)
end do
@@ -743,7 +743,7 @@
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_INTEGER, &
+ call MPI_Irecv(recvListPtr % ibuffer, recvListPtr % nlist, MPI_INTEGERKIND, &
recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
end if
recvListPtr => recvListPtr % next
@@ -755,7 +755,7 @@
allocate(sendListPtr % ibuffer(sendListPtr % nlist))
call packSendBuf1dInteger(nOwnedList, arrayIn, sendListPtr, 1, sendListPtr % nlist, &
sendListPtr % ibuffer, nPacked, lastPackedIdx)
- call MPI_Isend(sendListPtr % ibuffer, sendListPtr % nlist, MPI_INTEGER, &
+ 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
@@ -834,7 +834,7 @@
if (recvListPtr % procID /= dminfo % my_proc_id) then
d2 = dim1 * recvListPtr % nlist
allocate(recvListPtr % ibuffer(d2))
- call MPI_Irecv(recvListPtr % ibuffer, d2, MPI_INTEGER, &
+ call MPI_Irecv(recvListPtr % ibuffer, d2, MPI_INTEGERKIND, &
recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
end if
recvListPtr => recvListPtr % next
@@ -847,7 +847,7 @@
allocate(sendListPtr % ibuffer(d2))
call packSendBuf2dInteger(1, dim1, nOwnedList, arrayIn, sendListPtr, 1, d2, &
sendListPtr % ibuffer, nPacked, lastPackedIdx)
- call MPI_Isend(sendListPtr % ibuffer, d2, MPI_INTEGER, &
+ 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
@@ -1240,7 +1240,7 @@
n = (d1e-d1s+1) * (d2e-d2s+1)
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dInteger: Not enough space in buffer', &
+ write(0,*) 'packSendBuf3dInteger: Not enough space in buffer', &
' to fit a single slice.'
return
end if
@@ -1306,7 +1306,7 @@
n = de-ds+1
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dInteger: Not enough space in buffer', &
+ write(0,*) 'packSendBuf2dReal: Not enough space in buffer', &
' to fit a single slice.'
return
end if
@@ -1341,7 +1341,7 @@
n = (d1e-d1s+1) * (d2e-d2s+1)
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dInteger: Not enough space in buffer', &
+ write(0,*) 'packSendBuf3dReal: Not enough space in buffer', &
' to fit a single slice.'
return
end if
@@ -1455,96 +1455,6 @@
end subroutine unpackRecvBuf3dInteger
- subroutine unpackRecvBuf1dReal(nField, field, recvList, startUnpackIdx, nBuffer, buffer, nUnpacked, lastUnpackedIdx)
-
- implicit none
-
- integer, intent(in) :: nField, nBuffer, startUnpackIdx
- 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
-
- integer :: i
-
- nUnpacked = 0
- do i=startUnpackIdx, recvList % nlist
- nUnpacked = nUnpacked + 1
- if (nUnpacked > nBuffer) then
- nUnpacked = nUnpacked - 1
- lastUnpackedIdx = i - 1
- return
- end if
- field(recvList % list(i)) = buffer(nUnpacked)
- end do
- lastUnpackedIdx = recvList % nlist
-
- end subroutine unpackRecvBuf1dReal
-
-
- subroutine unpackRecvBuf2dReal(ds, de, nField, field, recvList, startUnpackIdx, nBuffer, buffer, nUnpacked, lastUnpackedIdx)
-
- implicit none
-
- integer, intent(in) :: ds, de, nField, nBuffer, startUnpackIdx
- 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
-
- integer :: i, n
-
- n = de-ds+1
-
- nUnpacked = 0
- do i=startUnpackIdx, recvList % nlist
- nUnpacked = nUnpacked + n
- if (nUnpacked > nBuffer) then
- nUnpacked = nUnpacked - n
- lastUnpackedIdx = i - 1
- return
- end if
- field(ds:de,recvList % list(i)) = buffer(nUnpacked-n+1:nUnpacked)
- end do
- lastUnpackedIdx = recvList % nlist
-
- end subroutine unpackRecvBuf2dReal
-
-
- subroutine unpackRecvBuf3dReal(d1s, d1e, d2s, d2e, nField, field, recvList, startUnpackIdx, nBuffer, buffer, &
- nUnpacked, lastUnpackedIdx)
-
- implicit none
-
- integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startUnpackIdx
- 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
-
- 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 unpackRecvBuf3dReal
-
-
subroutine dmpar_exch_halo_field1dInteger(dminfo, array, dim1, sendList, recvList)
implicit none
@@ -1735,6 +1645,96 @@
end subroutine dmpar_exch_halo_field3dInteger
+ subroutine unpackRecvBuf1dReal(nField, field, recvList, startUnpackIdx, nBuffer, buffer, nUnpacked, lastUnpackedIdx)
+
+ implicit none
+
+ integer, intent(in) :: nField, nBuffer, startUnpackIdx
+ 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
+
+ integer :: i
+
+ nUnpacked = 0
+ do i=startUnpackIdx, recvList % nlist
+ nUnpacked = nUnpacked + 1
+ if (nUnpacked > nBuffer) then
+ nUnpacked = nUnpacked - 1
+ lastUnpackedIdx = i - 1
+ return
+ end if
+ field(recvList % list(i)) = buffer(nUnpacked)
+ end do
+ lastUnpackedIdx = recvList % nlist
+
+ end subroutine unpackRecvBuf1dReal
+
+
+ subroutine unpackRecvBuf2dReal(ds, de, nField, field, recvList, startUnpackIdx, nBuffer, buffer, nUnpacked, lastUnpackedIdx)
+
+ implicit none
+
+ integer, intent(in) :: ds, de, nField, nBuffer, startUnpackIdx
+ 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
+
+ integer :: i, n
+
+ n = de-ds+1
+
+ nUnpacked = 0
+ do i=startUnpackIdx, recvList % nlist
+ nUnpacked = nUnpacked + n
+ if (nUnpacked > nBuffer) then
+ nUnpacked = nUnpacked - n
+ lastUnpackedIdx = i - 1
+ return
+ end if
+ field(ds:de,recvList % list(i)) = buffer(nUnpacked-n+1:nUnpacked)
+ end do
+ lastUnpackedIdx = recvList % nlist
+
+ end subroutine unpackRecvBuf2dReal
+
+
+ subroutine unpackRecvBuf3dReal(d1s, d1e, d2s, d2e, nField, field, recvList, startUnpackIdx, nBuffer, buffer, &
+ nUnpacked, lastUnpackedIdx)
+
+ implicit none
+
+ integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startUnpackIdx
+ 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
+
+ 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 unpackRecvBuf3dReal
+
+
subroutine dmpar_exch_halo_field1dReal(dminfo, array, dim1, sendList, recvList)
implicit none
</font>
</pre>