<p><b>mperego@fsu.edu</b> 2012-03-13 15:46:10 -0600 (Tue, 13 Mar 2012)</p><p>BRANCH COMMIT: update temporary greenland test<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/grid_tools/icesheet_legacy/src/module_read_topo.F
===================================================================
--- branches/land_ice_projects/grid_tools/icesheet_legacy/src/module_read_topo.F 2012-03-13 21:41:54 UTC (rev 1627)
+++ branches/land_ice_projects/grid_tools/icesheet_legacy/src/module_read_topo.F 2012-03-13 21:46:10 UTC (rev 1628)
@@ -28,8 +28,8 @@
! nferr = nf_open('topo/Greenland_5km_v1.1.nc', NF_SHARE, rd_ncid)
-! nferr = nf_open('topo/gis_5km.180511.nc', NF_SHARE, rd_ncid)
- nferr = nf_open('topo/gis_20km.180511.nc', NF_SHARE, rd_ncid)
+ nferr = nf_open('topo/gis_5km.180211.nc', NF_SHARE, rd_ncid)
+! nferr = nf_open('topo/gis_20km.180511.nc', NF_SHARE, rd_ncid)
write(6,*) ' nferr ncid ', nferr, rd_ncid
!
Modified: branches/land_ice_projects/grid_tools/icesheet_legacy/src/module_write_netcdf.F
===================================================================
--- branches/land_ice_projects/grid_tools/icesheet_legacy/src/module_write_netcdf.F 2012-03-13 21:41:54 UTC (rev 1627)
+++ branches/land_ice_projects/grid_tools/icesheet_legacy/src/module_write_netcdf.F 2012-03-13 21:46:10 UTC (rev 1628)
@@ -88,7 +88,7 @@
integer, intent(in) :: nVertLevels
integer, intent(in) :: vertexDegree
character (len=16) :: on_a_sphere
- double precision :: sphere_radius
+ real (kind = 8) :: sphere_radius
integer :: nferr
Modified: branches/land_ice_projects/grid_tools/icesheet_legacy/src/test_steady_state_greenland.F
===================================================================
--- branches/land_ice_projects/grid_tools/icesheet_legacy/src/test_steady_state_greenland.F 2012-03-13 21:41:54 UTC (rev 1627)
+++ branches/land_ice_projects/grid_tools/icesheet_legacy/src/test_steady_state_greenland.F 2012-03-13 21:46:10 UTC (rev 1628)
@@ -1,16 +1,13 @@
-
program test_steady_state_greenland
- use grid_types
- use configure
- use io_input
- use dmpar
- use timer
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_io_input
+ use mpas_dmpar
+ use mpas_timer
use read_topo
-
implicit none
-
type (dm_info), pointer :: dminfo
type (domain_type), pointer :: domain
@@ -24,29 +21,26 @@
integer, dimension(:), pointer :: recvEdgesArray
logical, dimension(:), pointer :: verticesMask
- real, dimension(:), pointer :: thckdata, elevdata, betadata, temperaturedata, u, v;
+ real, dimension(:), pointer :: thckdata, lowerSurfdata, betadata, temperaturedata, u, v, levelsRatio;
integer nVerticesSolve, nEdgesSolve, nCellsSolve, nVertices, nEdges, nCells, i, k, iCell, nLevels, nKeptVerteces
logical keep_vertex
-
-
-
allocate(dminfo)
allocate(list)
- call dmpar_init(dminfo)
+ call mpas_dmpar_init(dminfo)
! initialize lifev data structure
- call read_namelist(dminfo)
+ call mpas_read_namelist(dminfo)
- call allocate_domain(domain, dminfo)
+ call mpas_allocate_domain(domain, dminfo)
- call input_state_for_domain(domain)
+ call mpas_input_state_for_domain(domain)
!extract data from domain
nCellsSolve = domain % blocklist % mesh % nCellsSolve
@@ -56,46 +50,48 @@
nCells = domain % blocklist % mesh % nCells
nEdges = domain % blocklist % mesh % nEdges
nVertices = domain % blocklist % mesh % nVertices
- nLevels = domain % blocklist % mesh % nVertLevels
+ nLevels = 5;! domain % blocklist % mesh % nVertLevels
!get elevation and thickness on each cell
- call get_topo_fields(domain, thckdata, elevdata)
+ call get_topo_fields(domain, thckdata, lowerSurfdata)
! set beta and temperature
allocate(betadata(nCells))
allocate(temperaturedata(nCells*(nLevels+1)))
+ allocate(levelsRatio(nLevels))
+
+ do i = 1,nLevels
+ levelsRatio(i) = 1./real(nLevels);
+ end do
+
do i = 1,nCells
- betadata(i) = 1000
+ betadata(i) = 10000
do k = 1, nLevels+1
temperaturedata((i-1)*(nLevels+1)+k) = 270.963
end do
end do
-
-
!select vertices that should be kept depending on the thickness
!(discard a vertex if all the cells around it have zero thickness)
allocate(verticesMask(nVertices));
nKeptVerteces = 0;
do i = 1,nVertices
- keep_vertex = .false.
- do k = 1, 3
- iCell = domain % blocklist % mesh % cellsOnVertex % array(k,i)
- keep_vertex = keep_vertex.or.(thckdata(iCell) > 0.01)
- end do
- if (keep_vertex .and. (i <= nVerticesSolve )) then
- nKeptVerteces = nKeptVerteces + 1
- endif
- verticesMask(i) = keep_vertex
- end do
+ keep_vertex = .false.
+ do k = 1, 3
+ iCell = domain % blocklist % mesh % cellsOnVertex % array(k,i)
+ keep_vertex = keep_vertex.or.(thckdata(iCell) > 0.01)
+ end do
+ if (keep_vertex .and. (i <= nVerticesSolve )) then
+ nKeptVerteces = nKeptVerteces + 1
+ endif
+ verticesMask(i) = keep_vertex
+ end do
-write(6,*) 'nKeptVerteces', nKeptVerteces
-
- call init_lifev(dminfo % comm, nKeptVerteces)
+ call velocity_solver_init_mpi(dminfo % comm)
!build send and receive arrays using exchange_list
@@ -108,30 +104,55 @@
- !call lifeV and set the grid of the velocity solver
+ !call lifeV and set the grid of the velocity solver
+ call velocity_solver_set_grid_data(nCells, nEdges, nVertices, nLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, &
+ domain % blocklist % mesh % cellsOnEdge % array, domain % blocklist % mesh % cellsOnVertex % array, &
+ domain % blocklist % mesh % verticesOnCell % array, domain % blocklist % mesh % verticesOnEdge % array, &
+ domain % blocklist % mesh % xCell % array, domain % blocklist % mesh % yCell % array, domain % blocklist % mesh % zCell % array, &
+ sendCellsArray, recvCellsArray, sendEdgesArray, recvEdgesArray, sendVerticesArray, recvVerticesArray)
- call set_velocity_solver_grid(nCells, nEdges, nVertices, nCellsSolve, nEdgesSolve, nVerticesSolve, &
- domain % blocklist % mesh % cellsOnEdge % array, domain % blocklist % mesh % cellsOnVertex % array, &
- domain % blocklist % mesh % verticesOnCell % array, domain % blocklist % mesh % verticesOnEdge % array, &
- verticesMask, domain % blocklist % mesh % xCell % array, domain % blocklist % mesh % yCell % array, domain % blocklist % mesh % zCell % array, &
- thckdata, elevdata, &
- sendCellsArray, recvCellsArray, sendEdgesArray, recvEdgesArray, sendVerticesArray, recvVerticesArray)
+ call velocity_solver_compute_2d_grid(verticesMask);
+
+ !call velocity_solver_extrude_3d_grid(levelsRatio, lowerSurfdata, thckdata)
- ! call solve_lifev()
+
- call initialize_L1L2_solver(nLevels)
+ call velocity_solver_init_L1L2(levelsRatio)
+ !call velocity_solver_init_FO(levelsRatio)
+
allocate(u(nEdges*(nLevels+1)))
allocate(v(nEdges*(nLevels+1)))
- call velocity_solver_L1L2(elevdata, thckdata, betadata, temperaturedata, u, v)
+ !velocities u and v are the velocity projections along x and y axes
+ call velocity_solver_solve_L1L2(lowerSurfdata, thckdata, betadata, temperaturedata, u, v)
+ call velocity_solver_export_L1L2_velocity();
+ call velocity_solver_export_2d_data(lowerSurfdata, thckdata, betadata)
+
+ !call velocity_solver_solve_FO(lowerSurfdata, thckdata, betadata, temperaturedata, u, v)
+ !call velocity_solver_export_FO_velocity();
+
+ call velocity_solver_finalize();
+
+ deallocate(sendVerticesArray, recvVerticesArray, sendCellsArray, recvCellsArray, sendEdgesArray, recvEdgesArray)
+ deallocate(betadata, temperaturedata, verticesMask, u, v, thckdata, lowerSurfdata, levelsRatio);
+
+
+
- call deallocate_domain(domain)
+ call mpas_deallocate_domain(domain)
+ call mpas_dmpar_finalize(dminfo)
- call dmpar_finalize(dminfo)
+
+
+
+ deallocate(list)
+ deallocate(dminfo)
contains
+!----------------------------------------------------------------------
+
subroutine array_from_exchange_list(exlist, array)
implicit none
@@ -167,13 +188,15 @@
end subroutine array_from_exchange_list
- subroutine get_topo_fields(domain, thckdata, elevdata)
+!----------------------------------------------------------------------
+ subroutine get_topo_fields(domain, thckdata, lowerSurfdata)
+
implicit none
type (domain_type), pointer, intent (in) :: domain
- real, dimension(:), pointer, intent (out) :: thckdata, elevdata
+ real, dimension(:), pointer, intent (out) :: thckdata, lowerSurfdata
real(kind=4), allocatable, dimension(:) :: x, y
@@ -200,7 +223,7 @@
allocate(thckdata(nCells))
- allocate(elevdata(nCells))
+ allocate(lowerSurfdata(nCells))
!interpolate data from the structured mesh .nc to the current mesh
@@ -224,19 +247,15 @@
alpha_y = domain % blocklist % mesh % yCell % array(iCell) / delta_y + 1 - iy
thckdata(iCell) = (1.-alpha_x)*(1.-alpha_y)*thck(ix,iy) + alpha_x*(1.-alpha_y)*thck(ix_plus_1,iy) &
- + (1.- alpha_x)*alpha_y*thck(ix,iy_plus_1) + alpha_x*alpha_y*thck(ix_plus_1,iy_plus_1)
+ + (1.- alpha_x)*alpha_y*thck(ix,iy_plus_1) + alpha_x*alpha_y*thck(ix_plus_1,iy_plus_1) + 0.01
- elevdata(iCell) = thckdata(iCell) + (1.-alpha_x)*(1.-alpha_y)*topg(ix,iy) + alpha_x*(1.-alpha_y)*topg(ix_plus_1,iy) &
+ lowerSurfdata(iCell) = (1.-alpha_x)*(1.-alpha_y)*topg(ix,iy) + alpha_x*(1.-alpha_y)*topg(ix_plus_1,iy) &
+ (1.- alpha_x)*alpha_y*topg(ix,iy_plus_1) + alpha_x*alpha_y*topg(ix_plus_1, iy_plus_1)
enddo
+ deallocate(topg, x ,y)
+
end subroutine get_topo_fields
-
-
end program
-
-
-
-
</font>
</pre>