<p><b>mperego@fsu.edu</b> 2012-01-09 16:54:44 -0700 (Mon, 09 Jan 2012)</p><p></p><hr noshade><pre><font color="gray">Added: branches/land_ice/mpas/src/core_land_ice/test_steady_state_greenland.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/test_steady_state_greenland.F (rev 0)
+++ branches/land_ice/mpas/src/core_land_ice/test_steady_state_greenland.F 2012-01-09 23:54:44 UTC (rev 1333)
@@ -0,0 +1,235 @@
+
+
+program test_steady_state_greenland
+ use grid_types
+ use configure
+ use io_input
+ use dmpar
+ use timer
+ use read_topo
+
+
+ implicit none
+
+
+ type (dm_info), pointer :: dminfo
+ type (domain_type), pointer :: domain
+ type (exchange_list), pointer :: list
+
+ integer, dimension(:), pointer :: sendCellsArray
+ integer, dimension(:), pointer :: recvCellsArray
+ integer, dimension(:), pointer :: sendVerticesArray
+ integer, dimension(:), pointer :: recvVerticesArray
+ integer, dimension(:), pointer :: sendEdgesArray
+ integer, dimension(:), pointer :: recvEdgesArray
+ logical, dimension(:), pointer :: verticesMask
+
+ real, dimension(:), pointer :: thckdata, elevdata, betadata, temperaturedata, u, v;
+
+ integer nVerticesSolve, nEdgesSolve, nCellsSolve, nVertices, nEdges, nCells, i, k, iCell, nLevels
+ logical keep_vertex
+
+
+
+
+ allocate(dminfo)
+ allocate(list)
+
+
+ call dmpar_init(dminfo)
+
+ ! initialize lifev data structure
+ call init_lifev(dminfo % comm)
+
+ call read_namelist(dminfo)
+
+ call allocate_domain(domain, dminfo)
+
+
+ call input_state_for_domain(domain)
+
+ !extract data from domain
+ nCellsSolve = domain % blocklist % mesh % nCellsSolve
+ nEdgesSolve = domain % blocklist % mesh % nEdgesSolve
+ nVerticesSolve = domain % blocklist % mesh % nVerticesSolve
+
+ nCells = domain % blocklist % mesh % nCells
+ nEdges = domain % blocklist % mesh % nEdges
+ nVertices = domain % blocklist % mesh % nVertices
+ nLevels = 10 !domain % blocklist % mesh % nVertLevels
+
+ !get elevation and thickness on each cell
+
+ call get_topo_fields(domain, thckdata, elevdata)
+
+ ! set beta and temperature
+ allocate(betadata(nCells))
+ allocate(temperaturedata(nCells*(nLevels+1)))
+
+ do i = 1,nCells
+ betadata(i) = 1000
+ 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));
+
+ 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
+ verticesMask(i) = keep_vertex
+ end do
+
+
+ !build send and receive arrays using exchange_list
+
+ call array_from_exchange_list(domain % blocklist % parinfo % verticesToSend, sendVerticesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % verticesToRecv, recvVerticesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % cellsToSend, sendCellsArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % cellsToRecv, recvCellsArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % edgesToSend, sendEdgesArray)
+ call array_from_exchange_list(domain % blocklist % parinfo % edgesToRecv, recvEdgesArray)
+
+
+
+ !call lifeV and set the grid of the velocity solver
+
+ 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 solve_lifev()
+
+ call initialize_L1L2_solver(nLevels)
+
+ allocate(u(nEdges*(nLevels+1)))
+ allocate(v(nEdges*(nLevels+1)))
+
+ call velocity_solver_L1L2(elevdata, thckdata, betadata, temperaturedata, u, v)
+
+ call deallocate_domain(domain)
+
+ call dmpar_finalize(dminfo)
+
+ contains
+
+ subroutine array_from_exchange_list(exlist, array)
+
+ implicit none
+
+ type (exchange_list), pointer, intent(in) :: exlist
+ type (exchange_list), pointer :: listPtr
+
+ integer, dimension(:), pointer, intent(out) :: array
+ integer dataSize, offset, i
+
+ dataSize = 1 !in first position we will store the size of the vector
+ listPtr => exlist
+ do while (associated(listPtr))
+ dataSize = dataSize + (listPtr % nlist) + 2
+ listPtr => listPtr % next
+ end do
+
+ allocate(array(dataSize))
+
+ array(1) = dataSize;
+ offset = 2;
+ listPtr => exlist
+ do while (associated(listPtr))
+ array(offset) = listPtr % procID
+ offset = offset + 1
+ array(offset) = listPtr % nlist
+ do i=1,listPtr % nlist
+ array(i+offset) = listPtr % list(i) -1
+ end do
+ offset = offset + listPtr % nlist + 1
+ listPtr => listPtr % next
+ end do
+
+ end subroutine array_from_exchange_list
+
+ subroutine get_topo_fields(domain, thckdata, elevdata)
+
+ implicit none
+
+ type (domain_type), pointer, intent (in) :: domain
+
+ real, dimension(:), pointer, intent (out) :: thckdata, elevdata
+
+ real(kind=4), allocatable, dimension(:) :: x, y
+
+ real(kind=4), allocatable, dimension(:,:) :: thck, topg
+
+ real delta_x, delta_y, alpha_x, alpha_y
+
+ integer nCells, ix, iy, nx, ny, iCell, ix_plus_1, iy_plus_1
+
+ call read_topo_init( nx, ny)
+
+ allocate(thck(nx,ny))
+ thck = 0.0
+
+ allocate(topg(nx,ny))
+ topg = 0.0
+
+ allocate(x(nx))
+ allocate(y(ny))
+
+ call read_topo_fields(x,y,thck,topg)
+
+ nCells = domain % blocklist % mesh % nCells
+
+ allocate(thckdata(nCells))
+
+ allocate(elevdata(nCells))
+
+ !interpolate data from the structured mesh .nc to the current mesh
+
+ delta_x = (x(nx) -x(1))/(nx - 1);
+ delta_y = (y(ny) -y(1))/(ny - 1);
+
+ write(6,*) 'nx', nx, ny
+
+ do iCell=1,nCells
+
+ ix = floor( domain % blocklist % mesh % xCell % array(iCell) / delta_x )
+ iy = floor( domain % blocklist % mesh % yCell % array(iCell) / delta_y )
+
+ ix = min(nx,ix); iy=min(ny,iy)
+ ix = max(1,ix); iy=max(1,iy)
+
+ ix_plus_1 = min(ix+1, nx);
+ iy_plus_1 = min(iy+1, ny);
+
+ alpha_x = domain % blocklist % mesh % xCell % array(iCell) / delta_x + 1 - ix
+ 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)
+
+ elevdata(iCell) = thckdata(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
+
+ end subroutine get_topo_fields
+
+
+
+
+end program
+
+
+
+
</font>
</pre>