<p><b>duda</b> 2010-10-22 18:10:18 -0600 (Fri, 22 Oct 2010)</p><p>BRANCH COMMIT<br>
<br>
Update with respect to trunk.<br>
<br>
<br>
M namelist.input.sw<br>
M src/core_hyd_atmos/mpas_interface.F<br>
M src/core_sw/module_global_diagnostics.F<br>
M src/core_sw/mpas_interface.F<br>
A src/core_sw/module_advection.F<br>
M src/core_sw/Registry<br>
M src/core_sw/module_time_integration.F<br>
M src/core_sw/Makefile<br>
M src/driver/mpas.F<br>
M src/core_nhyd_atmos/mpas_interface.F<br>
M src/core_ocean/mpas_interface.F<br>
A src/core_ocean/module_advection.F<br>
M src/core_ocean/module_test_cases.F<br>
M src/core_ocean/Registry<br>
M src/core_ocean/module_time_integration.F<br>
M src/core_ocean/Makefile<br>
M src/framework/module_io_input.F<br>
M src/framework/module_block_decomp.F<br>
M src/framework/module_dmpar.F<br>
M namelist.input.ocean<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_nonhydrostatic/namelist.input.ocean
===================================================================
--- branches/atmos_nonhydrostatic/namelist.input.ocean 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/namelist.input.ocean 2010-10-23 00:10:18 UTC (rev 580)
@@ -42,6 +42,9 @@
config_vert_diffusion = 1.0e-4
/
&advection
- config_hor_tracer_adv = 'centered'
- config_vert_tracer_adv = 'centered'
+ config_vert_tracer_adv = 'upwind'
+ config_tracer_adv_order = 2
+ config_thickness_adv_order = 2
+ config_positive_definite = .false.
+ config_monotonic = .false.
/
Modified: branches/atmos_nonhydrostatic/namelist.input.sw
===================================================================
--- branches/atmos_nonhydrostatic/namelist.input.sw 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/namelist.input.sw 2010-10-23 00:10:18 UTC (rev 580)
@@ -6,6 +6,10 @@
config_output_interval = 500
config_stats_interval = 0
config_visc = 0.0
+ config_thickness_adv_order = 2
+ config_tracer_adv_order = 2
+ config_positive_definite = .false.
+ config_monotonic = .false.
/
&io
Modified: branches/atmos_nonhydrostatic/src/core_hyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_hyd_atmos/mpas_interface.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_hyd_atmos/mpas_interface.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -12,6 +12,20 @@
end subroutine mpas_setup_test_case
+subroutine mpas_init_domain(domain)
+! Initialize grid variables that are computed from the netcdf input file.
+
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ ! This is currently a stub.
+
+end subroutine mpas_init_domain
+
+
subroutine mpas_init(block, mesh, dt)
use grid_types
Modified: branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_nhyd_atmos/mpas_interface.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -18,6 +18,20 @@
end subroutine mpas_setup_test_case
+subroutine mpas_init_domain(domain)
+! Initialize grid variables that are computed from the netcdf input file.
+
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ ! This is currently a stub.
+
+end subroutine mpas_init_domain
+
+
subroutine mpas_init(block, mesh, dt)
use grid_types
Modified: branches/atmos_nonhydrostatic/src/core_ocean/Makefile
===================================================================
--- branches/atmos_nonhydrostatic/src/core_ocean/Makefile 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_ocean/Makefile 2010-10-23 00:10:18 UTC (rev 580)
@@ -1,6 +1,7 @@
.SUFFIXES: .F .o
OBJS = module_test_cases.o \
+ module_advection.o \
module_time_integration.o \
module_global_diagnostics.o \
mpas_interface.o
@@ -12,11 +13,13 @@
module_test_cases.o:
+module_advection.o:
+
module_time_integration.o:
module_global_diagnostics.o:
-mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o
+mpas_interface.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/atmos_nonhydrostatic/src/core_ocean/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_ocean/Registry 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_ocean/Registry 2010-10-23 00:10:18 UTC (rev 580)
@@ -29,8 +29,11 @@
namelist real vmix config_vmixTanhDiffMin 1.0e-5
namelist real vmix config_vmixTanhZMid -100
namelist real vmix config_vmixTanhZWidth 100
-namelist character advection config_hor_tracer_adv 'centered'
-namelist character advection config_vert_tracer_adv 'centered'
+namelist character advection config_vert_tracer_adv 'centered'
+namelist integer advection config_tracer_adv_order 2
+namelist integer advection config_thickness_adv_order 2
+namelist logical advection config_positive_definite false
+namelist logical advection config_monotonic false
#
# dim type name_in_file name_in_code
@@ -42,6 +45,8 @@
dim nVertices nVertices
dim TWO 2
dim R3 3
+dim FIFTEEN 15
+dim TWENTYONE 21
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
dim nVertLevelsP1 nVertLevels+1
@@ -99,6 +104,17 @@
var persistent real fVertex ( nVertices ) 0 iro fVertex mesh - -
var persistent real h_s ( nCells ) 0 iro h_s mesh - -
+# Space needed for advection
+var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 o deriv_two mesh - -
+var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
+
+# !! NOTE: the following arrays are needed to allow the use
+# !! of the module_advection.F w/o alteration
+# Space needed for deformation calculation weights
+var persistent real defc_a ( maxEdges nCells ) 0 - defc_a mesh - -
+var persistent real defc_b ( maxEdges nCells ) 0 - defc_b mesh - -
+var persistent real kdiff ( nVertLevels nCells Time ) 0 - kdiff mesh - -
+
# Arrays required for reconstruction of velocity field
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
@@ -112,6 +128,7 @@
# Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
+var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
# Prognostic variables: read from input, saved in restart, and written to output
Added: branches/atmos_nonhydrostatic/src/core_ocean/module_advection.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_ocean/module_advection.F (rev 0)
+++ branches/atmos_nonhydrostatic/src/core_ocean/module_advection.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -0,0 +1,934 @@
+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
+! type (grid_meta), 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
Modified: branches/atmos_nonhydrostatic/src/core_ocean/module_test_cases.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_ocean/module_test_cases.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_ocean/module_test_cases.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -25,15 +25,14 @@
type (block_type), pointer :: block_ptr
type (dm_info) :: dminfo
- ! mrp 100507: for diagnostic output
integer :: iTracer
real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
hZLevel, zMidZLevel, zTopZLevel, latCell,LonCell
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
+ real (kind=RKIND) :: centerx, centery
integer :: nCells, nEdges, nVertices, nVertLevels
- ! mrp 100507 end: for diagnostic output
if (config_test_case == 0) then
write(0,*) 'Using initial conditions supplied in input file'
@@ -146,6 +145,9 @@
! Set tracers, if not done in grid.nc file
!tracers = 0.0
+ centerx = 1.0e6
+ centery = 1.0e6
+ dist = 2.0e5
do iCell = 1,nCells
dist = sqrt( (latCell(iCell)-latCenter)**2 + (lonCell(iCell)-lonCenter)**2)
do iLevel = 1,nVertLevels
@@ -161,10 +163,15 @@
! tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = &
! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
- ! Tracer blob
- !if (dist.lt.pi/16) then
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
- !!else
+ !tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = block_ptr % mesh % lonCell % array(iCell)
+ !tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = block_ptr % mesh % latCell % array(iCell)
+
+ ! place tracer blob with radius dist at (centerx,centery)
+ !if ( sqrt( (xCell(iCell)-centerx)**2 &
+ ! + (yCell(iCell)-centery)**2) &
+ ! .lt.dist) then
+ ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
+ !else
! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 0.0
!endif
Modified: branches/atmos_nonhydrostatic/src/core_ocean/module_time_integration.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_ocean/module_time_integration.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_ocean/module_time_integration.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -679,11 +679,10 @@
type (state_type), intent(in) :: s
type (mesh_type), intent(in) :: grid
- integer :: iCell, iEdge, k, iTracer, cell1, cell2, upwindCell,&
- nEdges, nCells, nVertLevels
- real (kind=RKIND) :: h_tracer_eddy_diff2, h_tracer_eddy_diff4, invAreaCell1, invAreaCell2, tracer_turb_flux
- real (kind=RKIND) :: flux, tracer_edge, r
- real (kind=RKIND) :: dist
+ integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&
+ nEdges, nCells, nVertLevels, num_tracers
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
+ real (kind=RKIND) :: flux, tracer_edge, r, dist
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -694,16 +693,21 @@
type (dm_info) :: dminfo
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension(:), pointer :: &
- zTopZLevel
+ integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop, boundaryMask
real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer
real (kind=RKIND), dimension(:), allocatable:: vertDiffTop
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: coef_3rd_order
+
+
u => s % u % array
h => s % h % array
+ boundaryCell=> grid % boundaryCell % array
wTop => s % wTop % array
tracers => s % tracers % array
h_edge => s % h_edge % array
@@ -722,76 +726,159 @@
nEdges = grid % nEdges
nCells = grid % nCells
nVertLevels = grid % nVertLevels
+ num_tracers = s % num_tracers
+ deriv_two => grid % deriv_two % array
- h_tracer_eddy_diff2 = config_h_tracer_eddy_diff2
- h_tracer_eddy_diff4 = config_h_tracer_eddy_diff4
+ !
+ ! initialize tracer tendency (RHS of tracer equation) to zero.
+ !
+ tend_tr(:,:,:) = 0.0
!
! tracer tendency: horizontal advection term -div( h \phi u)
!
- tend_tr(:,:,:) = 0.0
- if (config_hor_tracer_adv.eq.'centered') then
+ 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
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
- do iTracer=1,s % num_tracers
- tracer_edge = 0.5 * ( tracers(iTracer,k,cell1) &
- + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) &
- * tracer_edge
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux
- end do
- end do
- end if
- end do
+ if (config_tracer_adv_order == 2) then
- elseif (config_hor_tracer_adv.eq.'upwind') then
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,nVertLevels
+ do iTracer=1,num_tracers
+ tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ end do
+ end do
+ end if
+ end do
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
- if (u(k,iEdge)>0.0) then
- upwindCell = cell1
- else
- upwindCell = cell2
- endif
- do iTracer=1,s % num_tracers
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) &
- * tracers(iTracer,k,upwindCell)
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux
- end do
- end do
- end if
- end do
+ else if (config_tracer_adv_order == 3) then
- endif
- do iCell=1,grid % nCellsSolve
- do k=1,grid % nVertLevelsSolve
- do iTracer=1,s % num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) / areaCell(iCell)
- end do
+ do iEdge=1,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,nVertLevels
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ do iTracer=1,num_tracers
+
+ !-- 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
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
+ end if
end do
- end do
+ else if (config_tracer_adv_order == 4) then
+
+ do iEdge=1,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,nVertLevels
+
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
+
+ do iTracer=1,num_tracers
+
+ !-- 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
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
+ end if
+ end do
+
+ endif ! if (config_tracer_adv_order == 2 )
+
+
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
- allocate(tracerTop(s % num_tracers,nVertLevels+1))
+ allocate(tracerTop(num_tracers,nVertLevels+1))
tracerTop(:,1)=0.0
tracerTop(:,nVertLevels+1)=0.0
do iCell=1,grid % nCellsSolve
if (config_vert_tracer_adv.eq.'centered') then
do k=2,nVertLevels
- do iTracer=1,s % num_tracers
+ do iTracer=1,num_tracers
tracerTop(iTracer,k) = ( tracers(iTracer,k-1,iCell) &
+tracers(iTracer,k ,iCell))/2.0
end do
@@ -804,7 +891,7 @@
else
upwindCell = k-1
endif
- do iTracer=1,s % num_tracers
+ do iTracer=1,num_tracers
tracerTop(iTracer,k) = tracers(iTracer,upwindCell,iCell)
end do
end do
@@ -812,7 +899,7 @@
endif
do k=1,nVertLevels
- do iTracer=1,s % num_tracers
+ do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- wTop(k+1,iCell)*tracerTop(iTracer,k+1))
@@ -825,7 +912,7 @@
!
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
!
- if ( h_tracer_eddy_diff2 > 0.0 ) then
+ if ( config_h_tracer_eddy_diff2 > 0.0 ) then
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
@@ -841,9 +928,9 @@
invAreaCell2 = 1.0/areaCell(cell2)
do k=1,grid % nVertLevels
- do iTracer=1,s % num_tracers
+ do iTracer=1,num_tracers
! \kappa_2 </font>
<font color="gray">abla \phi on edge
- tracer_turb_flux = h_tracer_eddy_diff2 &
+ tracer_turb_flux = config_h_tracer_eddy_diff2 &
*( tracers(iTracer,k,cell2) &
- tracers(iTracer,k,cell1))/dcEdge(iEdge)
@@ -865,7 +952,7 @@
! tracer tendency: del4 horizontal tracer diffusion, &
! div(h \kappa_4 </font>
<font color="black">abla [div(h </font>
<font color="gray">abla \phi)])
!
- if ( h_tracer_eddy_diff4 > 0.0 ) then
+ if ( config_h_tracer_eddy_diff4 > 0.0 ) then
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
@@ -874,7 +961,7 @@
boundaryMask = 1.0
where(boundaryEdge.eq.1) boundaryMask=0.0
- allocate(delsq_tracer(s % num_tracers,nVertLevels, nCells+1))
+ allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
delsq_tracer(:,:,:) = 0.
@@ -884,7 +971,7 @@
cell2 = grid % cellsOnEdge % array(2,iEdge)
do k=1,grid % nVertLevels
- do iTracer=1,s % num_tracers
+ do iTracer=1,num_tracers
delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ dvEdge(iEdge)*h_edge(k,iEdge) &
*(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) &
@@ -901,7 +988,7 @@
do iCell = 1, nCells
r = 1.0 / areaCell(iCell)
do k=1,nVertLevels
- do iTracer=1,s % num_tracers
+ do iTracer=1,num_tracers
delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
end do
end do
@@ -915,8 +1002,8 @@
invAreaCell2 = 1.0 / areaCell(cell2)
do k=1,grid % nVertLevels
- do iTracer=1,s % num_tracers
- tracer_turb_flux = h_tracer_eddy_diff4 &
+ do iTracer=1,num_tracers
+ 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
@@ -959,12 +1046,12 @@
call dmpar_abort(dminfo)
endif
- allocate(fluxVertTop(s % num_tracers,nVertLevels+1))
+ allocate(fluxVertTop(num_tracers,nVertLevels+1))
fluxVertTop(:,1) = 0.0
fluxVertTop(:,nVertLevels+1) = 0.0
do iCell=1,grid % nCellsSolve
do k=2,nVertLevels
- do iTracer=1,s % num_tracers
+ do iTracer=1,num_tracers
! compute \kappa_v d\phi/dz
fluxVertTop(iTracer,k) = vertDiffTop(k) &
* (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell) )&
@@ -974,7 +1061,7 @@
do k=1,nVertLevels
dist = zTop(k,iCell) - zTop(k+1,iCell)
- do iTracer=1,s % num_tracers
+ do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ h(k,iCell)*(fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1))/dist
enddo
@@ -987,7 +1074,7 @@
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
! 'iTracer,k, minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))'
-! do iTracer=1,s % num_tracers
+! do iTracer=1,num_tracers
! do k = 1,nVertLevels
! print '(2i5,20es12.4)', iTracer,k, &
! minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))
@@ -1032,8 +1119,11 @@
real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: r, h1, h2
@@ -1082,6 +1172,7 @@
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
hZLevel => grid % hZLevel % array
+ deriv_two => grid % deriv_two % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -1089,30 +1180,142 @@
nVertLevels = grid % nVertLevels
boundaryEdge => grid % boundaryEdge % array
+ boundaryCell => grid % boundaryCell % array
!
- ! Compute height on cell edges at velocity locations
+ ! Find those cells that have an edge on the boundary
!
+ boundaryCell(:,:) = 0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- elseif(cell1 <= nCells) then
- do k=1,nVertLevels
- h_edge(k,iEdge) = h(k,cell1)
- end do
- elseif(cell2 <= nCells) then
- do k=1,nVertLevels
- h_edge(k,iEdge) = h(k,cell2)
- end do
- end if
- end do
+ 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
!
Modified: branches/atmos_nonhydrostatic/src/core_ocean/mpas_interface.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_ocean/mpas_interface.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_ocean/mpas_interface.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -11,7 +11,19 @@
end subroutine mpas_setup_test_case
+subroutine mpas_init_domain(domain)
+! Initialize grid variables that are computed from the netcdf input file.
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ ! This is currently a stub.
+
+end subroutine mpas_init_domain
+
subroutine mpas_init(block, mesh, dt)
use grid_types
Modified: branches/atmos_nonhydrostatic/src/core_sw/Makefile
===================================================================
--- branches/atmos_nonhydrostatic/src/core_sw/Makefile 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_sw/Makefile 2010-10-23 00:10:18 UTC (rev 580)
@@ -1,21 +1,26 @@
.SUFFIXES: .F .o
-OBJS = module_test_cases.o \
- module_time_integration.o \
- module_global_diagnostics.o \
- mpas_interface.o
+OBJS = module_test_cases.o \
+ module_advection.o \
+ module_time_integration.o \
+ module_global_diagnostics.o \
+ mpas_interface.o
all: core_sw
core_sw: $(OBJS)
ar -ru libdycore.a $(OBJS)
-module_test_cases.o:
+module_test_cases.o:
-module_time_integration.o:
+module_advection.o:
-mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o
+module_time_integration.o:
+module_global_diagnostics.o:
+
+mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o module_advection.o
+
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/atmos_nonhydrostatic/src/core_sw/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_sw/Registry 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_sw/Registry 2010-10-23 00:10:18 UTC (rev 580)
@@ -8,6 +8,10 @@
namelist integer sw_model config_output_interval 500
namelist integer sw_model config_stats_interval 100
namelist real sw_model config_visc 0.0
+namelist integer sw_model config_thickness_adv_order 2
+namelist integer sw_model config_tracer_adv_order 2
+namelist logical sw_model config_positive_definite false
+namelist logical sw_model config_monotonic false
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -25,6 +29,8 @@
dim nVertices nVertices
dim TWO 2
dim R3 3
+dim FIFTEEN 15
+dim TWENTYONE 21
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
dim nTracers nTracers
@@ -83,12 +89,24 @@
var persistent real fCell ( nCells ) 0 iro fCell mesh - -
var persistent real h_s ( nCells ) 0 iro h_s mesh - -
+# Space needed for advection
+var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 o deriv_two mesh - -
+var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
+
+# !! NOTE: the following arrays are needed to allow the use
+# !! of the module_advection.F w/o alteration
+# Space needed for deformation calculation weights
+var persistent real defc_a ( maxEdges nCells ) 0 - defc_a mesh - -
+var persistent real defc_b ( maxEdges nCells ) 0 - defc_b mesh - -
+var persistent real kdiff ( nVertLevels nCells Time ) 0 - kdiff mesh - -
+
# Arrays required for reconstruction of velocity field
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
# Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
+var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
Added: branches/atmos_nonhydrostatic/src/core_sw/module_advection.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_sw/module_advection.F (rev 0)
+++ branches/atmos_nonhydrostatic/src/core_sw/module_advection.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -0,0 +1,933 @@
+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
Modified: branches/atmos_nonhydrostatic/src/core_sw/module_global_diagnostics.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_sw/module_global_diagnostics.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_sw/module_global_diagnostics.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -179,8 +179,8 @@
end do
keTend_CoriolisForce(iLevel,iEdge) = h_edge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
- iCell1 = cellsOnEdge(iEdge,1)
- iCell2 = cellsOnEdge(iEdge,2)
+ 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)))
@@ -197,8 +197,8 @@
end do
do iEdge = 1,nEdgesSolve
- iCell1 = cellsOnEdge(iEdge,1)
- iCell2 = cellsOnEdge(iEdge,2)
+ iCell1 = cellsOnEdge(1,iEdge)
+ iCell2 = cellsOnEdge(2,iEdge)
h_s_edge(iEdge) = 0.5*(h_s(iCell1) + h_s(iCell2))
end do
Modified: branches/atmos_nonhydrostatic/src/core_sw/module_time_integration.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_sw/module_time_integration.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_sw/module_time_integration.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -404,33 +404,158 @@
type (state_type), intent(in) :: s
type (mesh_type), intent(in) :: grid
- integer :: iCell, iEdge, k, iTracer, cell1, cell2
+ integer :: iCell, iEdge, k, iTracer, cell1, cell2, i
real (kind=RKIND) :: flux, tracer_edge
+ 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
- tend % tracers % array(:,:,:) = 0.0
+ 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
+ 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 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ 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 * (s % tracers % array(iTracer,k,cell1) + s % tracers % array(iTracer,k,cell2))
- flux = s % u % array(k,iEdge) * grid % dvEdge % array(iEdge) * s % h_edge % array(k,iEdge) * tracer_edge
- tend % tracers % array(iTracer,k,cell1) = tend % tracers % array(iTracer,k,cell1) - flux
- tend % tracers % array(iTracer,k,cell2) = tend % tracers % array(iTracer,k,cell2) + flux
+ 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
- do iCell=1,grid % nCellsSolve
- do k=1,grid % nVertLevelsSolve
- do iTracer=1,grid % nTracers
- tend % tracers % array(iTracer,k,iCell) = tend % tracers % array(iTracer,k,iCell) / grid % areaCell % array(iCell)
- 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
- 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 )
+
end subroutine compute_scalar_tend
@@ -458,11 +583,13 @@
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
+ 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
@@ -499,6 +626,7 @@
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
@@ -506,21 +634,141 @@
nVertLevels = grid % nVertLevels
boundaryEdge => grid % boundaryEdge % array
+ boundaryCell => grid % boundaryCell % array
!
- ! Compute height on cell edges at velocity locations
+ ! Find those cells that have an edge on the boundary
!
+ boundaryCell(:,:) = 0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end if
- end do
+ 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
!
Modified: branches/atmos_nonhydrostatic/src/core_sw/mpas_interface.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_sw/mpas_interface.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/core_sw/mpas_interface.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -12,6 +12,20 @@
end subroutine mpas_setup_test_case
+subroutine mpas_init_domain(domain)
+! Initialize grid variables that are computed from the netcdf input file.
+
+ use grid_types
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ ! This is currently a stub.
+
+end subroutine mpas_init_domain
+
+
subroutine mpas_init(block, mesh, dt)
use grid_types
Modified: branches/atmos_nonhydrostatic/src/driver/mpas.F
===================================================================
--- branches/atmos_nonhydrostatic/src/driver/mpas.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/driver/mpas.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -32,6 +32,8 @@
call input_state_for_domain(domain)
+ call mpas_init_domain(domain)
+
if (.not. config_do_restart) call mpas_setup_test_case(domain)
call timer_stop("initialize")
Modified: branches/atmos_nonhydrostatic/src/framework/module_block_decomp.F
===================================================================
--- branches/atmos_nonhydrostatic/src/framework/module_block_decomp.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/framework/module_block_decomp.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -147,7 +147,9 @@
do j=1,maxCells
if (cellsOnEdge(j,i) /= 0) exit
end do
-if (j > maxCells) write(0,*) 'Error in block_decomp_partitioned_edge_list: edge/vertex is not adjacent to any valid cells'
+ if (j > maxCells) &
+ write(0,*) 'Error in block_decomp_partitioned_edge_list: ',&
+ 'edge/vertex is not adjacent to any valid cells'
if (hash_search(h, cellsOnEdge(j,i))) then
lastEdge = lastEdge + 1
edgeIDList(lastEdge) = edgeIDListLocal(i)
@@ -155,14 +157,16 @@
ghostEdgeStart = ghostEdgeStart - 1
edgeIDList(ghostEdgeStart) = edgeIDListLocal(i)
end if
-if (ghostEdgeStart <= lastEdge) then
- write(0,*) 'block_decomp_partitioned_edge_list: Somehow we have more edges than we thought we should.'
-end if
+ if (ghostEdgeStart <= lastEdge) then
+ write(0,*) 'block_decomp_partitioned_edge_list: ',&
+ 'Somehow we have more edges than we thought we should.'
+ end if
end do
-if (ghostEdgeStart /= lastEdge + 1) then
- write(0,*) 'block_decomp_partitioned_edge_list: Somehow we didn''t have enough edges to fill edgeIDList.'
-end if
+ if (ghostEdgeStart /= lastEdge + 1) then
+ write(0,*) 'block_decomp_partitioned_edge_list:',&
+ ' Somehow we didn''t have enough edges to fill edgeIDList.'
+ end if
call hash_destroy(h)
@@ -202,10 +206,11 @@
do j=1,nEdgesOnCell(i)
if (.not. hash_search(h, edgesOnCell(j,i))) then
k = k + 1
-if (k > nEdges) then
- write(0,*) 'block_decomp_all_edges_in_block: Trying to add more edges than expected.'
- return
-end if
+ if (k > nEdges) then
+ write(0,*) 'block_decomp_all_edges_in_block: ',&
+ 'Trying to add more edges than expected.'
+ return
+ end if
edgeList(k) = edgesOnCell(j,i)
call hash_insert(h, edgesOnCell(j,i))
end if
@@ -214,9 +219,10 @@
call hash_destroy(h)
-if (k < nEdges) then
- write(0,*) 'block_decomp_all_edges_in_block: Listed fewer edges than expected.'
-end if
+ if (k < nEdges) then
+ write(0,*) 'block_decomp_all_edges_in_block: ',&
+ 'Listed fewer edges than expected.'
+ end if
end subroutine block_decomp_all_edges_in_block
@@ -236,6 +242,10 @@
call hash_init(h)
do i=1,local_graph_info % nVertices
+ call hash_insert(h, local_graph_info % vertexID(i))
+ end do
+
+ do i=1,local_graph_info % nVertices
do j=1,local_graph_info % nAdjacent(i)
if (local_graph_info % adjacencyList(j,i) /= 0) then
if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
@@ -259,8 +269,9 @@
call hash_init(h)
do i=1,local_graph_info % nVertices
-if (hash_search(h, local_graph_info % vertexID(i))) &
- write(0,*) 'block_decomp_add_halo: There appear to be duplicates in vertexID list.'
+ if (hash_search(h, local_graph_info % vertexID(i))) &
+ write(0,*) 'block_decomp_add_halo: ', &
+ 'There appear to be duplicates in vertexID list.'
call hash_insert(h, local_graph_info % vertexID(i))
local_graph_with_halo % vertexID(i) = local_graph_info % vertexID(i)
local_graph_with_halo % nAdjacent(i) = local_graph_info % nAdjacent(i)
@@ -268,8 +279,9 @@
end do
k = local_graph_with_halo % ghostStart
-if (hash_size(h) /= k-1) &
- write(0,*) 'block_decomp_add_halo: Somehow we don''t have the right number of non-ghost cells.'
+ if (hash_size(h) /= k-1) &
+ write(0,*) 'block_decomp_add_halo: ',&
+ 'Somehow we don''t have the right number of non-ghost cells.'
do i=1,local_graph_info % nVertices
do j=1,local_graph_info % nAdjacent(i)
if (local_graph_info % adjacencyList(j,i) /= 0) then
@@ -281,8 +293,9 @@
end if
end do
end do
-if (local_graph_with_halo % nVerticesTotal /= k-1) &
- write(0,*) 'block_decomp_add_halo: Somehow we don''t have the right number of total cells.'
+ if (local_graph_with_halo % nVerticesTotal /= k-1) &
+ write(0,*) 'block_decomp_add_halo: ',&
+ 'Somehow we don''t have the right number of total cells.'
call hash_destroy(h)
Modified: branches/atmos_nonhydrostatic/src/framework/module_dmpar.F
===================================================================
--- branches/atmos_nonhydrostatic/src/framework/module_dmpar.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/framework/module_dmpar.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -4,6 +4,7 @@
#ifdef _MPI
include 'mpif.h'
+ integer, parameter :: MPI_INTEGERKIND = MPI_INTEGER
#if (RKIND == 8)
integer, parameter :: MPI_REALKIND = MPI_DOUBLE_PRECISION
@@ -64,7 +65,8 @@
dminfo % nprocs = mpi_size
dminfo % my_proc_id = mpi_rank
- write(0,'(a,i5,a,i5,a)') 'task ', mpi_rank, ' of ', mpi_size, ' is running'
+ write(0,'(a,i5,a,i5,a)') 'task ', mpi_rank, ' of ', mpi_size, &
+ ' is running'
call open_streams(dminfo % my_proc_id)
@@ -139,7 +141,7 @@
#ifdef _MPI
integer :: mpi_ierr
- call MPI_Bcast(i, 1, MPI_INTEGER, IO_NODE, dminfo % comm, mpi_ierr)
+ call MPI_Bcast(i, 1, MPI_INTEGERKIND, IO_NODE, dminfo % comm, mpi_ierr)
#endif
end subroutine dmpar_bcast_int
@@ -156,7 +158,7 @@
#ifdef _MPI
integer :: mpi_ierr
- call MPI_Bcast(iarray, n, MPI_INTEGER, IO_NODE, dminfo % comm, mpi_ierr)
+ call MPI_Bcast(iarray, n, MPI_INTEGERKIND, IO_NODE, dminfo % comm, mpi_ierr)
#endif
end subroutine dmpar_bcast_ints
@@ -214,7 +216,7 @@
end if
end if
- call MPI_Bcast(itemp, 1, MPI_INTEGER, IO_NODE, dminfo % comm, mpi_ierr)
+ call MPI_Bcast(itemp, 1, MPI_INTEGERKIND, IO_NODE, dminfo % comm, mpi_ierr)
if (itemp == 1) then
l = .true.
@@ -253,7 +255,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(i, isum, 1, MPI_INTEGER, MPI_SUM, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(i, isum, 1, MPI_INTEGERKIND, MPI_SUM, dminfo % comm, mpi_ierr)
#else
isum = i
#endif
@@ -291,7 +293,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(i, imin, 1, MPI_INTEGER, MPI_MIN, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(i, imin, 1, MPI_INTEGERKIND, MPI_MIN, dminfo % comm, mpi_ierr)
#else
imin = i
#endif
@@ -329,7 +331,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(i, imax, 1, MPI_INTEGER, MPI_MAX, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(i, imax, 1, MPI_INTEGERKIND, MPI_MAX, dminfo % comm, mpi_ierr)
#else
imax = i
#endif
@@ -368,7 +370,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGER, MPI_SUM, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGERKIND, MPI_SUM, dminfo % comm, mpi_ierr)
#else
outArray = inArray
#endif
@@ -388,7 +390,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGER, MPI_MIN, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGERKIND, MPI_MIN, dminfo % comm, mpi_ierr)
#else
outArray = inArray
#endif
@@ -408,7 +410,7 @@
integer :: mpi_ierr
#ifdef _MPI
- call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGER, MPI_MAX, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(inArray, outArray, nElements, MPI_INTEGERKIND, MPI_MAX, dminfo % comm, mpi_ierr)
#else
outArray = inArray
#endif
@@ -489,7 +491,7 @@
#ifdef _MPI
integer :: mpi_ierr
- call MPI_Scatterv(inlist, counts, displs, MPI_INTEGER, outlist, noutlist, MPI_INTEGER, IO_NODE, dminfo % comm, mpi_ierr)
+ call MPI_Scatterv(inlist, counts, displs, MPI_INTEGERKIND, outlist, noutlist, MPI_INTEGERKIND, IO_NODE, dminfo % comm, mpi_ierr)
#endif
end subroutine dmpar_scatter_ints
@@ -532,13 +534,13 @@
#ifdef _MPI
else if (dminfo % my_proc_id == dminfo % nprocs - 1) then
- call MPI_Recv(global_start, 1, MPI_INTEGER, dminfo % my_proc_id - 1, 0, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call MPI_Recv(global_start, 1, MPI_INTEGERKIND, dminfo % my_proc_id - 1, 0, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
global_end = global_start + n - 1
else
- call MPI_Recv(global_start, 1, MPI_INTEGER, dminfo % my_proc_id - 1, 0, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
+ call MPI_Recv(global_start, 1, MPI_INTEGERKIND, dminfo % my_proc_id - 1, 0, dminfo % comm, MPI_STATUS_IGNORE, mpi_ierr)
global_end = global_start + n
- call MPI_Send(global_end, 1, MPI_INTEGER, dminfo % my_proc_id + 1, 0, dminfo % comm, mpi_ierr)
+ call MPI_Send(global_end, 1, MPI_INTEGERKIND, dminfo % my_proc_id + 1, 0, dminfo % comm, mpi_ierr)
global_end = global_end - 1
#endif
@@ -585,7 +587,7 @@
end do
call quicksort(nOwnedList, ownedListSorted)
- call MPI_Allreduce(nNeededList, totalSize, 1, MPI_INTEGER, MPI_MAX, dminfo % comm, mpi_ierr)
+ call MPI_Allreduce(nNeededList, totalSize, 1, MPI_INTEGERKIND, MPI_MAX, dminfo % comm, mpi_ierr)
allocate(ownerListIn(totalSize))
allocate(ownerListOut(totalSize))
@@ -634,12 +636,12 @@
end if
nMesgSend = nMesgRecv
- call MPI_Irecv(nMesgRecv, 1, MPI_INTEGER, recvNeighbor, i, dminfo % comm, mpi_rreq, mpi_ierr)
- call MPI_Isend(nMesgSend, 1, MPI_INTEGER, sendNeighbor, i, dminfo % comm, mpi_sreq, mpi_ierr)
+ call MPI_Irecv(nMesgRecv, 1, MPI_INTEGERKIND, recvNeighbor, i, dminfo % comm, mpi_rreq, mpi_ierr)
+ call MPI_Isend(nMesgSend, 1, MPI_INTEGERKIND, sendNeighbor, i, dminfo % comm, mpi_sreq, mpi_ierr)
call MPI_Wait(mpi_rreq, MPI_STATUS_IGNORE, mpi_ierr)
call MPI_Wait(mpi_sreq, MPI_STATUS_IGNORE, mpi_ierr)
- call MPI_Irecv(ownerListIn, nMesgRecv, MPI_INTEGER, recvNeighbor, i, dminfo % comm, mpi_rreq, mpi_ierr)
- call MPI_Isend(ownerListOut, nMesgSend, MPI_INTEGER, sendNeighbor, i, dminfo % comm, mpi_sreq, mpi_ierr)
+ call MPI_Irecv(ownerListIn, nMesgRecv, MPI_INTEGERKIND, recvNeighbor, i, dminfo % comm, mpi_rreq, mpi_ierr)
+ call MPI_Isend(ownerListOut, nMesgSend, MPI_INTEGERKIND, sendNeighbor, i, dminfo % comm, mpi_sreq, mpi_ierr)
call MPI_Wait(mpi_rreq, MPI_STATUS_IGNORE, mpi_ierr)
call MPI_Wait(mpi_sreq, MPI_STATUS_IGNORE, mpi_ierr)
end do
@@ -741,7 +743,7 @@
do while (associated(recvListPtr))
if (recvListPtr % procID /= dminfo % my_proc_id) then
allocate(recvListPtr % ibuffer(recvListPtr % nlist))
- call MPI_Irecv(recvListPtr % ibuffer, recvListPtr % nlist, MPI_INTEGER, &
+ call MPI_Irecv(recvListPtr % ibuffer, recvListPtr % nlist, MPI_INTEGERKIND, &
recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
end if
recvListPtr => recvListPtr % next
@@ -753,7 +755,7 @@
allocate(sendListPtr % ibuffer(sendListPtr % nlist))
call packSendBuf1dInteger(nOwnedList, arrayIn, sendListPtr, 1, sendListPtr % nlist, &
sendListPtr % ibuffer, nPacked, lastPackedIdx)
- call MPI_Isend(sendListPtr % ibuffer, sendListPtr % nlist, MPI_INTEGER, &
+ call MPI_Isend(sendListPtr % ibuffer, sendListPtr % nlist, MPI_INTEGERKIND, &
sendListPtr % procID, dminfo % my_proc_id, dminfo % comm, sendListPtr % reqID, mpi_ierr)
end if
sendListPtr => sendListPtr % next
@@ -781,7 +783,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
@@ -831,7 +834,7 @@
if (recvListPtr % procID /= dminfo % my_proc_id) then
d2 = dim1 * recvListPtr % nlist
allocate(recvListPtr % ibuffer(d2))
- call MPI_Irecv(recvListPtr % ibuffer, d2, MPI_INTEGER, &
+ call MPI_Irecv(recvListPtr % ibuffer, d2, MPI_INTEGERKIND, &
recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
end if
recvListPtr => recvListPtr % next
@@ -844,7 +847,7 @@
allocate(sendListPtr % ibuffer(d2))
call packSendBuf2dInteger(1, dim1, nOwnedList, arrayIn, sendListPtr, 1, d2, &
sendListPtr % ibuffer, nPacked, lastPackedIdx)
- call MPI_Isend(sendListPtr % ibuffer, d2, MPI_INTEGER, &
+ call MPI_Isend(sendListPtr % ibuffer, d2, MPI_INTEGERKIND, &
sendListPtr % procID, dminfo % my_proc_id, dminfo % comm, sendListPtr % reqID, mpi_ierr)
end if
sendListPtr => sendListPtr % next
@@ -873,7 +876,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
@@ -962,7 +966,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(1:nNeededList) = arrayIn(1:nOwnedList)
@@ -1054,7 +1059,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(:,1:nNeededList) = arrayIn(:,1:nOwnedList)
@@ -1146,7 +1152,8 @@
#else
if (nOwnedList /= nNeededList) then
- write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, arrayIn and arrayOut dims must match.'
+ write(0,*) 'Error in dmpar_alltoall_field: For non-dmpar, ', &
+ 'arrayIn and arrayOut dims must match.'
call dmpar_abort(dminfo)
else
arrayOut(:,:,1:nNeededList) = arrayIn(:,:,1:nOwnedList)
@@ -1198,7 +1205,8 @@
n = de-ds+1
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dInteger: Not enough space in buffer to fit a single slice.'
+ write(0,*) 'packSendBuf2dInteger: Not enough space in buffer', &
+ ' to fit a single slice.'
return
end if
@@ -1217,6 +1225,45 @@
end subroutine packSendBuf2dInteger
+ subroutine packSendBuf3dInteger(d1s, d1e, d2s, d2e, nField, field, sendList, startPackIdx, nBuffer, buffer, nPacked, lastPackedIdx)
+
+ implicit none
+
+ integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startPackIdx
+ integer, dimension(d1s:d1e,d2s:d2e,*), intent(in) :: field
+ type (exchange_list), intent(in) :: sendList
+ integer, dimension(nBuffer), intent(out) :: buffer
+ integer, intent(inout) :: nPacked, lastPackedIdx
+
+ integer :: i, j, k, n
+
+ n = (d1e-d1s+1) * (d2e-d2s+1)
+
+ if (n > nBuffer) then
+ write(0,*) 'packSendBuf3dInteger: Not enough space in buffer', &
+ ' to fit a single slice.'
+ return
+ end if
+
+ nPacked = 0
+ do i=startPackIdx, sendList % nlist
+ nPacked = nPacked + n
+ if (nPacked > nBuffer) then
+ nPacked = nPacked - n
+ lastPackedIdx = i - 1
+ return
+ end if
+ k = nPacked-n+1
+ do j=d2s,d2e
+ buffer(k:k+d1e-d1s) = field(d1s:d1e,j,sendList % list(i))
+ k = k + d1e-d1s+1
+ end do
+ end do
+ lastPackedIdx = sendList % nlist
+
+ end subroutine packSendBuf3dInteger
+
+
subroutine packSendBuf1dReal(nField, field, sendList, startPackIdx, nBuffer, buffer, nPacked, lastPackedIdx)
implicit none
@@ -1259,7 +1306,8 @@
n = de-ds+1
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dReal: Not enough space in buffer to fit a single slice.'
+ write(0,*) 'packSendBuf2dReal: Not enough space in buffer', &
+ ' to fit a single slice.'
return
end if
@@ -1293,7 +1341,8 @@
n = (d1e-d1s+1) * (d2e-d2s+1)
if (n > nBuffer) then
- write(0,*) 'packSendBuf2dReal: Not enough space in buffer to fit a single slice.'
+ write(0,*) 'packSendBuf3dReal: Not enough space in buffer', &
+ ' to fit a single slice.'
return
end if
@@ -1372,6 +1421,230 @@
end subroutine unpackRecvBuf2dInteger
+ subroutine unpackRecvBuf3dInteger(d1s, d1e, d2s, d2e, nField, field, recvList, startUnpackIdx, nBuffer, buffer, &
+ nUnpacked, lastUnpackedIdx)
+
+ implicit none
+
+ integer, intent(in) :: d1s, d1e, d2s, d2e, nField, nBuffer, startUnpackIdx
+ integer, dimension(d1s:d1e,d2s:d2e,*), intent(inout) :: field
+ type (exchange_list), intent(in) :: recvList
+ integer, dimension(nBuffer), intent(in) :: buffer
+ integer, intent(inout) :: nUnpacked, lastUnpackedIdx
+
+ integer :: i, j, k, n
+
+ n = (d1e-d1s+1) * (d2e-d2s+1)
+
+ nUnpacked = 0
+ do i=startUnpackIdx, recvList % nlist
+ nUnpacked = nUnpacked + n
+ if (nUnpacked > nBuffer) then
+ nUnpacked = nUnpacked - n
+ lastUnpackedIdx = i - 1
+ return
+ end if
+ k = nUnpacked-n+1
+ do j=d2s,d2e
+ field(d1s:d1e,j,recvList % list(i)) = buffer(k:k+d1e-d1s)
+ k = k + d1e-d1s+1
+ end do
+ end do
+ lastUnpackedIdx = recvList % nlist
+
+ end subroutine unpackRecvBuf3dInteger
+
+
+ subroutine dmpar_exch_halo_field1dInteger(dminfo, array, dim1, sendList, recvList)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: dim1
+ integer, dimension(*), intent(inout) :: array
+ type (exchange_list), pointer :: sendList, recvList
+
+ type (exchange_list), pointer :: sendListPtr, recvListPtr
+ integer :: lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ integer :: mpi_ierr
+
+#ifdef _MPI
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ allocate(recvListPtr % ibuffer(recvListPtr % nlist))
+ call MPI_Irecv(recvListPtr % ibuffer, recvListPtr % nlist, MPI_INTEGERKIND, &
+ recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ allocate(sendListPtr % ibuffer(sendListPtr % nlist))
+ call packSendBuf1dInteger(dim1, array, sendListPtr, 1, sendListPtr % nlist, sendListPtr % ibuffer, nPacked, lastPackedIdx)
+ call MPI_Isend(sendListPtr % ibuffer, sendListPtr % nlist, MPI_INTEGERKIND, &
+ sendListPtr % procID, dminfo % my_proc_id, dminfo % comm, sendListPtr % reqID, mpi_ierr)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(recvListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ call unpackRecvBuf1dInteger(dim1, array, recvListPtr, 1, recvListPtr % nlist, recvListPtr % ibuffer, nUnpacked, lastUnpackedIdx)
+ deallocate(recvListPtr % ibuffer)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(sendListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ deallocate(sendListPtr % ibuffer)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+#endif
+
+ end subroutine dmpar_exch_halo_field1dInteger
+
+
+ subroutine dmpar_exch_halo_field2dInteger(dminfo, array, dim1, dim2, sendList, recvList)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: dim1, dim2
+ integer, dimension(dim1,*), intent(inout) :: array
+ type (exchange_list), pointer :: sendList, recvList
+
+ type (exchange_list), pointer :: sendListPtr, recvListPtr
+ integer :: lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ integer :: mpi_ierr
+ integer :: d2
+
+#ifdef _MPI
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ d2 = dim1 * recvListPtr % nlist
+ allocate(recvListPtr % ibuffer(d2))
+ call MPI_Irecv(recvListPtr % ibuffer, d2, MPI_INTEGERKIND, &
+ recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ d2 = dim1 * sendListPtr % nlist
+ allocate(sendListPtr % ibuffer(d2))
+ call packSendBuf2dInteger(1, dim1, dim2, array, sendListPtr, 1, d2, sendListPtr % ibuffer, nPacked, lastPackedIdx)
+ call MPI_Isend(sendListPtr % ibuffer, d2, MPI_INTEGERKIND, &
+ sendListPtr % procID, dminfo % my_proc_id, dminfo % comm, sendListPtr % reqID, mpi_ierr)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(recvListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ d2 = dim1 * recvListPtr % nlist
+ call unpackRecvBuf2dInteger(1, dim1, dim2, array, recvListPtr, 1, d2, recvListPtr % ibuffer, nUnpacked, lastUnpackedIdx)
+ deallocate(recvListPtr % ibuffer)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(sendListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ deallocate(sendListPtr % ibuffer)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+#endif
+
+ end subroutine dmpar_exch_halo_field2dInteger
+
+
+ subroutine dmpar_exch_halo_field3dInteger(dminfo, array, dim1, dim2, dim3, sendList, recvList)
+
+ implicit none
+
+ type (dm_info), intent(in) :: dminfo
+ integer, intent(in) :: dim1, dim2, dim3
+ integer, dimension(dim1,dim2,*), intent(inout) :: array
+ type (exchange_list), pointer :: sendList, recvList
+
+ type (exchange_list), pointer :: sendListPtr, recvListPtr
+ integer :: lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ integer :: mpi_ierr
+ integer :: d3
+
+#ifdef _MPI
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ d3 = dim1 * dim2 * recvListPtr % nlist
+ allocate(recvListPtr % ibuffer(d3))
+ call MPI_Irecv(recvListPtr % ibuffer, d3, MPI_INTEGERKIND, &
+ recvListPtr % procID, recvListPtr % procID, dminfo % comm, recvListPtr % reqID, mpi_ierr)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ d3 = dim1 * dim2 * sendListPtr % nlist
+ allocate(sendListPtr % ibuffer(d3))
+ call packSendBuf3dInteger(1, dim1, 1, dim2, dim3, array, sendListPtr, 1, d3, &
+ sendListPtr % ibuffer, nPacked, lastPackedIdx)
+ call MPI_Isend(sendListPtr % ibuffer, d3, MPI_INTEGERKIND, &
+ sendListPtr % procID, dminfo % my_proc_id, dminfo % comm, sendListPtr % reqID, mpi_ierr)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+ recvListPtr => recvList
+ do while (associated(recvListPtr))
+ if (recvListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(recvListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ d3 = dim1 * dim2 * recvListPtr % nlist
+ call unpackRecvBuf3dInteger(1, dim1, 1, dim2, dim3, array, recvListPtr, 1, d3, &
+ recvListPtr % ibuffer, nUnpacked, lastUnpackedIdx)
+ deallocate(recvListPtr % ibuffer)
+ end if
+ recvListPtr => recvListPtr % next
+ end do
+
+ sendListPtr => sendList
+ do while (associated(sendListPtr))
+ if (sendListPtr % procID /= dminfo % my_proc_id) then
+ call MPI_Wait(sendListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ deallocate(sendListPtr % ibuffer)
+ end if
+ sendListPtr => sendListPtr % next
+ end do
+
+#endif
+
+ end subroutine dmpar_exch_halo_field3dInteger
+
+
subroutine unpackRecvBuf1dReal(nField, field, recvList, startUnpackIdx, nBuffer, buffer, nUnpacked, lastUnpackedIdx)
implicit none
@@ -1651,5 +1924,5 @@
end subroutine dmpar_exch_halo_field3dReal
-
+
end module dmpar
Modified: branches/atmos_nonhydrostatic/src/framework/module_io_input.F
===================================================================
--- branches/atmos_nonhydrostatic/src/framework/module_io_input.F 2010-10-22 22:57:52 UTC (rev 579)
+++ branches/atmos_nonhydrostatic/src/framework/module_io_input.F 2010-10-23 00:10:18 UTC (rev 580)
@@ -1127,7 +1127,8 @@
nferr = nf_get_att_real(input_obj % rd_ncid, NF_GLOBAL, attname, attvalue)
end if
if (nferr /= NF_NOERR) then
- write(0,*) 'Warning: Attribute '//trim(attname)//' not found in '//trim(input_obj % filename)
+ write(0,*) 'Warning: Attribute '//trim(attname)//&
+ ' not found in '//trim(input_obj % filename)
if (index(attname, 'sphere_radius') /= 0) then
write(0,*) ' Setting '//trim(attname)//' to 1.0'
attvalue = 1.0
@@ -1151,7 +1152,8 @@
nferr = nf_get_att_text(input_obj % rd_ncid, NF_GLOBAL, attname, attvalue)
if (nferr /= NF_NOERR) then
- write(0,*) 'Warning: Attribute '//trim(attname)//' not found in '//trim(input_obj % filename)
+ write(0,*) 'Warning: Attribute '//trim(attname)//&
+ ' not found in '//trim(input_obj % filename)
if (index(attname, 'on_a_sphere') /= 0) then
write(0,*) ' Setting '//trim(attname)//' to ''YES'''
attvalue = 'YES'
</font>
</pre>