<p><b>duda</b> 2009-12-09 14:29:29 -0700 (Wed, 09 Dec 2009)</p><p>Add 3rd- and 4th-order advection for theta and tracers.<br>
<br>
M src/module_test_cases.F<br>
M src/module_time_integration.F<br>
M src/module_grid_types.F<br>
A src/module_advection.F<br>
M src/module_sw_solver.F<br>
M src/Makefile<br>
M Registry/Registry<br>
</p><hr noshade><pre><font color="gray">Modified: branches/hyd_model/Registry/Registry
===================================================================
--- branches/hyd_model/Registry/Registry 2009-11-23 17:31:54 UTC (rev 79)
+++ branches/hyd_model/Registry/Registry 2009-12-09 21:29:29 UTC (rev 80)
@@ -10,10 +10,12 @@
namelist real sw_model config_v_mom_eddy_visc2 0.0
namelist real sw_model config_h_theta_eddy_visc2 0.0
namelist real sw_model config_v_theta_eddy_visc2 0.0
+namelist integer sw_model config_number_of_sub_steps 4
+namelist integer sw_model config_theta_adv_order 2
+namelist integer sw_model config_scalar_adv_order 2
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
-namelist integer sw_model config_number_of_sub_steps 4
#
# dim type name_in_file name_in_code
@@ -25,6 +27,8 @@
dim nVertices nVertices
dim TWO 2
dim THREE 3
+dim FIFTEEN 15
+dim TWENTYONE 21
dim nVertLevels nVertLevels
dim nTracers nTracers
dim nVertLevelsP1 nVertLevels+1
@@ -129,3 +133,9 @@
var real h_edge_old ( nVertLevels nEdges ) - h_edge_old
var real h_old ( nVertLevels nCells ) - h_old
var real pressure_old ( nVertLevelsP1 nCells ) - pressure_old
+var real tracers_old ( nTracers nVertLevels nCells ) - tracers_old
+
+# Space needed for advection
+var real deriv_two ( FIFTEEN TWO nEdges ) o deriv_two
+var integer advCells ( TWENTYONE nCells ) - advCells
+
Modified: branches/hyd_model/src/Makefile
===================================================================
--- branches/hyd_model/src/Makefile 2009-11-23 17:31:54 UTC (rev 79)
+++ branches/hyd_model/src/Makefile 2009-12-09 21:29:29 UTC (rev 80)
@@ -14,7 +14,8 @@
module_io_input.o \
module_io_output.o \
module_time_integration.o \
- streams.o
+ streams.o \
+ module_advection.o
all: swmodel
@@ -32,12 +33,15 @@
module_io_output.o: module_grid_types.o module_dmpar.o module_sort.o
-module_sw_solver.o: module_configure.o module_grid_types.o module_io_output.o module_time_integration.o module_dmpar.o module_timer.o
+module_sw_solver.o: module_configure.o module_grid_types.o module_io_output.o \
+ module_time_integration.o module_dmpar.o module_timer.o module_advection.o
module_time_integration.o: module_grid_types.o module_configure.o module_dmpar.o
module_test_cases.o: module_grid_types.o module_configure.o module_constants.o module_dmpar.o
+module_advection.o: module_configure.o module_grid_types.o module_constants.o
+
swmodel: $(OBJS)
$(FC) $(LDFLAGS) -o $@ $(OBJS) $(LIBS)
Added: branches/hyd_model/src/module_advection.F
===================================================================
--- branches/hyd_model/src/module_advection.F (rev 0)
+++ branches/hyd_model/src/module_advection.F 2009-12-09 21:29:29 UTC (rev 80)
@@ -0,0 +1,661 @@
+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 (grid_meta), 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 % 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
+
+!---
+
+ 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) <= 0) do_the_cell = .false.
+ end do
+
+
+ if ( .not. do_the_cell ) cycle
+
+! compute poynomial fit for this cell if all needed neighbors exist
+
+ 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
+
+ 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
+
+ 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
+ end do
+
+! fill second derivative stencil for rk advection
+
+ do i=1, grid % nEdgesOnCell % array (iCell)
+ iEdge = grid % EdgesOnCell % array (i,iCell)
+
+
+ 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
+ end do
+
+ end do ! end of loop over cells
+
+ if (debug) 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 = AMAX1(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
+
+end module advection
Modified: branches/hyd_model/src/module_grid_types.F
===================================================================
--- branches/hyd_model/src/module_grid_types.F 2009-11-23 17:31:54 UTC (rev 79)
+++ branches/hyd_model/src/module_grid_types.F 2009-12-09 21:29:29 UTC (rev 80)
@@ -3,7 +3,7 @@
use dmpar
integer, parameter :: nTimeLevs = 2
- integer, parameter :: storageFactor = 2 ! Additional storage used by time integration scheme
+ integer, parameter :: storageFactor = 1 ! Additional storage used by time integration scheme
! Derived type describing info for doing I/O specific to a field
Modified: branches/hyd_model/src/module_sw_solver.F
===================================================================
--- branches/hyd_model/src/module_sw_solver.F 2009-11-23 17:31:54 UTC (rev 79)
+++ branches/hyd_model/src/module_sw_solver.F 2009-12-09 21:29:29 UTC (rev 80)
@@ -6,6 +6,7 @@
use configure
use dmpar
use timer
+ use advection
integer :: output_frame
integer :: restart_frame
@@ -46,6 +47,9 @@
call compute_solver_constants(block_ptr % time_levs(1) % state, block_ptr % mesh)
call compute_state_diagnostics(block_ptr % time_levs(1) % state, block_ptr % mesh)
call compute_solve_diagnostics(dt, block_ptr % time_levs(1) % state, block_ptr % mesh)
+
+ call initialize_advection_rk( block_ptr % mesh )
+
if (.not. config_do_restart) block_ptr % time_levs(1) % state % xtime % scalar = 0.0
block_ptr => block_ptr % next
end do
Modified: branches/hyd_model/src/module_test_cases.F
===================================================================
--- branches/hyd_model/src/module_test_cases.F 2009-11-23 17:31:54 UTC (rev 79)
+++ branches/hyd_model/src/module_test_cases.F 2009-12-09 21:29:29 UTC (rev 80)
@@ -67,6 +67,8 @@
real (kind=RKIND), parameter :: omega_e = 7.29212e-05
real (kind=RKIND), parameter :: t0 = 288., delta_t = 4.8e+05, dtdz = 0.005, eta_t = 0.2
real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
+ real (kind=RKIND), parameter :: theta_c = pii/4.0
+ real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
real (kind=RKIND), dimension(:), pointer :: rdnu, rdnw, fnm, fnp, dbn, dnu, dnw
real (kind=RKIND), dimension(:), pointer :: surface_pressure
@@ -74,7 +76,7 @@
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
integer :: iCell, iEdge, vtx1, vtx2, ivtx, k, nz, nz1
- real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert
+ real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
real (kind=RKIND) :: ptop, p0, phi
@@ -336,10 +338,21 @@
end do
! scalars(2,:,:) = 1. ! transport test
- scalars(2,:,:) = theta ! transport test
- scalars(3,:,:) = theta + 100. ! transport test
- scalars(4,:,:) = theta + 200. ! transport test
- scalars(5,:,:) = theta + 300. ! transport test
+! scalars(2,:,:) = theta ! transport test
+ if (grid % nTracers >= 2) then
+ scalars(2,:,:) = 0.0
+ 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
+ do k=1,grid%nVertLevels
+ scalars(2,k,iCell) = (1.0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
+ end do
+ end if
+ end do
+ end if
+ if (grid % nTracers >= 3) scalars(3,:,:) = theta + 100. ! transport test
+ if (grid % nTracers >= 4) scalars(4,:,:) = theta + 200. ! transport test
+ if (grid % nTracers >= 5) scalars(5,:,:) = theta + 300. ! transport test
end subroutine hyd_test_case_1
Modified: branches/hyd_model/src/module_time_integration.F
===================================================================
--- branches/hyd_model/src/module_time_integration.F 2009-11-23 17:31:54 UTC (rev 79)
+++ branches/hyd_model/src/module_time_integration.F 2009-12-09 21:29:29 UTC (rev 80)
@@ -64,7 +64,6 @@
integer :: iCell, k
type (block_type), pointer :: block
- integer, parameter :: PROVIS = 2
integer, parameter :: TEND = 1
integer :: rk_step, number_of_sub_steps
@@ -166,7 +165,14 @@
block => domain % blocklist
do while (associated(block))
+ !
+ ! Note: Tracers in new time level shouldn't be overwritten, since their provisional values
+ ! from the previous RK step are needed to compute new scalar tendencies in advance_scalars.
+ ! A cleaner way of preserving tracers should be added in future.
+ !
+ block % mesh % tracers_old % array(:,:,:) = block % time_levs(2) % state % tracers % array(:,:,:)
call copy_state( block % time_levs(1) % state, block % time_levs(2) % state )
+ block % time_levs(2) % state % tracers % array(:,:,:) = block % mesh % tracers_old % array(:,:,:)
block => block % next
end do
@@ -307,14 +313,14 @@
block % mesh % areaCell % array (iCell)
tracer_min = min(tracer_min,block % time_levs(2) % state % tracers % array (2,k,iCell))
tracer_max = max(tracer_max,block % time_levs(2) % state % tracers % array (2,k,iCell))
- enddo
- enddo
+ end do
+ end do
block => block % next
- enddo
+ end do
write(0,*) ' mass in the domain = ',domain_mass
write(0,*) ' tracer mass in the domain = ',tracer_mass
write(0,*) ' tracer_min, tracer_max ',tracer_min, tracer_max
- endif
+ end if
end subroutine srk3
@@ -357,9 +363,9 @@
! do iq = 1, nTracers
do iq = 1, 1
grid % qtot % array(k,iCell) = grid % qtot % array(k,iCell) + s % tracers % array (iq, k, iCell)
- enddo
- enddo
- enddo
+ end do
+ end do
+ end do
do iEdge = 1, nEdges
do k = 1, nVertLevels
@@ -368,8 +374,8 @@
if (cell1 > 0 .and. cell2 > 0) then
grid % cqu % array(k,iEdge) = 1./( 1. + 0.5*(grid % qtot % array(k,cell1)+grid % qtot % array(k,cell2)) )
end if
- enddo
- enddo
+ end do
+ end do
end if
@@ -498,11 +504,13 @@
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, divergence, vorticity, ke, pv_edge, geopotential, theta, ww, &
h_diabatic, tend_theta
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND), dimension( grid % nVertLevels + 1 ) :: wdtn, wdun
real (kind=RKIND) :: theta_edge, theta_turb_flux, z1, z2, z3, z4, zm, z0, zp
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:), pointer :: rdnu, rdnw, fnm, fnp
@@ -529,6 +537,7 @@
areaCell => grid % areaCell % array
areaTriangle => grid % areaTriangle % array
fEdge => grid % fEdge % array
+ deriv_two => grid % deriv_two % array
vh => tend % vh % array
tend_u => tend % u % array
@@ -722,29 +731,94 @@
!
! horizontal advection for theta
!
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
- theta_edge = 0.5 * (theta(k,cell1) + theta(k,cell2))
+ if (config_theta_adv_order == 2) then
-! if(u(k,iEdge) .gt. 0.) then
-! theta_edge = theta_edge - (1./6.)*delsq_theta(k,cell1)
-! else
-! theta_edge = theta_edge - (1./6.)*delsq_theta(k,cell2)
-! end if
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+ do k=1,grid % nVertLevels
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) )
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ end do
+ end if
+ end do
- flux = dvEdge (iEdge) * h_edge(k,iEdge) * u(k,iEdge)* theta_edge
- tend_theta(k,cell1) = tend_theta(k,cell1) - flux
- tend_theta(k,cell2) = tend_theta(k,cell2) + flux
- end do
+ else if (config_theta_adv_order == 3) then
- end if
- end do
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,grid % nVertLevels
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+! 3rd order stencil
+ if( u(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
+ else
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
+ end if
+
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+
+ end do
+ end if
+ end do
+ else if (config_theta_adv_order == 4) then
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,grid % nVertLevels
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * theta(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * theta(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * theta(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * theta(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ flux = dvEdge(iEdge) * h_edge(k,iEdge) * u(k,iEdge) * ( &
+ 0.5*(theta(k,cell1) + theta(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
+
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ end do
+
+ end if
+
+ end do
+ end if
+
+
!
! vertical advection plus diabatic term
! Note: we are also dividing through by the cell area after the horizontal flux divergence
@@ -814,7 +888,8 @@
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
integer :: nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, dpsdt, surface_pressure
+ real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, dpsdt, &
+ surface_pressure
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, geopotential, alpha, theta, &
pressure, pressure_old, tend_theta, uhAvg, wwAvg, ww, u_old, &
@@ -1068,12 +1143,14 @@
type (grid_meta), intent(in) :: grid
real (kind=RKIND) :: dt
- integer :: iCell, iEdge, k, iTracer, cell1, cell2
- real (kind=RKIND) :: flux, tracer_edge
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+ real (kind=RKIND) :: flux, tracer_edge, d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: scalar_old, scalar_new, scalar_tend
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND), dimension(:,:), pointer :: uhAvg, h_old, h_new, wwAvg
- real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
+ integer, dimension(:,:), pointer :: cellsOnEdge
real (kind=RKIND), dimension( grid % nTracers, grid % nVertLevels + 1 ) :: wdtn
integer :: nVertLevels
@@ -1082,8 +1159,11 @@
scalar_old => s_old % tracers % array
scalar_new => s_new % tracers % array
+ deriv_two => grid % deriv_two % array
uhAvg => grid % uhAvg % array
dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
+ cellsOnEdge => grid % cellsOnEdge % array
scalar_tend => tend % tracers % array
h_old => s_old % h % array
h_new => s_new % h % array
@@ -1104,9 +1184,11 @@
!
! horizontal flux divergence, accumulate in scalar_tend
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (config_scalar_adv_order == 2) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
if (cell1 > 0 .and. cell2 > 0) then
do k=1,grid % nVertLevels
do iTracer=1,grid % nTracers
@@ -1117,8 +1199,86 @@
end do
end do
end if
- end do
+ end do
+ else if (config_scalar_adv_order == 3) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,grid % nVertLevels
+
+ do iTracer=1,grid % nTracers
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iTracer,k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ if (uhAvg(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
+ 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1) / 6. )
+ else
+ flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
+ 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell2) / 6. )
+ end if
+
+ scalar_tend(iTracer,k,cell1) = scalar_tend(iTracer,k,cell1) - flux/areaCell(cell1)
+ scalar_tend(iTracer,k,cell2) = scalar_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+
+ end do
+ end do
+ end if
+ end do
+
+ else if (config_scalar_adv_order == 4) then
+
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,grid % nVertLevels
+
+ do iTracer=1,grid % nTracers
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * scalar_new(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * scalar_new(iTracer,k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * scalar_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ flux = dvEdge(iEdge) * uhAvg(k,iEdge) * ( &
+ 0.5*(scalar_new(iTracer,k,cell1) + scalar_new(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
+
+ scalar_tend(iTracer,k,cell1) = scalar_tend(iTracer,k,cell1) - flux/areaCell(cell1)
+ scalar_tend(iTracer,k,cell2) = scalar_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ end do
+ end do
+ end if
+
+ end do
+ end if
+
+
!
! vertical flux divergence
!
@@ -1259,20 +1419,20 @@
if (cell1 > 0) then
do k=1,nVertLevels
divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- enddo
- endif
+ end do
+ end if
if(cell2 > 0) then
do k=1,nVertLevels
divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- enddo
+ end do
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
+ end do
+ end do
!
@@ -1340,8 +1500,8 @@
do k = 1,nVertLevels
gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
dvEdge(iEdge)
- enddo
- enddo
+ end do
+ end do
! tdr
!
@@ -1355,8 +1515,8 @@
if(iEdge > 0) then
do k=1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
- enddo
- endif
+ end do
+ end if
end do
end do
! tdr
@@ -1368,8 +1528,8 @@
do iEdge = 1,nEdges
do k = 1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * v(k,iEdge) * dt * gradPVt(k,iEdge)
- enddo
- enddo
+ end do
+ end do
! tdr
@@ -1384,10 +1544,10 @@
if( iCell > 0) then
do k = 1,nVertLevels
pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
- enddo
- endif
- enddo
- enddo
+ end do
+ end if
+ end do
+ end do
! tdr
! tdr
@@ -1401,9 +1561,9 @@
do k = 1,nVertLevels
gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
dcEdge(iEdge)
- enddo
- endif
- enddo
+ end do
+ end if
+ end do
! tdr
! Modify PV edge with upstream bias.
@@ -1411,8 +1571,8 @@
do iEdge = 1,nEdges
do k = 1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) - 0.5 * u(k,iEdge) * dt * gradPVn(k,iEdge)
- enddo
- enddo
+ end do
+ end do
end subroutine compute_solve_diagnostics
</font>
</pre>