<p><b>mhoffman@lanl.gov</b> 2012-05-17 11:02:17 -0600 (Thu, 17 May 2012)</p><p> BRANCH COMMIT -- land ice<br>
Major update to allow tracer advection using the Flux Corrected Transport routines in the ocean core. Currently only applies to temperature, but additional tracers should not require modification to this code. No vertical remapping of temperature is currently done. Not yet working on multiple processors.<br>
Added registry options to enable and adjust the FCT for tracers.<br>
Also added commented code to apply the FCT for thickness advection (for a future commit after a bit more development/testing).<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Makefile
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-05-16 20:31:21 UTC (rev 1916)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-05-17 17:02:17 UTC (rev 1917)
@@ -15,7 +15,17 @@
mpas_land_ice_vel.o \
mpas_land_ice_lifev.o \
mpas_land_ice_sia.o \
- mpas_land_ice_global_diagnostics.o
+ mpas_land_ice_global_diagnostics.o \
+ mpas_ocn_tracer_advection.o \
+ mpas_ocn_tracer_advection_std.o \
+ mpas_ocn_tracer_advection_std_hadv.o \
+ mpas_ocn_tracer_advection_std_vadv.o \
+ mpas_ocn_tracer_advection_std_vadv2.o \
+ mpas_ocn_tracer_advection_std_vadv3.o \
+ mpas_ocn_tracer_advection_std_vadv4.o \
+ mpas_ocn_tracer_advection_mono.o \
+ mpas_ocn_tracer_advection_helpers.o \
+ mpas_ocn_advection.o
all: core_land_ice
@@ -24,7 +34,7 @@
mpas_land_ice_time_integration.o: mpas_land_ice_time_integration_forwardeuler.o
-mpas_land_ice_time_integration_forwardeuler.o: mpas_land_ice_vel.o mpas_land_ice_tendency.o
+mpas_land_ice_time_integration_forwardeuler.o: mpas_land_ice_vel.o mpas_land_ice_tendency.o mpas_ocn_tracer_advection.o
mpas_land_ice_mask.o:
@@ -38,8 +48,32 @@
mpas_land_ice_global_diagnostics.o:
-mpas_land_ice_mpas_core.o: mpas_land_ice_global_diagnostics.o mpas_land_ice_time_integration.o
+mpas_land_ice_mpas_core.o: mpas_land_ice_global_diagnostics.o mpas_land_ice_time_integration.o mpas_ocn_advection.o
+
+# modules needed for FCT tracer advection from ocean core (eventually to be in operators?)
+mpas_ocn_tracer_advection.o: mpas_ocn_tracer_advection_std.o mpas_ocn_tracer_advection_mono.o
+
+mpas_ocn_tracer_advection_std.o: mpas_ocn_tracer_advection_std_hadv.o mpas_ocn_tracer_advection_std_vadv.o
+
+mpas_ocn_tracer_advection_std_hadv.o: mpas_ocn_tracer_advection_helpers.o
+
+mpas_ocn_tracer_advection_std_vadv.o: mpas_ocn_tracer_advection_std_vadv2.o mpas_ocn_tracer_advection_std_vadv3.o mpas_ocn_tracer_advection_std_vadv4.o
+
+mpas_ocn_tracer_advection_std_vadv2.o: mpas_ocn_tracer_advection_helpers.o
+
+mpas_ocn_tracer_advection_std_vadv3.o: mpas_ocn_tracer_advection_helpers.o
+
+mpas_ocn_tracer_advection_std_vadv4.o: mpas_ocn_tracer_advection_helpers.o
+
+mpas_ocn_tracer_advection_mono.o: mpas_ocn_tracer_advection_helpers.o
+
+mpas_ocn_tracer_advection_helpers.o:
+
+mpas_ocn_advection.o:
+
+
+
clean:
$(RM) *.o *.mod *.f90 libdycore.a
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-16 20:31:21 UTC (rev 1916)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-05-17 17:02:17 UTC (rev 1917)
@@ -19,7 +19,9 @@
namelist character land_ice_model config_time_integration ForwardEuler
% valid time integration is ForwardEuler
namelist character land_ice_model config_advection FO-Upwind
-% valid advection is FO-Upwind, None (currently only applies to thickness)
+% valid advection is FO-Upwind, None (FCT coming soon)
+namelist character land_ice_model config_tracer_advection FCT
+% valid options for tracer_advection are FCT (flux-corrected transport), 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
@@ -33,6 +35,15 @@
%namelist logical land_ice_model config_positive_definite false
%namelist logical land_ice_model config_monotonic false
+% Options needed for FCT tracer advection (from ocean core)
+namelist integer advection config_vert_tracer_adv_order 2
+% Currently, vert tracer advection is not implemented, so this values does not matter but needs to be declared
+namelist integer advection config_horiz_tracer_adv_order 2
+namelist real advection config_coef_3rd_order 0.25
+% coef_3rd_order only applies if previous config value is 3. 0.25 is supposed to be the optimal value.
+namelist logical advection config_monotonic true
+namelist logical advection config_check_monotonicity true
+
% I/O options:
namelist character io config_input_name land_ice_grid.nc
namelist character io config_output_name output.nc
@@ -58,6 +69,7 @@
dim nEdges nEdges
dim maxEdges maxEdges
dim maxEdges2 maxEdges2
+dim nAdvectionCells maxEdges2+0
dim nVertices nVertices
dim TWO 2
dim R3 3
@@ -132,17 +144,25 @@
%var persistent real meshScalingDel4 ( nEdges ) 0 ro meshScalingDel4 mesh - -
% Mesh variables needed for FCT advection
-% Space needed for advection
-% This requires calling something from mpas_ocn_advection module in block_init to calculate deriv_two
-%var persistent real deriv_two ( maxEdges2 TWO nEdges ) 0 - deriv_two mesh - -
-% Added for monotonic advection scheme
-%var persistent real adv_coefs ( nAdvectionCells nEdges ) 0 - adv_coefs mesh - -
-%var persistent real adv_coefs_2nd ( nAdvectionCells nEdges ) 0 - adv_coefs_2nd mesh - -
-%var persistent real adv_coefs_3rd ( nAdvectionCells nEdges ) 0 - adv_coefs_3rd mesh - -
-%var persistent integer advCellsForEdge ( nAdvectionCells nEdges ) 0 - advCellsForEdge mesh - -
-%var persistent integer nAdvCellsForEdge ( nEdges ) 0 - nAdvCellsForEdge mesh - -
-%var persistent integer highOrderAdvectionMask ( nVertLevels nEdges ) 0 - highOrderAdvectionMask mesh - -
-%var persistent integer lowOrderAdvectionMask ( nVertLevels nEdges ) 0 - lowOrderAdvectionMask mesh - -
+var persistent real deriv_two ( maxEdges2 TWO nEdges ) 0 - deriv_two mesh - -
+var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
+var persistent real adv_coefs ( nAdvectionCells nEdges ) 0 - adv_coefs mesh - -
+var persistent real adv_coefs_2nd ( nAdvectionCells nEdges ) 0 - adv_coefs_2nd mesh - -
+var persistent real adv_coefs_3rd ( nAdvectionCells nEdges ) 0 - adv_coefs_3rd mesh - -
+var persistent integer advCellsForEdge ( nAdvectionCells nEdges ) 0 - advCellsForEdge mesh - -
+var persistent integer nAdvCellsForEdge ( nEdges ) 0 - nAdvCellsForEdge mesh - -
+var persistent integer highOrderAdvectionMask ( nVertLevels nEdges ) 0 o highOrderAdvectionMask mesh - -
+var persistent integer lowOrderAdvectionMask ( nVertLevels nEdges ) 0 o lowOrderAdvectionMask mesh - -
+var persistent integer maxLevelCell ( nCells ) 0 - maxLevelCell mesh - -
+var persistent integer maxLevelEdgeTop ( nEdges ) 0 - maxLevelEdgeTop mesh - -
+var persistent integer maxLevelEdgeBot ( nEdges ) 0 - maxLevelEdgeBot mesh - -
+% !! NOTE: the following arrays are needed to allow the use
+% !! of the module_advection.F w/o alteration
+% Only needed to calculate deriv2, I think, so that calc could be moved out to a different module and then these variables would not be needed?
+% Space needed for deformation calculation weights
+var persistent real defc_a ( maxEdges nCells ) 0 - defc_a mesh - -
+var persistent real defc_b ( maxEdges nCells ) 0 - defc_b mesh - -
+var persistent real kdiff ( nVertLevels nCells Time ) 0 - kdiff mesh - -
% Arrays required for reconstruction of velocity field
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
@@ -150,8 +170,8 @@
% Boundary conditions: read from input, saved in restart and written to output (these are used by the ocean core but may be useful in this core eventually)
%var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
%var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
-%var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
+
% Mesh variables specific to land ice core
% Sigma coordinate: read from input, saved in restart and written to output
var persistent real layerThicknessFractions ( nVertLevels ) 0 iro layerThicknessFractions mesh - -
@@ -175,14 +195,14 @@
% 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 iro temperature state tracers dynamics
var persistent real normalVelocity ( nVertLevels nEdges Time ) 2 iro normalVelocity state - -
+%var persistent real sup_thickness ( nVertLevels nCells Time ) 2 iro sup_thickness state sup_thickness bob
-
% Tendency variables
-var persistent real tend_layerThickness ( nVertLevels nCells Time ) 1 - layerThickness tend - -
+var persistent real tend_layerThickness ( nVertLevels nCells Time ) 1 o layerThickness tend - -
var persistent real tend_thickness ( nCells Time ) 1 o thickness tend - -
-var persistent real tend_temperature ( nVertLevels nCells Time ) 1 - temperature tend tracers dynamics
+var persistent real tend_temperature ( nVertLevels nCells Time ) 1 o temperature tend tracers dynamics
+%var persistent real tend_sup_thickness ( nVertLevels nCells Time ) 1 o sup_thickness tend sup_thickness bob
-
% Diagnostic fields: only written to output
var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
var persistent real uReconstructY ( nVertLevels nCells Time ) 2 o uReconstructY state - -
@@ -190,15 +210,15 @@
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
var persistent real layerThickness ( nVertLevels nCells Time ) 2 o layerThickness state - -
+var persistent real layerThicknessEdge ( nVertLevels nEdges Time ) 2 o layerThicknessEdge state - -
var persistent real lowerSurface ( nCells Time ) 2 o lowerSurface state - -
var persistent real upperSurface ( nCells Time ) 2 o upperSurface state - -
var persistent real temperatureUpperBC ( nCells Time ) 2 o temperatureUpperBC state - -
var persistent real temperatureLowerBC ( nCells Time ) 2 o temperatureLowerBC state - -
-
% Other diagnostic variables: neither read nor written to any files
var persistent real layerThicknessVertex ( nVertLevels nVertices Time ) 2 - layerThicknessVertex state - -
-var persistent real layerThicknessEdge ( nVertLevels nEdges Time ) 2 - layerThicknessEdge state - -
+
% Masks: calculated diagnostically at each time step
var persistent integer cellMask ( nCells Time ) 2 o cellMask state - -
var persistent integer edgeMask ( nEdges Time ) 2 o edgeMask state - -
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-16 20:31:21 UTC (rev 1916)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -22,6 +22,7 @@
use mpas_grid_types
use mpas_constants
use land_ice_vel
+ use mpas_ocn_tracer_advection
implicit none
@@ -30,7 +31,7 @@
type (block_type), pointer :: block
- integer :: i, err
+ integer :: i, err, err_tmp
err = 0
@@ -53,6 +54,9 @@
call land_ice_vel_init(domain, err)
+ call mpas_ocn_tracer_advection_init(err_tmp)
+ err = ior(err,err_tmp)
+
block => domain % blocklist
do while (associated(block))
! Initialize blocks
@@ -149,6 +153,8 @@
use mpas_vector_reconstruction
use land_ice_vel
use land_ice_tendency, only: land_ice_diagnostic_solve
+ use mpas_ocn_tracer_advection
+ use ocn_advection
implicit none
@@ -168,6 +174,14 @@
! Call init routines ====
call land_ice_setup_vertical_coords(mesh)
+ ! Init for FCT tracer advection
+ mesh % maxLevelCell % array = mesh % nVertLevels ! Needed for FCT tracer advection
+ mesh % maxLevelEdgeTop % array = mesh % nVertLevels ! Needed for FCT tracer advection
+ mesh % maxLevelEdgeBot % array = mesh % nVertLevels ! Needed for FCT tracer advection
+ call ocn_initialize_advection_rk(mesh, err)
+ call mpas_ocn_tracer_advection_coefficients(mesh, err_tmp)
+ err = ior(err, err_tmp)
+
!call compute_mesh_scaling(mesh) ! may be useful for variable resolution meshes
call mpas_rbf_interp_initialize(mesh)
call mpas_init_reconstruct(mesh)
@@ -191,7 +205,7 @@
call land_ice_vel_solve(mesh, state, err)
call mpas_timer_stop("initial velocity calculation")
-
+
if(err == 1) then
print *, "An error has occurred in mpas_init_block. Aborting..."
call mpas_dmpar_abort(dminfo)
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-16 20:31:21 UTC (rev 1916)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -20,6 +20,7 @@
use mpas_constants
use mpas_dmpar
use land_ice_mask
+ use mpas_ocn_tracer_advection
implicit none
private
@@ -37,6 +38,7 @@
!
!--------------------------------------------------------------------
public :: land_ice_tendency_h, &
+ land_ice_tendency_tracers, &
land_ice_diagnostic_solve
!--------------------------------------------------------------------
@@ -65,7 +67,7 @@
!
!-----------------------------------------------------------------------
- subroutine land_ice_tendency_h(mesh, state, thickness_tend, dt, dminfo, err)
+ subroutine land_ice_tendency_h(mesh, state, layerThickness_tend, dt, dminfo, err)
!-----------------------------------------------------------------
!
@@ -91,8 +93,8 @@
!
!-----------------------------------------------------------------
- real (kind=RKIND), dimension(:), pointer, intent(inout) :: &
- thickness_tend !< Input/Output: thickness tendency
+ real (kind=RKIND), dimension(:,:), pointer, intent(inout) :: &
+ layerThickness_tend !< Input/Output: layer thickness tendency
!-----------------------------------------------------------------
!
@@ -108,33 +110,183 @@
!
!-----------------------------------------------------------------
+ ! 0 tendency
+ layerThickness_tend = 0.0
select case (config_advection)
- case ('FO-Upwind')
+ case ('FO-Upwind') !===================================================
!print *,'Using FO Upwind for thickness advection'
- 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
+ ! Alternate call to calculate thickness tendency for entire ice column (deprecated)
+ !call land_ice_tend_h_fo_upwind(mesh, state % normalVelocity % array, &
+ ! state % layerThickness % array, state % thickness % array, thickness_tend, dt, err)
+
+ call land_ice_tend_h_layers_fo_upwind(mesh, state % normalVelocity % array, &
+ state % layerThicknessEdge % array, layerThickness_tend, dt, err)
+
+ !case ('FCT') !===================================================
+ ! MJH TEMP TRYING IT WITH THICKNESS
+ ! stateOld % sup_thickness % array(1,:,:) = stateOld % layerThickness % array !!!!! MJH TEMP
+ ! block % tend % sup_thickness % array = 0.0
+ ! call mpas_ocn_tracer_advection_tend(stateOld % sup_thickness % array, uh , &
+ ! wTop, stateOld % sup_thickness % array(1,:,:), stateOld % sup_thickness % array(1,:,:), &
+ ! dt/SecondsInYear , mesh, 0.0*layerThickness_tend, block % tend % sup_thickness % array)
+ ! where (stateOld % sup_thickness % array .gt. 1.0e-9)
+ ! block % tend % sup_thickness % array = block % tend % sup_thickness % array / stateOld % sup_thickness % array
+ ! else where
+ ! block % tend % sup_thickness % array = 0.0
+ ! end where
+ ! assign to the thickness tend i actually use
+ ! layerThickness_tend = block % tend % sup_thickness % array(1,:,:)
+
+ case ('None') !===================================================
+ ! Do nothing
+ case default !===================================================
write(0,*) trim(config_advection), ' is not a valid thickness advection option.'
call mpas_dmpar_abort(dminfo)
- end select
+ end select !===================================================
+
! Add surface mass balance to tendency
- thickness_tend = thickness_tend + mesh % sfcMassBal % array * dt / SecondsInYear
+ ! NOTE: Need to decide where to put this for layer by layer thickness tendency, and how to deal with negative SMB that eliminates top layer or all ice (check for negative thickness?)
+ !thickness_tend = thickness_tend + mesh % sfcMassBal % array * dt / SecondsInYear
end subroutine land_ice_tendency_h
+
+
!***********************************************************************
!
+! subroutine land_ice_tendency_tracers
+!
+!> \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_tracers(mesh, state, layerThickness_tend, tracer_tendency, 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), dimension(:,:), pointer, intent(in) :: &
+ layerThickness_tend !< Input/Output: layer thickness 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) :: &
+ tracer_tendency !< Input/Output: tracers tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+ real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThickness, layerThicknessEdge
+ real (kind=RKIND), dimension(:,:), allocatable :: wTop, uh
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer :: iEdge, k
+
+ ! 0 tendency
+ tracer_tendency = 0.0
+
+ select case (config_tracer_advection)
+ case ('FCT') !===================================================
+ ! Assign pointers, etc.
+ layerThickness => state % layerThickness % array
+ layerThicknessEdge => state % layerThicknessEdge % array
+ normalVelocity => state % normalVelocity % array
+ tracers => state % tracers % array
+ allocate(wTop(mesh % nVertLevels + 1, mesh % nCells + 1))
+ allocate(uh(mesh % nVertLevels, mesh % nEdges + 1))
+
+ ! Setup 0 vertical velocity field - vertical velocity not needed with sigma coordinates
+ wTop = 0.0
+
+ ! prepare for FCT call
+
+ ! This is a kluge because the FCT code will divide by zero if thickness is zero \todo: deal with more elegantly!
+ where (layerThickness == 0.0)
+ layerThickness = 1.0e-12
+ end where
+
+ ! FCT code needs u*h on edges - this could be calculated elsewhere and saved, potentially (e.g. it's already calculated locally in the FO upwind thickness routine)
+ do iEdge = 1, mesh % nEdges
+ do k = 1, mesh % nVertLevels
+ uh(k, iEdge) = normalVelocity(k, iEdge) * layerThicknessEdge(k, iEdge)
+ end do
+ end do
+
+ !print *,'uh', maxval(uh), minval(uh)
+ !print *,'wTop', maxval(wTop), minval(wTop)
+ !print *,'layerThicknessOld', maxval(layerThicknessOld), minval(layerThicknessOld)
+ !print *,'layerThickness_tend', maxval(layerThickness_tend), minval(layerThickness_tend)
+ !print *,' temp min/max:', minval(stateOld % tracers % array(stateOld%index_temperature,:,1:mesh%nCells)), maxval(stateOld % tracers % array(stateOld%index_temperature,:,1:mesh%nCells))
+
+ ! Call the FCT! (this will likely move from the ocean core to operators eventually)
+ call mpas_ocn_tracer_advection_tend(tracers, uh , &
+ wTop, layerThickness, layerThickness, dt/SecondsInYear, &
+ mesh, layerThickness_tend, tracer_tendency)
+ !print *,'tracer_tendency', maxval(tracer_tendency), minval(tracer_tendency)
+
+ ! Set thickness back to 0 where needed. This is a kluge because the FCT code will divide by zero if thickness is zero \todo: deal with more elegantly!
+ where (layerThickness == 1.0e-12)
+ layerThickness = 0.0
+ end where
+
+
+ deallocate(wTop)
+ deallocate(uh)
+
+ case ('None') !===================================================
+ ! Do nothing
+ case default !===================================================
+ write(0,*) trim(config_tracer_advection), ' is not a valid tracer advection option.'
+ call mpas_dmpar_abort(dminfo)
+ end select !===================================================
+
+ end subroutine land_ice_tendency_tracers
+
+
+
+
+!***********************************************************************
+!
! subroutine land_ice_diagnostic_solve
!
!> \brief Computes diagnostic variables
@@ -181,17 +333,19 @@
integer :: iCell, iLevel, nCells, nVertLevels, nVertices
real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
lowerSurface, bedTopography, layerThicknessFractions
- real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge
integer, dimension(:), pointer :: vertexMask
- integer, dimension(:,:), pointer :: cellsOnVertex
- integer :: keep_vertex, i, k
+ integer, dimension(:,:), pointer :: cellsOnVertex, cellsOnEdge
+ integer :: keep_vertex, i, k, iEdge, nEdges, cell1, cell2
nCells = mesh % nCells
+ nEdges = mesh % nEdges
nVertices = mesh % nVertices
nVertLevels = mesh % nVertLevels
layerThicknessFractions => mesh % layerThicknessFractions % array
cellsOnVertex => mesh % cellsOnVertex % array
+ cellsOnEdge => mesh % cellsOnEdge % array
thickness => state % thickness % array
upperSurface => state % upperSurface % array
@@ -199,6 +353,7 @@
bedTopography => state % bedTopography % array
layerThickness => state % layerThickness % array
vertexMask => state % vertexMask % array
+ layerThicknessEdge => state % layerThicknessEdge % array
! no ice shelves -- lower surface same as bed topography
@@ -208,7 +363,7 @@
! as always, upper surface is the lower surface plus the thickness
upperSurface(:) = lowerSurface(:) + thickness(:)
- ! set up layerThickness from thickness using the layerThicknessFractions
+ ! set up layerThickness from thickness using the layerThicknessFractions (vertical remapping)
do iCell = 1, nCells
do iLevel = 1, nVertLevels
layerThickness(iLevel,iCell) = layerThicknessFractions(iLevel)*thickness(iCell)
@@ -219,7 +374,23 @@
! \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.
+ ! NOTE: This calculates FO upwind h edge
+ ! Both thickness and layerThickness should be updated by this time.
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1, nVertLevels
+ ! Calculate h on edges using first order
+ layerThicknessEdge(k,iEdge) = max(layerThickness(k, cell1), layerThickness(k, cell2))
+ end do
+ ! thickness_edge is not currently in registry and not currenly needed. If it is, uncomment the next line
+ !h_edge = max(thickness(cell1), thickness(cell2))
+ !!!h_edge = (thickness(k) + thickness(k) ) / 2.0 ! 2nd order
+ end do
+
+
+
! Calculate masks
call land_ice_calculate_mask(mesh, state, err)
@@ -352,7 +523,7 @@
!***********************************************************************
!
-! subroutine lice_tend_h_layers
+! subroutine land_ice_tend_h_layers_fo_upwind
!
!> \brief Computes tendency term from horizontal advection of thickness layers
!> \author Matt Hoffman
@@ -360,14 +531,14 @@
!> \version SVN:$Id$
!> \details
!> This routine computes the horizontal advection tendency for each
-!> thickness layer based on current state and user choices of forcings. Based on
+!> thickness layer using first-order upwinding. Based on
!> ocn_thick_hadv_tend in the ocean core. This is an alternative to lice_tend_h
-!> that caclculates the tendency for each layer, which would then need to be
+!> that calculates the tendency for each layer, which would then need to be
!> added up to calculate the change in thickness. The two methods yield identical
-!> results. Currently, this is not used but I'm leaving it here in case it is useful later.
+!> results.
!
!-----------------------------------------------------------------------
- subroutine lice_tend_h_layers(grid, normalVelocity, layerThickness, tend, dt, err)!{{{
+ subroutine land_ice_tend_h_layers_fo_upwind(grid, normalVelocity, layerThicknessEdge, tend, dt, err)!{{{
!-----------------------------------------------------------------
!
@@ -379,7 +550,7 @@
normalVelocity !< Input: velocity
real (kind=RKIND), dimension(:,:), intent(in) :: &
- layerThickness !< Input: thickness of each layer at edge
+ layerThicknessEdge !< Input: thickness of each layer at edge
type (mesh_type), intent(in) :: &
grid !< Input: grid information
@@ -413,7 +584,7 @@
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND) :: flux, VelSign, h_edge
+ real (kind=RKIND) :: flux
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, dcEdge
err = 0
@@ -426,22 +597,14 @@
dcEdge => grid % dcEdge % array
areaCell => grid % areaCell % array
- ! Zero the tendency before accumulating
- tend = 0.0
-
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1, nVertLevels
if ( (abs(normalVelocity(k, iEdge)) * dt/SecondsInYear) .gt. (0.5 * dcEdge(iEdge))) then
print *, 'CFL violation at level, edge', k, iEdge
- endif
- ! Calculate h on edges using first order
- VelSign = sign(1.0, normalVelocity(k, iEdge))
- h_edge = max(0.0, VelSign) * layerThickness(k, cell1) &
- - min(0.0, VelSign) * layerThickness(k, cell2)
- !h_edge = (layerThickness(k, cell1) +layerThickness(k, cell2) ) / 2.0
- flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * h_edge
+ endif
+ flux = normalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge)
tend(k,cell1) = tend(k,cell1) - flux / areaCell(cell1)
tend(k,cell2) = tend(k,cell2) + flux / areaCell(cell2)
end do
@@ -449,12 +612,8 @@
!--------------------------------------------------------------------
- end subroutine lice_tend_h_layers!}}}
+ end subroutine land_ice_tend_h_layers_fo_upwind!}}}
-
-
-
-
end module land_ice_tendency
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-16 20:31:21 UTC (rev 1916)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -46,18 +46,19 @@
type (mesh_type), pointer :: mesh
real (kind=RKIND), dimension(:), pointer :: thicknessOld, thicknessNew, thickness_tend, layerThicknessFractions
- real (kind=RKIND), dimension(:,:), pointer :: normalVelocityOld, normalVelocityNew, layerThicknessOld, layerThicknessNew !, layerThickness_tend
- real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ, &
- uReconstructZonal, uReconstructMeridional
+ real (kind=RKIND), dimension(:,:), pointer :: normalVelocityOld, normalVelocityNew, layerThicknessOld, layerThicknessNew, layerThickness_tend
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracer_tendency, tracersNew, tracersOld
integer, dimension(:), allocatable :: mask
integer :: blockVertexMaskChanged, procVertexMaskChanged, anyVertexMaskChanged
- integer :: nVertLevels, k
+ integer :: nVertLevels, k, iEdge, iTracer
integer :: err
+
+
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
! (time level 1 should not be modified.)
@@ -78,53 +79,66 @@
thicknessOld => stateOld % thickness % array
layerThicknessOld => stateOld % layerThickness % array
normalVelocityOld => stateOld % normalVelocity % array
+ tracersOld => stateOld % tracers % array
! State at time n+1 (advanced by dt by Forward Euler)
stateNew => block % state % time_levs(2) % state
thicknessNew => stateNew % thickness % array
layerThicknessNew => stateNew % layerThickness % array
normalVelocityNew => stateNew % normalVelocity % array
- uReconstructX => stateNew % uReconstructX % array
- uReconstructY => stateNew % uReconstructY % array
- uReconstructZ => stateNew % uReconstructZ % array
- uReconstructZonal => stateNew % uReconstructZonal % array
- uReconstructMeridional => stateNew % uReconstructMeridional % array
+ tracersNew => stateNew % tracers % array
! Tendencies
- !layerThickness_tend => block % tend % layerThickness % array
+ layerThickness_tend => block % tend % layerThickness % array
thickness_tend => block % tend % thickness % array
+ tracer_tendency => block % tend % tracers % array
allocate(mask(mesh%nCells + 1))
+
! === Implicit column physics (vertical temperature diffusion) ===========
!call ()
! === Calculate Tendencies ========================
call mpas_timer_start("calculate tendencies")
- ! Calculate thickness tendency using state at time n
- call land_ice_tendency_h(mesh, stateOld, thickness_tend, dt, dminfo, err)
+ ! Calculate thickness tendency using state at time n =========
+ call land_ice_tendency_h(mesh, stateOld, layerThickness_tend, dt, dminfo, err)
! 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, &
+ call mpas_dmpar_exch_halo_field2d_real(dminfo, layerthickness_tend, &
+ mesh % nVertLevels, mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
!print *,' Max abs thickness tend:', maxval(abs(thickness_tend(1:mesh % nCellsSolve)))
- ! (Calculate tracer tendencies)
- !call ()
+
+ ! Calculate tracer tendencies ==========
+ call land_ice_tendency_tracers(mesh, stateOld, layerThickness_tend, tracer_tendency, dt, dminfo, err)
+
+ ! update halos on tracer tend (perhaps move to land_ice_tendency_tracers after we've updated the halo update calls)
+ call mpas_dmpar_exch_halo_field3d_real(dminfo, tracer_tendency, &
+ size(tracer_tendency,dim=1), mesh % nVertLevels, mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
call mpas_timer_stop("calculate tendencies")
+
-
! === Compute new state for prognostic variables ==================================
! (once implicit column physics are added (i.e. temp diffusion), these calculations will need to be adjusted to apply to the new values as needed)
call mpas_timer_start("calc. new prognostic vars")
- ! Update thickness (and tracers eventually)
- thicknessNew = thicknessOld + thickness_tend * dt / SecondsInYear
+ ! Update thickness ======================
+
+ ! Commented out usage for advecting thickness as a column
+ !!!thicknessNew = thicknessOld + thickness_tend * dt / SecondsInYear
+ ! Commented out usage for using FCT for thickness
+ !!!stateNew % sup_thickness % array(1,:,:) = (stateOld % tracers % array(stateOld % index_temperature,:,:) * layerThicknessOld + tracer_tendency(stateOld % index_temperature, :, :) * dt / SecondsInYear) / (layerThicknessNew+1.0e-12)
- !Optionally print some information about the new state
+ layerThicknessNew = layerThicknessOld + layerThickness_tend * dt / SecondsInYear
+ thicknessNew = sum(layerThicknessNew, 1)
+
+ !Optionally print some information about the new thickness
!print *, 'thickness_tend maxval:', maxval(thickness_tend(1:mesh % nCellsSolve))
!print *, 'thicknessOld maxval:', maxval(thicknessOld(1:mesh % nCellsSolve))
print *, ' thicknessNew maxval:', maxval(thicknessNew(1:mesh % nCellsSolve))
@@ -143,8 +157,21 @@
end where
print *,' Cells with nonzero thickness:', sum(mask)
+
+ ! Calculate new tracer values =================
+ do iTracer = 1, size(tracersNew, 1)
+ where (layerThicknessNew .gt. 0.0)
+ tracersNew(iTracer,:,:) = (tracersOld(iTracer,:,:) * layerThicknessOld &
+ + tracer_tendency(iTracer,:,:) * dt / SecondsInYear) / (layerThicknessNew)
+ elsewhere
+ ! May or may not want to assign tracer values to non-ice cells
+ tracersNew(iTracer,:,:) = 0.0
+ end where
+ end do
+
call mpas_timer_stop("calc. new prognostic vars")
+
! === Calculate diagnostic variables for new state =====================
call mpas_timer_start("calc. diagnostic vars except vel")
call land_ice_diagnostic_solve(mesh, stateNew, err) ! perhaps velocity solve should move in here.
@@ -174,6 +201,9 @@
call mpas_timer_stop("calc. diagnostic vars except vel")
+
+ deallocate(mask)
+
block => block % next
end do
@@ -207,8 +237,6 @@
! ================================
- deallocate(mask)
-
end subroutine land_ice_time_integrator_ForwardEuler
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_advection.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_advection.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_advection.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_advection.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_advection.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_helpers.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_helpers.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_helpers.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection_helpers.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_helpers.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_mono.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_mono.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection_mono.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_mono.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection_std.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_hadv.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_hadv.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_hadv.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection_std_hadv.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_hadv.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection_std_vadv.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv2.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv2.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv2.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection_std_vadv2.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv2.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv3.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv3.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv3.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection_std_vadv3.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv3.F
___________________________________________________________________
Added: svn:special
+ *
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv4.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv4.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv4.F 2012-05-17 17:02:17 UTC (rev 1917)
@@ -0,0 +1 @@
+link ../core_ocean/mpas_ocn_tracer_advection_std_vadv4.F
\ No newline at end of file
Property changes on: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_ocn_tracer_advection_std_vadv4.F
___________________________________________________________________
Added: svn:special
+ *
</font>
</pre>