<p><b>gaw06e@fsu.edu</b> 2011-07-04 13:59:41 -0600 (Mon, 04 Jul 2011)</p><p>- Only mildly crippled now. Still has issues due to manually setting values for h_vertex and vorticity where boundaryVertex == 2.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/triangle_border_swm/src/core_sw/module_mpas_core.F
===================================================================
--- branches/ocean_projects/triangle_border_swm/src/core_sw/module_mpas_core.F 2011-06-29 20:21:36 UTC (rev 911)
+++ branches/ocean_projects/triangle_border_swm/src/core_sw/module_mpas_core.F 2011-07-04 19:59:41 UTC (rev 912)
@@ -53,7 +53,7 @@
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
- call alter_grid_for_triangle_borders(block % mesh)
+ !call alter_grid_for_triangle_borders(block % mesh)
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
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-29 20:21:36 UTC (rev 911)
+++ branches/ocean_projects/triangle_border_swm/src/core_sw/module_time_integration.F 2011-07-04 19:59:41 UTC (rev 912)
@@ -112,6 +112,7 @@
! BEGIN RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do rk_step = 1, 4
+! write(0,'(a18,i6)')'Starting RK4 step ',rk_step
! --- update halos for diagnostic variables
@@ -324,16 +325,25 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+! write(0,'(a18,i6,a14,i6,a14,i6)')'c_t,tend_h: iEdge ',iEdge,' cOE(1,iEdge) ',&
+! cellsOnEdge(1,iEdge),' cOE(2,iEdge) ',cellsOnEdge(2,iEdge)
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
+! write(0,'(a18,i6,a3,i6,a12,e14.7,a15,e14.7,a17,e14.7,a6,e14.7,a14,e14.7,a14,e14.7)')&
+! 'c_t,tend_h: iEdge ',iEdge,' k ',k,' u(k,iEdge) ',u(k,iEdge),' dvEdge(iEdge) ',&
+! dvEdge(iEdge),' h_edge(k,iEdge) ',h_edge(k,iEdge),' flux ',flux,' tend_h(k,c1) ',&
+! tend_h(k,cell1),' tend_h(k,c2) ',tend_h(k,cell2)
end do
end do
do iCell=1,grid % nCellsSolve
r = 1.0 / areaCell(iCell)
+! write(0,'(a18,i6,a17,e14.7,a21,e14.7)')'c_t,tend_h: iCell ',iCell,' areaCell(iCell) ',&
+! areaCell(iCell),' 1.0/areaCell(iCell) ',(1.0/areaCell(iCell))
do k=1,nVertLevels
tend_h(k,iCell) = tend_h(k,iCell) * r
+! write(0,'(a18,i6,a3,i6,a17,e14.7)')'c_t,tend_h: iCell ',iCell,' k ',k,' tend_h(k,iCell) ',tend_h(k,iCell)
end do
end do
@@ -347,13 +357,19 @@
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
-
+! write(0,'(a18,i6,a11,i6,a11,i6,a11,i6,a11,i6)')'c_t,tend_u: iEdge ',iEdge,' cOE(1,iE) ',&
+! cellsOnEdge(1,iEdge),' cOE(2,iE) ',cellsOnEdge(2,iEdge),' vOE(1,iE) ',&
+! verticesOnEdge(1,iEdge),' vOE(2,iE) ',verticesOnEdge(2,iEdge)
do k=1,nVertLevels
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(j,iEdge)
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
+ q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
+! write(0,'(a14,i6,a3,i6,a15,i6,a15,e14.7,a16,e14.7,a8,e14.7,a11,e14.7,a10,e14.7,a15,e14.7,a3,e14.7)')&
+! 'c_t,tend_u: k ',k,' j ',j,' eOE(j,iE)/eoe ',edgesOnEdge(j,iEdge),' pv_edge(k,iE) ',&
+! pv_edge(k,iEdge),' pv_edge(k,eoe) ',pv_edge(k,eoe),' workpv ',workpv,' wOE(j,iE) ',&
+! weightsOnEdge(j,iEdge),' u(k,eoe) ',u(k,eoe),' h_edge(k,eoe) ',h_edge(k,eoe), ' q ',q
end do
tend_u(k,iEdge) = &
@@ -361,6 +377,11 @@
- ( ke(k,cell2) - ke(k,cell1) + &
gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
) / dcEdge(iEdge)
+! write(0,'(a14,e14.7,a10,e14.7,a10,e14.7,a9,e14.7,a9,e14.7,a11,e14.7,a9,e14.7,a11,e14.7,a12,e14.7,a14,e14.7)')&
+! 'c_t,tend_u: q ',q,' ke(k,c2) ',ke(k,cell2),' ke(k,c1) ',ke(k,cell1),&
+! ' gravity ',gravity,' h(k,c2) ',h(k,cell2),' h_s(c2) ',h_s(cell2),&
+! ' h(k,c1) ',h(k,cell1),' h_s(c1) ',h_s(cell1),' dcEdge(iE) ',dcEdge(iEdge),&
+! ' tend_u(k,iE) ',tend_u(k,iEdge)
end do
end do
@@ -923,6 +944,10 @@
if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
do k=1,grid % nVertLevels
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+! write(0,'(a20,i6,a7,i6,a7,i6,a4,i6,a12,e14.7,a12,e14.7,a17,e14.7)')&
+! 'c_s_d,h_edge: iEdge ',iEdge,' cell1 ',cell1,' cell2 ',cell2,&
+! ' k ',k,' h(k,cell1) ',h(k,cell1),' h(k,cell2) ',h(k,cell2),&
+! ' h_edge(k,iEdge) ',h_edge(k,iEdge)
end do
end if
end do
@@ -1037,12 +1062,24 @@
do k=1,nVertLevels
circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
+! write(0,'(a18,i6,a4,i6,a12,e14.7,a9,e14.7,a19,e14.7,a12,e14.7,a17,e14.7)')&
+! 'c_s_d,circ: iEdge ',iEdge,' k ',k,' dcEdge(iE) ',dcEdge(iEdge),&
+! ' u(k,iE) ',u(k,iEdge),' circ(k,vOE(1,iE)) ',circulation(k,verticesOnEdge(1,iEdge)),&
+! ' circ(k,vOE(2,iE)) ',circulation(k,verticesOnEdge(2,iEdge))
end do
end do
do iVertex=1,nVertices
- do k=1,nVertLevels
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
- end do
+
+ do k=1,nVertLevels
+ if (boundaryVertex(k,iVertex).ne.2) then
+ vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+! write(0,'(a15,i6,a12,e14.7,a11,e14.7,a12,e14.7)')&
+! 'c_s_d,vort: iV ',iVertex,' circ(k,iV) ',circulation(k,iVertex),' areaT(iV) ',&
+! areaTriangle(iVertex),' vort(k,iV) ',vorticity(k,iVertex)
+ else
+ vorticity(k,iVertex) = 0.0
+ end if
+ end do
end do
@@ -1122,15 +1159,34 @@
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
do iVertex = 1,nVertices
- do k=1,nVertLevels
- h_vertex(k,iVertex) = 0.0
- do i=1,grid % vertexDegree
- h_vertex(k,iVertex) = h_vertex(k,iVertex) + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
- end do
- h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)
-
- pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
- end do
+! write(0,'(a20,i6,a22,i6)')'c_s_d,pv_vertex: iV ',iVertex,' boundaryVertex(1,iV) ',&
+! boundaryVertex(1,iVertex)
+ do k=1,nVertLevels
+ if (boundaryVertex(k,iVertex).eq.0) then
+ h_vertex(k,iVertex) = 0.0
+ do i=1,grid % vertexDegree
+ h_vertex(k,iVertex) = h_vertex(k,iVertex) + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
+! write(0,'(a19,i6,a3,i6,a16,e14.7,a12,e14.7,a16,e14.7)')'c_s_d,pv_vertex: k ',k,' i ',i,&
+! ' h(k,cOV(i,iV)) ',h(k,cellsOnVertex(i,iVertex)),' kAOV(i,iV) ',&
+! kiteAreasOnVertex(i,iVertex),' h_vertex(k,iV) ',h_vertex(k,iVertex)
+ end do
+ h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)
+! write(0,'(a24,e14.7,a16,e14.7)')'c_s_d,pv_vertex: aT(iV) ',areaTriangle(iVertex),&
+! ' h_vertex(k,iV) ',h_vertex(k,iVertex)
+ if (h_vertex(k,iVertex).eq.0.0) then
+ h_vertex(k,iVertex)=700.0
+ write(0,'(a7)')'WHOOPS!' ! this hasn't happened yet, in my experience, but was left for debug
+ end if
+ elseif(boundaryVertex(k,iVertex).eq.2) then
+ h_vertex(k,iVertex) = 700.0
+! write(0,'(a19,i6,a16,e14.7)')'c_s_d,pv_vertex: k ',k,' h_vertex(k,iV) ',h_vertex(k,iVertex)
+ end if
+
+ pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
+! write(0,'(a29,e14.7,a17,e14.7,a16,e14.7,a17,e14.7)'),'c_s_d,pv_vertex: fVertex(iV) ',&
+! fVertex(iVertex),' vorticity(k,iV) ',vorticity(k,iVertex),' h_vertex(k,iV) ',&
+! h_vertex(k,iVertex),' pv_vertex(k,iV) ',pv_vertex(k,iVertex)
+ end do
end do
@@ -1142,6 +1198,10 @@
do k = 1,nVertLevels
gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
dvEdge(iEdge)
+! write(0,'(a18,i6,a3,i6,a24,e14.7,a24,e14.7,a12,e14.7,a15,e14.7)')'c_s_d,gradPVt: iE ',&
+! iEdge,' k ',k,' pv_vertex(k,vOE(2,iE)) ',pv_vertex(k,verticesOnEdge(2,iEdge)),&
+! ' pv_vertex(k,vOE(1,iE)) ',pv_vertex(k,verticesOnEdge(1,iEdge)),' dvEdge(iE) ',&
+! dvEdge(iEdge),' gradPVt(k,iE) ',gradPVt(k,iEdge)
enddo
enddo
</font>
</pre>