<p><b>gaw06e@fsu.edu</b> 2011-06-04 17:59:53 -0600 (Sat, 04 Jun 2011)</p><p>- removed hardcoded if statements that altered areaCell/dvEdge values now that alter_grid_for_triangle_borders is implemented<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F 2011-06-04 23:37:20 UTC (rev 874)
+++ branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F 2011-06-04 23:59:53 UTC (rev 875)
@@ -274,7 +274,6 @@
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
real (kind=RKIND) :: ke_edge
- real (kind=RKIND) :: thisdvEdge
h => s % h % array
u => s % u % array
@@ -325,24 +324,14 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
do k=1,nVertLevels
- !flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- flux = u(k,iEdge) * thisdvEdge * h_edge(k,iEdge)
+ 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
end do
end do
do iCell=1,grid % nCellsSolve
- if (boundaryCell(1,iCell).eq.1) then
- r = 1.0 / ( areaCell(iCell) * 2.0 )
- else
- r = 1.0 / areaCell(iCell)
- end if
+ r = 1.0 / areaCell(iCell)
do k=1,nVertLevels
tend_h(k,iCell) = tend_h(k,iCell) * r
end do
@@ -412,15 +401,9 @@
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
do k=1,nVertLevels
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- ! -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / thisdvEdge
+ -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
end do
@@ -480,24 +463,15 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
do k=1,nVertLevels
delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
- + delsq_u(k,iEdge)*thisdvEdge
+ + delsq_u(k,iEdge)*dvEdge(iEdge)
delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
- - delsq_u(k,iEdge)*thisdvEdge
+ - delsq_u(k,iEdge)*dvEdge(iEdge)
end do
end do
do iCell = 1,nCells
- if (boundaryCell(1,iCell).eq.1) then
- r = 1.0 / ( areaCell(iCell) * 2.0 )
- else
- r = 1.0 / areaCell(iCell)
- end if
+ r = 1.0 / areaCell(iCell)
do k = 1,nVertLevels
delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
end do
@@ -510,17 +484,12 @@
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
do k=1,nVertLevels
u_diffusion = ( delsq_divergence(k,cell2) &
- delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
-( delsq_vorticity(k,vertex2) &
- - delsq_vorticity(k,vertex1) ) / thisdvEdge
+ - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
u_diffusion = config_h_mom_eddy_visc4 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
@@ -592,8 +561,6 @@
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND), dimension(:,:), pointer :: u, h_edge
- real (kind=RKIND) :: thisdvEdge
-
u => s % u % array
h_edge => s % h_edge % array
dcEdge => grid % dcEdge % array
@@ -619,26 +586,13 @@
do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- if (boundaryCell(1,cell1).eq.1) then
- invAreaCell1 = 1.0 / ( areaCell(cell1) * 2.0 )
- else
- invAreaCell1 = 1.0 / areaCell(cell1)
- end if
- if (boundaryCell(1,cell2).eq.1) then
- invAreaCell2 = 1.0 / ( areaCell(cell2) * 2.0 )
- else
- invAreaCell2 = 1.0 / areaCell(cell2)
- end if
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * thisdvEdge * h_edge(k,iEdge) * tracer_edge
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1
tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2
end do
@@ -769,21 +723,8 @@
do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
- if (boundaryCell(1,cell1).eq.1) then
- invAreaCell1 = 1.0 / ( areaCell(cell1) * 2.0 )
- else
- invAreaCell1 = 1.0 / areaCell(cell1)
- end if
- if (boundaryCell(1,cell2).eq.1) then
- invAreaCell2 = 1.0 / ( areaCell(cell2) * 2.0 )
- else
- invAreaCell2 = 1.0 / areaCell(cell2)
- end if
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
do k=1,grid % nVertLevels
do iTracer=1, grid % nTracers
! \kappa_2 </font>
<font color="gray">abla \phi on edge
@@ -791,7 +732,7 @@
*( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
! div(h \kappa_2 </font>
<font color="gray">abla \phi) at cell center
- flux = thisdvEdge * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) - flux * invAreaCell2
end do
@@ -824,29 +765,19 @@
do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
do k=1,grid % nVertLevels
do iTracer=1, grid % nTracers
delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
- + thisdvEdge * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+ + dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
- - thisdvEdge * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
+ - dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
end do
end do
end do
do iCell = 1, grid % nCells
- if (boundaryCell(1,cell1).eq.1) then
- r = 1.0 / ( areaCell(cell1) * 2.0 )
- else
- r = 1.0 / areaCell(cell1)
- end if
- !r = 1.0 / grid % areaCell % array(iCell)
+ r = 1.0 / areaCell(cell1)
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
@@ -858,28 +789,12 @@
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
- if (boundaryCell(1,cell1).eq.1) then
- invAreaCell1 = 1.0 / ( areaCell(cell1) * 2.0 )
- else
- invAreaCell1 = 1.0 / areaCell(cell1)
- end if
- if (boundaryCell(1,cell2).eq.1) then
- invAreaCell2 = 1.0 / ( areaCell(cell2) * 2.0 )
- else
- invAreaCell2 = 1.0 / areaCell(cell2)
- end if
- !invAreaCell1 = 1.0 / grid % areaCell % array(cell1)
- !invAreaCell2 = 1.0 / grid % areaCell % array(cell2)
-
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) / dcEdge(iEdge)
- flux = thisdvEdge * tracer_turb_flux
+ flux = dvEdge(iEdge) * tracer_turb_flux
tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1 * boundaryMask(k,iEdge)
tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2 * boundaryMask(k,iEdge)
end do
@@ -926,7 +841,6 @@
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
- real (kind=RKIND) :: thisdvEdge
h => s % h % array
u => s % u % array
@@ -1139,29 +1053,19 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
if (cell1 <= nCells) then
do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*thisdvEdge
+ 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)*thisdvEdge
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
end if
end do
do iCell = 1,nCells
- if (boundaryCell(1,iCell).eq.1) then
- r = 1.0 / ( areaCell(iCell) * 2.0 )
- else
- r = 1.0 / areaCell(iCell)
- end if
- !r = 1.0 / areaCell(iCell)
+ r = 1.0 / areaCell(iCell)
do k = 1,nVertLevels
divergence(k,iCell) = divergence(k,iCell) * r
enddo
@@ -1172,20 +1076,11 @@
!
ke(:,:) = 0.0
do iCell=1,nCells
- if (boundaryCell(1,iCell).eq.1) then
- r = 1.0 / ( areaCell(iCell) * 2.0 )
- else
- r = 1.0 / areaCell(iCell)
- end if
+ r = 1.0 / areaCell(iCell)
do i=1,nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
do k=1,nVertLevels
- ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * thisdvEdge * u(k,iEdge)**2.0
+ ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
end do
end do
do k=1,nVertLevels
@@ -1244,14 +1139,9 @@
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
do iEdge = 1,nEdges
- if (boundaryEdge(1,iEdge).eq.1) then
- thisdvEdge = dvEdge(iEdge) / 2
- else
- thisdvEdge = dvEdge(iEdge)
- end if
do k = 1,nVertLevels
gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
- thisdvEdge
+ dvEdge(iEdge)
enddo
enddo
@@ -1289,11 +1179,7 @@
do iVertex = 1, nVertices
do i=1,grid % vertexDegree
iCell = cellsOnVertex(i,iVertex)
- if (boundaryCell(1,iCell).eq.1) then
- r = 1.0 / ( areaCell(iCell) * 2.0 )
- else
- r = 1.0 / areaCell(iCell)
- end if
+ r = 1.0 / areaCell(iCell)
if (iCell <= nCells) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) * r
</font>
</pre>