<p><b>xylar@lanl.gov</b> 2010-05-24 13:42:21 -0600 (Mon, 24 May 2010)</p><p>BRANCH COMMIT<br>
<br>
Merging changes from the trunk<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/IBinterp/mpas/Makefile
===================================================================
--- branches/ocean_projects/IBinterp/mpas/Makefile 2010-05-24 19:22:28 UTC (rev 304)
+++ branches/ocean_projects/IBinterp/mpas/Makefile 2010-05-24 19:42:21 UTC (rev 305)
@@ -1,7 +1,7 @@
#MODEL_FORMULATION = -DNCAR_FORMULATION
MODEL_FORMULATION = -DLANL_FORMULATION
-#EXPAND_LEVELS = -DEXPAND_LEVELS=26
+EXPAND_LEVELS = -DEXPAND_LEVELS=26
FILE_OFFSET = -DOFFSET64BIT
#########################
Modified: branches/ocean_projects/IBinterp/mpas/src/core_hyd_atmos/module_advection.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_hyd_atmos/module_advection.F 2010-05-24 19:22:28 UTC (rev 304)
+++ branches/ocean_projects/IBinterp/mpas/src/core_hyd_atmos/module_advection.F 2010-05-24 19:42:21 UTC (rev 305)
@@ -25,6 +25,7 @@
! 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
@@ -102,50 +103,62 @@
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
+ 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. )
+ 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
+ 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) )
- 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
- 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
- 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
-! 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
+ else ! On an x-y plane
- do i=1,n-1
- xp(i) = cos(thetat(i)) * dl_sphere(i)
- yp(i) = sin(thetat(i)) * dl_sphere(i)
- end do
+ do i=1,n-1
+ xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
+ yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+ end do
+ end if
+
+
ma = n-1
mw = grid % nEdgesOnCell % array (iCell)
@@ -244,20 +257,25 @@
yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/a
zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- call arc_bisect( xv1, yv1, zv1, &
- xv2, yv2, zv2, &
- xec, yec, zec )
+ 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
+ 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
- thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
+ 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
@@ -266,31 +284,40 @@
iEdge = grid % EdgesOnCell % array (i,iCell)
- if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
+ 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
+ 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
+ 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(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)
+ 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)
+ deriv_two(j,2,iEdge) = deriv_two(j,1,iEdge)
end do
end if
end do
Modified: branches/ocean_projects/IBinterp/mpas/src/framework/module_grid_types.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/framework/module_grid_types.F 2010-05-24 19:22:28 UTC (rev 304)
+++ branches/ocean_projects/IBinterp/mpas/src/framework/module_grid_types.F 2010-05-24 19:42:21 UTC (rev 305)
@@ -70,6 +70,9 @@
#include "field_dimensions.inc"
+ logical :: on_a_sphere
+ real (kind=RKIND) :: sphere_radius
+
#include "time_invariant_fields.inc"
end type grid_meta
Modified: branches/ocean_projects/IBinterp/mpas/src/framework/module_io_input.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/framework/module_io_input.F 2010-05-24 19:22:28 UTC (rev 304)
+++ branches/ocean_projects/IBinterp/mpas/src/framework/module_io_input.F 2010-05-24 19:42:21 UTC (rev 305)
@@ -49,6 +49,9 @@
integer :: i, j, k
type (io_input_object) :: input_obj
#include "dim_decls.inc"
+
+ character (len=16) :: c_on_a_sphere
+ real (kind=RKIND) :: r_sphere_radius
integer :: readCellStart, readCellEnd, nReadCells
integer :: readEdgeStart, readEdgeEnd, nReadEdges
@@ -729,6 +732,18 @@
#include "dim_dummy_args.inc"
)
+ !
+ ! Read attributes
+ !
+ call io_input_get_att_text(input_obj, 'on_a_sphere', c_on_a_sphere)
+ call io_input_get_att_real(input_obj, 'sphere_radius', r_sphere_radius)
+ if (index(c_on_a_sphere, 'YES') /= 0) then
+ domain % blocklist % mesh % on_a_sphere = .true.
+ else
+ domain % blocklist % mesh % on_a_sphere = .false.
+ end if
+ domain % blocklist % mesh % sphere_radius = r_sphere_radius
+
if (.not. config_do_restart) then
input_obj % time = 1
else
@@ -1100,7 +1115,59 @@
end subroutine io_input_get_dimension
+
+ subroutine io_input_get_att_real(input_obj, attname, attvalue)
+
+ implicit none
+ type (io_input_object), intent(in) :: input_obj
+ character (len=*), intent(in) :: attname
+ real (kind=RKIND), intent(out) :: attvalue
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ if (RKIND == 8) then
+ nferr = nf_get_att_double(input_obj % rd_ncid, NF_GLOBAL, attname, attvalue)
+ else
+ 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)
+ if (index(attname, 'sphere_radius') /= 0) then
+ write(0,*) ' Setting '//trim(attname)//' to 1.0'
+ attvalue = 1.0
+ end if
+ end if
+
+ end subroutine io_input_get_att_real
+
+
+ subroutine io_input_get_att_text(input_obj, attname, attvalue)
+
+ implicit none
+
+ type (io_input_object), intent(in) :: input_obj
+ character (len=*), intent(in) :: attname
+ character (len=*), intent(out) :: attvalue
+
+ include 'netcdf.inc'
+
+ integer :: nferr
+
+ 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)
+ if (index(attname, 'on_a_sphere') /= 0) then
+ write(0,*) ' Setting '//trim(attname)//' to ''YES'''
+ attvalue = 'YES'
+ end if
+ end if
+
+ end subroutine io_input_get_att_text
+
+
subroutine io_input_field0dReal(input_obj, field)
implicit none
Modified: branches/ocean_projects/IBinterp/mpas/src/framework/module_io_output.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/framework/module_io_output.F 2010-05-24 19:22:28 UTC (rev 304)
+++ branches/ocean_projects/IBinterp/mpas/src/framework/module_io_output.F 2010-05-24 19:42:21 UTC (rev 305)
@@ -90,6 +90,7 @@
! although in future, work needs to be done to write model state
! from many distributed blocks
call io_output_init(output_obj, domain % dminfo, &
+ block_ptr % mesh, &
#include "output_dim_actual_args.inc"
)
@@ -326,6 +327,7 @@
subroutine io_output_init( output_obj, &
dminfo, &
+ mesh, &
#include "dim_dummy_args.inc"
)
@@ -335,6 +337,7 @@
type (io_output_object), intent(inout) :: output_obj
type (dm_info), intent(in) :: dminfo
+ type (grid_meta), intent(in) :: mesh
#include "dim_dummy_decls.inc"
integer :: nferr
@@ -348,6 +351,17 @@
#endif
#include "netcdf_def_dims_vars.inc"
+
+ if (mesh % on_a_sphere) then
+ nferr = nf_put_att_text(output_obj % wr_ncid, NF_GLOBAL, 'on_a_sphere', 16, 'YES ')
+ else
+ nferr = nf_put_att_text(output_obj % wr_ncid, NF_GLOBAL, 'on_a_sphere', 16, 'NO ')
+ end if
+ if (RKIND == 8) then
+ nferr = nf_put_att_double(output_obj % wr_ncid, NF_GLOBAL, 'sphere_radius', NF_DOUBLE, 1, mesh % sphere_radius)
+ else
+ nferr = nf_put_att_real(output_obj % wr_ncid, NF_GLOBAL, 'sphere_radius', NF_FLOAT, 1, mesh % sphere_radius)
+ end if
nferr = nf_enddef(output_obj % wr_ncid)
end if
</font>
</pre>