<p><b>duda</b> 2009-08-03 12:33:16 -0600 (Mon, 03 Aug 2009)</p><p>Compute u field from streamfunction sampled at vertices to ensure <br>
zero divergence in initial conditions for test cases 1, 2, 5.<br>
The initial conditions for test case 6 still need to be corrected,<br>
since the initial height field appears to be incorrect.<br>
<br>
M module_test_cases.F<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/module_test_cases.F
===================================================================
--- trunk/swmodel/module_test_cases.F 2009-08-02 01:28:52 UTC (rev 14)
+++ trunk/swmodel/module_test_cases.F 2009-08-03 18:33:16 UTC (rev 15)
@@ -108,8 +108,9 @@
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
real (kind=RKIND), parameter :: alpha = pii/4.0
- integer :: iCell, iEdge
+ integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: r, u, v
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex
!
! Scale all distances and areas from a unit sphere to one with radius a
@@ -132,19 +133,22 @@
!
! Initialize wind field
!
+ allocate(psiVertex(grid % nVertices))
+ do iVtx=1,grid % nVertices
+ psiVertex(iVtx) = -a * u0 * ( &
+ sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
+ cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
+ )
+ end do
do iEdge=1,grid % nEdges
- u = cos(grid%latEdge%array(iEdge)) * cos(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(grid%lonEdge%array(iEdge)) * sin(alpha)
- v = -sin(grid%lonEdge%array(iEdge)) * sin(alpha)
- state % u % array(1,iEdge) = u * cos(grid%angleEdge%array(iEdge)) + v * sin(grid%angleEdge%array(iEdge))
+ state % u % array(1,iEdge) = -1.0 * ( &
+ psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
+ psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
+ ) / grid%dvEdge%array(iEdge)
end do
+ deallocate(psiVertex)
!
- ! Scale wind field from 1.0 m/s at equator for solid-body rotation
- !
- state % u % array = state % u % array * u0
-
- !
! Initialize cosine bell at (theta_c, lambda_c)
!
do iCell=1,grid % nCells
@@ -180,6 +184,7 @@
integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: u, v
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex
!
@@ -199,23 +204,27 @@
grid % areaCell % array = grid % areaCell % array * a**2.0
grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
+
!
! Initialize wind field
!
+ allocate(psiVertex(grid % nVertices))
+ do iVtx=1,grid % nVertices
+ psiVertex(iVtx) = -a * u0 * ( &
+ sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
+ cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
+ )
+ end do
do iEdge=1,grid % nEdges
- u = cos(grid%latEdge%array(iEdge)) * cos(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(grid%lonEdge%array(iEdge)) * sin(alpha)
- v = -sin(grid%lonEdge%array(iEdge)) * sin(alpha)
- state % u % array(1,iEdge) = u * cos(grid%angleEdge%array(iEdge)) + v * sin(grid%angleEdge%array(iEdge))
+ state % u % array(1,iEdge) = -1.0 * ( &
+ psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
+ psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
+ ) / grid%dvEdge%array(iEdge)
end do
+ deallocate(psiVertex)
!
- ! Scale wind field from 1.0 m/s at equator for solid-body rotation
- !
- state % u % array = state % u % array * u0
-
- !
! Generate rotated Coriolis field
!
do iEdge=1,grid % nEdges
@@ -270,6 +279,7 @@
integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: r, u, v
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex
!
@@ -293,19 +303,22 @@
!
! Initialize wind field
!
+ allocate(psiVertex(grid % nVertices))
+ do iVtx=1,grid % nVertices
+ psiVertex(iVtx) = -a * u0 * ( &
+ sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
+ cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
+ )
+ end do
do iEdge=1,grid % nEdges
- u = cos(grid%latEdge%array(iEdge)) * cos(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(grid%lonEdge%array(iEdge)) * sin(alpha)
- v = -sin(grid%lonEdge%array(iEdge)) * sin(alpha)
- state % u % array(1,iEdge) = u * cos(grid%angleEdge%array(iEdge)) + v * sin(grid%angleEdge%array(iEdge))
+ state % u % array(1,iEdge) = -1.0 * ( &
+ psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
+ psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
+ ) / grid%dvEdge%array(iEdge)
end do
+ deallocate(psiVertex)
!
- ! Scale wind field from 1.0 m/s at equator for solid-body rotation
- !
- state % u % array = state % u % array * u0
-
- !
! Generate rotated Coriolis field
!
do iEdge=1,grid % nEdges
@@ -379,8 +392,9 @@
real (kind=RKIND), parameter :: K = 7.848e-6
real (kind=RKIND), parameter :: R = 4.0
- integer :: iCell, iEdge
+ integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: u, v
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex
!
@@ -404,25 +418,27 @@
!
! Initialize wind field
!
+ allocate(psiVertex(grid % nVertices))
+ do iVtx=1,grid % nVertices
+ psiVertex(iVtx) = -a * a * w * sin(grid%latVertex%array(iVtx)) + &
+ a *a * K * (cos(grid%latVertex%array(iVtx))**R) * sin(grid%latVertex%array(iVtx)) * cos(R * grid%lonVertex%array(iVtx))
+ end do
do iEdge=1,grid % nEdges
- u = a * w * cos(grid%latEdge%array(iEdge)) + &
- a * K * cos(grid%latEdge%array(iEdge))**(R-1.0) * &
- (R * sin(grid%latEdge%array(iEdge))**2.0 - cos(grid%latEdge%array(iEdge))**2.0) * &
- cos(R*grid%lonEdge%array(iEdge))
- v = -a * K * R * cos(grid%latEdge%array(iEdge))**(R-1.0) * sin(grid%latEdge%array(iEdge)) * &
- sin(R*grid%lonEdge%array(iEdge))
- state % u % array(1,iEdge) = u * cos(grid%angleEdge%array(iEdge)) + v * sin(grid%angleEdge%array(iEdge))
+ state % u % array(1,iEdge) = -1.0 * ( &
+ psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
+ psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
+ ) / grid%dvEdge%array(iEdge)
end do
+ deallocate(psiVertex)
!
! Initialize height field (actually, fluid thickness field)
!
do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gravity * h0 + a**2.0*AA(grid%latCell%array(iCell)) + &
- a**2.0*BB(grid%latCell%array(iCell))*cos(R*grid%lonCell%array(iCell)) + &
- a**2.0*CC(grid%latCell%array(iCell))*cos(2.0*R*grid%lonCell%array(iCell))&
- ) / &
- gravity
+ state % h % array(1,iCell) = (gravity * h0 + a*a*AA(grid%latCell%array(iCell)) + &
+ a*a*BB(grid%latCell%array(iCell)) * cos(R*grid%lonCell%array(iCell)) + &
+ a*a*CC(grid%latCell%array(iCell)) * cos(2.0*R*grid%lonCell%array(iCell)) &
+ ) / gravity
end do
end subroutine sw_test_case_6
@@ -461,7 +477,7 @@
real (kind=RKIND), intent(in) :: theta
AA = 0.5 * w * (2.0 * omega + w) * cos(theta)**2.0 + &
- 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*R**2.0*cos(theta)**2.0)
+ 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*(R*cos(theta))**2.0)
end function AA
@@ -479,7 +495,7 @@
real (kind=RKIND), intent(in) :: theta
- BB = 2.0*(omega + w)*K / ((R+1.0)*(R+2.0)) * cos(theta)**R * ((R**2.0 + 2.0*R + 2.0) - (R+1.0)**2.0*cos(theta)**2.0)
+ BB = (2.0*(omega + w)*K / ((R+1.0)*(R+2.0))) * cos(theta)**R * ((R**2.0 + 2.0*R + 2.0) - ((R+1.0)*cos(theta))**2.0)
end function BB
</font>
</pre>