<p><b>xylar@lanl.gov</b> 2012-01-20 17:14:16 -0700 (Fri, 20 Jan 2012)</p><p>BRANCH COMMIT<br>
<br>
Added missing variable initializations to land ice test case 1<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_test_cases.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_test_cases.F        2012-01-20 23:54:35 UTC (rev 1406)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_test_cases.F        2012-01-21 00:14:16 UTC (rev 1407)
@@ -64,14 +64,14 @@
end subroutine setup_land_ice_test_case
- subroutine land_ice_test_case_1(grid, state)
+ subroutine land_ice_test_case_1(mesh, state)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup land ice test case 1: Isothermal dome
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
- type (mesh_type), intent(inout) :: grid
+ type (mesh_type), intent(inout) :: mesh
type (state_type), intent(inout) :: state
real (kind=RKIND), parameter :: r0 = 60000.0*sqrt(0.125) !meters
@@ -79,24 +79,47 @@
real (kind=RKIND), parameter :: x0 = 30000.0 !meters
real (kind=RKIND), parameter :: y0 = 30000.0 !meters
- integer :: iCell
+ integer :: iCell, iLevel, nCells, nVertLevels, index_temperature
real (kind=RKIND) :: r
- real (kind=RKIND), dimension(:), pointer :: x, y, thickness
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, thickness, upperSurface, &
+ lowerSurface, bedTopography, layerThicknessFractions
+ real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThickness
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- if(grid % on_a_sphere) then
- print *, "Land ice test case 1 only supports planar grids."
+
+ if(mesh % on_a_sphere) then
+ print *, "Land ice test case 1 only supports planar meshes."
stop
end if
+ nCells = mesh % nCells
+ nVertLevels = mesh % nVertLevels
- x => grid % xCell % array
- y => grid % yCell % array
+ xCell => mesh % xCell % array
+ yCell => mesh % yCell % array
+ layerThicknessFractions => mesh % layerThicknessFractions % array
thickness => state % thickness % array
+ upperSurface => state % upperSurface % array
+ lowerSurface => state % lowerSurface % array
+ bedTopography => state % bedTopography % array
+ normalVelocity => state % normalVelocity % array
+ layerThickness => state % layerThickness % array
+ tracers => state % tracers % array
+ index_temperature = state % index_temperature
- do iCell = 1,grid % nCells
- r = sqrt((x(iCell)-x0)**2 + (y(iCell)-y0)**2)
+ ! zero velocity everywhere
+ normalVelocity(:,:) = 0.0
+
+ ! flat bed at sea level
+ bedTopography(:) = 0.0
+
+ ! constant, arbitrary temperature
+ tracers(index_temperature,:,:) = -20.0 ! degrees C
+
+ do iCell = 1, nCells
+ r = sqrt((xCell(iCell)-x0)**2 + (yCell(iCell)-y0)**2)
if(r < r0) then
thickness(iCell) = h0*sqrt(1 - (r/r0)**2)
else
@@ -104,7 +127,20 @@
end if
end do
+ ! no ice shelves -- lower surface same as bed topography
+ lowerSurface(:) = bedTopography(:)
+ ! as always, upper surface is the lower surface plus the thickness
+ upperSurface(:) = lowerSurface(:) + thickness(:)
+
+ ! set up layerThickness from thickness using the layerThicknessFractions
+ do iCell = 1, nCells
+ do iLevel = 1, nVertLevels
+ layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel)*thickness(iCell)
+ end do
+ end do
+
+
end subroutine land_ice_test_case_1
!!! subroutine sw_test_case_1(grid, state)
</font>
</pre>