<p><b>mhoffman@lanl.gov</b> 2012-05-07 16:33:23 -0600 (Mon, 07 May 2012)</p><p>BRANCH COMMIT - land ice<br>
<br>
Added a thickness advection config option called 'None' that holds thickness constant (will be useful for spinning up temperature once the core includes temperature calculations). It required a little reorganization.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-05-07 21:41:10 UTC (rev 1876)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-05-07 22:33:23 UTC (rev 1877)
@@ -15,8 +15,8 @@
% valid dycores are SIA, L1L2, FO, Stokes. All but SIA require compiling with LifeV libraries.
namelist character land_ice_model config_time_integration ForwardEuler
% valid time integration is ForwardEuler
-% %namelist character land_ice_model config_advection FO-Upwind
-%% valid time integration is FO-Upwind - this could be added once other advection schemes are added.
+namelist character land_ice_model config_advection FO-Upwind
+% valid advection is FO-Upwind, None (currently only applies to thickness)
namelist real land_ice_model config_ice_density 900.0
% The following options are used by the ocean core and may be useful in the future
namelist integer land_ice_model config_stats_interval 100
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-05-07 21:41:10 UTC (rev 1876)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-05-07 22:33:23 UTC (rev 1877)
@@ -148,7 +148,7 @@
use mpas_rbf_interpolation
use mpas_vector_reconstruction
use land_ice_vel
- use land_ice_tendency, only: lice_diagnostic_solve
+ use land_ice_tendency, only: land_ice_diagnostic_solve
implicit none
@@ -180,7 +180,7 @@
! Initialize state ====
call mpas_timer_start("initial state calculation")
- call lice_diagnostic_solve(mesh, state, err)
+ call land_ice_diagnostic_solve(mesh, state, err)
! Set the initial flag to 1 so LifeV will calculate its mesh information the first time
state % anyVertexMaskChanged % scalar = 1
call mpas_timer_stop("initial state calculation")
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-05-07 21:41:10 UTC (rev 1876)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-05-07 22:33:23 UTC (rev 1877)
@@ -35,9 +35,8 @@
! Public member functions
!
!--------------------------------------------------------------------
- public :: lice_tend_h, &
- lice_tend_h_layers, &
- lice_diagnostic_solve
+ public :: land_ice_tendency_h, &
+ land_ice_diagnostic_solve
!--------------------------------------------------------------------
!
@@ -50,11 +49,199 @@
contains
+!***********************************************************************
+!
+! subroutine land_ice_tendency_h
+!
+!> \brief Computes tendency term from horizontal advection of thickness
+!> \author Matt Hoffman
+!> \date 16 April 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal advection tendency for
+!> thickness based on current state and user choices of forcings. Based on
+!> ocn_thick_hadv_tend in the ocean core.
+!
+!-----------------------------------------------------------------------
+ subroutine land_ice_tendency_h(mesh, state, thickness_tend, dt, dminfo, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: grid information
+
+ type (state_type), intent(in) :: &
+ state !< Input: state to use to calculate tendency
+
+ real (kind=RKIND), intent(in) :: &
+ dt !< Input: dt
+
+ type (dm_info), pointer, intent(in) :: &
+ dminfo !< Input: domain info
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:), pointer, intent(inout) :: &
+ thickness_tend !< Input/Output: thickness tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+
+ select case (config_advection)
+ case ('FO-Upwind')
+ call land_ice_tend_h_fo_upwind(mesh, state % normalVelocity % array, &
+ state % layerThickness % array, state % thickness % array, thickness_tend, dt, err)
+ ! Commented lines below calculate the thickness tendency per layer instead of for the whole column (last 2 lines should be moved to following 'section' of code)
+ !call lice_tend_h_layers(mesh, normalVelocity, layerThickness, layerThickness_tend, dt, err)
+ !layerThickness = layerThickness + layerThickness_tend * dt / SecondsInYear
+ !thickness = sum(layerThickness, 1)
+
+ case ('None')
+ thickness_tend = 0.0
+ case default
+ write(0,*) trim(config_advection), ' is not a valid thickness advection option.'
+ call mpas_dmpar_abort(dminfo)
+ end select
+
+ end subroutine land_ice_tendency_h
+
+
+
!***********************************************************************
!
-! subroutine lice_tend_h
+! subroutine land_ice_diagnostic_solve
!
+!> \brief Computes diagnostic variables
+!> \author Matt Hoffman
+!> \date 19 April 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the diagnostic variables for land ice
+!
+!-----------------------------------------------------------------------
+ subroutine land_ice_diagnostic_solve(mesh, state, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(inout) :: &
+ state !< Input/Output: state for which to update diagnostic variables
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+ integer :: iCell, iLevel, nCells, nVertLevels, nVertices
+ real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
+ lowerSurface, bedTopography, layerThicknessFractions
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+ integer, dimension(:), pointer :: vertexMask
+ integer, dimension(:,:), pointer :: cellsOnVertex
+ integer :: keep_vertex, i, k
+
+
+ nCells = mesh % nCells
+ nVertices = mesh % nVertices
+ nVertLevels = mesh % nVertLevels
+ layerThicknessFractions => mesh % layerThicknessFractions % array
+ cellsOnVertex => mesh % cellsOnVertex % array
+
+ thickness => state % thickness % array
+ upperSurface => state % upperSurface % array
+ lowerSurface => state % lowerSurface % array
+ bedTopography => state % bedTopography % array
+ layerThickness => state % layerThickness % array
+ vertexMask => state % vertexMask % array
+
+
+ ! no ice shelves -- lower surface same as bed topography
+ ! Eventually this should consider floatation!
+ 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
+
+
+ !\todo Calculate h_edge that can be used by SIA solve and FO advection scheme. ocn_diagnostic_solve in mpas_ocn_tendency.F has 2, 3, & 4th order calculations for h_edge that can be used. Right now I am using 2nd order h_edge in SIA solve and 1st order h_edge in FO advection.
+
+ !\todo Calculate masks. Right now sia driver and Forward Euler time integration have need of a mask and both do it on their own. We need to define how the masks should work and implement them here.
+ ! simple vertex mask needed by LifeV
+ do i = 1,nVertices
+ keep_vertex = 0
+ do k = 1, 3
+ iCell = cellsOnVertex(k,i)
+ if (thickness(iCell) .gt. 0.01) then
+ keep_vertex = 1
+ endif
+ end do
+ vertexMask(i) = keep_vertex
+ end do
+
+
+
+ end subroutine land_ice_diagnostic_solve
+
+
+
+
+ ! private subroutines================================================
+
+
+
+
+!***********************************************************************
+!
+! subroutine land_ice_tend_h_fo_upwind
+!
!> \brief Computes tendency term from horizontal advection of thickness
!> \author Matt Hoffman
!> \date 16 April 2012
@@ -66,7 +253,7 @@
!
!-----------------------------------------------------------------------
- subroutine lice_tend_h(grid, normalVelocity, layerThickness, thickness, thickness_tend, dt, err)!{{{
+ subroutine land_ice_tend_h_fo_upwind(grid, normalVelocity, layerThickness, thickness, thickness_tend, dt, err)!{{{
!-----------------------------------------------------------------
!
@@ -160,7 +347,7 @@
end do
!--------------------------------------------------------------------
- end subroutine lice_tend_h!}}}
+ end subroutine land_ice_tend_h_fo_upwind !}}}
@@ -268,108 +455,7 @@
-!***********************************************************************
-!
-! subroutine lice_diagnostic_solve
-!
-!> \brief Computes diagnostic variables
-!> \author Matt Hoffman
-!> \date 19 April 2012
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the diagnostic variables for land ice
-!
-!-----------------------------------------------------------------------
- subroutine lice_diagnostic_solve(mesh, state, err)!{{{
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
- type (mesh_type), intent(in) :: &
- mesh !< Input: grid information
-
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
-
- type (state_type), intent(inout) :: &
- state !< Input/Output: state for which to update diagnostic variables
-
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
- integer :: iCell, iLevel, nCells, nVertLevels, nVertices
- real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
- lowerSurface, bedTopography, layerThicknessFractions
- real (kind=RKIND), dimension(:,:), pointer :: layerThickness
- integer, dimension(:), pointer :: vertexMask
- integer, dimension(:,:), pointer :: cellsOnVertex
- integer :: keep_vertex, i, k
-
-
- nCells = mesh % nCells
- nVertices = mesh % nVertices
- nVertLevels = mesh % nVertLevels
- layerThicknessFractions => mesh % layerThicknessFractions % array
- cellsOnVertex => mesh % cellsOnVertex % array
-
- thickness => state % thickness % array
- upperSurface => state % upperSurface % array
- lowerSurface => state % lowerSurface % array
- bedTopography => state % bedTopography % array
- layerThickness => state % layerThickness % array
- vertexMask => state % vertexMask % array
-
-
- ! no ice shelves -- lower surface same as bed topography
- ! Eventually this should consider floatation!
- 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
-
-
- !\todo Calculate h_edge that can be used by SIA solve and FO advection scheme. ocn_diagnostic_solve in mpas_ocn_tendency.F has 2, 3, & 4th order calculations for h_edge that can be used. Right now I am using 2nd order h_edge in SIA solve and 1st order h_edge in FO advection.
-
- !\todo Calculate masks. Right now sia driver and Forward Euler time integration have need of a mask and both do it on their own. We need to define how the masks should work and implement them here.
- ! simple vertex mask needed by LifeV
- do i = 1,nVertices
- keep_vertex = 0
- do k = 1, 3
- iCell = cellsOnVertex(k,i)
- if (thickness(iCell) .gt. 0.01) then
- keep_vertex = 1
- endif
- end do
- vertexMask(i) = keep_vertex
- end do
-
-
-
- end subroutine lice_diagnostic_solve
-
-
end module land_ice_tendency
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-05-07 21:41:10 UTC (rev 1876)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-05-07 22:33:23 UTC (rev 1877)
@@ -52,11 +52,11 @@
case ('ForwardEuler')
call land_ice_time_integrator_forwardeuler(domain, dt)
case ('RK4')
- write(*,*) trim(config_time_integration), ' is not currently supported.'
- !call mpas_dmpar_abort(dminfo)
+ write(0,*) trim(config_time_integration), ' is not currently supported.'
+ call mpas_dmpar_abort(domain % dminfo)
case default
- write(*,*) trim(config_time_integration), ' is not a valid land ice time integration option.'
- !call mpas_dmpar_abort(dminfo)
+ write(0,*) trim(config_time_integration), ' is not a valid land ice time integration option.'
+ call mpas_dmpar_abort(domain % dminfo)
end select
block => domain % blocklist
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-05-07 21:41:10 UTC (rev 1876)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-05-07 22:33:23 UTC (rev 1877)
@@ -104,14 +104,9 @@
! === Calculate Tendencies ========================
call mpas_timer_start("calculate tendencies")
! Calculate thickness tendency using state at time n
- call lice_tend_h(mesh, normalVelocityOld, layerThicknessOld, thicknessOld, thickness_tend, dt, err)
+ call land_ice_tendency_h(mesh, stateOld, thickness_tend, dt, dminfo, err)
- ! Commented lines below calculate the thickness tendency per layer instead of for the whole column (last 2 lines should be moved to following 'section' of code)
- !call lice_tend_h_layers(mesh, normalVelocity, layerThickness, layerThickness_tend, dt, err)
- !layerThickness = layerThickness + layerThickness_tend * dt / SecondsInYear
- !thickness = sum(layerThickness, 1)
-
- ! update halos on thickness tend
+ ! update halos on thickness tend (perhaps move to land_ice_tendency_h after we've updated the halo update calls)
call mpas_dmpar_exch_halo_field1d_real(dminfo, thickness_tend, &
mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
@@ -150,7 +145,7 @@
! === Calculate diagnostic variables for new state =====================
call mpas_timer_start("calc. diagnostic vars except vel")
- call lice_diagnostic_solve(mesh, stateNew, err) ! perhaps velocity solve should move in here.
+ call land_ice_diagnostic_solve(mesh, stateNew, err) ! perhaps velocity solve should move in here.
! Determine if the vertex mask changed during this time step for this block (needed for LifeV)
if ( sum(stateNew % vertexMask % array - stateOld % vertexMask % array) .ne. 0 ) then
</font>
</pre>