<p><b>mhoffman@lanl.gov</b> 2012-04-19 15:35:38 -0600 (Thu, 19 Apr 2012)</p><p>BRANCH COMMIT<br>
Created a diagnostic solve subroutine that is called from the time stepper to calculate e.g. upperSurface. <br>
<br>
The velocity solve could be called from that subroutine as well, but I've left it separate since it is a very special diagnostic variable. However should we introduce diagnostic variables that are derived from velocity, we'll probably want to move it in there (though the init call will need to be different if a velocity field was supplied with input).<br>
<br>
I also reverted the velocity solve calls to not take 'normalVelocity' as an additional argument but simply take the current state. normalVelocity is a part of that state and we want it consistent with the state passed in, so there is no need for a separate input field for normalVelocity.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F 2012-04-18 22:42:58 UTC (rev 1796)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_lifev.F 2012-04-19 21:35:38 UTC (rev 1797)
@@ -289,7 +289,7 @@
!***********************************************************************
- subroutine land_ice_lifev_solve(mesh, state, normalVelocity, err)
+ subroutine land_ice_lifev_solve(mesh, state, err)
!-----------------------------------------------------------------
!
@@ -300,19 +300,14 @@
type (mesh_type), intent(in) :: &
mesh !< Input: mesh information
- type (state_type), intent(in) :: &
- state !< Input: state information
- ! Note: This state is meant to be unmodified (input only).
-
!-----------------------------------------------------------------
!
! input/output variables
!
!-----------------------------------------------------------------
- real (kind=RKIND), intent(inout), dimension(:,:) :: &
- normalVelocity !< Input/Output: the result of the velocity solve
- ! Note: this velocity should be used for both input or output - it may be different than the velocity included in state, and the velocity included in state should not be touched. If land_ice_vel_solve() returns anything else besides normalVelocity, those variables should also be treated in this way. This allow the time integration scheme to decide how to treat multiple time levels.
+ type (state_type), intent(inout) :: &
+ state !< Input: state information
!-----------------------------------------------------------------
!
@@ -329,6 +324,8 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:), pointer :: &
thickness, lowerSurface, beta
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ normalVelocity
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers
@@ -336,6 +333,7 @@
err = 0
+ normalVelocity => state % normalVelocity % array
thickness => state % thickness % array
lowerSurface => state % lowerSurface % array
beta => state % beta % array
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-04-18 22:42:58 UTC (rev 1796)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-04-19 21:35:38 UTC (rev 1797)
@@ -55,12 +55,6 @@
! Assign initial time stamp
block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
- ! Calculate an initial (diagnostic) velocity field consistent with the initial thickness field
- ! \todo: skip this step if a velocity field was supplied in the I.C. input file
- call land_ice_vel_solve(block % mesh, block % state % time_levs(1) % state, &
- block % state % time_levs(1) % state % normalVelocity % array, err)
- ! (halo update of velocity field)
-
! Make sure all time levels have a copy of the initial state
do i=2,nTimeLevs
call mpas_copy_state(block % state % time_levs(i) % state, block % state % time_levs(1) % state)
@@ -144,8 +138,8 @@
use mpas_grid_types
use mpas_rbf_interpolation
use mpas_vector_reconstruction
- use land_ice_time_integration
- use land_ice_vel, only: land_ice_vel_block_init
+ use land_ice_vel
+ use land_ice_tendency, only: lice_diagnostic_solve
implicit none
@@ -156,57 +150,30 @@
type (dm_info) :: dminfo
integer :: err, err_tmp
- integer :: iCell, iLevel, nCells, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
- lowerSurface, bedTopography, layerThicknessFractions
- real (kind=RKIND), dimension(:,:), pointer :: layerThickness
-
+ type (state_type), pointer :: state
+
err = 0
- nCells = mesh % nCells
- nVertLevels = mesh % nVertLevels
- layerThicknessFractions => mesh % layerThicknessFractions % array
- thickness => block % state % time_levs(1) % state % thickness % array
- upperSurface => block % state % time_levs(1) % state % upperSurface % array
- lowerSurface => block % state % time_levs(1) % state % lowerSurface % array
- bedTopography => block % state % time_levs(1) % state % bedTopography % array
- layerThickness => block % state % time_levs(1) % state % layerThickness % array
-
+ state => block % state % time_levs(1) % state ! initial state
-
+ ! Call init routines ====
call compute_mesh_scaling(mesh)
call mpas_rbf_interp_initialize(mesh)
call mpas_init_reconstruct(mesh)
- call mpas_reconstruct(mesh, block % state % time_levs(1) % state % normalVelocity % array, &
- block % state % time_levs(1) % state % uReconstructX % array, &
- block % state % time_levs(1) % state % uReconstructY % array, &
- block % state % time_levs(1) % state % uReconstructZ % array, &
- block % state % time_levs(1) % state % uReconstructZonal % array, &
- block % state % time_levs(1) % state % uReconstructMeridional % array &
- )
-
+ call land_ice_vel_block_init(mesh, err_tmp)
+ err = ior(err, err_tmp)
- ! Initialize state ==== eventually this should be moved to its own subroutine/module
- ! no ice shelves -- lower surface same as bed topography
- ! Eventually this should consider floatation!
- lowerSurface(:) = bedTopography(:)
+ ! Initialize state ====
+ call lice_diagnostic_solve(mesh, state, err)
- ! as always, upper surface is the lower surface plus the thickness
- upperSurface(:) = lowerSurface(:) + thickness(:)
+ ! Calculate an initial (diagnostic) velocity field consistent with the initial thickness field
+ ! \todo: skip this step if a velocity field was supplied in the I.C. input file
+ call land_ice_vel_solve(mesh, state, err)
+ ! (halo update of velocity field)
- ! 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 initalize state ====
-
- call land_ice_vel_block_init(mesh, err_tmp)
- err = ior(err, err_tmp)
-
+
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_sia.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-04-18 22:42:58 UTC (rev 1796)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_sia.F 2012-04-19 21:35:38 UTC (rev 1797)
@@ -261,7 +261,7 @@
!***********************************************************************
- subroutine land_ice_sia_solve(mesh, state, normalVelocity, err)
+ subroutine land_ice_sia_solve(mesh, state, err)
use mpas_constants, only: gravity
!-----------------------------------------------------------------
@@ -273,19 +273,14 @@
type (mesh_type), intent(in) :: &
mesh !< Input: mesh information
- type (state_type), intent(in) :: &
- state !< Input: state information
- ! Note: This state is meant to be unmodified (input only).
-
!-----------------------------------------------------------------
!
! input/output variables
!
!-----------------------------------------------------------------
- real (kind=RKIND), intent(inout), dimension(:,:) :: &
- normalVelocity !< Input/Output: the result of the velocity solve
- ! Note: this velocity should be used for both input or output - it may be different than the velocity included in state, and the velocity included in state should not be touched. If land_ice_vel_solve() returns anything else besides normalVelocity, those variables should also be treated in this way. This allow the time integration scheme to decide how to treat multiple time levels.
+ type (state_type), intent(inout) :: &
+ state !< Input: state information
!-----------------------------------------------------------------
!
@@ -303,7 +298,7 @@
real (kind=RKIND), dimension(:), pointer :: thickness, layerThicknessFractions
real (kind=RKIND), dimension(:), pointer :: dcEdge
- real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness, normalVelocity
integer, dimension(:,:), pointer :: cellsOnEdge
real (kind=RKIND), dimension(mesh%nVertLevels) :: layerFractionAboveBed
integer :: nVertLevels, nCells, nEdges, iCell, iLevel, iEdge
@@ -322,6 +317,7 @@
cellsOnEdge => mesh % cellsOnEdge % array
layerThicknessFractions => mesh % layerThicknessFractions % array
+ normalVelocity => state % normalVelocity % array
thickness => state % thickness % array
layerThickness => state % layerThickness % array
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-04-18 22:42:58 UTC (rev 1796)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_tendency.F 2012-04-19 21:35:38 UTC (rev 1797)
@@ -36,7 +36,8 @@
!
!--------------------------------------------------------------------
public :: lice_tend_h, &
- lice_tend_h_layers
+ lice_tend_h_layers, &
+ lice_diagnostic_solve
!--------------------------------------------------------------------
!
@@ -266,4 +267,86 @@
+
+!***********************************************************************
+!
+! 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
+ real (kind=RKIND), dimension(:), pointer :: thickness, upperSurface, &
+ lowerSurface, bedTopography, layerThicknessFractions
+ real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+
+
+ nCells = mesh % nCells
+ nVertLevels = mesh % nVertLevels
+ layerThicknessFractions => mesh % layerThicknessFractions % array
+
+ thickness => state % thickness % array
+ upperSurface => state % upperSurface % array
+ lowerSurface => state % lowerSurface % array
+ bedTopography => state % bedTopography % array
+ layerThickness => state % layerThickness % 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
+
+
+ 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-04-18 22:42:58 UTC (rev 1796)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration.F 2012-04-19 21:35:38 UTC (rev 1797)
@@ -14,7 +14,6 @@
module land_ice_time_integration
- use mpas_vector_reconstruction
use mpas_grid_types
use mpas_configure
use mpas_constants
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-04-18 22:42:58 UTC (rev 1796)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_time_integration_forwardeuler.F 2012-04-19 21:35:38 UTC (rev 1797)
@@ -14,7 +14,6 @@
module land_ice_time_integration_forwardeuler
- use mpas_vector_reconstruction
use mpas_grid_types
use mpas_configure
use mpas_constants
@@ -91,7 +90,10 @@
allocate(mask(mesh%nCells + 1))
+! === Implicit column physics (vertical temperature diffusion) ===========
+ !call ()
+
! === Calculate Tendencies ========================
! Calculate thickness tendency using state at time n
call lice_tend_h(mesh, normalVelocityOld, layerThicknessOld, thicknessOld, thickness_tend, dt, err)
@@ -110,7 +112,7 @@
! === 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)
! Update thickness (and tracers eventually)
thicknessNew = thicknessOld + thickness_tend * dt / SecondsInYear
@@ -134,25 +136,18 @@
! === Calculate diagnostic variables for new state =====================
- ! Remap layer thicknesses (could be moved to diagnostic solve)
- do k = 1, nVertLevels
- layerThicknessNew(k, :) = thicknessNew * layerThicknessFractions(k)
- end do
- ! call land_ice_diagnostic_solve() ! perhaps layer thickness remapping and velocity solve should move in here.
+ call lice_diagnostic_solve(mesh, stateNew, err) ! perhaps velocity solve should move in here.
! Compute the diagnostic velocity for this time step using the updated (new) state
! Pass in the old velocity to be used as an initial guess for velocity solvers that need one.
normalVelocityNew = normalVelocityOld
- call land_ice_vel_solve(mesh, stateNew, normalVelocityNew, err)
+ call land_ice_vel_solve(mesh, stateNew, err)
! update halos on velocity
! call ()
- ! Calculate reconstructed velocities for this time step (could be moved to diagnostic solve or vel_solve)
- call mpas_reconstruct(mesh, normalVelocityNew, &
- uReconstructX, uReconstructY, uReconstructZ, &
- uReconstructZonal, uReconstructMeridional)
+
! === Cleanup & Misc. =============================
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F 2012-04-18 22:42:58 UTC (rev 1796)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_vel.F 2012-04-19 21:35:38 UTC (rev 1797)
@@ -238,9 +238,10 @@
!***********************************************************************
- subroutine land_ice_vel_solve(mesh, state, normalVelocity, err)
+ subroutine land_ice_vel_solve(mesh, state, err)
use land_ice_SIA
+ use mpas_vector_reconstruction
!-----------------------------------------------------------------
!
@@ -251,19 +252,14 @@
type (mesh_type), intent(in) :: &
mesh !< Input: mesh information
- type (state_type), intent(in) :: &
- state !< Input: state information
- ! Note: This state is meant to be unmodified (input only).
-
!-----------------------------------------------------------------
!
! input/output variables
!
!-----------------------------------------------------------------
- real (kind=RKIND), intent(inout), dimension(:,:) :: &
- normalVelocity !< Input/Output: the result of the velocity solve
- ! Note: this velocity should be used for both input or output - it may be different than the velocity included in state, and the velocity included in state should not be touched. If land_ice_vel_solve() returns anything else besides normalVelocity, those variables should also be treated in this way. This allow the time integration scheme to decide how to treat multiple time levels.
+ type (state_type), intent(inout) :: &
+ state !< Input: state information
!-----------------------------------------------------------------
!
@@ -283,15 +279,24 @@
select case (config_dycore)
case ('L1L2')
- call land_ice_lifev_solve(mesh, state, normalVelocity, err)
+ call land_ice_lifev_solve(mesh, state, err)
case ('SIA')
- call land_ice_sia_solve(mesh, state, normalVelocity, err)
+ call land_ice_sia_solve(mesh, state, err)
case default
write(*,*) trim(config_dycore), ' is not a valid land ice dycore option.'
err = 1
return
end select
+ ! Calculate reconstructed velocities for this time step
+ call mpas_reconstruct(mesh, state % normalVelocity % array,&
+ state % uReconstructX % array, &
+ state % uReconstructY % array, &
+ state % uReconstructZ % array, &
+ state % uReconstructZonal % array, &
+ state % uReconstructMeridional % array &
+ )
+
!--------------------------------------------------------------------
end subroutine land_ice_vel_solve
</font>
</pre>