<p><b>dwj07@fsu.edu</b> 2012-04-17 11:03:18 -0600 (Tue, 17 Apr 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Cleaning up some of the deriv_two computation.<br>
Adding documentation and changing variable names.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/advection_cleanup/src/core_ocean/mpas_ocn_advection.F
===================================================================
--- branches/ocean_projects/advection_cleanup/src/core_ocean/mpas_ocn_advection.F 2012-04-17 17:02:45 UTC (rev 1785)
+++ branches/ocean_projects/advection_cleanup/src/core_ocean/mpas_ocn_advection.F 2012-04-17 17:03:18 UTC (rev 1786)
@@ -31,22 +31,18 @@
! local variables
- real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
+ real (kind=RKIND), dimension(2, grid % nEdges) :: edgeAngleToEast
real (kind=RKIND), dimension(grid % nEdges) :: xe, ye
- real (kind=RKIND), dimension(grid % nCells) :: theta_abs
- real (kind=RKIND), dimension(25) :: xc, yc, zc ! cell center coordinates
- real (kind=RKIND), dimension(25) :: thetav, thetat, dl_sphere
- real (kind=RKIND) :: xm, ym, zm, dl, xec, yec, zec
- real (kind=RKIND) :: thetae_tmp, xe_tmp, ye_tmp
- real (kind=RKIND) :: xv1, xv2, yv1, yv2, zv1, zv2
- integer :: i, j, k, ip1, ip2, m, n, ip1a, ii
+ real (kind=RKIND), dimension(25) :: xCellScaled, yCellScaled, zCellScaled
+ real (kind=RKIND), dimension(25) :: angleToEast, wedgeAngle, sphericalDistance
+ real (kind=RKIND) :: referenceAngleToEast
+ real (kind=RKIND) :: xEdgeMidpoint, yEdgeMidpoint, zEdgeMidpoint
+ real (kind=RKIND) :: edgeAngleToEast_tmp
+ real (kind=RKIND) :: xVertexScaled1, xVertexScaled2, yVertexScaled1, yVertexScaled2, zVertexScaled1, zVertexScaled2
+ integer :: i, j, k, ip1, ip2, nAdvectionCells
integer :: iCell, iEdge
- real (kind=RKIND) :: pii
- real (kind=RKIND) :: x0, y0, x1, y1, x2, y2, x3, y3, x4, y4, x5, y5
- real (kind=RKIND) :: pdx1, pdx2, pdx3, pdy1, pdy2, pdy3, dx1, dx2, dy1, dy2
- real (kind=RKIND) :: angv1, angv2, dl1, dl2
- real (kind=RKIND), dimension(25) :: dxe, dye, x2v, y2v, xp, yp
+ real (kind=RKIND), dimension(25) :: xCellPlanar, yCellPlanar
real (kind=RKIND) :: amatrix(25,25), bmatrix(25,25), wmatrix(25,25)
real (kind=RKIND) :: length_scale
@@ -63,7 +59,6 @@
logical, parameter :: reset_poly = .true.
real (kind=RKIND) :: rcell, cos2t, costsint, sin2t
- real (kind=RKIND), dimension(grid%maxEdges) :: angle_2d
!---
err = 0
@@ -74,116 +69,109 @@
return
end if
- pii = 2.*asin(1.0)
-
-! advCells => grid % advCells % array
- allocate(advCells(grid % maxEdges2))
+ allocate(advCells(grid % maxEdges2+1))
deriv_two => grid % deriv_two % array
deriv_two(:,:,:) = 0.
- do iCell = 1, grid % nCells ! is this correct? - we need first halo cell also...
+ ! Loop over cells
+ do iCell = 1, grid % nCells
cell_list(1) = iCell
do i=2, grid % nEdgesOnCell % array(iCell)+1
cell_list(i) = grid % CellsOnCell % array(i-1,iCell)
end do
- n = grid % nEdgesOnCell % array(iCell) + 1
+ nAdvectionCells = grid % nEdgesOnCell % array(iCell) + 1
if ( polynomial_order > 2 ) then
do i=2,grid % nEdgesOnCell % array(iCell) + 1
do j=1,grid % nEdgesOnCell % array ( cell_list(i) )
cell_add = grid % CellsOnCell % array (j,cell_list(i))
add_the_cell = .true.
- do k=1,n
+ do k=1,nAdvectionCells
if ( cell_add == cell_list(k) ) add_the_cell = .false.
end do
if (add_the_cell) then
- n = n+1
- cell_list(n) = cell_add
+ nAdvectionCells = nAdvectionCells+1
+ cell_list(nAdvectionCells) = cell_add
end if
end do
end do
end if
- advCells(1) = n
+ advCells(1) = nAdvectionCells
-! check to see if we are reaching outside the halo
-
+ ! Check if any cell is outside the 2 halo of the current block
do_the_cell = .true.
- do i=1,n
+ do i=1,nAdvectionCells
if (cell_list(i) > grid % nCells) do_the_cell = .false.
end do
-
+ ! Check to see if the cell should be skipped
if ( .not. do_the_cell ) cycle
-
-! compute poynomial fit for this cell if all needed neighbors exist
if ( grid % on_a_sphere ) then
- do i=1,n
+ ! Scale cell centers to a unit sphere
+ do i=1,nAdvectionCells
advCells(i+1) = cell_list(i)
- xc(i) = grid % xCell % array(advCells(i+1))/a
- yc(i) = grid % yCell % array(advCells(i+1))/a
- zc(i) = grid % zCell % array(advCells(i+1))/a
+ xCellScaled(i) = grid % xCell % array(advCells(i+1))/grid % sphere_radius
+ yCellScaled(i) = grid % yCell % array(advCells(i+1))/grid % sphere_radius
+ zCellScaled(i) = grid % zCell % array(advCells(i+1))/grid % sphere_radius
end do
- theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
- xc(2), yc(2), zc(2), &
+ ! Compute the angle between the x axis, and the vector connecting the
+ ! iCell and cellsOnCell(1, iCell)
+ referenceAngleToEast = pii/2. - sphere_angle( xCellScaled(1), yCellScaled(1), zCellScaled(1), &
+ xCellScaled(2), yCellScaled(2), zCellScaled(2), &
0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
-! angles from cell center to neighbor centers (thetav)
-
- do i=1,n-1
+ ! wedgeAngle is the angle between two vectors connecting iCell to consecutive entries in cellsonCell
+ ! and sphericalDistance is the distance on the sphere between iCell and cellsOnCell(:, iCell)
+ do i=1,nAdvectionCells-1
ip2 = i+2
- if (ip2 > n) ip2 = 2
+ if (ip2 > nAdvectionCells) ip2 = 2
- thetav(i) = sphere_angle( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1), &
- xc(ip2), yc(ip2), zc(ip2) )
+ wedgeAngle(i) = sphere_angle( xCellScaled(1), yCellScaled(1), zCellScaled(1), &
+ xCellScaled(i+1), yCellScaled(i+1), zCellScaled(i+1), &
+ xCellScaled(ip2), yCellScaled(ip2), zCellScaled(ip2) )
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1) )
+ sphericalDistance(i) = grid % sphere_radius*arc_length( xCellScaled(1), yCellScaled(1), zCellScaled(1), &
+ xCellScaled(i+1), yCellScaled(i+1), zCellScaled(i+1) )
end do
- length_scale = 1.
- do i=1,n-1
- dl_sphere(i) = dl_sphere(i)/length_scale
+ ! Using the referenceAngleToXAxis variable, and the wedgeAngle
+ ! build the angleToXAxis variable, which contains the angle between
+ ! each iCell -> cellsOnCell(:, iCell) vector and the local east direction (x axis in plane)
+ angleToEast(1) = referenceAngleToEast
+ do i=2,nAdvectionCells-1
+ angleToEast(i) = angleToEast(i-1) + wedgeAngle(i-1)
end do
-
-! thetat(1) = 0. ! this defines the x direction, cell center 1 ->
- thetat(1) = theta_abs(iCell) ! this defines the x direction, longitude line
- do i=2,n-1
- thetat(i) = thetat(i-1) + thetav(i-1)
- end do
- do i=1,n-1
- xp(i) = cos(thetat(i)) * dl_sphere(i)
- yp(i) = sin(thetat(i)) * dl_sphere(i)
+ do i=1,nAdvectionCells-1
+ xCellPlanar(i) = cos(angleToEast(i)) * sphericalDistance(i)
+ yCellPlanar(i) = sin(angleToEast(i)) * sphericalDistance(i)
end do
else ! On an x-y plane
- do i=1,n-1
+ ! angleEdge is the angle the edge makes with local east
+ do i=1,nAdvectionCells-1
- angle_2d(i) = grid%angleEdge%array(grid % EdgesOnCell % array(i,iCell))
+ angleToEast(i) = grid%angleEdge%array(grid % EdgesOnCell % array(i,iCell))
iEdge = grid % EdgesOnCell % array(i,iCell)
if ( iCell .ne. grid % CellsOnEdge % array(1,iEdge)) &
- angle_2d(i) = angle_2d(i) - pii
+ angleToEast(i) = angleToEast(i) - pii
-! xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
-! yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+ xCellPlanar(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * cos(angleToEast(i))
+ yCellPlanar(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * sin(angleToEast(i))
- xp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * cos(angle_2d(i))
- yp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * sin(angle_2d(i))
-
end do
end if
-
- ma = n-1
+ ! Build polynomial matrices to solve for second derivative
+ ma = nAdvectionCells-1
mw = grid % nEdgesOnCell % array (iCell)
bmatrix = 0.
@@ -198,11 +186,11 @@
wmatrix(1,1) = 1.
do i=2,ma
amatrix(i,1) = 1.
- amatrix(i,2) = xp(i-1)
- amatrix(i,3) = yp(i-1)
- amatrix(i,4) = xp(i-1)**2
- amatrix(i,5) = xp(i-1) * yp(i-1)
- amatrix(i,6) = yp(i-1)**2
+ amatrix(i,2) = xCellPlanar(i-1)
+ amatrix(i,3) = yCellPlanar(i-1)
+ amatrix(i,4) = xCellPlanar(i-1)**2
+ amatrix(i,5) = xCellPlanar(i-1) * yCellPlanar(i-1)
+ amatrix(i,6) = yCellPlanar(i-1)**2
wmatrix(i,i) = 1.
end do
@@ -215,17 +203,17 @@
wmatrix(1,1) = 1.
do i=2,ma
amatrix(i,1) = 1.
- amatrix(i,2) = xp(i-1)
- amatrix(i,3) = yp(i-1)
+ amatrix(i,2) = xCellPlanar(i-1)
+ amatrix(i,3) = yCellPlanar(i-1)
- amatrix(i,4) = xp(i-1)**2
- amatrix(i,5) = xp(i-1) * yp(i-1)
- amatrix(i,6) = yp(i-1)**2
+ amatrix(i,4) = xCellPlanar(i-1)**2
+ amatrix(i,5) = xCellPlanar(i-1) * yCellPlanar(i-1)
+ amatrix(i,6) = yCellPlanar(i-1)**2
- amatrix(i,7) = xp(i-1)**3
- amatrix(i,8) = yp(i-1) * (xp(i-1)**2)
- amatrix(i,9) = xp(i-1) * (yp(i-1)**2)
- amatrix(i,10) = yp(i-1)**3
+ amatrix(i,7) = xCellPlanar(i-1)**3
+ amatrix(i,8) = yCellPlanar(i-1) * (xCellPlanar(i-1)**2)
+ amatrix(i,9) = xCellPlanar(i-1) * (yCellPlanar(i-1)**2)
+ amatrix(i,10) = yCellPlanar(i-1)**3
wmatrix(i,i) = 1.
@@ -239,23 +227,23 @@
wmatrix(1,1) = 1.
do i=2,ma
amatrix(i,1) = 1.
- amatrix(i,2) = xp(i-1)
- amatrix(i,3) = yp(i-1)
+ amatrix(i,2) = xCellPlanar(i-1)
+ amatrix(i,3) = yCellPlanar(i-1)
- amatrix(i,4) = xp(i-1)**2
- amatrix(i,5) = xp(i-1) * yp(i-1)
- amatrix(i,6) = yp(i-1)**2
+ amatrix(i,4) = xCellPlanar(i-1)**2
+ amatrix(i,5) = xCellPlanar(i-1) * yCellPlanar(i-1)
+ amatrix(i,6) = yCellPlanar(i-1)**2
- amatrix(i,7) = xp(i-1)**3
- amatrix(i,8) = yp(i-1) * (xp(i-1)**2)
- amatrix(i,9) = xp(i-1) * (yp(i-1)**2)
- amatrix(i,10) = yp(i-1)**3
+ amatrix(i,7) = xCellPlanar(i-1)**3
+ amatrix(i,8) = yCellPlanar(i-1) * (xCellPlanar(i-1)**2)
+ amatrix(i,9) = xCellPlanar(i-1) * (yCellPlanar(i-1)**2)
+ amatrix(i,10) = yCellPlanar(i-1)**3
- amatrix(i,11) = xp(i-1)**4
- amatrix(i,12) = yp(i-1) * (xp(i-1)**3)
- amatrix(i,13) = (xp(i-1)**2)*(yp(i-1)**2)
- amatrix(i,14) = xp(i-1) * (yp(i-1)**3)
- amatrix(i,15) = yp(i-1)**4
+ amatrix(i,11) = xCellPlanar(i-1)**4
+ amatrix(i,12) = yCellPlanar(i-1) * (xCellPlanar(i-1)**3)
+ amatrix(i,13) = (xCellPlanar(i-1)**2)*(yCellPlanar(i-1)**2)
+ amatrix(i,14) = xCellPlanar(i-1) * (yCellPlanar(i-1)**3)
+ amatrix(i,15) = yCellPlanar(i-1)**4
wmatrix(i,i) = 1.
@@ -267,45 +255,43 @@
end if
+ ! Solve matrix equation for second derivative matrix
call ocn_poly_fit_2( amatrix, bmatrix, wmatrix, ma, na, 25 )
+ ! Compute edge location, as midpoint between two vertices.
+ ! Also compute angle the edge's normal makes with local east.
do i=1,grid % nEdgesOnCell % array (iCell)
ip1 = i+1
- if (ip1 > n-1) ip1 = 1
+ if (ip1 > nAdvectionCells-1) ip1 = 1
iEdge = grid % EdgesOnCell % array (i,iCell)
- xv1 = grid % xVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- yv1 = grid % yVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- zv1 = grid % zVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- xv2 = grid % xVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/a
+ xVertexScaled1 = grid % xVertex % array(grid % verticesOnEdge % array (1,iedge))/grid % sphere_radius
+ yVertexScaled1 = grid % yVertex % array(grid % verticesOnEdge % array (1,iedge))/grid % sphere_radius
+ zVertexScaled1 = grid % zVertex % array(grid % verticesOnEdge % array (1,iedge))/grid % sphere_radius
+ xVertexScaled2 = grid % xVertex % array(grid % verticesOnEdge % array (2,iedge))/grid % sphere_radius
+ yVertexScaled2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/grid % sphere_radius
+ zVertexScaled2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/grid % sphere_radius
if ( grid % on_a_sphere ) then
- call ocn_arc_bisect( xv1, yv1, zv1, &
- xv2, yv2, zv2, &
- xec, yec, zec )
+ call ocn_arc_bisect( xVertexScaled1, yVertexScaled1, zVertexScaled1, &
+ xVertexScaled2, yVertexScaled2, zVertexScaled2, &
+ xEdgeMidpoint, yEdgeMidpoint, zEdgeMidpoint )
- thetae_tmp = sphere_angle( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1), &
- xec, yec, zec )
- thetae_tmp = thetae_tmp + thetat(i)
+ edgeAngleToEast_tmp = sphere_angle( xCellScaled(1), yCellScaled(1), zCellScaled(1), &
+ xCellScaled(i+1), yCellScaled(i+1), zCellScaled(i+1), &
+ xEdgeMidpoint, yEdgeMidpoint, zEdgeMidpoint )
+ edgeAngleToEast_tmp = edgeAngleToEast_tmp + angleToEast(i)
if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
- thetae(1,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ edgeAngleToEast(1,grid % EdgesOnCell % array (i,iCell)) = edgeAngleToEast_tmp
else
- thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ edgeAngleToEast(2,grid % EdgesOnCell % array (i,iCell)) = edgeAngleToEast_tmp
end if
-! else
-!
-! xe(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (xv1 + xv2)
-! ye(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (yv1 + yv2)
end if
end do
-! fill second derivative stencil for rk advection
-
+ ! Build deriv_two array using bmatrix and the edge's angle to east.
do i=1, grid % nEdgesOnCell % array (iCell)
iEdge = grid % EdgesOnCell % array (i,iCell)
@@ -313,26 +299,26 @@
if ( grid % on_a_sphere ) then
if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
- cos2t = cos(thetae(1,grid % EdgesOnCell % array (i,iCell)))
- sin2t = sin(thetae(1,grid % EdgesOnCell % array (i,iCell)))
+ cos2t = cos(edgeAngleToEast(1,grid % EdgesOnCell % array (i,iCell)))
+ sin2t = sin(edgeAngleToEast(1,grid % EdgesOnCell % array (i,iCell)))
costsint = cos2t*sin2t
cos2t = cos2t**2
sin2t = sin2t**2
- do j=1,n
+ do j=1,nAdvectionCells
deriv_two(j,1,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ 2.*costsint*bmatrix(5,j) &
+ 2.*sin2t*bmatrix(6,j)
end do
else
- cos2t = cos(thetae(2,grid % EdgesOnCell % array (i,iCell)))
- sin2t = sin(thetae(2,grid % EdgesOnCell % array (i,iCell)))
+ cos2t = cos(edgeAngleToEast(2,grid % EdgesOnCell % array (i,iCell)))
+ sin2t = sin(edgeAngleToEast(2,grid % EdgesOnCell % array (i,iCell)))
costsint = cos2t*sin2t
cos2t = cos2t**2
sin2t = sin2t**2
- do j=1,n
+ do j=1,nAdvectionCells
deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ 2.*costsint*bmatrix(5,j) &
+ 2.*sin2t*bmatrix(6,j)
@@ -341,27 +327,20 @@
else
- cos2t = cos(angle_2d(i))
- sin2t = sin(angle_2d(i))
+ cos2t = cos(angleToEast(i))
+ sin2t = sin(angleToEast(i))
costsint = cos2t*sin2t
cos2t = cos2t**2
sin2t = sin2t**2
-! do j=1,n
-!
-! deriv_two(j,1,iEdge) = 2.*xe(iEdge)*xe(iEdge)*bmatrix(4,j) &
-! + 2.*xe(iEdge)*ye(iEdge)*bmatrix(5,j) &
-! + 2.*ye(iEdge)*ye(iEdge)*bmatrix(6,j)
-! end do
-
if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
- do j=1,n
+ do j=1,nAdvectionCells
deriv_two(j,1,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ 2.*costsint*bmatrix(5,j) &
+ 2.*sin2t*bmatrix(6,j)
end do
else
- do j=1,n
+ do j=1,nAdvectionCells
deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ 2.*costsint*bmatrix(5,j) &
+ 2.*sin2t*bmatrix(6,j)
@@ -371,29 +350,10 @@
end if
end do
- end do ! end of loop over cells
+ end do ! end iCell loop
- if (debug) stop
+ deallocate(advCells)
-
-! write(0,*) ' check for deriv2 coefficients, iEdge 4 '
-!
-! iEdge = 4
-! j = 1
-! iCell = grid % cellsOnEdge % array(1,iEdge)
-! write(0,*) ' j, icell, coef ',j,iCell,deriv_two(j,1,iEdge)
-! do j=2,7
-! write(0,*) ' j, icell, coef ',j,grid % CellsOnCell % array(j-1,iCell),deriv_two(j,1,iEdge)
-! end do
-!
-! j = 1
-! iCell = grid % cellsOnEdge % array(2,iEdge)
-! write(0,*) ' j, icell, coef ',j,iCell,deriv_two(j,2,iEdge)
-! do j=2,7
-! write(0,*) ' j, icell, coef ',j,grid % CellsOnCell % array(j-1,iCell),deriv_two(j,2,iEdge)
-! end do
-! stop
-
end subroutine ocn_initialize_advection_rk!}}}
</font>
</pre>