<p><b>mpetersen@lanl.gov</b> 2010-06-16 15:51:04 -0600 (Wed, 16 Jun 2010)</p><p>Put if constructs around computations specific to isopycnal or zlevel grids, like MontPot and pZLevel. For zlevel grid, removed computation of tend_h for k>1, which is zero anyway. These changes made a x3 quad grid 7% faster.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F 2010-06-14 16:44:45 UTC (rev 356)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/module_test_cases.F 2010-06-16 21:51:04 UTC (rev 357)
@@ -26,7 +26,7 @@
type (dm_info) :: dminfo
! mrp 100507: for diagnostic output
- integer :: iTracer, cell1, cell2, cell
+ integer :: iTracer, cell1, cell2, cell, k
real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
hZLevel, zMidZLevel, zTopZLevel, latCell, lonCell
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
@@ -130,9 +130,12 @@
nVertLevels = block % mesh % nVertLevels
vertexDegree = block % mesh % vertexDegree
- ! mrp 100608 For now, use a flat bottom
-! maxLevelCell = nVertLevels
+ ! Isopycnal grid uses all vertical cells
+ if (config_vert_grid_type.eq.'isopycnal') then
+ maxLevelCell = nVertLevels
+ endif
+ ! mrp insert topography mesa for testing
! do iCell = 1,nCells
! if (latCell(iCell)> 0.0/180.*3.14.and.&
! latCell(iCell)< 30.0/180.*3.14.and. &
@@ -176,15 +179,6 @@
end do
end do
- if (config_vert_grid_type.eq.'zlevel') then
- ! These should eventually be in an input file. For now
- ! I just read them in from h(:,1).
- hZLevel = h(:,1)
- zTopZLevel(1) = 0.0
- do iLevel = 1,nVertLevels
- zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
- zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
- enddo
if (config_vert_grid_type.eq.'isopycnal') then
print *, ' Using isopycnal coordinates'
elseif (config_vert_grid_type.eq.'zlevel') then
@@ -195,10 +189,22 @@
call dmpar_abort(dminfo)
endif
+
+ if (config_vert_grid_type.eq.'zlevel') then
+
+ ! These should eventually be in an input file. For now
+ ! I just read them in from h(:,1).
+ hZLevel = h(:,1)
+ zTopZLevel(1) = 0.0
+ do iLevel = 1,nVertLevels
+ zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
+ zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
+ enddo
+
! Set tracers, if not done in grid.nc file
- tracers = -2.0
+! tracers = -2.0
do iCell = 1,nCells
- tracers(index_tracer2,1,iCell) = maxLevelCell(iCell)
+! tracers(index_tracer2,1,iCell) = maxLevelCell(iCell)
! tracers(index_tracer1,1,iCell) = latCell(iCell)
! tracers(index_tracer2,1,iCell) = lonCell(iCell)
do iLevel = 1,maxLevelCell(iCell)
@@ -207,10 +213,10 @@
! tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
! for x3, 25 layer test
- tracers(index_temperature,iLevel,iCell) = 10.0 ! temperature
- tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
+! tracers(index_temperature,iLevel,iCell) = 10.0 ! temperature
+! tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
- tracers(index_tracer1,iLevel,iCell) = 1.0
+! tracers(index_tracer1,iLevel,iCell) = 1.0
! tracers(index_tracer1,iLevel,iCell) = maxLevelCell(iCell)
! tracers(index_tracer2,iLevel,iCell) = &
! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
@@ -224,7 +230,6 @@
endif
-
! print some diagnostics
print '(10a)', 'ilevel',&
' rho ',&
@@ -252,23 +257,22 @@
block => block % next
end do
-! --- update halos for maxLevel variables
-! mrp 100610 Integer halo updates currently cause a seg fault.
+ ! 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_field2dInteger(domain % dminfo, maxLevelCell, &
-! block % mesh % nVertLevels, block % mesh % nCells, &
-! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-! call dmpar_exch_halo_field2dInteger(domain % dminfo, maxLevelEdge, &
-! block % mesh % nVertLevels, block % mesh % nEdges, &
-! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-! call dmpar_exch_halo_field2dInteger(domain % dminfo, maxLevelVertex, &
-! block % mesh % nVertLevels, block % mesh % nVertices, &
-! block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+ maxLevelCell, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+ maxLevelEdge, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+ maxLevelVertex, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
block => block % next
end do
Modified: branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F 2010-06-14 16:44:45 UTC (rev 356)
+++ branches/ocean_projects/topography_mrp/src/core_ocean/module_time_integration.F 2010-06-16 21:51:04 UTC (rev 357)
@@ -250,7 +250,8 @@
type (grid_state), intent(in) :: s
type (grid_meta), 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
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
@@ -328,32 +329,39 @@
!
! Initialize tendencies to zero
!
- tend_h(:,:) = 0.0
- tend_u(:,:) = 0.0
+ tend_h = 0.0
+ tend_u = 0.0
!
! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
!
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
! for explanation of divergence operator.
+ !
+ ! 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,maxLevelCell(cell1)
+ 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,maxLevelCell(cell2)
+ 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
@@ -361,22 +369,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,maxLevelCell(iCell)
- tend_h(k,iCell) = tend_h(k,iCell) &
- - (wTop(k,iCell) - wTop(k+1,iCell))
- end do
-
end do
endif ! coordinate type
@@ -408,12 +404,6 @@
deallocate(w_dudzTopEdge)
endif
- if (1==2) then !if (isNaN(sum(tend_u))) then
- write(0,*) '1 Abort: NaN detected in vertical advection'
- call dmpar_abort(dminfo)
- endif
-
-
!
! velocity tendency: pressure gradient
!
@@ -439,10 +429,6 @@
enddo
endif
- if (1==2) then !if (isNaN(sum(tend_u))) then
- write(0,*) '2 Abort: NaN detected in p grad'
- call dmpar_abort(dminfo)
- endif
!
! velocity tendency: -q(h u^\perp) - </font>
<font color="black">abla K
! +</font>
<font color="black">u_h(</font>
<font color="black">abla \delta + {\bf k}\times </font>
<font color="gray">abla \xi)
@@ -475,10 +461,6 @@
end do
end do
- if (1==2) then !if (isNaN(sum(tend_u))) then
- write(0,*) '3 Abort: NaN detected in hor dif'
- call dmpar_abort(dminfo)
- endif
!
! velocity tendency: forcing and bottom drag
!
@@ -520,11 +502,6 @@
enddo
- if (1==2) then !if (isNaN(sum(tend_u))) then
- write(0,*) '4 Abort: NaN detected in forcing'
- call dmpar_abort(dminfo)
- endif
-
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
@@ -569,11 +546,6 @@
end do
deallocate(fluxVertTop, vertViscTop)
- if (1==2) then !if (isNaN(sum(tend_u))) then
- write(0,*) '3 Abort: NaN detected in vertical mix'
- call dmpar_abort(dminfo)
- endif
-
end subroutine compute_tend
@@ -691,11 +663,6 @@
end do
end do
- if (1==2) then !if (isNaN(sum(tend_tr))) then
- write(0,*) '1 Abort: NaN detected tend_tr'
- call dmpar_abort(dminfo)
- endif
-
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
@@ -737,11 +704,6 @@
enddo ! iCell
deallocate(tracerTop)
- if (1==2) then !if (isNaN(sum(tend_tr))) then
- write(0,*) ' 2 Abort: NaN detected tend_tr'
- call dmpar_abort(dminfo)
- endif
-
!
! tracer tendency: horizontal tracer diffusion
! div(h \kappa_h </font>
<font color="gray">abla\phi )
@@ -798,11 +760,6 @@
enddo
deallocate(tr_flux, tr_div)
- if (1==2) then !if (isNaN(sum(tend_tr))) then
- write(0,*) '3 Abort: NaN detected tend_tr'
- call dmpar_abort(dminfo)
- endif
-
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
@@ -852,11 +809,6 @@
enddo ! iCell loop
deallocate(fluxVertTop, vertDiffTop)
- if (1==2) then !if (isNaN(sum(tend_tr))) then
- write(0,*) '4 Abort: NaN detected in tend tr vert diff'
- call dmpar_abort(dminfo)
- endif
-
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
! 'iTracer,k, minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))'
@@ -896,7 +848,7 @@
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, &
@@ -958,8 +910,10 @@
fEdge => grid % fEdge % array
hZLevel => grid % hZLevel % array
maxLevelCell => grid % maxLevelCell % array
- maxLevelEdge => grid % maxLevelEdge % array
- maxLevelVertex => grid % maxLevelVertex % array
+ maxLevelEdge => grid % maxLevelEdge % array
+ maxLevelVertex => grid % maxLevelVertex % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -1225,73 +1179,121 @@
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
+ 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,maxLevelCell(iCell)
- delta_p = rho(k,iCell)*gravity*hZLevel(k)
- pZLevel(k,iCell) = pZLevel(k-1,iCell) + 0.5*delta_p
- end do
+ 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,maxLevelCell(iCell)
+ delta_p = rho(k,iCell)*gravity*hZLevel(k)
+ pZLevel(k,iCell) = pZLevel(k-1,iCell) + 0.5*delta_p
+ end do
- end do
+ end do
+ endif
+
! compute vertical velocity
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
</font>
</pre>