<p><b>dwj07@fsu.edu</b> 2011-10-17 09:44:46 -0600 (Mon, 17 Oct 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Performing renaming of core_sw<br>
</p><hr noshade><pre><font color="gray">Modified: branches/source_renaming/src/core_sw/Makefile
===================================================================
--- branches/source_renaming/src/core_sw/Makefile        2011-10-17 15:40:02 UTC (rev 1086)
+++ branches/source_renaming/src/core_sw/Makefile        2011-10-17 15:44:46 UTC (rev 1087)
@@ -1,25 +1,25 @@
.SUFFIXES: .F .o
-OBJS =         module_mpas_core.o \
- module_test_cases.o \
-        module_advection.o \
-        module_time_integration.o \
-        module_global_diagnostics.o
+OBJS =         mpas_sw_mpas_core.o \
+ mpas_sw_test_cases.o \
+        mpas_sw_advection.o \
+        mpas_sw_time_integration.o \
+        mpas_sw_global_diagnostics.o
all: core_sw
core_sw: $(OBJS)
        ar -ru libdycore.a $(OBJS)
-module_test_cases.o:
+mpas_sw_test_cases.o:
-module_advection.o:
+mpas_sw_advection.o:
-module_time_integration.o:
+mpas_sw_time_integration.o:
-module_global_diagnostics.o:
+mpas_sw_global_diagnostics.o:
-module_mpas_core.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o module_advection.o
+mpas_sw_mpas_core.o: mpas_sw_global_diagnostics.o mpas_sw_test_cases.o mpas_sw_time_integration.o mpas_sw_advection.o
clean:
        $(RM) *.o *.mod *.f90 libdycore.a
Deleted: branches/source_renaming/src/core_sw/module_advection.F
===================================================================
--- branches/source_renaming/src/core_sw/module_advection.F        2011-10-17 15:40:02 UTC (rev 1086)
+++ branches/source_renaming/src/core_sw/module_advection.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -1,933 +0,0 @@
-module advection
-
- use grid_types
- use configure
- use constants
-
-
- contains
-
-
- subroutine initialize_advection_rk( grid )
-
-!
-! compute the cell coefficients for the polynomial fit.
-! this is performed during setup for model integration.
-! WCS, 31 August 2009
-!
- implicit none
-
- type (mesh_type), intent(in) :: grid
-
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- integer, dimension(:,:), pointer :: advCells
-
-! local variables
-
- real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
- 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
- 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) :: amatrix(25,25), bmatrix(25,25), wmatrix(25,25)
- real (kind=RKIND) :: length_scale
- integer :: ma,na, cell_add, mw, nn
- integer, dimension(25) :: cell_list
-
-
- integer :: cell1, cell2
- integer, parameter :: polynomial_order = 2
-! logical, parameter :: debug = .true.
- logical, parameter :: debug = .false.
-! logical, parameter :: least_squares = .false.
- logical, parameter :: least_squares = .true.
- logical :: add_the_cell, do_the_cell
-
- logical, parameter :: reset_poly = .true.
-
- real (kind=RKIND) :: rcell, cos2t, costsint, sin2t
- real (kind=RKIND), dimension(grid%maxEdges) :: angle_2d
-
-!---
-
- pii = 2.*asin(1.0)
-
- advCells => grid % advCells % array
- deriv_two => grid % deriv_two % array
- deriv_two(:,:,:) = 0.
-
- do iCell = 1, grid % nCells ! is this correct? - we need first halo cell also...
-
- 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
-
- 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
- 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
- end if
- end do
- end do
- end if
-
- advCells(1,iCell) = n
-
-! check to see if we are reaching outside the halo
-
- do_the_cell = .true.
- do i=1,n
- if (cell_list(i) > grid % nCells) do_the_cell = .false.
- end do
-
-
- 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
- advCells(i+1,iCell) = cell_list(i)
- xc(i) = grid % xCell % array(advCells(i+1,iCell))/a
- yc(i) = grid % yCell % array(advCells(i+1,iCell))/a
- zc(i) = grid % zCell % array(advCells(i+1,iCell))/a
- end do
-
- theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
- xc(2), yc(2), zc(2), &
- 0., 0., 1. )
-
-! angles from cell center to neighbor centers (thetav)
-
- do i=1,n-1
-
- ip2 = i+2
- if (ip2 > n) 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) )
-
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1) )
- end do
-
- length_scale = 1.
- do i=1,n-1
- dl_sphere(i) = dl_sphere(i)/length_scale
- 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)
- end do
-
- else ! On an x-y plane
-
- do i=1,n-1
-
- angle_2d(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
-
-! xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
-! yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
-
- 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
- mw = grid % nEdgesOnCell % array (iCell)
-
- bmatrix = 0.
- amatrix = 0.
- wmatrix = 0.
-
- if (polynomial_order == 2) then
- na = 6
- ma = ma+1
-
- amatrix(1,1) = 1.
- 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
-
- wmatrix(i,i) = 1.
- end do
-
- else if (polynomial_order == 3) then
- na = 10
- ma = ma+1
-
- amatrix(1,1) = 1.
- 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,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
-
- wmatrix(i,i) = 1.
-
- end do
-
- else
- na = 15
- ma = ma+1
-
- amatrix(1,1) = 1.
- 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,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,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
-
- wmatrix(i,i) = 1.
-
- end do
-
- do i=1,mw
- wmatrix(i,i) = 1.
- end do
-
- end if
-
- call poly_fit_2( amatrix, bmatrix, wmatrix, ma, na, 25 )
-
- do i=1,grid % nEdgesOnCell % array (iCell)
- ip1 = i+1
- if (ip1 > n-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
-
- if ( grid % on_a_sphere ) then
- call arc_bisect( xv1, yv1, zv1, &
- xv2, yv2, zv2, &
- xec, yec, zec )
-
- 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)
- if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
- thetae(1,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
- else
- thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_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
-
- do i=1, grid % nEdgesOnCell % array (iCell)
- iEdge = grid % EdgesOnCell % array (i,iCell)
-
-
- 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)))
- costsint = cos2t*sin2t
- cos2t = cos2t**2
- sin2t = sin2t**2
-
- do j=1,n
- 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)))
- costsint = cos2t*sin2t
- cos2t = cos2t**2
- sin2t = sin2t**2
-
- do j=1,n
- deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
- + 2.*costsint*bmatrix(5,j) &
- + 2.*sin2t*bmatrix(6,j)
- end do
- end if
-
- else
-
- cos2t = cos(angle_2d(i))
- sin2t = sin(angle_2d(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
- 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
- deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
- + 2.*costsint*bmatrix(5,j) &
- + 2.*sin2t*bmatrix(6,j)
- end do
- end if
-
- end if
- end do
-
- end do ! end of loop over cells
-
- if (debug) stop
-
-
-! 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 initialize_advection_rk
-
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! FUNCTION SPHERE_ANGLE
- !
- ! Computes the angle between arcs AB and AC, given points A, B, and C
- ! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
-
- implicit none
-
- real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz
-
- real (kind=RKIND) :: a, b, c ! Side lengths of spherical triangle ABC
-
- real (kind=RKIND) :: ABx, ABy, ABz ! The components of the vector AB
- real (kind=RKIND) :: mAB ! The magnitude of AB
- real (kind=RKIND) :: ACx, ACy, ACz ! The components of the vector AC
- real (kind=RKIND) :: mAC ! The magnitude of AC
-
- real (kind=RKIND) :: Dx ! The i-components of the cross product AB x AC
- real (kind=RKIND) :: Dy ! The j-components of the cross product AB x AC
- real (kind=RKIND) :: Dz ! The k-components of the cross product AB x AC
-
- real (kind=RKIND) :: s ! Semiperimeter of the triangle
- real (kind=RKIND) :: sin_angle
-
- a = acos(max(min(bx*cx + by*cy + bz*cz,1.0),-1.0)) ! Eqn. (3)
- b = acos(max(min(ax*cx + ay*cy + az*cz,1.0),-1.0)) ! Eqn. (2)
- c = acos(max(min(ax*bx + ay*by + az*bz,1.0),-1.0)) ! Eqn. (1)
-
- ABx = bx - ax
- ABy = by - ay
- ABz = bz - az
-
- ACx = cx - ax
- ACy = cy - ay
- ACz = cz - az
-
- Dx = (ABy * ACz) - (ABz * ACy)
- Dy = -((ABx * ACz) - (ABz * ACx))
- Dz = (ABx * ACy) - (ABy * ACx)
-
- s = 0.5*(a + b + c)
-! sin_angle = sqrt((sin(s-b)*sin(s-c))/(sin(b)*sin(c))) ! Eqn. (28)
- sin_angle = sqrt(min(1.,max(0.,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
-
- if ((Dx*ax + Dy*ay + Dz*az) >= 0.0) then
- sphere_angle = 2.0 * asin(max(min(sin_angle,1.0),-1.0))
- else
- sphere_angle = -2.0 * asin(max(min(sin_angle,1.0),-1.0))
- end if
-
- end function sphere_angle
-
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! FUNCTION PLANE_ANGLE
- !
- ! Computes the angle between vectors AB and AC, given points A, B, and C, and
- ! a vector (u,v,w) normal to the plane.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
-
- implicit none
-
- real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w
-
- real (kind=RKIND) :: ABx, ABy, ABz ! The components of the vector AB
- real (kind=RKIND) :: mAB ! The magnitude of AB
- real (kind=RKIND) :: ACx, ACy, ACz ! The components of the vector AC
- real (kind=RKIND) :: mAC ! The magnitude of AC
-
- real (kind=RKIND) :: Dx ! The i-components of the cross product AB x AC
- real (kind=RKIND) :: Dy ! The j-components of the cross product AB x AC
- real (kind=RKIND) :: Dz ! The k-components of the cross product AB x AC
-
- real (kind=RKIND) :: cos_angle
-
- ABx = bx - ax
- ABy = by - ay
- ABz = bz - az
- mAB = sqrt(ABx**2.0 + ABy**2.0 + ABz**2.0)
-
- ACx = cx - ax
- ACy = cy - ay
- ACz = cz - az
- mAC = sqrt(ACx**2.0 + ACy**2.0 + ACz**2.0)
-
-
- Dx = (ABy * ACz) - (ABz * ACy)
- Dy = -((ABx * ACz) - (ABz * ACx))
- Dz = (ABx * ACy) - (ABy * ACx)
-
- cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
-
- if ((Dx*u + Dy*v + Dz*w) >= 0.0) then
- plane_angle = acos(max(min(cos_angle,1.0),-1.0))
- else
- plane_angle = -acos(max(min(cos_angle,1.0),-1.0))
- end if
-
- end function plane_angle
-
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! FUNCTION ARC_LENGTH
- !
- ! Returns the length of the great circle arc from A=(ax, ay, az) to
- ! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
- ! same sphere centered at the origin.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real function arc_length(ax, ay, az, bx, by, bz)
-
- implicit none
-
- real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz
-
- real (kind=RKIND) :: r, c
- real (kind=RKIND) :: cx, cy, cz
-
- cx = bx - ax
- cy = by - ay
- cz = bz - az
-
-! r = ax*ax + ay*ay + az*az
-! c = cx*cx + cy*cy + cz*cz
-!
-! arc_length = sqrt(r) * acos(1.0 - c/(2.0*r))
-
- r = sqrt(ax*ax + ay*ay + az*az)
- c = sqrt(cx*cx + cy*cy + cz*cz)
-! arc_length = sqrt(r) * 2.0 * asin(c/(2.0*r))
- arc_length = r * 2.0 * asin(c/(2.0*r))
-
- end function arc_length
-
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! SUBROUTINE ARC_BISECT
- !
- ! Returns the point C=(cx, cy, cz) that bisects the great circle arc from
- ! A=(ax, ay, az) to B=(bx, by, bz). It is assumed that A and B lie on the
- ! surface of a sphere centered at the origin.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine arc_bisect(ax, ay, az, bx, by, bz, cx, cy, cz)
-
- implicit none
-
- real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz
- real (kind=RKIND), intent(out) :: cx, cy, cz
-
- real (kind=RKIND) :: r ! Radius of the sphere
- real (kind=RKIND) :: d
-
- r = sqrt(ax*ax + ay*ay + az*az)
-
- cx = 0.5*(ax + bx)
- cy = 0.5*(ay + by)
- cz = 0.5*(az + bz)
-
- if (cx == 0. .and. cy == 0. .and. cz == 0.) then
- write(0,*) 'Error: arc_bisect: A and B are diametrically opposite'
- else
- d = sqrt(cx*cx + cy*cy + cz*cz)
- cx = r * cx / d
- cy = r * cy / d
- cz = r * cz / d
- end if
-
- end subroutine arc_bisect
-
-
- subroutine poly_fit_2(a_in,b_out,weights_in,m,n,ne)
-
- implicit none
-
- integer, intent(in) :: m,n,ne
- real (kind=RKIND), dimension(ne,ne), intent(in) :: a_in, weights_in
- real (kind=RKIND), dimension(ne,ne), intent(out) :: b_out
-
- ! local storage
-
- real (kind=RKIND), dimension(m,n) :: a
- real (kind=RKIND), dimension(n,m) :: b
- real (kind=RKIND), dimension(m,m) :: w,wt,h
- real (kind=RKIND), dimension(n,m) :: at, ath
- real (kind=RKIND), dimension(n,n) :: ata, ata_inv, atha, atha_inv
- integer, dimension(n) :: indx
- integer :: i,j
-
- if ( (ne<n) .or. (ne<m) ) then
- write(6,*) ' error in poly_fit_2 inversion ',m,n,ne
- stop
- end if
-
-! a(1:m,1:n) = a_in(1:n,1:m)
- a(1:m,1:n) = a_in(1:m,1:n)
- w(1:m,1:m) = weights_in(1:m,1:m)
- b_out(:,:) = 0.
-
- wt = transpose(w)
- h = matmul(wt,w)
- at = transpose(a)
- ath = matmul(at,h)
- atha = matmul(ath,a)
-
- ata = matmul(at,a)
-
-! if (m == n) then
-! call migs(a,n,b,indx)
-! else
-
- call migs(atha,n,atha_inv,indx)
-
- b = matmul(atha_inv,ath)
-
-! call migs(ata,n,ata_inv,indx)
-! b = matmul(ata_inv,at)
-! end if
- b_out(1:n,1:m) = b(1:n,1:m)
-
-! do i=1,n
-! write(6,*) ' i, indx ',i,indx(i)
-! end do
-!
-! write(6,*) ' '
-
- end subroutine poly_fit_2
-
-
-! Updated 10/24/2001.
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! !
-! Please Note: !
-! !
-! (1) This computer program is written by Tao Pang in conjunction with !
-! his book, "An Introduction to Computational Physics," published !
-! by Cambridge University Press in 1997. !
-! !
-! (2) No warranties, express or implied, are made for this program. !
-! !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-SUBROUTINE MIGS (A,N,X,INDX)
-!
-! Subroutine to invert matrix A(N,N) with the inverse stored
-! in X(N,N) in the output. Copyright (c) Tao Pang 2001.
-!
- IMPLICIT NONE
- INTEGER, INTENT (IN) :: N
- INTEGER :: I,J,K
- INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
- REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
- REAL (kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
- REAL (kind=RKIND), DIMENSION (N,N) :: B
-!
- DO I = 1, N
- DO J = 1, N
- B(I,J) = 0.0
- END DO
- END DO
- DO I = 1, N
- B(I,I) = 1.0
- END DO
-!
- CALL ELGS (A,N,INDX)
-!
- DO I = 1, N-1
- DO J = I+1, N
- DO K = 1, N
- B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K)
- END DO
- END DO
- END DO
-!
- DO I = 1, N
- X(N,I) = B(INDX(N),I)/A(INDX(N),N)
- DO J = N-1, 1, -1
- X(J,I) = B(INDX(J),I)
- DO K = J+1, N
- X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
- END DO
- X(J,I) = X(J,I)/A(INDX(J),J)
- END DO
- END DO
-END SUBROUTINE MIGS
-
-
-SUBROUTINE ELGS (A,N,INDX)
-!
-! Subroutine to perform the partial-pivoting Gaussian elimination.
-! A(N,N) is the original matrix in the input and transformed matrix
-! plus the pivoting element ratios below the diagonal in the output.
-! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001.
-!
- IMPLICIT NONE
- INTEGER, INTENT (IN) :: N
- INTEGER :: I,J,K,ITMP
- INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
- REAL (kind=RKIND) :: C1,PI,PI1,PJ
- REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
- REAL (kind=RKIND), DIMENSION (N) :: C
-!
-! Initialize the index
-!
- DO I = 1, N
- INDX(I) = I
- END DO
-!
-! Find the rescaling factors, one from each row
-!
- DO I = 1, N
- C1= 0.0
- DO J = 1, N
- C1 = MAX(C1,ABS(A(I,J)))
- END DO
- C(I) = C1
- END DO
-!
-! Search the pivoting (largest) element from each column
-!
- DO J = 1, N-1
- PI1 = 0.0
- DO I = J, N
- PI = ABS(A(INDX(I),J))/C(INDX(I))
- IF (PI.GT.PI1) THEN
- PI1 = PI
- K = I
- ENDIF
- END DO
-!
-! Interchange the rows via INDX(N) to record pivoting order
-!
- ITMP = INDX(J)
- INDX(J) = INDX(K)
- INDX(K) = ITMP
- DO I = J+1, N
- PJ = A(INDX(I),J)/A(INDX(J),J)
-!
-! Record pivoting ratios below the diagonal
-!
- A(INDX(I),J) = PJ
-!
-! Modify other elements accordingly
-!
- DO K = J+1, N
- A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
- END DO
- END DO
- END DO
-!
-END SUBROUTINE ELGS
-
-!-------------------------------------------------------------
-
- subroutine initialize_deformation_weights( grid )
-
-!
-! compute the cell coefficients for the deformation calculations
-! WCS, 13 July 2010
-!
- implicit none
-
- type (mesh_type), intent(in) :: grid
-
- real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
- integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
-
-! local variables
-
- real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
- 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
- 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, xpt, ypt
-
- real (kind=RKIND) :: length_scale
- integer :: ma,na, cell_add, mw, nn
- integer, dimension(25) :: cell_list
-
- integer :: cell1, cell2, iv
- logical :: do_the_cell
- real (kind=RKIND) :: area_cell, sint2, cost2, sint_cost, sumw1, sumw2, xptt, area_cellt
-
- logical, parameter :: debug = .false.
-
- if (debug) write(0,*) ' in def weight calc '
-
- defc_a => grid % defc_a % array
- defc_b => grid % defc_b % array
- cellsOnEdge => grid % cellsOnEdge % array
- edgesOnCell => grid % edgesOnCell % array
-
- defc_a(:,:) = 0.
- defc_b(:,:) = 0.
-
- pii = 2.*asin(1.0)
-
- if (debug) write(0,*) ' beginning cell loop '
-
- do iCell = 1, grid % nCells
-
- if (debug) write(0,*) ' cell loop ', iCell
-
- 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
-
-! check to see if we are reaching outside the halo
-
- if (debug) write(0,*) ' points ', n
-
- do_the_cell = .true.
- do i=1,n
- if (cell_list(i) > grid % nCells) do_the_cell = .false.
- end do
-
-
- if (.not. do_the_cell) cycle
-
-
-! compute poynomial fit for this cell if all needed neighbors exist
- if (grid % on_a_sphere) then
-
- xc(1) = grid % xCell % array(iCell)/a
- yc(1) = grid % yCell % array(iCell)/a
- zc(1) = grid % zCell % array(iCell)/a
-
-
- do i=2,n
- iv = grid % verticesOnCell % array(i-1,iCell)
- xc(i) = grid % xVertex % array(iv)/a
- yc(i) = grid % yVertex % array(iv)/a
- zc(i) = grid % zVertex % array(iv)/a
- end do
-
- theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
- xc(2), yc(2), zc(2), &
- 0., 0., 1. )
-
-! angles from cell center to neighbor centers (thetav)
-
- do i=1,n-1
-
- ip2 = i+2
- if (ip2 > n) 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) )
-
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
- xc(i+1), yc(i+1), zc(i+1) )
- end do
-
- length_scale = 1.
- do i=1,n-1
- dl_sphere(i) = dl_sphere(i)/length_scale
- 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)
- end do
-
- else ! On an x-y plane
-
- xp(1) = grid % xCell % array(iCell)
- yp(1) = grid % yCell % array(iCell)
-
-
- do i=2,n
- iv = grid % verticesOnCell % array(i-1,iCell)
- xp(i) = grid % xVertex % array(iv)
- yp(i) = grid % yVertex % array(iv)
- end do
-
- end if
-
-! thetat(1) = 0.
- thetat(1) = theta_abs(iCell)
- do i=2,n-1
- ip1 = i+1
- if (ip1 == n) ip1 = 1
- thetat(i) = plane_angle( 0.,0.,0., &
- xp(i)-xp(i-1), yp(i)-yp(i-1), 0., &
- xp(ip1)-xp(i), yp(ip1)-yp(i), 0., &
- 0., 0., 1.)
- thetat(i) = thetat(i) + thetat(i-1)
- end do
-
- area_cell = 0.
- area_cellt = 0.
- do i=1,n-1
- ip1 = i+1
- if (ip1 == n) ip1 = 1
- dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
- area_cell = area_cell + 0.25*(xp(i)+xp(ip1))*(yp(ip1)-yp(i)) - 0.25*(yp(i)+yp(ip1))*(xp(ip1)-xp(i))
- area_cellt = area_cellt + (0.25*(xp(i)+xp(ip1))*cos(thetat(i)) + 0.25*(yp(i)+yp(ip1))*sin(thetat(i)))*dl
- end do
- if (debug) write(0,*) ' area_cell, area_cellt ',area_cell, area_cellt,area_cell-area_cellt
-
- do i=1,n-1
- ip1 = i+1
- if (ip1 == n) ip1 = 1
- dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
- sint2 = (sin(thetat(i)))**2
- cost2 = (cos(thetat(i)))**2
- sint_cost = sin(thetat(i))*cos(thetat(i))
- defc_a(i,iCell) = dl*(cost2 - sint2)/area_cell
- defc_b(i,iCell) = dl*2.*sint_cost/area_cell
- if (cellsOnEdge(1,EdgesOnCell(i,iCell)) /= iCell) then
- defc_a(i,iCell) = - defc_a(i,iCell)
- defc_b(i,iCell) = - defc_b(i,iCell)
- end if
-
- end do
-
- end do
-
- if (debug) write(0,*) ' exiting def weight calc '
-
- end subroutine initialize_deformation_weights
-
-end module advection
Deleted: branches/source_renaming/src/core_sw/module_global_diagnostics.F
===================================================================
--- branches/source_renaming/src/core_sw/module_global_diagnostics.F        2011-10-17 15:40:02 UTC (rev 1086)
+++ branches/source_renaming/src/core_sw/module_global_diagnostics.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -1,384 +0,0 @@
-module global_diagnostics
-
- use grid_types
- use configure
- use constants
- use dmpar
-
- implicit none
- save
- public
-
- contains
-
- subroutine computeGlobalDiagnostics(dminfo, state, grid, timeIndex, dt)
-
- ! Note: this routine assumes that there is only one block per processor. No looping
- ! is preformed over blocks.
- ! dminfo is the domain info needed for global communication
- ! state contains the state variables needed to compute global diagnostics
- ! grid conains the meta data about the grid
- ! timeIndex is the current time step counter
- ! dt is the duration of each time step
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! INSTRUCTIONS !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! To add a new Diagnostic as a Global Stat, follow these steps.
- ! 1. Define the array to integrate, and the variable for the value above.
- ! 2. Allocate the array with the correct dimensions.
- ! 3. Fill the array with the data to be integrated.
- ! eg. GlobalFluidThickness = Sum(h dA)/Sum(dA), See below for array filling
- ! 4. Call Function to compute Global Stat that you want.
- ! 5. Finish computing the global stat/integral
- ! 6. Write out your global stat to the file
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (dm_info), intent(in) :: dminfo
- type (state_type), intent(inout) :: state
- type (mesh_type), intent(in) :: grid
- integer, intent(in) :: timeIndex
- real (kind=RKIND), intent(in) :: dt
-
- integer :: nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, nCellsGlobal, nEdgesGlobal, nVerticesGlobal, iTracer
- integer :: nCells
-
- ! Step 1
- ! 1. Define the array to integrate, and the variable for the value to be stored in after the integration
- real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, h_s, fCell, fEdge
- real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, pv_edge, pv_vertex, pv_cell, h_vertex, weightsOnEdge
-
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
-
- real (kind=RKIND), dimension(:), allocatable :: volumeWeightedPotentialEnergyReservoir, averageThickness
- real (kind=RKIND), dimension(:), allocatable :: potentialEnstrophyReservior, areaEdge, h_s_edge
-
- real (kind=RKIND), dimension(:,:), allocatable :: cellVolume, cellArea, volumeWeightedPotentialVorticity
- real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnstrophy, vertexVolume, volumeWeightedKineticEnergy
- real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnergy, volumeWeightedPotentialEnergyTopography
- real (kind=RKIND), dimension(:,:), allocatable :: keTend_CoriolisForce, keTend_PressureGradient
- real (kind=RKIND), dimension(:,:), allocatable ::peTend_DivThickness, refAreaWeightedSurfaceHeight, refAreaWeightedSurfaceHeight_edge
-
- real (kind=RKIND) :: sumCellVolume, sumCellArea, sumVertexVolume, sumrefAreaWeightedSurfaceHeight
-
- real (kind=RKIND) :: globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, globalEnergy
- real (kind=RKIND) :: globalCoriolisEnergyTendency, globalKEPETendency, globalPotentialEnstrophyReservoir
- real (kind=RKIND) :: globalKineticEnergy, globalPotentialEnergy, globalPotentialEnergyReservoir
- real (kind=RKIND) :: globalKineticEnergyTendency, globalPotentialEnergyTendency
- real (kind=RKIND) :: global_temp, workpv, q
- real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
-
- integer :: elementIndex, variableIndex, nVariables, nSums, nMaxes, nMins
- integer :: timeLevel, eoe, iLevel, iCell, iEdge, iVertex
- integer :: fileID, iCell1, iCell2, j
-
- integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgesOnEdge
- integer, dimension(:), pointer :: nEdgesOnEdge
-
- cellsOnEdge => grid % cellsOnEdge % array
- edgesOnCell => grid % edgesOnCell % array
-
- nVertLevels = grid % nVertLevels
- nCellsSolve = grid % nCellsSolve
- nEdgesSolve = grid % nEdgesSolve
- nVerticesSolve = grid % nVerticesSolve
- nCells = grid % nCells
-
- h_s => grid % h_s % array
- areaCell => grid % areaCell % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaTriangle => grid % areaTriangle % array
- fCell => grid % fCell % array
- fEdge => grid % fEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
-
- allocate(areaEdge(1:nEdgesSolve))
- areaEdge = dcEdge(1:nEdgesSolve)*dvEdge(1:nEdgesSolve)
- weightsOnEdge => grid % weightsOnEdge % array
-
- h => state % h % array
- u => state % u % array
- v => state % v % array
- tracers => state % tracers % array
- h_edge => state % h_edge % array
- h_vertex => state % h_vertex % array
- pv_edge => state % pv_edge % array
- pv_vertex => state % pv_vertex % array
- pv_cell => state % pv_cell % array
-
- ! Step 2
- ! 2. Allocate the array with the correct dimensions.
- allocate(cellVolume(nVertLevels,nCellsSolve))
- allocate(cellArea(nVertLevels,nCellsSolve))
- allocate(refAreaWeightedSurfaceHeight(nVertLevels,nCellsSolve))
- allocate(refAreaWeightedSurfaceHeight_edge(nVertLevels,nEdgesSolve))
- allocate(volumeWeightedPotentialVorticity(nVertLevels,nVerticesSolve))
- allocate(volumeWeightedPotentialEnstrophy(nVertLevels,nVerticesSolve))
- allocate(potentialEnstrophyReservior(nCellsSolve))
- allocate(vertexVolume(nVertLevels,nVerticesSolve))
- allocate(volumeWeightedKineticEnergy(nVertLevels,nEdgesSolve))
- allocate(volumeWeightedPotentialEnergy(nVertLevels,nCellsSolve))
- allocate(volumeWeightedPotentialEnergyTopography(nVertLevels,nCellsSolve))
- allocate(volumeWeightedPotentialEnergyReservoir(nCellsSolve))
- allocate(keTend_CoriolisForce(nVertLevels,nEdgesSolve))
- allocate(keTend_PressureGradient(nVertLevels,nEdgesSolve))
- allocate(peTend_DivThickness(nVertLevels,nCells))
-
- allocate(averageThickness(nCellsSolve))
-
- allocate(h_s_edge(nEdgesSOlve))
-
-
- cellVolume = 0
- refAreaWeightedSurfaceHeight = 0
- refAreaWeightedSurfaceHeight_edge = 0
- vertexVolume = 0
- cellArea = 0
- averageThickness = 0
- volumeWeightedPotentialVorticity = 0
- volumeWeightedPotentialEnstrophy = 0
- volumeWeightedKineticEnergy = 0
- volumeWeightedPotentialEnergy = 0
- volumeWeightedPotentialEnergyTopography = 0
- volumeWeightedPotentialEnergyReservoir = 0
- keTend_PressureGradient = 0
- peTend_DivThickness = 0
- keTend_CoriolisForce = 0
- h_s_edge = 0
-
- ! Build Arrays for Global Integrals
- ! Step 3
- ! 3. Fill the array with the data to be integrated.
- ! eg. GlobalFluidThickness = Sum(h dA)/Sum(dA), See below for array filling
- do iLevel = 1,nVertLevels
- ! eg. GlobalFluidThickness top (Sum( h dA)) = Sum(cellVolume)
- cellVolume(iLevel,:) = h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
- ! eg. GlobalFluidThickness bot (Sum(dA)) = Sum(cellArea)
- cellArea(iLevel,:) = areaCell(1:nCellsSolve)
- volumeWeightedPotentialVorticity(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
- *h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
- volumeWeightedPotentialEnstrophy(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
- *pv_vertex(iLevel,1:nVerticesSolve)*h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
- vertexVolume(iLevel,:) = h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
- volumeWeightedKineticEnergy(iLevel,:) = u(iLevel,1:nEdgesSolve)*u(iLevel,1:nEdgesSolve) &
- *h_edge(iLevel,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
- volumeWeightedPotentialEnergy(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
- volumeWeightedPotentialEnergyTopography(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h_s(1:nCellsSolve)*areaCell(1:nCellsSolve)
- refAreaWeightedSurfaceHeight(iLevel,:) = areaCell(1:nCellsSolve)*(h(iLevel,1:nCellsSolve)+h_s(1:nCellsSolve))
-
- do iEdge = 1,nEdgesSolve
- q = 0.0
- do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- workpv = 0.5 * (pv_edge(iLevel,iEdge) + pv_edge(iLevel,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * h_edge(iLevel,eoe)
- end do
- keTend_CoriolisForce(iLevel,iEdge) = h_edge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
-
- iCell1 = cellsOnEdge(1,iEdge)
- iCell2 = cellsOnEdge(2,iEdge)
-
- refAreaWeightedSurfaceHeight_edge(iLevel,iEdge) = areaEdge(iEdge)*(h_edge(iLevel,iEdge) + 0.5*(h_s(iCell1) + h_s(iCell2)))
-
- keTend_PressureGradient(iLevel,iEdge) = areaEdge(iEdge)*h_edge(iLevel,iEdge)*u(iLevel,iEdge) &
- *gravity*(h(iLevel,iCell2)+h_s(iCell2) - h(iLevel,iCell1)-h_s(iCell1))/dcEdge(iEdge)
- peTend_DivThickness(iLevel,iCell1) = peTend_DivThickness(iLevel,iCell1) &
- + h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
- peTend_DivThickness(iLevel,iCell2) = peTend_DivThickness(iLevel,iCell2) &
- - h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
- end do
-
- peTend_DivThickness(iLevel,:) = peTend_DivThickness(iLevel,1:nCells)*gravity &
- *(h(iLevel,1:nCells)+h_s(1:nCells))
- end do
-
- do iEdge = 1,nEdgesSolve
- iCell1 = cellsOnEdge(1,iEdge)
- iCell2 = cellsOnEdge(2,iEdge)
-
- h_s_edge(iEdge) = 0.5*(h_s(iCell1) + h_s(iCell2))
- end do
-
- ! Step 4
- ! 4. Call Function to compute Global Stat that you want.
- ! Computing Kinetic and Potential Energy Tendency Terms
- call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, keTend_PressureGradient, globalKineticEnergyTendency)
- call computeGlobalSum(dminfo, nVertLevels, nCells, peTend_DivThickness, globalPotentialEnergyTendency)
-
- ! Computing top and bottom of global mass integral
- call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, cellVolume, sumCellVolume)
- call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, cellArea, sumCellArea)
-
- globalKineticEnergyTendency = globalKineticEnergyTendency / sumCellVolume
- globalPotentialEnergyTendency = globalPotentialEnergyTendency / sumCellVolume
-
- ! Step 5
- ! 5. Finish computing the global stat/integral
- globalFluidThickness = sumCellVolume/sumCellArea
-
- ! Compute Average Sea Surface Height for Potential Energy and Enstrophy
- ! Reservoir computations
- call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, refAreaWeightedSurfaceHeight, sumrefAreaWeightedSurfaceHeight)
-
- averageThickness(:) = (sumrefAreaWeightedSurfaceHeight/sumCellArea)-h_s(1:nCellsSolve)
-
- ! Compute Volume Weighted Averages of Potential Vorticity and Potential Enstrophy
- call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
- call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialEnstrophy, globalPotentialEnstrophy)
- call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, vertexVolume, sumVertexVolume)
-
- globalPotentialVorticity = globalPotentialVorticity/sumVertexVolume
- globalPotentialEnstrophy = globalPotentialEnstrophy/sumVertexVolume
-
- ! Compte Potential Enstrophy Reservior
- potentialEnstrophyReservior(:) = areaCell(:)*fCell(:)*fCell(:)/averageThickness
- call computeGlobalSum(dminfo, 1, nCellsSolve, potentialEnstrophyReservior, globalPotentialEnstrophyReservoir)
- globalPotentialEnstrophyReservoir = globalPotentialEnstrophyReservoir/sumCellVolume
-
- globalPotentialEnstrophy = globalPotentialEnstrophy - globalPotentialEnstrophyReservoir
-
- ! Compute Kinetic and Potential Energy terms to be combined into total energy
- call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, volumeWeightedKineticEnergy, globalKineticEnergy)
- call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergy, globalPotentialEnergy)
- call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyTopography, global_temp)
-
- globalKineticEnergy = globalKineticEnergy/sumCellVolume
- globalPotentialEnergy = (globalPotentialEnergy + global_temp)/sumCellVolume
-
- ! Compute Potential energy reservoir to be subtracted from potential energy term
- volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*averageThickness*gravity*0.5
- call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, globalPotentialEnergyReservoir)
- volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*h_s(1:nCellsSolve)*gravity
- call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, global_temp)
-
- globalPotentialEnergyReservoir = (globalPotentialEnergyReservoir + global_temp)/sumCellVolume
-
- globalPotentialEnergy = globalPotentialEnergy - globalPotentialEnergyReservoir
- globalEnergy = globalKineticEnergy + globalPotentialEnergy
-
- ! Compute Coriolis energy tendency term
- call computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, keTend_CoriolisForce, globalCoriolisEnergyTendency)
- globalCoriolisEnergyTendency = globalCoriolisEnergyTendency/sumCellVolume
-
- ! Step 6
- ! 6. Write out your global stat to the file
- if (dminfo % my_proc_id == IO_NODE) then
- fileID = getFreeUnit()
-
- if (timeIndex/config_stats_interval == 1) then
- open(fileID, file='GlobalIntegrals.txt',STATUS='unknown')
- else
- open(fileID, file='GlobalIntegrals.txt',POSITION='append')
- endif
- write(fileID,'(1i0, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
- globalEnergy, globalCoriolisEnergyTendency, globalKineticEnergyTendency+globalPotentialEnergyTendency, &
- globalKineticEnergy, globalPotentialEnergy
- close(fileID)
- end if
-
- deallocate(areaEdge)
- end subroutine computeGlobalDiagnostics
-
- integer function getFreeUnit()
- implicit none
-
- integer :: index
- logical :: isOpened
-
- getFreeUnit = 0
- do index = 1,99
- if((index /= 5) .and. (index /= 6)) then
- inquire(unit = index, opened = isOpened)
- if( .not. isOpened) then
- getFreeUnit = index
- return
- end if
- end if
- end do
- end function getFreeUnit
-
- subroutine computeGlobalSum(dminfo, nVertLevels, nElements, field, globalSum)
-
- implicit none
-
- type (dm_info), intent(in) :: dminfo
- integer, intent(in) :: nVertLevels, nElements
- real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
- real (kind=RKIND), intent(out) :: globalSum
-
- real (kind=RKIND) :: localSum
-
- localSum = sum(field)
- call dmpar_sum_real(dminfo, localSum, globalSum)
-
- end subroutine computeGlobalSum
-
- subroutine computeGlobalMin(dminfo, nVertLevels, nElements, field, globalMin)
-
- implicit none
-
- type (dm_info), intent(in) :: dminfo
- integer, intent(in) :: nVertLevels, nElements
- real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
- real (kind=RKIND), intent(out) :: globalMin
-
- real (kind=RKIND) :: localMin
-
- localMin = minval(field)
- call dmpar_min_real(dminfo, localMin, globalMin)
-
- end subroutine computeGlobalMin
-
- subroutine computeGlobalMax(dminfo, nVertLevels, nElements, field, globalMax)
-
- implicit none
-
- type (dm_info), intent(in) :: dminfo
- integer, intent(in) :: nVertLevels, nElements
- real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
- real (kind=RKIND), intent(out) :: globalMax
-
- real (kind=RKIND) :: localMax
-
- localMax = maxval(field)
- call dmpar_max_real(dminfo, localMax, globalMax)
-
- end subroutine computeGlobalMax
-
- subroutine computeGlobalVertSumHorizMin(dminfo, nVertLevels, nElements, field, globalMin)
-
- implicit none
-
- type (dm_info), intent(in) :: dminfo
- integer, intent(in) :: nVertLevels, nElements
- real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
- real (kind=RKIND), intent(out) :: globalMin
-
- real (kind=RKIND) :: localMin
-
- localMin = minval(sum(field,1))
- call dmpar_min_real(dminfo, localMin, globalMin)
-
- end subroutine computeGlobalVertSumHorizMin
-
- subroutine computeGlobalVertSumHorizMax(dminfo, nVertLevels, nElements, field, globalMax)
-
- implicit none
-
- type (dm_info), intent(in) :: dminfo
- integer, intent(in) :: nVertLevels, nElements
- real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
- real (kind=RKIND), intent(out) :: globalMax
-
- real (kind=RKIND) :: localMax
-
- localMax = maxval(sum(field,1))
- call dmpar_max_real(dminfo, localMax, globalMax)
-
- end subroutine computeGlobalVertSumHorizMax
-
-end module global_diagnostics
Deleted: branches/source_renaming/src/core_sw/module_mpas_core.F
===================================================================
--- branches/source_renaming/src/core_sw/module_mpas_core.F        2011-10-17 15:40:02 UTC (rev 1086)
+++ branches/source_renaming/src/core_sw/module_mpas_core.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -1,382 +0,0 @@
-module mpas_core
-
- use mpas_framework
- use mpas_timekeeping
-
- type (io_output_object) :: restart_obj
- integer :: restart_frame
- integer :: current_outfile_frames
-
- type (MPAS_Clock_type) :: clock
-
- integer, parameter :: outputAlarmID = 1
- integer, parameter :: restartAlarmID = 2
- !integer, parameter :: statsAlarmID = 3
-
- contains
-
- subroutine mpas_core_init(domain, startTimeStamp)
-
- use configure
- use grid_types
- use test_cases
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- character(len=*), intent(out) :: startTimeStamp
-
- real (kind=RKIND) :: dt
- type (block_type), pointer :: block
-
-
- if (.not. config_do_restart) call setup_sw_test_case(domain)
-
- !
- ! Initialize core
- !
- dt = config_dt
-
- call simulation_clock_init(domain, dt, startTimeStamp)
-
- block => domain % blocklist
- do while (associated(block))
- call mpas_init_block(block, block % mesh, dt)
- block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
- block => block % next
- end do
-
- restart_frame = 1
- current_outfile_frames = 0
-
- end subroutine mpas_core_init
-
-
- subroutine simulation_clock_init(domain, dt, startTimeStamp)
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- real (kind=RKIND), intent(in) :: dt
- character(len=*), intent(out) :: startTimeStamp
-
- type (MPAS_Time_Type) :: startTime, stopTime, alarmStartTime
- type (MPAS_TimeInterval_type) :: runDuration, timeStep, alarmTimeStep
- integer :: ierr
-
- call MPAS_setTime(curr_time=startTime, dateTimeString=config_start_time, ierr=ierr)
- call MPAS_setTimeInterval(timeStep, dt=dt, ierr=ierr)
-
- if (trim(config_run_duration) /= "none") then
- call MPAS_setTimeInterval(runDuration, timeString=config_run_duration, ierr=ierr)
- call MPAS_createClock(clock, startTime=startTime, timeStep=timeStep, runDuration=runDuration, ierr=ierr)
-
- if (trim(config_stop_time) /= "none") then
- call MPAS_setTime(curr_time=stopTime, dateTimeString=config_stop_time, ierr=ierr)
- if(startTime + runduration /= stopTime) then
- write(0,*) 'Warning: config_run_duration and config_stop_time are inconsitent: using config_run_duration.'
- end if
- end if
- else if (trim(config_stop_time) /= "none") then
- call MPAS_setTime(curr_time=stopTime, dateTimeString=config_stop_time, ierr=ierr)
- call MPAS_createClock(clock, startTime=startTime, timeStep=timeStep, stopTime=stopTime, ierr=ierr)
- else
- write(0,*) 'Error: Neither config_run_duration nor config_stop_time were specified.'
- call dmpar_abort(domain % dminfo)
- end if
-
- ! set output alarm
- call MPAS_setTimeInterval(alarmTimeStep, timeString=config_output_interval, ierr=ierr)
- alarmStartTime = startTime + alarmTimeStep
- call MPAS_addClockAlarm(clock, outputAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
-
- ! set restart alarm, if necessary
- if (trim(config_restart_interval) /= "none") then
- call MPAS_setTimeInterval(alarmTimeStep, timeString=config_restart_interval, ierr=ierr)
- alarmStartTime = startTime + alarmTimeStep
- call MPAS_addClockAlarm(clock, restartAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
- end if
-
- !TODO: use this code if we desire to convert config_stats_interval to alarms
- !(must also change config_stats_interval type to character)
- ! set stats alarm, if necessary
- !if (trim(config_stats_interval) /= "none") then
- ! call MPAS_setTimeInterval(alarmTimeStep, timeString=config_stats_interval, ierr=ierr)
- ! alarmStartTime = startTime + alarmTimeStep
- ! call MPAS_addClockAlarm(clock, statsAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
- !end if
-
- call MPAS_getTime(curr_time=startTime, dateTimeString=startTimeStamp, ierr=ierr)
-
- end subroutine simulation_clock_init
-
-
- subroutine mpas_init_block(block, mesh, dt)
-
- use grid_types
- use time_integration
- use RBF_interpolation
- use vector_reconstruction
-
- implicit none
-
- type (block_type), intent(inout) :: block
- type (mesh_type), intent(inout) :: mesh
- real (kind=RKIND), intent(in) :: dt
-
-
- call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
- call compute_mesh_scaling(mesh)
-
- call rbfInterp_initialize(mesh)
- call init_reconstruct(mesh)
- call reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
- block % state % time_levs(1) % state % uReconstructX % array, &
- block % state % time_levs(1) % state % uReconstructY % array, &
- block % state % time_levs(1) % state % uReconstructZ % array, &
- block % state % time_levs(1) % state % uReconstructZonal % array, &
- block % state % time_levs(1) % state % uReconstructMeridional % array &
- )
-
-
- end subroutine mpas_init_block
-
-
- subroutine mpas_core_run(domain, output_obj, output_frame)
-
- use grid_types
- use io_output
- use timer
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- type (io_output_object), intent(inout) :: output_obj
- integer, intent(inout) :: output_frame
-
- integer :: itimestep
- real (kind=RKIND) :: dt
- type (block_type), pointer :: block_ptr
-
- type (MPAS_Time_Type) :: currTime
- character(len=32) :: timeStamp
- integer :: ierr
-
- ! Eventually, dt should be domain specific
- dt = config_dt
-
- currTime = MPAS_getClockTime(clock, MPAS_NOW, ierr)
- call MPAS_getTime(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
- write(0,*) 'Initial timestep ', timeStamp
-
- call write_output_frame(output_obj, output_frame, domain)
-
- ! During integration, time level 1 stores the model state at the beginning of the
- ! time step, and time level 2 stores the state advanced dt in time by timestep(...)
- itimestep = 0
- do while (.not. MPAS_isClockStopTime(clock))
-
- itimestep = itimestep + 1
- call MPAS_advanceClock(clock)
-
- currTime = MPAS_getClockTime(clock, MPAS_NOW, ierr)
- call MPAS_getTime(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
- write(0,*) 'Doing timestep ', timeStamp
-
- call timer_start("time integration")
- call mpas_timestep(domain, itimestep, dt, timeStamp)
- call timer_stop("time integration")
-
- ! Move time level 2 fields back into time level 1 for next time step
- call shift_time_levels_state(domain % blocklist % state)
-
- !TODO: MPAS_getClockRingingAlarms is probably faster than multiple MPAS_isAlarmRinging...
-
- if (MPAS_isAlarmRinging(clock, outputAlarmID, ierr=ierr)) then
- call MPAS_resetClockAlarm(clock, outputAlarmID, ierr=ierr)
- if(output_frame == 1) call output_state_init(output_obj, domain, "OUTPUT", trim(timeStamp)) ! output_frame will always be > 1 here unless it is reset after the output file is finalized
- call write_output_frame(output_obj, output_frame, domain)
- end if
-
- if (MPAS_isAlarmRinging(clock, restartAlarmID, ierr=ierr)) then
- call MPAS_resetClockAlarm(clock, restartAlarmID, ierr=ierr)
- if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
- call output_state_for_domain(restart_obj, domain, restart_frame)
- restart_frame = restart_frame + 1
- end if
-
- end do
-
- end subroutine mpas_core_run
-
-
- subroutine write_output_frame(output_obj, output_frame, domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields for a domain and write model state to output file
- !
- ! Input/Output: domain - contains model state; diagnostic field are computed
- ! before returning
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use grid_types
- use io_output
-
- implicit none
-
- type (io_output_object), intent(inout) :: output_obj
- integer, intent(inout) :: output_frame
- type (domain_type), intent(inout) :: domain
-
- integer :: i, j, k
- integer :: eoe
- type (block_type), pointer :: block_ptr
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call compute_output_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % mesh)
- block_ptr => block_ptr % next
- end do
-
- call output_state_for_domain(output_obj, domain, output_frame)
- output_frame = output_frame + 1
-
- ! if the maximum number of frames per outfile has been reached, finalize outfile and reset frame
- if (config_frames_per_outfile > 0) then
- current_outfile_frames = current_outfile_frames + 1
- if(current_outfile_frames >= config_frames_per_outfile) then
- current_outfile_frames = 0
- call output_state_finalize(output_obj, domain % dminfo)
- output_frame = 1
- end if
- end if
-
- end subroutine write_output_frame
-
-
- subroutine compute_output_diagnostics(state, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields for a domain
- !
- ! Input: state - contains model prognostic fields
- ! grid - contains grid metadata
- !
- ! Output: state - upon returning, diagnostic fields will have be computed
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use grid_types
-
- implicit none
-
- type (state_type), intent(inout) :: state
- type (mesh_type), intent(in) :: grid
-
- integer :: i, eoe
- integer :: iEdge, k
-
- end subroutine compute_output_diagnostics
-
-
- subroutine mpas_timestep(domain, itimestep, dt, timeStamp)
-
- use grid_types
- use time_integration
- use timer
- use global_diagnostics
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- integer, intent(in) :: itimestep
- real (kind=RKIND), intent(in) :: dt
- character(len=*), intent(in) :: timeStamp
-
- type (block_type), pointer :: block_ptr
- integer :: ierr
-
- call timestep(domain, dt, timeStamp)
-
- if(config_stats_interval .gt. 0) then
- if(mod(itimestep, config_stats_interval) == 0) then
- block_ptr => domain % blocklist
- if(associated(block_ptr % next)) then
- write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
- 'that there is only one block per processor.'
- end if
-
- call timer_start("global_diagnostics")
- call computeGlobalDiagnostics(domain % dminfo, &
- block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
- itimestep, dt)
- call timer_stop("global_diagnostics")
- end if
- end if
-
- !TODO: replace the above code block with this if we desire to convert config_stats_interval to use alarms
- !if (MPAS_isAlarmRinging(clock, statsAlarmID, ierr=ierr)) then
- ! call MPAS_resetClockAlarm(clock, statsAlarmID, ierr=ierr)
-
- ! block_ptr => domain % blocklist
- ! if(associated(block_ptr % next)) then
- ! write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
- ! 'that there is only one block per processor.'
- ! end if
-
- ! call timer_start("global_diagnostics")
- ! call computeGlobalDiagnostics(domain % dminfo, &
- ! block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
- ! timeStamp, dt)
- ! call timer_stop("global_diagnostics")
- !end if
-
- end subroutine mpas_timestep
-
-
- subroutine mpas_core_finalize(domain)
-
- use grid_types
-
- implicit none
-
- integer :: ierr
-
- type (domain_type), intent(inout) :: domain
-
- if (restart_frame > 1) call output_state_finalize(restart_obj, domain % dminfo)
-
- call MPAS_destroyClock(clock, ierr)
-
- end subroutine mpas_core_finalize
-
-
- subroutine compute_mesh_scaling(mesh)
-
- use grid_types
-
- implicit none
-
- type (mesh_type), intent(inout) :: mesh
-
- integer :: iEdge, cell1, cell2
- real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4
-
- meshDensity => mesh % meshDensity % array
- meshScalingDel2 => mesh % meshScalingDel2 % array
- meshScalingDel4 => mesh % meshScalingDel4 % array
-
- !
- ! Compute the scaling factors to be used in the del2 and del4 dissipation
- !
- meshScalingDel2(:) = 1.0
- meshScalingDel4(:) = 1.0
- if (config_h_ScaleWithMesh) then
- do iEdge=1,mesh%nEdges
- cell1 = mesh % cellsOnEdge % array(1,iEdge)
- cell2 = mesh % cellsOnEdge % array(2,iEdge)
- meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/12.0)
- meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/6.0)
- end do
- end if
-
- end subroutine compute_mesh_scaling
-
-end module mpas_core
Deleted: branches/source_renaming/src/core_sw/module_test_cases.F
===================================================================
--- branches/source_renaming/src/core_sw/module_test_cases.F        2011-10-17 15:40:02 UTC (rev 1086)
+++ branches/source_renaming/src/core_sw/module_test_cases.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -1,527 +0,0 @@
-module test_cases
-
- use grid_types
- use configure
- use constants
-
-
- contains
-
-
- subroutine setup_sw_test_case(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Configure grid metadata and model state for the shallow water test case
- ! specified in the namelist
- !
- ! Output: block - a subset (not necessarily proper) of the model domain to be
- ! initialized
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: i
- type (block_type), pointer :: block_ptr
-
- if (config_test_case == 0) then
- write(0,*) 'Using initial conditions supplied in input file'
-
- else if (config_test_case == 1) then
- write(0,*) 'Setting up shallow water test case 1'
- write(0,*) ' -- Advection of Cosine Bell over the Pole'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- do i=2,nTimeLevs
- call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 2) then
- write(0,*) 'Setting up shallow water test case 2'
- write(0,*) ' -- Setup shallow water test case 2: Global Steady State Nonlinear Zonal Geostrophic Flow'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_2(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- do i=2,nTimeLevs
- call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 5) then
- write(0,*) 'Setting up shallow water test case 5'
- write(0,*) ' -- Setup shallow water test case 5: Zonal Flow over an Isolated Mountain'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_5(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- do i=2,nTimeLevs
- call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 6) then
- write(0,*) 'Setting up shallow water test case 6'
- write(0,*) ' -- Rossby-Haurwitz Wave'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_6(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- do i=2,nTimeLevs
- call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- else
- write(0,*) 'Only test case 1, 2, 5, and 6 are currently supported.'
- stop
- end if
-
- end subroutine setup_sw_test_case
-
-
- subroutine sw_test_case_1(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 1: Advection of Cosine Bell over the Pole
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
- real (kind=RKIND), parameter :: h0 = 1000.0
- real (kind=RKIND), parameter :: theta_c = 0.0
- real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: alpha = pii/4.0
-
- 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
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- 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
- 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 cosine bell at (theta_c, lambda_c)
- !
- do iCell=1,grid % nCells
- r = sphere_distance(theta_c, lambda_c, grid % latCell % array(iCell), grid % lonCell % array(iCell), a)
- if (r < a/3.0) then
- state % h % array(1,iCell) = (h0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
- else
- state % h % array(1,iCell) = 0.0
- end if
- end do
-
- end subroutine sw_test_case_1
-
-
- subroutine sw_test_case_2(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 2: Global Steady State Nonlinear Zonal
- ! Geostrophic Flow
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
- real (kind=RKIND), parameter :: gh0 = 29400.0
- real (kind=RKIND), parameter :: alpha = 0.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- 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
- 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)
-
- !
- ! Generate rotated Coriolis field
- !
- do iEdge=1,grid % nEdges
- grid % fEdge % array(iEdge) = 2.0 * omega * &
- ( -cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(alpha) &
- )
- end do
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 2.0 * omega * &
- (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) &
- )
- end do
-
- !
- ! Initialize height field (actually, fluid thickness field)
- !
- do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
- (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
- sin(grid%latCell%array(iCell)) * cos(alpha) &
- )**2.0 &
- ) / &
- gravity
- end do
-
- end subroutine sw_test_case_2
-
-
- subroutine sw_test_case_5(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 5: Zonal Flow over an Isolated Mountain
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 20.
- real (kind=RKIND), parameter :: gh0 = 5960.0*gravity
- real (kind=RKIND), parameter :: hs0 = 2000.
- real (kind=RKIND), parameter :: theta_c = pii/6.0
- real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: rr = pii/9.0
- real (kind=RKIND), parameter :: alpha = 0.0
-
- 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
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- 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
- 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)
-
- !
- ! Generate rotated Coriolis field
- !
- do iEdge=1,grid % nEdges
- grid % fEdge % array(iEdge) = 2.0 * omega * &
- (-cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(alpha) &
- )
- end do
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 2.0 * omega * &
- (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) &
- )
- end do
-
- !
- ! Initialize mountain
- !
- do iCell=1,grid % nCells
- if (grid % lonCell % array(iCell) < 0.0) grid % lonCell % array(iCell) = grid % lonCell % array(iCell) + 2.0 * pii
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- grid % h_s % array(iCell) = hs0 * (1.0 - r/rr)
- end do
-
- !
- ! Initialize tracer fields
- !
- do iCell=1,grid % nCells
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- state % tracers % array(1,1,iCell) = 1.0 - r/rr
- end do
- if (grid%nTracers > 1) then
- do iCell=1,grid % nCells
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + &
- (grid % latCell % array(iCell) - theta_c - pii/6.0)**2.0 &
- ) &
- )
- state % tracers % array(2,1,iCell) = 1.0 - r/rr
- end do
- end if
-
- !
- ! Initialize height field (actually, fluid thickness field)
- !
- do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
- (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
- sin(grid%latCell%array(iCell)) * cos(alpha) &
- )**2.0 &
- ) / &
- gravity
- state % h % array(1,iCell) = state % h % array(1,iCell) - grid % h_s % array(iCell)
- end do
-
- end subroutine sw_test_case_5
-
-
- subroutine sw_test_case_6(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 6: Rossby-Haurwitz Wave
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: h0 = 8000.0
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- 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 * 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
- 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*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
-
-
- real function sphere_distance(lat1, lon1, lat2, lon2, radius)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
- ! sphere with given radius.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
- real (kind=RKIND) :: arg1
-
- arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
- cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
- sphere_distance = 2.*radius*asin(arg1)
-
- end function sphere_distance
-
-
- real function AA(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! A, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- 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))
-
- end function AA
-
-
- real function BB(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! B, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- 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)*cos(theta))**2.0)
-
- end function BB
-
-
- real function CC(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! C, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- real (kind=RKIND), intent(in) :: theta
-
- CC = 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 - R - 2.0)
-
- end function CC
-
-end module test_cases
Deleted: branches/source_renaming/src/core_sw/module_time_integration.F
===================================================================
--- branches/source_renaming/src/core_sw/module_time_integration.F        2011-10-17 15:40:02 UTC (rev 1086)
+++ branches/source_renaming/src/core_sw/module_time_integration.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -1,1287 +0,0 @@
-module time_integration
-
- use vector_reconstruction
- use grid_types
- use configure
- use constants
- use dmpar
-
-
- contains
-
-
- subroutine timestep(domain, dt, timeStamp)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Advance model state forward in time by the specified time step
- !
- ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
- ! plus grid meta-data
- ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
- ! model state advanced forward in time by dt seconds
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- real (kind=RKIND), intent(in) :: dt
- character(len=*), intent(in) :: timeStamp
-
- type (block_type), pointer :: block
-
- if (trim(config_time_integration) == 'RK4') then
- call rk4(domain, dt)
- else
- write(0,*) 'Unknown time integration option '//trim(config_time_integration)
- write(0,*) 'Currently, only ''RK4'' is supported.'
- stop
- end if
-
- block => domain % blocklist
- do while (associated(block))
- block % state % time_levs(2) % state % xtime % scalar = timeStamp
- block => block % next
- end do
-
- end subroutine timestep
-
-
- subroutine rk4(domain, dt)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Advance model state forward in time by the specified time step using
- ! 4th order Runge-Kutta
- !
- ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
- ! plus grid meta-data
- ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
- ! model state advanced forward in time by dt seconds
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- real (kind=RKIND), intent(in) :: dt
-
- integer :: iCell, k
- type (block_type), pointer :: block
- type (state_type) :: provis
-
- integer :: rk_step
-
- real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
-
- block => domain % blocklist
- call allocate_state(provis, &
- block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
- block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels, &
- block % mesh % nTracers)
-
- !
- ! Initialize time_levs(2) with state at current time
- ! Initialize first RK state
- ! Couple tracers time_levs(2) with h in time-levels
- ! Initialize RK weights
- !
- block => domain % blocklist
- do while (associated(block))
-
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
- do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % h % array(k,iCell)
- end do
- end do
-
- call copy_state(provis, block % state % time_levs(1) % state)
-
- block => block % next
- end do
-
- rk_weights(1) = dt/6.
- rk_weights(2) = dt/3.
- rk_weights(3) = dt/3.
- rk_weights(4) = dt/6.
-
- rk_substep_weights(1) = dt/2.
- rk_substep_weights(2) = dt/2.
- rk_substep_weights(3) = dt
- rk_substep_weights(4) = 0.
-
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! BEGIN RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do rk_step = 1, 4
-
-! --- update halos for diagnostic variables
-
- block => domain % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_edge % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
- if (config_h_mom_eddy_visc4 > 0.0) then
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nVertices, &
- block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- end if
-
- block => block % next
- end do
-
-! --- compute tendencies
-
- block => domain % blocklist
- do while (associated(block))
- call compute_tend(block % tend, provis, block % mesh)
- call compute_scalar_tend(block % tend, provis, block % mesh)
- call enforce_boundaryEdge(block % tend, block % mesh)
- block => block % next
- end do
-
-! --- update halos for prognostic variables
-
- block => domain % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field3dReal(domain % dminfo, block % tend % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
-
-! --- compute next substep state
-
- if (rk_step < 4) then
- block => domain % blocklist
- do while (associated(block))
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
- provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- provis % tracers % array(:,k,iCell) = ( &
- block % state % time_levs(1) % state % h % array(k,iCell) * &
- block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
- ) / provis % h % array(k,iCell)
- end do
- end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- end if
- call compute_solve_diagnostics(dt, provis, block % mesh)
- block => block % next
- end do
- end if
-
-!--- accumulate update (for RK4)
-
- block => domain % blocklist
- do while (associated(block))
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
- + rk_weights(rk_step) * block % tend % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
- + rk_weights(rk_step) * block % tend % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
- end do
- end do
- block => block % next
- end do
-
- end do
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! END RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
- !
- ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
- !
- block => domain % blocklist
- do while (associated(block))
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- / block % state % time_levs(2) % state % h % array(k,iCell)
- end do
- end do
-
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- end if
-
- call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
-
- call reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
- block % state % time_levs(2) % state % uReconstructX % array, &
- block % state % time_levs(2) % state % uReconstructY % array, &
- block % state % time_levs(2) % state % uReconstructZ % array, &
- block % state % time_levs(2) % state % uReconstructZonal % array, &
- block % state % time_levs(2) % state % uReconstructMeridional % array &
- )
-
- block => block % next
- end do
-
- call deallocate_state(provis)
-
- end subroutine rk4
-
-
- subroutine compute_tend(tend, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute height and normal wind tendencies, as well as diagnostic variables
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (mesh_type), intent(in) :: grid
-
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
-
- integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- meshScalingDel2, meshScalingDel4
- real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, divergence, h_vertex
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r, u_diffusion
-
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
- real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
- real (kind=RKIND) :: ke_edge
-
-
- h => s % h % array
- u => s % u % array
- v => s % v % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- pv_edge => s % pv_edge % array
- vh => s % vh % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- fVertex => grid % fVertex % array
- fEdge => grid % fEdge % array
-
- tend_h => tend % h % array
- tend_u => tend % u % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- u_src => grid % u_src % array
-
- meshScalingDel2 => grid % meshScalingDel2 % array
- meshScalingDel4 => grid % meshScalingDel4 % array
-
-
- !
- ! Compute height tendency for each cell
- !
- tend_h(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = 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
- end do
- end do
- do iCell=1,grid % nCellsSolve
- do k=1,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
- end do
- end do
-
-#ifdef LANL_FORMULATION
- !
- ! Compute u (normal) velocity tendency for each edge (cell face)
- !
- tend_u(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = 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)
- end do
-
- tend_u(k,iEdge) = &
- q &
- - ( ke(k,cell2) - ke(k,cell1) + &
- gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ) / dcEdge(iEdge)
- end do
- end do
-
-
-#endif
-
-#ifdef NCAR_FORMULATION
- !
- ! Compute u (normal) velocity tendency for each edge (cell face)
- !
- tend_u(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,nVertLevels
- vorticity_abs = fEdge(iEdge) + (circulation(k,vertex1) + circulation(k,vertex2)) / &
- (areaTriangle(vertex1) + areaTriangle(vertex2))
-
- workpv = 2.0 * vorticity_abs / (h(k,cell1) + h(k,cell2))
-
- tend_u(k,iEdge) = workpv * vh(k,iEdge) - &
- (ke(k,cell2) - ke(k,cell1) + &
- gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ) / &
- dcEdge(iEdge)
- end do
- end do
-#endif
-
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
- ! only valid for visc == constant
- if (config_h_mom_eddy_visc2 > 0.0) then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
-
- do k=1,nVertLevels
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
- end do
- end do
- end if
-
- !
- ! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="red">abla^4 u
- ! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- ! applied recursively.
- ! strictly only valid for h_mom_eddy_visc4 == constant
- !
- if (config_h_mom_eddy_visc4 > 0.0) then
- allocate(delsq_divergence(nVertLevels, nCells+1))
- allocate(delsq_u(nVertLevels, nEdges+1))
- allocate(delsq_circulation(nVertLevels, nVertices+1))
- allocate(delsq_vorticity(nVertLevels, nVertices+1))
-
- delsq_u(:,:) = 0.0
-
- ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
-
- do k=1,nVertLevels
-
- delsq_u(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
-
- end do
- end do
-
- ! vorticity using </font>
<font color="red">abla^2 u
- delsq_circulation(:,:) = 0.0
- do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
- delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
- - dcEdge(iEdge) * delsq_u(k,iEdge)
- delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
- + dcEdge(iEdge) * delsq_u(k,iEdge)
- end do
- end do
- do iVertex=1,nVertices
- r = 1.0 / areaTriangle(iVertex)
- do k=1,nVertLevels
- delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
- end do
- end do
-
- ! Divergence using </font>
<font color="red">abla^2 u
- delsq_divergence(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
- delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
- + delsq_u(k,iEdge)*dvEdge(iEdge)
- delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
- - delsq_u(k,iEdge)*dvEdge(iEdge)
- end do
- end do
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
- delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
- end do
- end do
-
- ! Compute - \kappa </font>
<font color="red">abla^4 u
- ! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="red">abla^2 u) )
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
-
- do k=1,nVertLevels
-
- u_diffusion = ( delsq_divergence(k,cell2) &
- - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
- -( delsq_vorticity(k,vertex2) &
- - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
-
- u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
- tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
-
- end do
- end do
-
- deallocate(delsq_divergence)
- deallocate(delsq_u)
- deallocate(delsq_circulation)
- deallocate(delsq_vorticity)
-
- end if
-
- ! Compute u (velocity) tendency from wind stress (u_src)
- if(config_wind_stress) then
- do iEdge=1,grid % nEdges
- tend_u(1,iEdge) = tend_u(1,iEdge) &
- + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
- end do
- endif
-
- if (config_bottom_drag) then
- do iEdge=1,grid % nEdges
- ! bottom drag is the same as POP:
- ! -c |u| u where c is unitless and 1.0e-3.
- ! see POP Reference guide, section 3.4.4.
- ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
- + ke(1,cellsOnEdge(2,iEdge)))
-
- tend_u(1,iEdge) = tend_u(1,iEdge) &
- - 1.0e-3*u(1,iEdge) &
- *sqrt(2.0*ke_edge)/h_edge(1,iEdge)
- end do
- endif
-
- end subroutine compute_tend
-
-
- subroutine compute_scalar_tend(tend, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed scalar tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (mesh_type), intent(in) :: grid
-
- integer :: iCell, iEdge, k, iTracer, cell1, cell2, i
- real (kind=RKIND) :: flux, tracer_edge, r
- real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
- integer, dimension(:,:), pointer :: boundaryEdge
- real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
- real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
-
- real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
- integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
- real (kind=RKIND), dimension(:,:), pointer :: u, h_edge
-
- u => s % u % array
- h_edge => s % h_edge % array
- dcEdge => grid % dcEdge % array
- deriv_two => grid % deriv_two % array
- dvEdge => grid % dvEdge % array
- tracers => s % tracers % array
- cellsOnEdge => grid % cellsOnEdge % array
- boundaryCell=> grid % boundaryCell % array
- boundaryEdge=> grid % boundaryEdge % array
- areaCell => grid % areaCell % array
- tracer_tend => tend % tracers % array
-
- coef_3rd_order = 0.
- if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
- if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
-
-
- tracer_tend(:,:,:) = 0.0
-
- if (config_tracer_adv_order == 2) then
-
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- 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) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- end do
- end do
- end if
- end do
-
- else if (config_tracer_adv_order == 3) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,grid % nTracers
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- !-- else u <= 0:
- else
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- end if
-
- !-- update tendency
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
- end do
-
- else if (config_tracer_adv_order == 4) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if an edge is not on the outer-most ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,grid % nTracers
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- !-- update tendency
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
- end do
-
- endif ! if (config_tracer_adv_order == 2 )
-
- !
- ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="red">abla \phi)
- !
- if ( config_h_tracer_eddy_diff2 > 0.0 ) then
-
- !
- ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
- !
- allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
- boundaryMask = 1.0
- where(boundaryEdge.eq.1) boundaryMask=0.0
-
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- 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="red">abla \phi on edge
- tracer_turb_flux = config_h_tracer_eddy_diff2 &
- *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
-
- ! div(h \kappa_2 </font>
<font color="red">abla \phi) at cell center
- 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
- end do
-
- end do
-
- deallocate(boundaryMask)
-
- end if
-
- !
- ! tracer tendency: del4 horizontal tracer diffusion, &
- ! div(h \kappa_4 </font>
<font color="black">abla [div(h </font>
<font color="red">abla \phi)])
- !
- if ( config_h_tracer_eddy_diff4 > 0.0 ) then
-
- !
- ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
- !
- allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
- boundaryMask = 1.0
- where(boundaryEdge.eq.1) boundaryMask=0.0
-
- allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
-
- delsq_tracer(:,:,:) = 0.
-
- ! first del2: div(h </font>
<font color="red">abla \phi) at cell center
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,grid % nVertLevels
- do iTracer=1, grid % nTracers
- delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
- + 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) &
- - 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
- r = 1.0 / grid % areaCell % array(iCell)
- do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
- end do
- end do
- end do
-
- ! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- invAreaCell1 = 1.0 / grid % areaCell % array(cell1)
- invAreaCell2 = 1.0 / grid % areaCell % array(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 = 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
- enddo
-
- end do
-
- deallocate(delsq_tracer)
- deallocate(boundaryMask)
-
- end if
-
- end subroutine compute_scalar_tend
-
-
- subroutine compute_solve_diagnostics(dt, s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields used in the tendency computations
- !
- ! Input: grid - grid metadata
- !
- ! Output: s - computed diagnostics
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: dt
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
-
-
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, workpv
-
- integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- h_vertex, vorticity_cell
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r, h1, h2
- real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
-
- h => s % h % array
- u => s % u % array
- v => s % v % array
- vh => s % vh % array
- h_edge => s % h_edge % array
- h_vertex => s % h_vertex % array
- tend_h => s % h % array
- tend_u => s % u % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- pv_edge => s % pv_edge % array
- pv_vertex => s % pv_vertex % array
- pv_cell => s % pv_cell % array
- vorticity_cell => s % vorticity_cell % array
- gradPVn => s % gradPVn % array
- gradPVt => s % gradPVt % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- fVertex => grid % fVertex % array
- fEdge => grid % fEdge % array
- deriv_two => grid % deriv_two % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- boundaryEdge => grid % boundaryEdge % array
- boundaryCell => grid % boundaryCell % array
-
- !
- ! Find those cells that have an edge on the boundary
- !
- boundaryCell(:,:) = 0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if(boundaryEdge(k,iEdge).eq.1) then
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- boundaryCell(k,cell1) = 1
- boundaryCell(k,cell2) = 1
- endif
- enddo
- enddo
-
- !
- ! Compute height on cell edges at velocity locations
- ! Namelist options control the order of accuracy of the reconstructed h_edge value
- !
-
- coef_3rd_order = 0.
- if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
- if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
-
- if (config_thickness_adv_order == 2) then
-
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- 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))
- end do
- end if
- end do
-
- else if (config_thickness_adv_order == 3) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- !-- else u <= 0:
- else
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end if
-
- end do ! do k
- end if ! if (cell1 <=
- end do ! do iEdge
-
- else if (config_thickness_adv_order == 4) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- end do ! do k
- end if ! if (cell1 <=
- end do ! do iEdge
-
- endif ! if(config_thickness_adv_order == 2)
-
- !
- ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
- ! used to when reading for edges that do not exist
- !
- u(:,nEdges+1) = 0.0
-
- !
- ! Compute circulation and relative vorticity at each vertex
- !
- circulation(:,:) = 0.0
- do iEdge=1,nEdges
- 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)
- end do
- end do
- do iVertex=1,nVertices
- do k=1,nVertLevels
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
- end do
- end do
-
-
- !
- ! Compute the divergence at each cell center
- !
- divergence(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- 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)*dvEdge(iEdge)
- enddo
- end if
- end do
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
- divergence(k,iCell) = divergence(k,iCell) * r
- enddo
- enddo
-
- !
- ! Compute kinetic energy in each cell
- !
- ke(:,:) = 0.0
- do iCell=1,nCells
- do i=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i,iCell)
- do k=1,nVertLevels
- 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
- ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- end do
- end do
-
- !
- ! Compute v (tangential) velocities
- !
- v(:,:) = 0.0
- do iEdge = 1,nEdges
- do i=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(i,iEdge)
- do k = 1,nVertLevels
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
- end do
- end do
- end do
-
-#ifdef NCAR_FORMULATION
- !
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
- vh(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- do j=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- do k=1,nVertLevels
- vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * u(k,eoe) * h_edge(k,eoe)
- end do
- end do
- end do
-#endif
-
-
- !
- ! Compute height at vertices, pv at vertices, and average pv to edge locations
- ! ( 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
- end do
-
-
- !
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
- !
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
- dvEdge(iEdge)
- enddo
- enddo
-
- !
- ! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells )
- !
- pv_edge(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
- do k=1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
- end do
- end do
- end do
-
-
- !
- ! Modify PV edge with upstream bias.
- !
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
- enddo
- enddo
-
-
- !
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
- !
- pv_cell(:,:) = 0.0
- vorticity_cell(:,:) = 0.0
- do iVertex = 1, nVertices
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- if (iCell <= nCells) then
- do k = 1,nVertLevels
- pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
- vorticity_cell(k,iCell) = vorticity_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) / areaCell(iCell)
- enddo
- endif
- enddo
- enddo
-
-
- !
- ! Compute gradient of PV in normal direction
- ! ( this computes gradPVn for all edges bounding real cells )
- !
- gradPVn(:,:) = 0.0
- do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
- do k = 1,nVertLevels
- gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
- dcEdge(iEdge)
- enddo
- endif
- enddo
-
- ! Modify PV edge with upstream bias.
- !
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
- enddo
- enddo
-
- !
- ! set pv_edge = fEdge / h_edge at boundary points
- !
- ! if (maxval(boundaryEdge).ge.0) then
- ! do iEdge = 1,nEdges
- ! cell1 = cellsOnEdge(1,iEdge)
- ! cell2 = cellsOnEdge(2,iEdge)
- ! do k = 1,nVertLevels
- ! if(boundaryEdge(k,iEdge).eq.1) then
- ! v(k,iEdge) = 0.0
- ! if(cell1.gt.0) then
- ! h1 = h(k,cell1)
- ! pv_edge(k,iEdge) = fEdge(iEdge) / h1
- ! h_edge(k,iEdge) = h1
- ! else
- ! h2 = h(k,cell2)
- ! pv_edge(k,iEdge) = fEdge(iEdge) / h2
- ! h_edge(k,iEdge) = h2
- ! endif
- ! endif
- ! enddo
- ! enddo
- ! endif
-
-
- end subroutine compute_solve_diagnostics
-
-
- subroutine enforce_boundaryEdge(tend, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Enforce any boundary conditions on the normal velocity at each edge
- !
- ! Input: grid - grid metadata
- !
- ! Output: tend_u set to zero at boundaryEdge == 1 locations
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (mesh_type), intent(in) :: grid
-
- integer, dimension(:,:), pointer :: boundaryEdge
- real (kind=RKIND), dimension(:,:), pointer :: tend_u
- integer :: nCells, nEdges, nVertices, nVertLevels
- integer :: iEdge, k
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- boundaryEdge => grid % boundaryEdge % array
- tend_u => tend % u % array
-
- if(maxval(boundaryEdge).le.0) return
-
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
-
- if(boundaryEdge(k,iEdge).eq.1) then
- tend_u(k,iEdge) = 0.0
- endif
-
- enddo
- enddo
-
- end subroutine enforce_boundaryEdge
-
-
-end module time_integration
Copied: branches/source_renaming/src/core_sw/mpas_sw_advection.F (from rev 1084, branches/source_renaming/src/core_sw/module_advection.F)
===================================================================
--- branches/source_renaming/src/core_sw/mpas_sw_advection.F         (rev 0)
+++ branches/source_renaming/src/core_sw/mpas_sw_advection.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -0,0 +1,933 @@
+module advection
+
+ use grid_types
+ use configure
+ use constants
+
+
+ contains
+
+
+ subroutine sw_initialize_advection_rk( grid )
+
+!
+! compute the cell coefficients for the polynomial fit.
+! this is performed during setup for model integration.
+! WCS, 31 August 2009
+!
+ implicit none
+
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ integer, dimension(:,:), pointer :: advCells
+
+! local variables
+
+ real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
+ 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
+ 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) :: amatrix(25,25), bmatrix(25,25), wmatrix(25,25)
+ real (kind=RKIND) :: length_scale
+ integer :: ma,na, cell_add, mw, nn
+ integer, dimension(25) :: cell_list
+
+
+ integer :: cell1, cell2
+ integer, parameter :: polynomial_order = 2
+! logical, parameter :: debug = .true.
+ logical, parameter :: debug = .false.
+! logical, parameter :: least_squares = .false.
+ logical, parameter :: least_squares = .true.
+ logical :: add_the_cell, do_the_cell
+
+ logical, parameter :: reset_poly = .true.
+
+ real (kind=RKIND) :: rcell, cos2t, costsint, sin2t
+ real (kind=RKIND), dimension(grid%maxEdges) :: angle_2d
+
+!---
+
+ pii = 2.*asin(1.0)
+
+ advCells => grid % advCells % array
+ deriv_two => grid % deriv_two % array
+ deriv_two(:,:,:) = 0.
+
+ do iCell = 1, grid % nCells ! is this correct? - we need first halo cell also...
+
+ 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
+
+ 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
+ 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
+ end if
+ end do
+ end do
+ end if
+
+ advCells(1,iCell) = n
+
+! check to see if we are reaching outside the halo
+
+ do_the_cell = .true.
+ do i=1,n
+ if (cell_list(i) > grid % nCells) do_the_cell = .false.
+ end do
+
+
+ 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
+ advCells(i+1,iCell) = cell_list(i)
+ xc(i) = grid % xCell % array(advCells(i+1,iCell))/a
+ yc(i) = grid % yCell % array(advCells(i+1,iCell))/a
+ zc(i) = grid % zCell % array(advCells(i+1,iCell))/a
+ end do
+
+ theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
+ xc(2), yc(2), zc(2), &
+ 0., 0., 1. )
+
+! angles from cell center to neighbor centers (thetav)
+
+ do i=1,n-1
+
+ ip2 = i+2
+ if (ip2 > n) 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) )
+
+ dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1) )
+ end do
+
+ length_scale = 1.
+ do i=1,n-1
+ dl_sphere(i) = dl_sphere(i)/length_scale
+ 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)
+ end do
+
+ else ! On an x-y plane
+
+ do i=1,n-1
+
+ angle_2d(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
+
+! xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
+! yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+
+ 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
+ mw = grid % nEdgesOnCell % array (iCell)
+
+ bmatrix = 0.
+ amatrix = 0.
+ wmatrix = 0.
+
+ if (polynomial_order == 2) then
+ na = 6
+ ma = ma+1
+
+ amatrix(1,1) = 1.
+ 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
+
+ wmatrix(i,i) = 1.
+ end do
+
+ else if (polynomial_order == 3) then
+ na = 10
+ ma = ma+1
+
+ amatrix(1,1) = 1.
+ 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,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
+
+ wmatrix(i,i) = 1.
+
+ end do
+
+ else
+ na = 15
+ ma = ma+1
+
+ amatrix(1,1) = 1.
+ 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,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,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
+
+ wmatrix(i,i) = 1.
+
+ end do
+
+ do i=1,mw
+ wmatrix(i,i) = 1.
+ end do
+
+ end if
+
+ call sw_poly_fit_2( amatrix, bmatrix, wmatrix, ma, na, 25 )
+
+ do i=1,grid % nEdgesOnCell % array (iCell)
+ ip1 = i+1
+ if (ip1 > n-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
+
+ if ( grid % on_a_sphere ) then
+ call sw_arc_bisect( xv1, yv1, zv1, &
+ xv2, yv2, zv2, &
+ xec, yec, zec )
+
+ 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)
+ if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
+ thetae(1,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ else
+ thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_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
+
+ do i=1, grid % nEdgesOnCell % array (iCell)
+ iEdge = grid % EdgesOnCell % array (i,iCell)
+
+
+ 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)))
+ costsint = cos2t*sin2t
+ cos2t = cos2t**2
+ sin2t = sin2t**2
+
+ do j=1,n
+ 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)))
+ costsint = cos2t*sin2t
+ cos2t = cos2t**2
+ sin2t = sin2t**2
+
+ do j=1,n
+ deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ + 2.*costsint*bmatrix(5,j) &
+ + 2.*sin2t*bmatrix(6,j)
+ end do
+ end if
+
+ else
+
+ cos2t = cos(angle_2d(i))
+ sin2t = sin(angle_2d(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
+ 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
+ deriv_two(j,2,iEdge) = 2.*cos2t*bmatrix(4,j) &
+ + 2.*costsint*bmatrix(5,j) &
+ + 2.*sin2t*bmatrix(6,j)
+ end do
+ end if
+
+ end if
+ end do
+
+ end do ! end of loop over cells
+
+ if (debug) stop
+
+
+! 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 sw_initialize_advection_rk
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! FUNCTION SPHERE_ANGLE
+ !
+ ! Computes the angle between arcs AB and AC, given points A, B, and C
+ ! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz
+
+ real (kind=RKIND) :: a, b, c ! Side lengths of spherical triangle ABC
+
+ real (kind=RKIND) :: ABx, ABy, ABz ! The components of the vector AB
+ real (kind=RKIND) :: mAB ! The magnitude of AB
+ real (kind=RKIND) :: ACx, ACy, ACz ! The components of the vector AC
+ real (kind=RKIND) :: mAC ! The magnitude of AC
+
+ real (kind=RKIND) :: Dx ! The i-components of the cross product AB x AC
+ real (kind=RKIND) :: Dy ! The j-components of the cross product AB x AC
+ real (kind=RKIND) :: Dz ! The k-components of the cross product AB x AC
+
+ real (kind=RKIND) :: s ! Semiperimeter of the triangle
+ real (kind=RKIND) :: sin_angle
+
+ a = acos(max(min(bx*cx + by*cy + bz*cz,1.0),-1.0)) ! Eqn. (3)
+ b = acos(max(min(ax*cx + ay*cy + az*cz,1.0),-1.0)) ! Eqn. (2)
+ c = acos(max(min(ax*bx + ay*by + az*bz,1.0),-1.0)) ! Eqn. (1)
+
+ ABx = bx - ax
+ ABy = by - ay
+ ABz = bz - az
+
+ ACx = cx - ax
+ ACy = cy - ay
+ ACz = cz - az
+
+ Dx = (ABy * ACz) - (ABz * ACy)
+ Dy = -((ABx * ACz) - (ABz * ACx))
+ Dz = (ABx * ACy) - (ABy * ACx)
+
+ s = 0.5*(a + b + c)
+! sin_angle = sqrt((sin(s-b)*sin(s-c))/(sin(b)*sin(c))) ! Eqn. (28)
+ sin_angle = sqrt(min(1.,max(0.,(sin(s-b)*sin(s-c))/(sin(b)*sin(c))))) ! Eqn. (28)
+
+ if ((Dx*ax + Dy*ay + Dz*az) >= 0.0) then
+ sphere_angle = 2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ else
+ sphere_angle = -2.0 * asin(max(min(sin_angle,1.0),-1.0))
+ end if
+
+ end function sphere_angle
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! FUNCTION PLANE_ANGLE
+ !
+ ! Computes the angle between vectors AB and AC, given points A, B, and C, and
+ ! a vector (u,v,w) normal to the plane.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w
+
+ real (kind=RKIND) :: ABx, ABy, ABz ! The components of the vector AB
+ real (kind=RKIND) :: mAB ! The magnitude of AB
+ real (kind=RKIND) :: ACx, ACy, ACz ! The components of the vector AC
+ real (kind=RKIND) :: mAC ! The magnitude of AC
+
+ real (kind=RKIND) :: Dx ! The i-components of the cross product AB x AC
+ real (kind=RKIND) :: Dy ! The j-components of the cross product AB x AC
+ real (kind=RKIND) :: Dz ! The k-components of the cross product AB x AC
+
+ real (kind=RKIND) :: cos_angle
+
+ ABx = bx - ax
+ ABy = by - ay
+ ABz = bz - az
+ mAB = sqrt(ABx**2.0 + ABy**2.0 + ABz**2.0)
+
+ ACx = cx - ax
+ ACy = cy - ay
+ ACz = cz - az
+ mAC = sqrt(ACx**2.0 + ACy**2.0 + ACz**2.0)
+
+
+ Dx = (ABy * ACz) - (ABz * ACy)
+ Dy = -((ABx * ACz) - (ABz * ACx))
+ Dz = (ABx * ACy) - (ABy * ACx)
+
+ cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
+
+ if ((Dx*u + Dy*v + Dz*w) >= 0.0) then
+ plane_angle = acos(max(min(cos_angle,1.0),-1.0))
+ else
+ plane_angle = -acos(max(min(cos_angle,1.0),-1.0))
+ end if
+
+ end function plane_angle
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! FUNCTION ARC_LENGTH
+ !
+ ! Returns the length of the great circle arc from A=(ax, ay, az) to
+ ! B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
+ ! same sphere centered at the origin.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ real function arc_length(ax, ay, az, bx, by, bz)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz
+
+ real (kind=RKIND) :: r, c
+ real (kind=RKIND) :: cx, cy, cz
+
+ cx = bx - ax
+ cy = by - ay
+ cz = bz - az
+
+! r = ax*ax + ay*ay + az*az
+! c = cx*cx + cy*cy + cz*cz
+!
+! arc_length = sqrt(r) * acos(1.0 - c/(2.0*r))
+
+ r = sqrt(ax*ax + ay*ay + az*az)
+ c = sqrt(cx*cx + cy*cy + cz*cz)
+! arc_length = sqrt(r) * 2.0 * asin(c/(2.0*r))
+ arc_length = r * 2.0 * asin(c/(2.0*r))
+
+ end function arc_length
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTine sw_arc_bisect
+ !
+ ! Returns the point C=(cx, cy, cz) that bisects the great circle arc from
+ ! A=(ax, ay, az) to B=(bx, by, bz). It is assumed that A and B lie on the
+ ! surface of a sphere centered at the origin.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine sw_arc_bisect(ax, ay, az, bx, by, bz, cx, cy, cz)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz
+ real (kind=RKIND), intent(out) :: cx, cy, cz
+
+ real (kind=RKIND) :: r ! Radius of the sphere
+ real (kind=RKIND) :: d
+
+ r = sqrt(ax*ax + ay*ay + az*az)
+
+ cx = 0.5*(ax + bx)
+ cy = 0.5*(ay + by)
+ cz = 0.5*(az + bz)
+
+ if (cx == 0. .and. cy == 0. .and. cz == 0.) then
+ write(0,*) 'Error: arc_bisect: A and B are diametrically opposite'
+ else
+ d = sqrt(cx*cx + cy*cy + cz*cz)
+ cx = r * cx / d
+ cy = r * cy / d
+ cz = r * cz / d
+ end if
+
+ end subroutine sw_arc_bisect
+
+
+ subroutine sw_poly_fit_2(a_in,b_out,weights_in,m,n,ne)
+
+ implicit none
+
+ integer, intent(in) :: m,n,ne
+ real (kind=RKIND), dimension(ne,ne), intent(in) :: a_in, weights_in
+ real (kind=RKIND), dimension(ne,ne), intent(out) :: b_out
+
+ ! local storage
+
+ real (kind=RKIND), dimension(m,n) :: a
+ real (kind=RKIND), dimension(n,m) :: b
+ real (kind=RKIND), dimension(m,m) :: w,wt,h
+ real (kind=RKIND), dimension(n,m) :: at, ath
+ real (kind=RKIND), dimension(n,n) :: ata, ata_inv, atha, atha_inv
+ integer, dimension(n) :: indx
+ integer :: i,j
+
+ if ( (ne<n) .or. (ne<m) ) then
+ write(6,*) ' error in poly_fit_2 inversion ',m,n,ne
+ stop
+ end if
+
+! a(1:m,1:n) = a_in(1:n,1:m)
+ a(1:m,1:n) = a_in(1:m,1:n)
+ w(1:m,1:m) = weights_in(1:m,1:m)
+ b_out(:,:) = 0.
+
+ wt = transpose(w)
+ h = matmul(wt,w)
+ at = transpose(a)
+ ath = matmul(at,h)
+ atha = matmul(ath,a)
+
+ ata = matmul(at,a)
+
+! if (m == n) then
+! call sw_migs(a,n,b,indx)
+! else
+
+ call sw_migs(atha,n,atha_inv,indx)
+
+ b = matmul(atha_inv,ath)
+
+! call sw_migs(ata,n,ata_inv,indx)
+! b = matmul(ata_inv,at)
+! end if
+ b_out(1:n,1:m) = b(1:n,1:m)
+
+! do i=1,n
+! write(6,*) ' i, indx ',i,indx(i)
+! end do
+!
+! write(6,*) ' '
+
+ end subroutine sw_poly_fit_2
+
+
+! Updated 10/24/2001.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! !
+! Please Note: !
+! !
+! (1) This computer program is written by Tao Pang in conjunction with !
+! his book, "An Introduction to Computational Physics," published !
+! by Cambridge University Press in 1997. !
+! !
+! (2) No warranties, express or implied, are made for this program. !
+! !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+SUBROUTine sw_migs (A,N,X,INDX)
+!
+! Subroutine to invert matrix A(N,N) with the inverse stored
+! in X(N,N) in the output. Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ INTEGER, INTENT (IN) :: N
+ INTEGER :: I,J,K
+ INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
+ REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
+ REAL (kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
+ REAL (kind=RKIND), DIMENSION (N,N) :: B
+!
+ DO I = 1, N
+ DO J = 1, N
+ B(I,J) = 0.0
+ END DO
+ END DO
+ DO I = 1, N
+ B(I,I) = 1.0
+ END DO
+!
+ call sw_elgs (A,N,INDX)
+!
+ DO I = 1, N-1
+ DO J = I+1, N
+ DO K = 1, N
+ B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K)
+ END DO
+ END DO
+ END DO
+!
+ DO I = 1, N
+ X(N,I) = B(INDX(N),I)/A(INDX(N),N)
+ DO J = N-1, 1, -1
+ X(J,I) = B(INDX(J),I)
+ DO K = J+1, N
+ X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
+ END DO
+ X(J,I) = X(J,I)/A(INDX(J),J)
+ END DO
+ END DO
+END SUBROUTine sw_migs
+
+
+SUBROUTine sw_elgs (A,N,INDX)
+!
+! Subroutine to perform the partial-pivoting Gaussian elimination.
+! A(N,N) is the original matrix in the input and transformed matrix
+! plus the pivoting element ratios below the diagonal in the output.
+! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ INTEGER, INTENT (IN) :: N
+ INTEGER :: I,J,K,ITMP
+ INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
+ REAL (kind=RKIND) :: C1,PI,PI1,PJ
+ REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
+ REAL (kind=RKIND), DIMENSION (N) :: C
+!
+! Initialize the index
+!
+ DO I = 1, N
+ INDX(I) = I
+ END DO
+!
+! Find the rescaling factors, one from each row
+!
+ DO I = 1, N
+ C1= 0.0
+ DO J = 1, N
+ C1 = MAX(C1,ABS(A(I,J)))
+ END DO
+ C(I) = C1
+ END DO
+!
+! Search the pivoting (largest) element from each column
+!
+ DO J = 1, N-1
+ PI1 = 0.0
+ DO I = J, N
+ PI = ABS(A(INDX(I),J))/C(INDX(I))
+ IF (PI.GT.PI1) THEN
+ PI1 = PI
+ K = I
+ ENDIF
+ END DO
+!
+! Interchange the rows via INDX(N) to record pivoting order
+!
+ ITMP = INDX(J)
+ INDX(J) = INDX(K)
+ INDX(K) = ITMP
+ DO I = J+1, N
+ PJ = A(INDX(I),J)/A(INDX(J),J)
+!
+! Record pivoting ratios below the diagonal
+!
+ A(INDX(I),J) = PJ
+!
+! Modify other elements accordingly
+!
+ DO K = J+1, N
+ A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
+ END DO
+ END DO
+ END DO
+!
+END SUBROUTine sw_elgs
+
+!-------------------------------------------------------------
+
+ subroutine sw_initialize_deformation_weights( grid )
+
+!
+! compute the cell coefficients for the deformation calculations
+! WCS, 13 July 2010
+!
+ implicit none
+
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
+
+! local variables
+
+ real (kind=RKIND), dimension(2, grid % nEdges) :: thetae
+ 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
+ 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, xpt, ypt
+
+ real (kind=RKIND) :: length_scale
+ integer :: ma,na, cell_add, mw, nn
+ integer, dimension(25) :: cell_list
+
+ integer :: cell1, cell2, iv
+ logical :: do_the_cell
+ real (kind=RKIND) :: area_cell, sint2, cost2, sint_cost, sumw1, sumw2, xptt, area_cellt
+
+ logical, parameter :: debug = .false.
+
+ if (debug) write(0,*) ' in def weight calc '
+
+ defc_a => grid % defc_a % array
+ defc_b => grid % defc_b % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+
+ defc_a(:,:) = 0.
+ defc_b(:,:) = 0.
+
+ pii = 2.*asin(1.0)
+
+ if (debug) write(0,*) ' beginning cell loop '
+
+ do iCell = 1, grid % nCells
+
+ if (debug) write(0,*) ' cell loop ', iCell
+
+ 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
+
+! check to see if we are reaching outside the halo
+
+ if (debug) write(0,*) ' points ', n
+
+ do_the_cell = .true.
+ do i=1,n
+ if (cell_list(i) > grid % nCells) do_the_cell = .false.
+ end do
+
+
+ if (.not. do_the_cell) cycle
+
+
+! compute poynomial fit for this cell if all needed neighbors exist
+ if (grid % on_a_sphere) then
+
+ xc(1) = grid % xCell % array(iCell)/a
+ yc(1) = grid % yCell % array(iCell)/a
+ zc(1) = grid % zCell % array(iCell)/a
+
+
+ do i=2,n
+ iv = grid % verticesOnCell % array(i-1,iCell)
+ xc(i) = grid % xVertex % array(iv)/a
+ yc(i) = grid % yVertex % array(iv)/a
+ zc(i) = grid % zVertex % array(iv)/a
+ end do
+
+ theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
+ xc(2), yc(2), zc(2), &
+ 0., 0., 1. )
+
+! angles from cell center to neighbor centers (thetav)
+
+ do i=1,n-1
+
+ ip2 = i+2
+ if (ip2 > n) 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) )
+
+ dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1) )
+ end do
+
+ length_scale = 1.
+ do i=1,n-1
+ dl_sphere(i) = dl_sphere(i)/length_scale
+ 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)
+ end do
+
+ else ! On an x-y plane
+
+ xp(1) = grid % xCell % array(iCell)
+ yp(1) = grid % yCell % array(iCell)
+
+
+ do i=2,n
+ iv = grid % verticesOnCell % array(i-1,iCell)
+ xp(i) = grid % xVertex % array(iv)
+ yp(i) = grid % yVertex % array(iv)
+ end do
+
+ end if
+
+! thetat(1) = 0.
+ thetat(1) = theta_abs(iCell)
+ do i=2,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ thetat(i) = plane_angle( 0.,0.,0., &
+ xp(i)-xp(i-1), yp(i)-yp(i-1), 0., &
+ xp(ip1)-xp(i), yp(ip1)-yp(i), 0., &
+ 0., 0., 1.)
+ thetat(i) = thetat(i) + thetat(i-1)
+ end do
+
+ area_cell = 0.
+ area_cellt = 0.
+ do i=1,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
+ area_cell = area_cell + 0.25*(xp(i)+xp(ip1))*(yp(ip1)-yp(i)) - 0.25*(yp(i)+yp(ip1))*(xp(ip1)-xp(i))
+ area_cellt = area_cellt + (0.25*(xp(i)+xp(ip1))*cos(thetat(i)) + 0.25*(yp(i)+yp(ip1))*sin(thetat(i)))*dl
+ end do
+ if (debug) write(0,*) ' area_cell, area_cellt ',area_cell, area_cellt,area_cell-area_cellt
+
+ do i=1,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
+ sint2 = (sin(thetat(i)))**2
+ cost2 = (cos(thetat(i)))**2
+ sint_cost = sin(thetat(i))*cos(thetat(i))
+ defc_a(i,iCell) = dl*(cost2 - sint2)/area_cell
+ defc_b(i,iCell) = dl*2.*sint_cost/area_cell
+ if (cellsOnEdge(1,EdgesOnCell(i,iCell)) /= iCell) then
+ defc_a(i,iCell) = - defc_a(i,iCell)
+ defc_b(i,iCell) = - defc_b(i,iCell)
+ end if
+
+ end do
+
+ end do
+
+ if (debug) write(0,*) ' exiting def weight calc '
+
+ end subroutine sw_initialize_deformation_weights
+
+end module advection
Copied: branches/source_renaming/src/core_sw/mpas_sw_global_diagnostics.F (from rev 1084, branches/source_renaming/src/core_sw/module_global_diagnostics.F)
===================================================================
--- branches/source_renaming/src/core_sw/mpas_sw_global_diagnostics.F         (rev 0)
+++ branches/source_renaming/src/core_sw/mpas_sw_global_diagnostics.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -0,0 +1,384 @@
+module global_diagnostics
+
+ use grid_types
+ use configure
+ use constants
+ use dmpar
+
+ implicit none
+ save
+ public
+
+ contains
+
+ subroutine sw_compute_global_diagnostics(dminfo, state, grid, timeIndex, dt)
+
+ ! Note: this routine assumes that there is only one block per processor. No looping
+ ! is preformed over blocks.
+ ! dminfo is the domain info needed for global communication
+ ! state contains the state variables needed to compute global diagnostics
+ ! grid conains the meta data about the grid
+ ! timeIndex is the current time step counter
+ ! dt is the duration of each time step
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! INSTRUCTIONS !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! To add a new Diagnostic as a Global Stat, follow these steps.
+ ! 1. Define the array to integrate, and the variable for the value above.
+ ! 2. Allocate the array with the correct dimensions.
+ ! 3. Fill the array with the data to be integrated.
+ ! eg. GlobalFluidThickness = Sum(h dA)/Sum(dA), See below for array filling
+ ! 4. Call Function to compute Global Stat that you want.
+ ! 5. Finish computing the global stat/integral
+ ! 6. Write out your global stat to the file
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+ integer, intent(in) :: timeIndex
+ real (kind=RKIND), intent(in) :: dt
+
+ integer :: nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, nCellsGlobal, nEdgesGlobal, nVerticesGlobal, iTracer
+ integer :: nCells
+
+ ! Step 1
+ ! 1. Define the array to integrate, and the variable for the value to be stored in after the integration
+ real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, h_s, fCell, fEdge
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, pv_edge, pv_vertex, pv_cell, h_vertex, weightsOnEdge
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
+ real (kind=RKIND), dimension(:), allocatable :: volumeWeightedPotentialEnergyReservoir, averageThickness
+ real (kind=RKIND), dimension(:), allocatable :: potentialEnstrophyReservior, areaEdge, h_s_edge
+
+ real (kind=RKIND), dimension(:,:), allocatable :: cellVolume, cellArea, volumeWeightedPotentialVorticity
+ real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnstrophy, vertexVolume, volumeWeightedKineticEnergy
+ real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnergy, volumeWeightedPotentialEnergyTopography
+ real (kind=RKIND), dimension(:,:), allocatable :: keTend_CoriolisForce, keTend_PressureGradient
+ real (kind=RKIND), dimension(:,:), allocatable ::peTend_DivThickness, refAreaWeightedSurfaceHeight, refAreaWeightedSurfaceHeight_edge
+
+ real (kind=RKIND) :: sumCellVolume, sumCellArea, sumVertexVolume, sumrefAreaWeightedSurfaceHeight
+
+ real (kind=RKIND) :: globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, globalEnergy
+ real (kind=RKIND) :: globalCoriolisEnergyTendency, globalKEPETendency, globalPotentialEnstrophyReservoir
+ real (kind=RKIND) :: globalKineticEnergy, globalPotentialEnergy, globalPotentialEnergyReservoir
+ real (kind=RKIND) :: globalKineticEnergyTendency, globalPotentialEnergyTendency
+ real (kind=RKIND) :: global_temp, workpv, q
+ real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
+
+ integer :: elementIndex, variableIndex, nVariables, nSums, nMaxes, nMins
+ integer :: timeLevel, eoe, iLevel, iCell, iEdge, iVertex
+ integer :: fileID, iCell1, iCell2, j
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgesOnEdge
+ integer, dimension(:), pointer :: nEdgesOnEdge
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+
+ nVertLevels = grid % nVertLevels
+ nCellsSolve = grid % nCellsSolve
+ nEdgesSolve = grid % nEdgesSolve
+ nVerticesSolve = grid % nVerticesSolve
+ nCells = grid % nCells
+
+ h_s => grid % h_s % array
+ areaCell => grid % areaCell % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaTriangle => grid % areaTriangle % array
+ fCell => grid % fCell % array
+ fEdge => grid % fEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+
+ allocate(areaEdge(1:nEdgesSolve))
+ areaEdge = dcEdge(1:nEdgesSolve)*dvEdge(1:nEdgesSolve)
+ weightsOnEdge => grid % weightsOnEdge % array
+
+ h => state % h % array
+ u => state % u % array
+ v => state % v % array
+ tracers => state % tracers % array
+ h_edge => state % h_edge % array
+ h_vertex => state % h_vertex % array
+ pv_edge => state % pv_edge % array
+ pv_vertex => state % pv_vertex % array
+ pv_cell => state % pv_cell % array
+
+ ! Step 2
+ ! 2. Allocate the array with the correct dimensions.
+ allocate(cellVolume(nVertLevels,nCellsSolve))
+ allocate(cellArea(nVertLevels,nCellsSolve))
+ allocate(refAreaWeightedSurfaceHeight(nVertLevels,nCellsSolve))
+ allocate(refAreaWeightedSurfaceHeight_edge(nVertLevels,nEdgesSolve))
+ allocate(volumeWeightedPotentialVorticity(nVertLevels,nVerticesSolve))
+ allocate(volumeWeightedPotentialEnstrophy(nVertLevels,nVerticesSolve))
+ allocate(potentialEnstrophyReservior(nCellsSolve))
+ allocate(vertexVolume(nVertLevels,nVerticesSolve))
+ allocate(volumeWeightedKineticEnergy(nVertLevels,nEdgesSolve))
+ allocate(volumeWeightedPotentialEnergy(nVertLevels,nCellsSolve))
+ allocate(volumeWeightedPotentialEnergyTopography(nVertLevels,nCellsSolve))
+ allocate(volumeWeightedPotentialEnergyReservoir(nCellsSolve))
+ allocate(keTend_CoriolisForce(nVertLevels,nEdgesSolve))
+ allocate(keTend_PressureGradient(nVertLevels,nEdgesSolve))
+ allocate(peTend_DivThickness(nVertLevels,nCells))
+
+ allocate(averageThickness(nCellsSolve))
+
+ allocate(h_s_edge(nEdgesSOlve))
+
+
+ cellVolume = 0
+ refAreaWeightedSurfaceHeight = 0
+ refAreaWeightedSurfaceHeight_edge = 0
+ vertexVolume = 0
+ cellArea = 0
+ averageThickness = 0
+ volumeWeightedPotentialVorticity = 0
+ volumeWeightedPotentialEnstrophy = 0
+ volumeWeightedKineticEnergy = 0
+ volumeWeightedPotentialEnergy = 0
+ volumeWeightedPotentialEnergyTopography = 0
+ volumeWeightedPotentialEnergyReservoir = 0
+ keTend_PressureGradient = 0
+ peTend_DivThickness = 0
+ keTend_CoriolisForce = 0
+ h_s_edge = 0
+
+ ! Build Arrays for Global Integrals
+ ! Step 3
+ ! 3. Fill the array with the data to be integrated.
+ ! eg. GlobalFluidThickness = Sum(h dA)/Sum(dA), See below for array filling
+ do iLevel = 1,nVertLevels
+ ! eg. GlobalFluidThickness top (Sum( h dA)) = Sum(cellVolume)
+ cellVolume(iLevel,:) = h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
+ ! eg. GlobalFluidThickness bot (Sum(dA)) = Sum(cellArea)
+ cellArea(iLevel,:) = areaCell(1:nCellsSolve)
+ volumeWeightedPotentialVorticity(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
+ *h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+ volumeWeightedPotentialEnstrophy(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
+ *pv_vertex(iLevel,1:nVerticesSolve)*h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+ vertexVolume(iLevel,:) = h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+ volumeWeightedKineticEnergy(iLevel,:) = u(iLevel,1:nEdgesSolve)*u(iLevel,1:nEdgesSolve) &
+ *h_edge(iLevel,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
+ volumeWeightedPotentialEnergy(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
+ volumeWeightedPotentialEnergyTopography(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h_s(1:nCellsSolve)*areaCell(1:nCellsSolve)
+ refAreaWeightedSurfaceHeight(iLevel,:) = areaCell(1:nCellsSolve)*(h(iLevel,1:nCellsSolve)+h_s(1:nCellsSolve))
+
+ do iEdge = 1,nEdgesSolve
+ q = 0.0
+ do j = 1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ workpv = 0.5 * (pv_edge(iLevel,iEdge) + pv_edge(iLevel,eoe))
+ q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * h_edge(iLevel,eoe)
+ end do
+ keTend_CoriolisForce(iLevel,iEdge) = h_edge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
+
+ iCell1 = cellsOnEdge(1,iEdge)
+ iCell2 = cellsOnEdge(2,iEdge)
+
+ refAreaWeightedSurfaceHeight_edge(iLevel,iEdge) = areaEdge(iEdge)*(h_edge(iLevel,iEdge) + 0.5*(h_s(iCell1) + h_s(iCell2)))
+
+ keTend_PressureGradient(iLevel,iEdge) = areaEdge(iEdge)*h_edge(iLevel,iEdge)*u(iLevel,iEdge) &
+ *gravity*(h(iLevel,iCell2)+h_s(iCell2) - h(iLevel,iCell1)-h_s(iCell1))/dcEdge(iEdge)
+ peTend_DivThickness(iLevel,iCell1) = peTend_DivThickness(iLevel,iCell1) &
+ + h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
+ peTend_DivThickness(iLevel,iCell2) = peTend_DivThickness(iLevel,iCell2) &
+ - h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
+ end do
+
+ peTend_DivThickness(iLevel,:) = peTend_DivThickness(iLevel,1:nCells)*gravity &
+ *(h(iLevel,1:nCells)+h_s(1:nCells))
+ end do
+
+ do iEdge = 1,nEdgesSolve
+ iCell1 = cellsOnEdge(1,iEdge)
+ iCell2 = cellsOnEdge(2,iEdge)
+
+ h_s_edge(iEdge) = 0.5*(h_s(iCell1) + h_s(iCell2))
+ end do
+
+ ! Step 4
+ ! 4. Call Function to compute Global Stat that you want.
+ ! Computing Kinetic and Potential Energy Tendency Terms
+ call sw_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, keTend_PressureGradient, globalKineticEnergyTendency)
+ call sw_compute_global_sum(dminfo, nVertLevels, nCells, peTend_DivThickness, globalPotentialEnergyTendency)
+
+ ! Computing top and bottom of global mass integral
+ call sw_compute_global_sum(dminfo, nVertLevels, nCellsSolve, cellVolume, sumCellVolume)
+ call sw_compute_global_sum(dminfo, nVertLevels, nCellsSolve, cellArea, sumCellArea)
+
+ globalKineticEnergyTendency = globalKineticEnergyTendency / sumCellVolume
+ globalPotentialEnergyTendency = globalPotentialEnergyTendency / sumCellVolume
+
+ ! Step 5
+ ! 5. Finish computing the global stat/integral
+ globalFluidThickness = sumCellVolume/sumCellArea
+
+ ! Compute Average Sea Surface Height for Potential Energy and Enstrophy
+ ! Reservoir computations
+ call sw_compute_global_sum(dminfo, nVertLevels, nCellsSolve, refAreaWeightedSurfaceHeight, sumrefAreaWeightedSurfaceHeight)
+
+ averageThickness(:) = (sumrefAreaWeightedSurfaceHeight/sumCellArea)-h_s(1:nCellsSolve)
+
+ ! Compute Volume Weighted Averages of Potential Vorticity and Potential Enstrophy
+ call sw_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
+ call sw_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialEnstrophy, globalPotentialEnstrophy)
+ call sw_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, vertexVolume, sumVertexVolume)
+
+ globalPotentialVorticity = globalPotentialVorticity/sumVertexVolume
+ globalPotentialEnstrophy = globalPotentialEnstrophy/sumVertexVolume
+
+ ! Compte Potential Enstrophy Reservior
+ potentialEnstrophyReservior(:) = areaCell(:)*fCell(:)*fCell(:)/averageThickness
+ call sw_compute_global_sum(dminfo, 1, nCellsSolve, potentialEnstrophyReservior, globalPotentialEnstrophyReservoir)
+ globalPotentialEnstrophyReservoir = globalPotentialEnstrophyReservoir/sumCellVolume
+
+ globalPotentialEnstrophy = globalPotentialEnstrophy - globalPotentialEnstrophyReservoir
+
+ ! Compute Kinetic and Potential Energy terms to be combined into total energy
+ call sw_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, volumeWeightedKineticEnergy, globalKineticEnergy)
+ call sw_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergy, globalPotentialEnergy)
+ call sw_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyTopography, global_temp)
+
+ globalKineticEnergy = globalKineticEnergy/sumCellVolume
+ globalPotentialEnergy = (globalPotentialEnergy + global_temp)/sumCellVolume
+
+ ! Compute Potential energy reservoir to be subtracted from potential energy term
+ volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*averageThickness*gravity*0.5
+ call sw_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, globalPotentialEnergyReservoir)
+ volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*h_s(1:nCellsSolve)*gravity
+ call sw_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, global_temp)
+
+ globalPotentialEnergyReservoir = (globalPotentialEnergyReservoir + global_temp)/sumCellVolume
+
+ globalPotentialEnergy = globalPotentialEnergy - globalPotentialEnergyReservoir
+ globalEnergy = globalKineticEnergy + globalPotentialEnergy
+
+ ! Compute Coriolis energy tendency term
+ call sw_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, keTend_CoriolisForce, globalCoriolisEnergyTendency)
+ globalCoriolisEnergyTendency = globalCoriolisEnergyTendency/sumCellVolume
+
+ ! Step 6
+ ! 6. Write out your global stat to the file
+ if (dminfo % my_proc_id == IO_NODE) then
+ fileID = getFreeUnit()
+
+ if (timeIndex/config_stats_interval == 1) then
+ open(fileID, file='GlobalIntegrals.txt',STATUS='unknown')
+ else
+ open(fileID, file='GlobalIntegrals.txt',POSITION='append')
+ endif
+ write(fileID,'(1i0, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
+ globalEnergy, globalCoriolisEnergyTendency, globalKineticEnergyTendency+globalPotentialEnergyTendency, &
+ globalKineticEnergy, globalPotentialEnergy
+ close(fileID)
+ end if
+
+ deallocate(areaEdge)
+ end subroutine sw_compute_global_diagnostics
+
+ integer function getFreeUnit()
+ implicit none
+
+ integer :: index
+ logical :: isOpened
+
+ getFreeUnit = 0
+ do index = 1,99
+ if((index /= 5) .and. (index /= 6)) then
+ inquire(unit = index, opened = isOpened)
+ if( .not. isOpened) then
+ getFreeUnit = index
+ return
+ end if
+ end if
+ end do
+ end function getFreeUnit
+
+ subroutine sw_compute_global_sum(dminfo, nVertLevels, nElements, field, globalSum)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalSum
+
+ real (kind=RKIND) :: localSum
+
+ localSum = sum(field)
+ call dmpar_sum_real(dminfo, localSum, globalSum)
+
+ end subroutine sw_compute_global_sum
+
+ subroutine sw_compute_global_min(dminfo, nVertLevels, nElements, field, globalMin)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMin
+
+ real (kind=RKIND) :: localMin
+
+ localMin = minval(field)
+ call dmpar_min_real(dminfo, localMin, globalMin)
+
+ end subroutine sw_compute_global_min
+
+ subroutine sw_compute_global_max(dminfo, nVertLevels, nElements, field, globalMax)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMax
+
+ real (kind=RKIND) :: localMax
+
+ localMax = maxval(field)
+ call dmpar_max_real(dminfo, localMax, globalMax)
+
+ end subroutine sw_compute_global_max
+
+ subroutine computeGlobalVertSumHorizMin(dminfo, nVertLevels, nElements, field, globalMin)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMin
+
+ real (kind=RKIND) :: localMin
+
+ localMin = minval(sum(field,1))
+ call dmpar_min_real(dminfo, localMin, globalMin)
+
+ end subroutine computeGlobalVertSumHorizMin
+
+ subroutine sw_compute_global_vert_sum_horiz_max(dminfo, nVertLevels, nElements, field, globalMax)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: nVertLevels, nElements
+ real (kind=RKIND), dimension(nVertLevels, nElements), intent(in) :: field
+ real (kind=RKIND), intent(out) :: globalMax
+
+ real (kind=RKIND) :: localMax
+
+ localMax = maxval(sum(field,1))
+ call dmpar_max_real(dminfo, localMax, globalMax)
+
+ end subroutine sw_compute_global_vert_sum_horiz_max
+
+end module global_diagnostics
Copied: branches/source_renaming/src/core_sw/mpas_sw_mpas_core.F (from rev 1084, branches/source_renaming/src/core_sw/module_mpas_core.F)
===================================================================
--- branches/source_renaming/src/core_sw/mpas_sw_mpas_core.F         (rev 0)
+++ branches/source_renaming/src/core_sw/mpas_sw_mpas_core.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -0,0 +1,382 @@
+module mpas_core
+
+ use mpas_framework
+ use mpas_timekeeping
+
+ type (io_output_object) :: restart_obj
+ integer :: restart_frame
+ integer :: current_outfile_frames
+
+ type (MPAS_Clock_type) :: clock
+
+ integer, parameter :: outputAlarmID = 1
+ integer, parameter :: restartAlarmID = 2
+ !integer, parameter :: statsAlarmID = 3
+
+ contains
+
+ subroutine mpas_core_init(domain, startTimeStamp)
+
+ use configure
+ use grid_types
+ use test_cases
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ character(len=*), intent(out) :: startTimeStamp
+
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block
+
+
+ if (.not. config_do_restart) call setup_sw_test_case(domain)
+
+ !
+ ! Initialize core
+ !
+ dt = config_dt
+
+ call simulation_clock_init(domain, dt, startTimeStamp)
+
+ block => domain % blocklist
+ do while (associated(block))
+ call mpas_init_block(block, block % mesh, dt)
+ block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
+ block => block % next
+ end do
+
+ restart_frame = 1
+ current_outfile_frames = 0
+
+ end subroutine mpas_core_init
+
+
+ subroutine simulation_clock_init(domain, dt, startTimeStamp)
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ real (kind=RKIND), intent(in) :: dt
+ character(len=*), intent(out) :: startTimeStamp
+
+ type (MPAS_Time_Type) :: startTime, stopTime, alarmStartTime
+ type (MPAS_TimeInterval_type) :: runDuration, timeStep, alarmTimeStep
+ integer :: ierr
+
+ call MPAS_setTime(curr_time=startTime, dateTimeString=config_start_time, ierr=ierr)
+ call MPAS_setTimeInterval(timeStep, dt=dt, ierr=ierr)
+
+ if (trim(config_run_duration) /= "none") then
+ call MPAS_setTimeInterval(runDuration, timeString=config_run_duration, ierr=ierr)
+ call MPAS_createClock(clock, startTime=startTime, timeStep=timeStep, runDuration=runDuration, ierr=ierr)
+
+ if (trim(config_stop_time) /= "none") then
+ call MPAS_setTime(curr_time=stopTime, dateTimeString=config_stop_time, ierr=ierr)
+ if(startTime + runduration /= stopTime) then
+ write(0,*) 'Warning: config_run_duration and config_stop_time are inconsitent: using config_run_duration.'
+ end if
+ end if
+ else if (trim(config_stop_time) /= "none") then
+ call MPAS_setTime(curr_time=stopTime, dateTimeString=config_stop_time, ierr=ierr)
+ call MPAS_createClock(clock, startTime=startTime, timeStep=timeStep, stopTime=stopTime, ierr=ierr)
+ else
+ write(0,*) 'Error: Neither config_run_duration nor config_stop_time were specified.'
+ call dmpar_abort(domain % dminfo)
+ end if
+
+ ! set output alarm
+ call MPAS_setTimeInterval(alarmTimeStep, timeString=config_output_interval, ierr=ierr)
+ alarmStartTime = startTime + alarmTimeStep
+ call MPAS_addClockAlarm(clock, outputAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
+
+ ! set restart alarm, if necessary
+ if (trim(config_restart_interval) /= "none") then
+ call MPAS_setTimeInterval(alarmTimeStep, timeString=config_restart_interval, ierr=ierr)
+ alarmStartTime = startTime + alarmTimeStep
+ call MPAS_addClockAlarm(clock, restartAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
+ end if
+
+ !TODO: use this code if we desire to convert config_stats_interval to alarms
+ !(must also change config_stats_interval type to character)
+ ! set stats alarm, if necessary
+ !if (trim(config_stats_interval) /= "none") then
+ ! call MPAS_setTimeInterval(alarmTimeStep, timeString=config_stats_interval, ierr=ierr)
+ ! alarmStartTime = startTime + alarmTimeStep
+ ! call MPAS_addClockAlarm(clock, statsAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
+ !end if
+
+ call MPAS_getTime(curr_time=startTime, dateTimeString=startTimeStamp, ierr=ierr)
+
+ end subroutine simulation_clock_init
+
+
+ subroutine mpas_init_block(block, mesh, dt)
+
+ use grid_types
+ use time_integration
+ use RBF_interpolation
+ use vector_reconstruction
+
+ implicit none
+
+ type (block_type), intent(inout) :: block
+ type (mesh_type), intent(inout) :: mesh
+ real (kind=RKIND), intent(in) :: dt
+
+
+ call sw_compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
+ call compute_mesh_scaling(mesh)
+
+ call rbfInterp_initialize(mesh)
+ call init_reconstruct(mesh)
+ call reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
+ block % state % time_levs(1) % state % uReconstructX % array, &
+ block % state % time_levs(1) % state % uReconstructY % array, &
+ block % state % time_levs(1) % state % uReconstructZ % array, &
+ block % state % time_levs(1) % state % uReconstructZonal % array, &
+ block % state % time_levs(1) % state % uReconstructMeridional % array &
+ )
+
+
+ end subroutine mpas_init_block
+
+
+ subroutine mpas_core_run(domain, output_obj, output_frame)
+
+ use grid_types
+ use io_output
+ use timer
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ type (io_output_object), intent(inout) :: output_obj
+ integer, intent(inout) :: output_frame
+
+ integer :: itimestep
+ real (kind=RKIND) :: dt
+ type (block_type), pointer :: block_ptr
+
+ type (MPAS_Time_Type) :: currTime
+ character(len=32) :: timeStamp
+ integer :: ierr
+
+ ! Eventually, dt should be domain specific
+ dt = config_dt
+
+ currTime = MPAS_getClockTime(clock, MPAS_NOW, ierr)
+ call MPAS_getTime(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
+ write(0,*) 'Initial timestep ', timeStamp
+
+ call write_output_frame(output_obj, output_frame, domain)
+
+ ! During integration, time level 1 stores the model state at the beginning of the
+ ! time step, and time level 2 stores the state advanced dt in time by timestep(...)
+ itimestep = 0
+ do while (.not. MPAS_isClockStopTime(clock))
+
+ itimestep = itimestep + 1
+ call MPAS_advanceClock(clock)
+
+ currTime = MPAS_getClockTime(clock, MPAS_NOW, ierr)
+ call MPAS_getTime(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
+ write(0,*) 'Doing timestep ', timeStamp
+
+ call timer_start("time integration")
+ call mpas_timestep(domain, itimestep, dt, timeStamp)
+ call timer_stop("time integration")
+
+ ! Move time level 2 fields back into time level 1 for next time step
+ call shift_time_levels_state(domain % blocklist % state)
+
+ !TODO: MPAS_getClockRingingAlarms is probably faster than multiple MPAS_isAlarmRinging...
+
+ if (MPAS_isAlarmRinging(clock, outputAlarmID, ierr=ierr)) then
+ call MPAS_resetClockAlarm(clock, outputAlarmID, ierr=ierr)
+ if(output_frame == 1) call output_state_init(output_obj, domain, "OUTPUT", trim(timeStamp)) ! output_frame will always be > 1 here unless it is reset after the output file is finalized
+ call write_output_frame(output_obj, output_frame, domain)
+ end if
+
+ if (MPAS_isAlarmRinging(clock, restartAlarmID, ierr=ierr)) then
+ call MPAS_resetClockAlarm(clock, restartAlarmID, ierr=ierr)
+ if (restart_frame == 1) call output_state_init(restart_obj, domain, "RESTART")
+ call output_state_for_domain(restart_obj, domain, restart_frame)
+ restart_frame = restart_frame + 1
+ end if
+
+ end do
+
+ end subroutine mpas_core_run
+
+
+ subroutine write_output_frame(output_obj, output_frame, domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields for a domain and write model state to output file
+ !
+ ! Input/Output: domain - contains model state; diagnostic field are computed
+ ! before returning
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use grid_types
+ use io_output
+
+ implicit none
+
+ type (io_output_object), intent(inout) :: output_obj
+ integer, intent(inout) :: output_frame
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, j, k
+ integer :: eoe
+ type (block_type), pointer :: block_ptr
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call compute_output_diagnostics(block_ptr % state % time_levs(1) % state, block_ptr % mesh)
+ block_ptr => block_ptr % next
+ end do
+
+ call output_state_for_domain(output_obj, domain, output_frame)
+ output_frame = output_frame + 1
+
+ ! if the maximum number of frames per outfile has been reached, finalize outfile and reset frame
+ if (config_frames_per_outfile > 0) then
+ current_outfile_frames = current_outfile_frames + 1
+ if(current_outfile_frames >= config_frames_per_outfile) then
+ current_outfile_frames = 0
+ call output_state_finalize(output_obj, domain % dminfo)
+ output_frame = 1
+ end if
+ end if
+
+ end subroutine write_output_frame
+
+
+ subroutine compute_output_diagnostics(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields for a domain
+ !
+ ! Input: state - contains model prognostic fields
+ ! grid - contains grid metadata
+ !
+ ! Output: state - upon returning, diagnostic fields will have be computed
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ use grid_types
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, eoe
+ integer :: iEdge, k
+
+ end subroutine compute_output_diagnostics
+
+
+ subroutine mpas_timestep(domain, itimestep, dt, timeStamp)
+
+ use grid_types
+ use time_integration
+ use timer
+ use global_diagnostics
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ integer, intent(in) :: itimestep
+ real (kind=RKIND), intent(in) :: dt
+ character(len=*), intent(in) :: timeStamp
+
+ type (block_type), pointer :: block_ptr
+ integer :: ierr
+
+ call sw_timestep(domain, dt, timeStamp)
+
+ if(config_stats_interval .gt. 0) then
+ if(mod(itimestep, config_stats_interval) == 0) then
+ block_ptr => domain % blocklist
+ if(associated(block_ptr % next)) then
+ write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
+ 'that there is only one block per processor.'
+ end if
+
+ call timer_start("global_diagnostics")
+ call sw_compute_global_diagnostics(domain % dminfo, &
+ block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
+ itimestep, dt)
+ call timer_stop("global_diagnostics")
+ end if
+ end if
+
+ !TODO: replace the above code block with this if we desire to convert config_stats_interval to use alarms
+ !if (MPAS_isAlarmRinging(clock, statsAlarmID, ierr=ierr)) then
+ ! call MPAS_resetClockAlarm(clock, statsAlarmID, ierr=ierr)
+
+ ! block_ptr => domain % blocklist
+ ! if(associated(block_ptr % next)) then
+ ! write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
+ ! 'that there is only one block per processor.'
+ ! end if
+
+ ! call timer_start("global_diagnostics")
+ ! call sw_compute_global_diagnostics(domain % dminfo, &
+ ! block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
+ ! timeStamp, dt)
+ ! call timer_stop("global_diagnostics")
+ !end if
+
+ end subroutine mpas_timestep
+
+
+ subroutine mpas_core_finalize(domain)
+
+ use grid_types
+
+ implicit none
+
+ integer :: ierr
+
+ type (domain_type), intent(inout) :: domain
+
+ if (restart_frame > 1) call output_state_finalize(restart_obj, domain % dminfo)
+
+ call MPAS_destroyClock(clock, ierr)
+
+ end subroutine mpas_core_finalize
+
+
+ subroutine compute_mesh_scaling(mesh)
+
+ use grid_types
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer :: iEdge, cell1, cell2
+ real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4
+
+ meshDensity => mesh % meshDensity % array
+ meshScalingDel2 => mesh % meshScalingDel2 % array
+ meshScalingDel4 => mesh % meshScalingDel4 % array
+
+ !
+ ! Compute the scaling factors to be used in the del2 and del4 dissipation
+ !
+ meshScalingDel2(:) = 1.0
+ meshScalingDel4(:) = 1.0
+ if (config_h_ScaleWithMesh) then
+ do iEdge=1,mesh%nEdges
+ cell1 = mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = mesh % cellsOnEdge % array(2,iEdge)
+ meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/12.0)
+ meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/6.0)
+ end do
+ end if
+
+ end subroutine compute_mesh_scaling
+
+end module mpas_core
Copied: branches/source_renaming/src/core_sw/mpas_sw_test_cases.F (from rev 1084, branches/source_renaming/src/core_sw/module_test_cases.F)
===================================================================
--- branches/source_renaming/src/core_sw/mpas_sw_test_cases.F         (rev 0)
+++ branches/source_renaming/src/core_sw/mpas_sw_test_cases.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -0,0 +1,527 @@
+module test_cases
+
+ use grid_types
+ use configure
+ use constants
+
+
+ contains
+
+
+ subroutine setup_sw_test_case(domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Configure grid metadata and model state for the shallow water test case
+ ! specified in the namelist
+ !
+ ! Output: block - a subset (not necessarily proper) of the model domain to be
+ ! initialized
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i
+ type (block_type), pointer :: block_ptr
+
+ if (config_test_case == 0) then
+ write(0,*) 'Using initial conditions supplied in input file'
+
+ else if (config_test_case == 1) then
+ write(0,*) 'Setting up shallow water test case 1'
+ write(0,*) ' -- Advection of Cosine Bell over the Pole'
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call sw_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
+ end do
+
+ block_ptr => block_ptr % next
+ end do
+
+ else if (config_test_case == 2) then
+ write(0,*) 'Setting up shallow water test case 2'
+ write(0,*) ' -- Setup shallow water test case 2: Global Steady State Nonlinear Zonal Geostrophic Flow'
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call sw_test_case_2(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
+ end do
+
+ block_ptr => block_ptr % next
+ end do
+
+ else if (config_test_case == 5) then
+ write(0,*) 'Setting up shallow water test case 5'
+ write(0,*) ' -- Setup shallow water test case 5: Zonal Flow over an Isolated Mountain'
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call sw_test_case_5(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
+ end do
+
+ block_ptr => block_ptr % next
+ end do
+
+ else if (config_test_case == 6) then
+ write(0,*) 'Setting up shallow water test case 6'
+ write(0,*) ' -- Rossby-Haurwitz Wave'
+
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ call sw_test_case_6(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
+ end do
+
+ block_ptr => block_ptr % next
+ end do
+
+ else
+ write(0,*) 'Only test case 1, 2, 5, and 6 are currently supported.'
+ stop
+ end if
+
+ end subroutine setup_sw_test_case
+
+
+ subroutine sw_test_case_1(grid, state)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Setup shallow water test case 1: Advection of Cosine Bell over the Pole
+ !
+ ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
+ ! Approximations to the Shallow Water Equations in Spherical
+ ! Geometry" J. of Comp. Phys., 102, pp. 211--224
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+ type (state_type), intent(inout) :: state
+
+ real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
+ real (kind=RKIND), parameter :: h0 = 1000.0
+ real (kind=RKIND), parameter :: theta_c = 0.0
+ real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
+ real (kind=RKIND), parameter :: alpha = pii/4.0
+
+ 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
+ !
+ grid % xCell % array = grid % xCell % array * a
+ grid % yCell % array = grid % yCell % array * a
+ grid % zCell % array = grid % zCell % array * a
+ grid % xVertex % array = grid % xVertex % array * a
+ grid % yVertex % array = grid % yVertex % array * a
+ grid % zVertex % array = grid % zVertex % array * a
+ grid % xEdge % array = grid % xEdge % array * a
+ grid % yEdge % array = grid % yEdge % array * a
+ grid % zEdge % array = grid % zEdge % array * a
+ grid % dvEdge % array = grid % dvEdge % array * a
+ grid % dcEdge % array = grid % dcEdge % array * a
+ 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
+ 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 cosine bell at (theta_c, lambda_c)
+ !
+ do iCell=1,grid % nCells
+ r = sphere_distance(theta_c, lambda_c, grid % latCell % array(iCell), grid % lonCell % array(iCell), a)
+ if (r < a/3.0) then
+ state % h % array(1,iCell) = (h0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
+ else
+ state % h % array(1,iCell) = 0.0
+ end if
+ end do
+
+ end subroutine sw_test_case_1
+
+
+ subroutine sw_test_case_2(grid, state)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Setup shallow water test case 2: Global Steady State Nonlinear Zonal
+ ! Geostrophic Flow
+ !
+ ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
+ ! Approximations to the Shallow Water Equations in Spherical
+ ! Geometry" J. of Comp. Phys., 102, pp. 211--224
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+ type (state_type), intent(inout) :: state
+
+ real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
+ real (kind=RKIND), parameter :: gh0 = 29400.0
+ real (kind=RKIND), parameter :: alpha = 0.0
+
+ integer :: iCell, iEdge, iVtx
+ real (kind=RKIND) :: u, v
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex
+
+
+ !
+ ! Scale all distances and areas from a unit sphere to one with radius a
+ !
+ grid % xCell % array = grid % xCell % array * a
+ grid % yCell % array = grid % yCell % array * a
+ grid % zCell % array = grid % zCell % array * a
+ grid % xVertex % array = grid % xVertex % array * a
+ grid % yVertex % array = grid % yVertex % array * a
+ grid % zVertex % array = grid % zVertex % array * a
+ grid % xEdge % array = grid % xEdge % array * a
+ grid % yEdge % array = grid % yEdge % array * a
+ grid % zEdge % array = grid % zEdge % array * a
+ grid % dvEdge % array = grid % dvEdge % array * a
+ grid % dcEdge % array = grid % dcEdge % array * a
+ 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
+ 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)
+
+ !
+ ! Generate rotated Coriolis field
+ !
+ do iEdge=1,grid % nEdges
+ grid % fEdge % array(iEdge) = 2.0 * omega * &
+ ( -cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha) + &
+ sin(grid%latEdge%array(iEdge)) * cos(alpha) &
+ )
+ end do
+ do iVtx=1,grid % nVertices
+ grid % fVertex % array(iVtx) = 2.0 * omega * &
+ (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
+ sin(grid%latVertex%array(iVtx)) * cos(alpha) &
+ )
+ end do
+
+ !
+ ! Initialize height field (actually, fluid thickness field)
+ !
+ do iCell=1,grid % nCells
+ state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
+ (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
+ sin(grid%latCell%array(iCell)) * cos(alpha) &
+ )**2.0 &
+ ) / &
+ gravity
+ end do
+
+ end subroutine sw_test_case_2
+
+
+ subroutine sw_test_case_5(grid, state)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Setup shallow water test case 5: Zonal Flow over an Isolated Mountain
+ !
+ ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
+ ! Approximations to the Shallow Water Equations in Spherical
+ ! Geometry" J. of Comp. Phys., 102, pp. 211--224
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+ type (state_type), intent(inout) :: state
+
+ real (kind=RKIND), parameter :: u0 = 20.
+ real (kind=RKIND), parameter :: gh0 = 5960.0*gravity
+ real (kind=RKIND), parameter :: hs0 = 2000.
+ real (kind=RKIND), parameter :: theta_c = pii/6.0
+ real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
+ real (kind=RKIND), parameter :: rr = pii/9.0
+ real (kind=RKIND), parameter :: alpha = 0.0
+
+ 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
+ !
+ grid % xCell % array = grid % xCell % array * a
+ grid % yCell % array = grid % yCell % array * a
+ grid % zCell % array = grid % zCell % array * a
+ grid % xVertex % array = grid % xVertex % array * a
+ grid % yVertex % array = grid % yVertex % array * a
+ grid % zVertex % array = grid % zVertex % array * a
+ grid % xEdge % array = grid % xEdge % array * a
+ grid % yEdge % array = grid % yEdge % array * a
+ grid % zEdge % array = grid % zEdge % array * a
+ grid % dvEdge % array = grid % dvEdge % array * a
+ grid % dcEdge % array = grid % dcEdge % array * a
+ 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
+ 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)
+
+ !
+ ! Generate rotated Coriolis field
+ !
+ do iEdge=1,grid % nEdges
+ grid % fEdge % array(iEdge) = 2.0 * omega * &
+ (-cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha) + &
+ sin(grid%latEdge%array(iEdge)) * cos(alpha) &
+ )
+ end do
+ do iVtx=1,grid % nVertices
+ grid % fVertex % array(iVtx) = 2.0 * omega * &
+ (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
+ sin(grid%latVertex%array(iVtx)) * cos(alpha) &
+ )
+ end do
+
+ !
+ ! Initialize mountain
+ !
+ do iCell=1,grid % nCells
+ if (grid % lonCell % array(iCell) < 0.0) grid % lonCell % array(iCell) = grid % lonCell % array(iCell) + 2.0 * pii
+ r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
+ grid % h_s % array(iCell) = hs0 * (1.0 - r/rr)
+ end do
+
+ !
+ ! Initialize tracer fields
+ !
+ do iCell=1,grid % nCells
+ r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
+ state % tracers % array(1,1,iCell) = 1.0 - r/rr
+ end do
+ if (grid%nTracers > 1) then
+ do iCell=1,grid % nCells
+ r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + &
+ (grid % latCell % array(iCell) - theta_c - pii/6.0)**2.0 &
+ ) &
+ )
+ state % tracers % array(2,1,iCell) = 1.0 - r/rr
+ end do
+ end if
+
+ !
+ ! Initialize height field (actually, fluid thickness field)
+ !
+ do iCell=1,grid % nCells
+ state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
+ (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
+ sin(grid%latCell%array(iCell)) * cos(alpha) &
+ )**2.0 &
+ ) / &
+ gravity
+ state % h % array(1,iCell) = state % h % array(1,iCell) - grid % h_s % array(iCell)
+ end do
+
+ end subroutine sw_test_case_5
+
+
+ subroutine sw_test_case_6(grid, state)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Setup shallow water test case 6: Rossby-Haurwitz Wave
+ !
+ ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
+ ! Approximations to the Shallow Water Equations in Spherical
+ ! Geometry" J. of Comp. Phys., 102, pp. 211--224
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+ type (state_type), intent(inout) :: state
+
+ real (kind=RKIND), parameter :: h0 = 8000.0
+ real (kind=RKIND), parameter :: w = 7.848e-6
+ real (kind=RKIND), parameter :: K = 7.848e-6
+ real (kind=RKIND), parameter :: R = 4.0
+
+ integer :: iCell, iEdge, iVtx
+ real (kind=RKIND) :: u, v
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex
+
+
+ !
+ ! Scale all distances and areas from a unit sphere to one with radius a
+ !
+ grid % xCell % array = grid % xCell % array * a
+ grid % yCell % array = grid % yCell % array * a
+ grid % zCell % array = grid % zCell % array * a
+ grid % xVertex % array = grid % xVertex % array * a
+ grid % yVertex % array = grid % yVertex % array * a
+ grid % zVertex % array = grid % zVertex % array * a
+ grid % xEdge % array = grid % xEdge % array * a
+ grid % yEdge % array = grid % yEdge % array * a
+ grid % zEdge % array = grid % zEdge % array * a
+ grid % dvEdge % array = grid % dvEdge % array * a
+ grid % dcEdge % array = grid % dcEdge % array * a
+ 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 * 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
+ 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*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
+
+
+ real function sphere_distance(lat1, lon1, lat2, lon2, radius)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
+ ! sphere with given radius.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
+
+ real (kind=RKIND) :: arg1
+
+ arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
+ cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
+ sphere_distance = 2.*radius*asin(arg1)
+
+ end function sphere_distance
+
+
+ real function AA(theta)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! A, used in height field computation for Rossby-Haurwitz wave
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), parameter :: w = 7.848e-6
+ real (kind=RKIND), parameter :: K = 7.848e-6
+ real (kind=RKIND), parameter :: R = 4.0
+
+ 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))
+
+ end function AA
+
+
+ real function BB(theta)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! B, used in height field computation for Rossby-Haurwitz wave
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), parameter :: w = 7.848e-6
+ real (kind=RKIND), parameter :: K = 7.848e-6
+ real (kind=RKIND), parameter :: R = 4.0
+
+ 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)*cos(theta))**2.0)
+
+ end function BB
+
+
+ real function CC(theta)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! C, used in height field computation for Rossby-Haurwitz wave
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), parameter :: w = 7.848e-6
+ real (kind=RKIND), parameter :: K = 7.848e-6
+ real (kind=RKIND), parameter :: R = 4.0
+
+ real (kind=RKIND), intent(in) :: theta
+
+ CC = 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 - R - 2.0)
+
+ end function CC
+
+end module test_cases
Copied: branches/source_renaming/src/core_sw/mpas_sw_time_integration.F (from rev 1084, branches/source_renaming/src/core_sw/module_time_integration.F)
===================================================================
--- branches/source_renaming/src/core_sw/mpas_sw_time_integration.F         (rev 0)
+++ branches/source_renaming/src/core_sw/mpas_sw_time_integration.F        2011-10-17 15:44:46 UTC (rev 1087)
@@ -0,0 +1,1287 @@
+module time_integration
+
+ use vector_reconstruction
+ use grid_types
+ use configure
+ use constants
+ use dmpar
+
+
+ contains
+
+
+ subroutine sw_timestep(domain, dt, timeStamp)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Advance model state forward in time by the specified time step
+ !
+ ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
+ ! plus grid meta-data
+ ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
+ ! model state advanced forward in time by dt seconds
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ real (kind=RKIND), intent(in) :: dt
+ character(len=*), intent(in) :: timeStamp
+
+ type (block_type), pointer :: block
+
+ if (trim(config_time_integration) == 'RK4') then
+ call sw_rk4(domain, dt)
+ else
+ write(0,*) 'Unknown time integration option '//trim(config_time_integration)
+ write(0,*) 'Currently, only ''RK4'' is supported.'
+ stop
+ end if
+
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(2) % state % xtime % scalar = timeStamp
+ block => block % next
+ end do
+
+ end subroutine sw_timestep
+
+
+ subroutine sw_rk4(domain, dt)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Advance model state forward in time by the specified time step using
+ ! 4th order Runge-Kutta
+ !
+ ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
+ ! plus grid meta-data
+ ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
+ ! model state advanced forward in time by dt seconds
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ real (kind=RKIND), intent(in) :: dt
+
+ integer :: iCell, k
+ type (block_type), pointer :: block
+ type (state_type) :: provis
+
+ integer :: rk_step
+
+ real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
+
+ block => domain % blocklist
+ call allocate_state(provis, &
+ block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
+ block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels, &
+ block % mesh % nTracers)
+
+ !
+ ! Initialize time_levs(2) with state at current time
+ ! Initialize first RK state
+ ! Couple tracers time_levs(2) with h in time-levels
+ ! Initialize RK weights
+ !
+ block => domain % blocklist
+ do while (associated(block))
+
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ do iCell=1,block % mesh % nCells ! couple tracers to h
+ do k=1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ * block % state % time_levs(1) % state % h % array(k,iCell)
+ end do
+ end do
+
+ call copy_state(provis, block % state % time_levs(1) % state)
+
+ block => block % next
+ end do
+
+ rk_weights(1) = dt/6.
+ rk_weights(2) = dt/3.
+ rk_weights(3) = dt/3.
+ rk_weights(4) = dt/6.
+
+ rk_substep_weights(1) = dt/2.
+ rk_substep_weights(2) = dt/2.
+ rk_substep_weights(3) = dt
+ rk_substep_weights(4) = 0.
+
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! BEGIN RK loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do rk_step = 1, 4
+
+! --- update halos for diagnostic variables
+
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ end if
+
+ block => block % next
+ end do
+
+! --- compute tendencies
+
+ block => domain % blocklist
+ do while (associated(block))
+ call sw_compute_tend(block % tend, provis, block % mesh)
+ call sw_compute_scalar_tend(block % tend, provis, block % mesh)
+ call sw_enforce_boundary_edge(block % tend, block % mesh)
+ block => block % next
+ end do
+
+! --- update halos for prognostic variables
+
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % h % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(domain % dminfo, block % tend % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+
+! --- compute next substep state
+
+ if (rk_step < 4) then
+ block => domain % blocklist
+ do while (associated(block))
+ provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
+ provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ provis % tracers % array(:,k,iCell) = ( &
+ block % state % time_levs(1) % state % h % array(k,iCell) * &
+ block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
+ ) / provis % h % array(k,iCell)
+ end do
+ end do
+ if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ end if
+ call sw_compute_solve_diagnostics(dt, provis, block % mesh)
+ block => block % next
+ end do
+ end if
+
+!--- accumulate update (for RK4)
+
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
+ + rk_weights(rk_step) * block % tend % u % array(:,:)
+ block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
+ + rk_weights(rk_step) * block % tend % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
+ end do
+ end do
+ block => block % next
+ end do
+
+ end do
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! END RK loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+ !
+ ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
+ !
+ block => domain % blocklist
+ do while (associated(block))
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ / block % state % time_levs(2) % state % h % array(k,iCell)
+ end do
+ end do
+
+ if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ end if
+
+ call sw_compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+
+ call reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
+ block % state % time_levs(2) % state % uReconstructX % array, &
+ block % state % time_levs(2) % state % uReconstructY % array, &
+ block % state % time_levs(2) % state % uReconstructZ % array, &
+ block % state % time_levs(2) % state % uReconstructZonal % array, &
+ block % state % time_levs(2) % state % uReconstructMeridional % array &
+ )
+
+ block => block % next
+ end do
+
+ call deallocate_state(provis)
+
+ end subroutine sw_rk4
+
+
+ subroutine sw_compute_tend(tend, s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute height and normal wind tendencies, as well as diagnostic variables
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
+ real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
+
+ integer :: nCells, nEdges, nVertices, nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ meshScalingDel2, meshScalingDel4
+ real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
+ circulation, vorticity, ke, pv_edge, divergence, h_vertex
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ real (kind=RKIND) :: r, u_diffusion
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+ real (kind=RKIND) :: ke_edge
+
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ pv_edge => s % pv_edge % array
+ vh => s % vh % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+
+ tend_h => tend % h % array
+ tend_u => tend % u % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ u_src => grid % u_src % array
+
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
+
+ !
+ ! Compute height tendency for each cell
+ !
+ tend_h(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = 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
+ end do
+ end do
+ do iCell=1,grid % nCellsSolve
+ do k=1,nVertLevels
+ tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ end do
+ end do
+
+#ifdef LANL_FORMULATION
+ !
+ ! Compute u (normal) velocity tendency for each edge (cell face)
+ !
+ tend_u(:,:) = 0.0
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = 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)
+ end do
+
+ tend_u(k,iEdge) = &
+ q &
+ - ( ke(k,cell2) - ke(k,cell1) + &
+ gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
+ ) / dcEdge(iEdge)
+ end do
+ end do
+
+
+#endif
+
+#ifdef NCAR_FORMULATION
+ !
+ ! Compute u (normal) velocity tendency for each edge (cell face)
+ !
+ tend_u(:,:) = 0.0
+ do iEdge=1,grid % nEdgesSolve
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+ vorticity_abs = fEdge(iEdge) + (circulation(k,vertex1) + circulation(k,vertex2)) / &
+ (areaTriangle(vertex1) + areaTriangle(vertex2))
+
+ workpv = 2.0 * vorticity_abs / (h(k,cell1) + h(k,cell2))
+
+ tend_u(k,iEdge) = workpv * vh(k,iEdge) - &
+ (ke(k,cell2) - ke(k,cell1) + &
+ gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
+ ) / &
+ dcEdge(iEdge)
+ end do
+ end do
+#endif
+
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="blue">abla vorticity
+ ! only valid for visc == constant
+ if (config_h_mom_eddy_visc2 > 0.0) then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+ end do
+ end do
+ end if
+
+ !
+ ! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="blue">abla^4 u
+ ! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ ! applied recursively.
+ ! strictly only valid for h_mom_eddy_visc4 == constant
+ !
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ allocate(delsq_divergence(nVertLevels, nCells+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
+
+ delsq_u(:,:) = 0.0
+
+ ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+
+ delsq_u(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+
+ end do
+ end do
+
+ ! vorticity using </font>
<font color="blue">abla^2 u
+ delsq_circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
+ - dcEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
+ + dcEdge(iEdge) * delsq_u(k,iEdge)
+ end do
+ end do
+ do iVertex=1,nVertices
+ r = 1.0 / areaTriangle(iVertex)
+ do k=1,nVertLevels
+ delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
+ end do
+ end do
+
+ ! Divergence using </font>
<font color="blue">abla^2 u
+ delsq_divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+ + delsq_u(k,iEdge)*dvEdge(iEdge)
+ delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
+ - delsq_u(k,iEdge)*dvEdge(iEdge)
+ end do
+ end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,nVertLevels
+ delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
+ end do
+ end do
+
+ ! Compute - \kappa </font>
<font color="blue">abla^4 u
+ ! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="blue">abla^2 u) )
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k=1,nVertLevels
+
+ u_diffusion = ( delsq_divergence(k,cell2) &
+ - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( delsq_vorticity(k,vertex2) &
+ - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
+
+ end do
+ end do
+
+ deallocate(delsq_divergence)
+ deallocate(delsq_u)
+ deallocate(delsq_circulation)
+ deallocate(delsq_vorticity)
+
+ end if
+
+ ! Compute u (velocity) tendency from wind stress (u_src)
+ if(config_wind_stress) then
+ do iEdge=1,grid % nEdges
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+ end do
+ endif
+
+ if (config_bottom_drag) then
+ do iEdge=1,grid % nEdges
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
+ ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
+ + ke(1,cellsOnEdge(2,iEdge)))
+
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ - 1.0e-3*u(1,iEdge) &
+ *sqrt(2.0*ke_edge)/h_edge(1,iEdge)
+ end do
+ endif
+
+ end subroutine sw_compute_tend
+
+
+ subroutine sw_compute_scalar_tend(tend, s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed scalar tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer :: iCell, iEdge, k, iTracer, cell1, cell2, i
+ real (kind=RKIND) :: flux, tracer_edge, r
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
+ integer, dimension(:,:), pointer :: boundaryEdge
+ real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
+
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
+ integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND), dimension(:,:), pointer :: u, h_edge
+
+ u => s % u % array
+ h_edge => s % h_edge % array
+ dcEdge => grid % dcEdge % array
+ deriv_two => grid % deriv_two % array
+ dvEdge => grid % dvEdge % array
+ tracers => s % tracers % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ boundaryCell=> grid % boundaryCell % array
+ boundaryEdge=> grid % boundaryEdge % array
+ areaCell => grid % areaCell % array
+ tracer_tend => tend % tracers % array
+
+ coef_3rd_order = 0.
+ if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
+ if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
+
+ tracer_tend(:,:,:) = 0.0
+
+ if (config_tracer_adv_order == 2) then
+
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ 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) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ end do
+ end do
+ end if
+ end do
+
+ else if (config_tracer_adv_order == 3) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ !-- if a cell not on the most outside ring of the halo
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+
+ do k=1,grid % nVertLevels
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ do iTracer=1,grid % nTracers
+
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
+
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ endif
+
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ !-- else u <= 0:
+ else
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ end if
+
+ !-- update tendency
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
+ end if
+ end do
+
+ else if (config_tracer_adv_order == 4) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ !-- if an edge is not on the outer-most ring of the halo
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+
+ do k=1,grid % nVertLevels
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ do iTracer=1,grid % nTracers
+
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
+
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ endif
+
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
+
+ !-- update tendency
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
+ end if
+ end do
+
+ endif ! if (config_tracer_adv_order == 2 )
+
+ !
+ ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="blue">abla \phi)
+ !
+ if ( config_h_tracer_eddy_diff2 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ 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="blue">abla \phi on edge
+ tracer_turb_flux = config_h_tracer_eddy_diff2 &
+ *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
+
+ ! div(h \kappa_2 </font>
<font color="blue">abla \phi) at cell center
+ 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
+ end do
+
+ end do
+
+ deallocate(boundaryMask)
+
+ end if
+
+ !
+ ! tracer tendency: del4 horizontal tracer diffusion, &
+ ! div(h \kappa_4 </font>
<font color="black">abla [div(h </font>
<font color="blue">abla \phi)])
+ !
+ if ( config_h_tracer_eddy_diff4 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
+
+ delsq_tracer(:,:,:) = 0.
+
+ ! first del2: div(h </font>
<font color="blue">abla \phi) at cell center
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1, grid % nTracers
+ delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ + 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) &
+ - 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
+ r = 1.0 / grid % areaCell % array(iCell)
+ do k=1,grid % nVertLevels
+ do iTracer=1,grid % nTracers
+ delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
+ end do
+ end do
+ end do
+
+ ! second del2: div(h </font>
<font color="blue">abla [delsq_tracer]) at cell center
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ invAreaCell1 = 1.0 / grid % areaCell % array(cell1)
+ invAreaCell2 = 1.0 / grid % areaCell % array(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 = 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
+ enddo
+
+ end do
+
+ deallocate(delsq_tracer)
+ deallocate(boundaryMask)
+
+ end if
+
+ end subroutine sw_compute_scalar_tend
+
+
+ subroutine sw_compute_solve_diagnostics(dt, s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields used in the tendency computations
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: dt
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
+ real (kind=RKIND) :: flux, vorticity_abs, workpv
+
+ integer :: nCells, nEdges, nVertices, nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
+ circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
+ h_vertex, vorticity_cell
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ real (kind=RKIND) :: r, h1, h2
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: coef_3rd_order
+
+ h => s % h % array
+ u => s % u % array
+ v => s % v % array
+ vh => s % vh % array
+ h_edge => s % h_edge % array
+ h_vertex => s % h_vertex % array
+ tend_h => s % h % array
+ tend_u => s % u % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ pv_edge => s % pv_edge % array
+ pv_vertex => s % pv_vertex % array
+ pv_cell => s % pv_cell % array
+ vorticity_cell => s % vorticity_cell % array
+ gradPVn => s % gradPVn % array
+ gradPVt => s % gradPVt % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ deriv_two => grid % deriv_two % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ boundaryEdge => grid % boundaryEdge % array
+ boundaryCell => grid % boundaryCell % array
+
+ !
+ ! Find those cells that have an edge on the boundary
+ !
+ boundaryCell(:,:) = 0
+ do iEdge=1,nEdges
+ do k=1,nVertLevels
+ if(boundaryEdge(k,iEdge).eq.1) then
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ boundaryCell(k,cell1) = 1
+ boundaryCell(k,cell2) = 1
+ endif
+ enddo
+ enddo
+
+ !
+ ! Compute height on cell edges at velocity locations
+ ! Namelist options control the order of accuracy of the reconstructed h_edge value
+ !
+
+ coef_3rd_order = 0.
+ if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
+ if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
+ if (config_thickness_adv_order == 2) then
+
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ 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))
+ end do
+ end if
+ end do
+
+ else if (config_thickness_adv_order == 3) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ !-- if a cell not on the most outside ring of the halo
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+
+ do k=1,grid % nVertLevels
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+ end do
+
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ endif
+
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ !-- else u <= 0:
+ else
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ end if
+
+ end do ! do k
+ end if ! if (cell1 <=
+ end do ! do iEdge
+
+ else if (config_thickness_adv_order == 4) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ !-- if a cell not on the most outside ring of the halo
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+
+ do k=1,grid % nVertLevels
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+ end do
+
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ endif
+
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+ end do ! do k
+ end if ! if (cell1 <=
+ end do ! do iEdge
+
+ endif ! if(config_thickness_adv_order == 2)
+
+ !
+ ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
+ ! used to when reading for edges that do not exist
+ !
+ u(:,nEdges+1) = 0.0
+
+ !
+ ! Compute circulation and relative vorticity at each vertex
+ !
+ circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ 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)
+ end do
+ end do
+ do iVertex=1,nVertices
+ do k=1,nVertLevels
+ vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+ end do
+ end do
+
+
+ !
+ ! Compute the divergence at each cell center
+ !
+ divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= nCells) then
+ do k=1,nVertLevels
+ 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)*dvEdge(iEdge)
+ enddo
+ end if
+ end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,nVertLevels
+ divergence(k,iCell) = divergence(k,iCell) * r
+ enddo
+ enddo
+
+ !
+ ! Compute kinetic energy in each cell
+ !
+ ke(:,:) = 0.0
+ do iCell=1,nCells
+ do i=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i,iCell)
+ do k=1,nVertLevels
+ 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
+ ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
+ end do
+ end do
+
+ !
+ ! Compute v (tangential) velocities
+ !
+ v(:,:) = 0.0
+ do iEdge = 1,nEdges
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ do k = 1,nVertLevels
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
+ end do
+ end do
+
+#ifdef NCAR_FORMULATION
+ !
+ ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
+ !
+ vh(:,:) = 0.0
+ do iEdge=1,grid % nEdgesSolve
+ do j=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ do k=1,nVertLevels
+ vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * u(k,eoe) * h_edge(k,eoe)
+ end do
+ end do
+ end do
+#endif
+
+
+ !
+ ! Compute height at vertices, pv at vertices, and average pv to edge locations
+ ! ( 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
+ end do
+
+
+ !
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
+ !
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+ gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
+ dvEdge(iEdge)
+ enddo
+ enddo
+
+ !
+ ! Compute pv at the edges
+ ! ( this computes pv_edge at all edges bounding real cells )
+ !
+ pv_edge(:,:) = 0.0
+ do iVertex = 1,nVertices
+ do i=1,grid % vertexDegree
+ iEdge = edgesOnVertex(i,iVertex)
+ do k=1,nVertLevels
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
+ end do
+ end do
+ end do
+
+
+ !
+ ! Modify PV edge with upstream bias.
+ !
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ enddo
+ enddo
+
+
+ !
+ ! Compute pv at cell centers
+ ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
+ !
+ pv_cell(:,:) = 0.0
+ vorticity_cell(:,:) = 0.0
+ do iVertex = 1, nVertices
+ do i=1,grid % vertexDegree
+ iCell = cellsOnVertex(i,iVertex)
+ if (iCell <= nCells) then
+ do k = 1,nVertLevels
+ pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
+ vorticity_cell(k,iCell) = vorticity_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) / areaCell(iCell)
+ enddo
+ endif
+ enddo
+ enddo
+
+
+ !
+ ! Compute gradient of PV in normal direction
+ ! ( this computes gradPVn for all edges bounding real cells )
+ !
+ gradPVn(:,:) = 0.0
+ do iEdge = 1,nEdges
+ if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
+ do k = 1,nVertLevels
+ gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
+ dcEdge(iEdge)
+ enddo
+ endif
+ enddo
+
+ ! Modify PV edge with upstream bias.
+ !
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ enddo
+ enddo
+
+ !
+ ! set pv_edge = fEdge / h_edge at boundary points
+ !
+ ! if (maxval(boundaryEdge).ge.0) then
+ ! do iEdge = 1,nEdges
+ ! cell1 = cellsOnEdge(1,iEdge)
+ ! cell2 = cellsOnEdge(2,iEdge)
+ ! do k = 1,nVertLevels
+ ! if(boundaryEdge(k,iEdge).eq.1) then
+ ! v(k,iEdge) = 0.0
+ ! if(cell1.gt.0) then
+ ! h1 = h(k,cell1)
+ ! pv_edge(k,iEdge) = fEdge(iEdge) / h1
+ ! h_edge(k,iEdge) = h1
+ ! else
+ ! h2 = h(k,cell2)
+ ! pv_edge(k,iEdge) = fEdge(iEdge) / h2
+ ! h_edge(k,iEdge) = h2
+ ! endif
+ ! endif
+ ! enddo
+ ! enddo
+ ! endif
+
+
+ end subroutine sw_compute_solve_diagnostics
+
+
+ subroutine sw_enforce_boundary_edge(tend, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Enforce any boundary conditions on the normal velocity at each edge
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: tend_u set to zero at boundaryEdge == 1 locations
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (mesh_type), intent(in) :: grid
+
+ integer, dimension(:,:), pointer :: boundaryEdge
+ real (kind=RKIND), dimension(:,:), pointer :: tend_u
+ integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: iEdge, k
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ boundaryEdge => grid % boundaryEdge % array
+ tend_u => tend % u % array
+
+ if(maxval(boundaryEdge).le.0) return
+
+ do iEdge = 1,nEdges
+ do k = 1,nVertLevels
+
+ if(boundaryEdge(k,iEdge).eq.1) then
+ tend_u(k,iEdge) = 0.0
+ endif
+
+ enddo
+ enddo
+
+ end subroutine sw_enforce_boundary_edge
+
+
+end module time_integration
</font>
</pre>