<p><b>mhoffman@lanl.gov</b> 2013-02-07 09:52:42 -0700 (Thu, 07 Feb 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
Adding a first-order accurate vertical remapping scheme for thickness and tracers, based on the one in CISM. For the tests I have done using both the FO tracer advection scheme in my last commit and Doug's FCT tracer advection, energy is conserved on 1 and 4 processors with advancing boundaries and an initial temperature field that varies both horizontally and vertically. I have not tested with SMB/BMB.<br>
<br>
With this commit tracer advection can be considered functional. Note that vertical diffusion still needs to be added to get full temperature evolution.<br>
<br>
I have adjusted Registry so that FO tracer advection is the default (I think further testing is needed to understand why FCT tracer advection works for advancing cells before encouraging its use in the land ice core), and made temperature an output variable (it had been turned off for the ice2sea runs).<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 2013-02-07 16:06:05 UTC (rev 2453)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2013-02-07 16:52:42 UTC (rev 2454)
@@ -22,8 +22,8 @@
% valid time integration is ForwardEuler
namelist character land_ice_model config_thickness_advection FO-Upwind
% valid advection is FO-Upwind, None (FCT coming soon)
-namelist character land_ice_model config_tracer_advection None
-% valid options for tracer_advection are FCT (flux-corrected transport), None
+namelist character land_ice_model config_tracer_advection FO
+% valid options for tracer_advection are FCT (flux-corrected transport), FO (first-order upwind), None
namelist real land_ice_model config_ice_density 900.0
namelist real land_ice_model config_ocean_density 1000.0
namelist real land_ice_model config_sealevel 0.0
@@ -210,14 +210,14 @@
var persistent real thickness ( nCells Time ) 2 iro thickness state - -
var persistent real bedTopography ( nCells Time ) 2 ir bedTopography state - -
% bedTopography is really a mesh variable now, but would be a state variable if isostasy is included.
-var persistent real temperature ( nVertLevels nCells Time ) 2 ir temperature state tracers dynamics
+var persistent real temperature ( nVertLevels nCells Time ) 2 iro temperature state tracers dynamics
var persistent real normalVelocity ( nVertLevels nEdges Time ) 2 ir normalVelocity state - -
-%var persistent real sup_thickness ( nVertLevels nCells Time ) 2 iro sup_thickness state sup_thickness bob
+%var persistent real sup_layerThickness ( nVertLevels nCells Time ) 2 iro sup_layerThickness state sup_layerThickness bob
% Tendency variables
var persistent real tend_layerThickness ( nVertLevels nCells Time ) 1 - layerThickness tend - -
var persistent real tend_temperature ( nVertLevels nCells Time ) 1 - temperature tend tracers dynamics
-%var persistent real tend_sup_thickness ( nVertLevels nCells Time ) 1 o sup_thickness tend sup_thickness bob
+%var persistent real tend_sup_layerThickness ( nVertLevels nCells Time ) 1 o sup_layerThickness tend sup_layerThickness bob
% Diagnostic fields: only written to output
var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
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 2013-02-07 16:06:05 UTC (rev 2453)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2013-02-07 16:52:42 UTC (rev 2454)
@@ -439,16 +439,9 @@
! as always, upper surface is the lower surface plus the thickness
upperSurface(:) = lowerSurface(:) + thickness(:)
- ! set up layerThickness from thickness using the layerThicknessFractions (vertical remapping)
- do iCell = 1, nCells
- thisThk = sum(layerThickness(:,iCell))
- do iLevel = 1, nVertLevels
- layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel) * thisThk
- end do
- ! Check for conservation of mass. Put any residual in the top layer.
- layerThickness(1,iCell) = layerThickness(1,iCell) + (thisThk - sum(layerThickness(:,iCell)) )
- end do
+ call vertical_remap(layerThickness, thickness, state % tracers % array, mesh, err)
+
end subroutine land_ice_diagnostic_solve
@@ -1143,6 +1136,157 @@
end subroutine tracer_advection_tend_fo
+
+!***********************************************************************
+!
+! subroutine vertical_remap
+!
+!> \brief Vertical remapping of thickness and tracers
+!> \author Matt Hoffman
+!> \date 06 February 2013
+!> \version SVN:$Id$
+!> \details
+!> This routine performs vertical remapping of thickness and tracers from one vertical
+!> coordinate system to another, as is required for our sigma coordinate system.
+!> The remapping is first-order accurate.
+!> This uses code from the CISM glissade_transport.F90 module written by Bill Lipscomb.
+!> I have altered the array structures to work with MPAS. Indexing/looping order is a bit
+!> of a hodgepodge at the moment and should be optimized.
+!
+!-----------------------------------------------------------------------
+ subroutine vertical_remap(layerThickness, thickness, tracers, mesh, err)
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:), intent(in) :: &
+ thickness !< Input:
+
+ type (mesh_type), intent(in) :: &
+ mesh !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ layerThickness !< Input:
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tracers !< Input:
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ ! pointers to mesh arrays
+ real (kind=RKIND), dimension(:), pointer :: layerThicknessFractions, layerInterfaceSigma
+ ! local arrays
+ real (kind=RKIND), dimension(:), allocatable :: recipThickness
+ real (kind=RKIND), dimension(:,:), allocatable :: layerInterfaceSigma_Input
+ real (kind=RKIND), dimension(:,:,:), allocatable :: hTsum
+ ! counters, mesh variables, index variables
+ integer :: nTracers, nCells, nVertLevels
+ integer :: iCell, k, k1, k2, nt
+ ! stuff for making calculations
+ real(kind=RKIND) :: thisThk, zhi, zlo, hOverlap
+
+ nCells = mesh % nCells
+ nVertLevels = mesh % nVertLevels
+ nTracers = size(tracers, 1)
+
+ layerThicknessFractions => mesh % layerThicknessFractions % array
+ layerInterfaceSigma => mesh % layerInterfaceSigma % array
+
+ allocate(recipThickness(nCells+1))
+ allocate(layerInterfaceSigma_Input(nVertLevels+1, nCells+1))
+ allocate(hTsum(nCells+1, nTracers, nVertLevels))
+
+ ! *** Calculate reciprocal thickness so we don't divide by 0
+ where (thickness > 0.0_RKIND)
+ recipThickness = 1.0_RKIND / thickness
+ elsewhere
+ recipThickness = 0.0_RKIND
+ end where
+
+ ! *** Calculate vertical sigma coordinates of each layer interface for the input non-sigma state and desired new sigma-based state (we already have that as mesh % layerInterfaceSigma)
+ layerInterfaceSigma_Input(1,:) = 0.0_RKIND
+ do k = 2, nVertLevels
+ layerInterfaceSigma_Input(k,:) = layerInterfaceSigma_Input(k-1,:) + layerThickness(k-1,:) * recipThickness(:)
+ end do
+ layerInterfaceSigma_Input(nVertLevels+1,:) = 1.0_RKIND
+
+ ! *** Compute new layer thicknesses (layerInterfaceSigma coordinates)
+ do iCell = 1, nCells
+ thisThk = sum(layerThickness(:,iCell))
+ do k = 1, nVertLevels
+ layerThickness(k,iCell) = layerThicknessFractions(k) * thisThk
+ end do
+ ! Check for conservation of mass. Put any residual in the top layer.
+ layerThickness(1,iCell) = layerThickness(1,iCell) + (thisThk - sum(layerThickness(:,iCell)) )
+ end do
+ ! TODO This conservation check may make layerThicknesses inconsistent with the sigma levels (which are used below in tracer remapping.
+
+ !-----------------------------------------------------------------
+ ! Compute sum of h*T for each new layer (k2) by integrating
+ ! over the regions of overlap with old layers (k1).
+ ! Note: It might be worth trying a more efficient
+ ! search algorithm if the number of layers is large.
+ ! This algorithm scales as nlyr^2.
+ ! Also, may want to rearrange loop order if there are many tracers.
+ !-----------------------------------------------------------------
+
+ do k2 = 1, nVertLevels
+ hTsum(:,:,k2) = 0.d0
+ do k1 = 1, nVertLevels
+ do nt = 1, nTracers
+ do iCell = 1, nCells
+ zhi = min (layerInterfaceSigma_Input(k1+1,iCell), layerInterfaceSigma(k2+1))
+ zlo = max (layerInterfaceSigma_Input(k1,iCell), layerInterfaceSigma(k2))
+ hOverlap = max (zhi-zlo, 0.d0) * thickness(iCell)
+ hTsum(iCell,nt,k2) = htsum(iCell,nt,k2) &
+ + tracers(nt,k1,iCell) * hOverlap
+ enddo ! iCell
+ enddo ! nt
+ enddo ! k1
+ enddo ! k2
+
+ !-----------------------------------------------------------------
+ ! Compute tracer values in new layers
+ !-----------------------------------------------------------------
+
+ do k = 1, nVertLevels
+ do nt = 1, nTracers
+ do iCell = 1, nCells
+ if (layerThickness(k, iCell) > 0.0_RKIND) then
+ tracers(nt,k,iCell) = hTsum(iCell,nt,k) / layerThickness(k, iCell)
+ else
+ tracers(nt,k,iCell) = 0.0_RKIND
+ endif
+ enddo ! iCell
+ enddo ! nt
+ enddo ! k
+
+ deallocate(recipThickness)
+ deallocate(layerInterfaceSigma_Input)
+ deallocate(hTsum)
+
+ end subroutine vertical_remap
+
end module land_ice_tendency
</font>
</pre>