<p><b>mhoffman@lanl.gov</b> 2012-06-08 13:00:37 -0600 (Fri, 08 Jun 2012)</p><p>BRANCH COMMIT -- land ice<br>
<br>
Created option for using annual forcing from the input file. <br>
<br>
See Registry for invocation: config_forcing_frequency = 'Annual'<br>
The forcing is setup in a new module: land_ice_annual_forcing<br>
<br>
* This commit will break any existing grids, because a land ice grid will now need new fields: betaTimeSeries, sfcMassBalTimeSeries, sfcAirTempTimeSeries, basalHeatFluxTimeSeries<br>
* These fields need to be present even if time-series are not supplied (in which case they would have a single time level).<br>
* The scripts in grid_tools and test_cases have been modified to work with the requirements, so any new grids generated with them will work with this revision of the land ice core.<br>
* Using annual forcing with large datasets may create memory problems since the entire forcing time-series is read into memory during init and held for the entire simulation.<br>
[On my laptop with 4 processors/8GB total memory, adding a single 100-yr forcing field for 5km Greenland (346,000 cells) created only a minor performance hit after file input was complete. However, adding 2 200-yr forcing fields created a large performance hit (the time spent on time integration was 2-4x greater despite not touching that part of the code).]<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/grid_tools/add_land_ice_variables_to_mpas_grid.py
===================================================================
--- branches/land_ice_projects/grid_tools/add_land_ice_variables_to_mpas_grid.py 2012-06-07 17:25:03 UTC (rev 1973)
+++ branches/land_ice_projects/grid_tools/add_land_ice_variables_to_mpas_grid.py 2012-06-08 19:00:37 UTC (rev 1974)
@@ -65,6 +65,11 @@
# Create nVertLevelsPlus2 dimension
# fileout.createDimension('nVertLevelsPlus2', filein.dimensions['nVertLevels'] + 2)
+# Create the dimensions needed for time-dependent forcings
+fileout.createDimension('nBetaTimeSlices', 1)
+fileout.createDimension('nSfcMassBalTimeSlices', 1)
+fileout.createDimension('nSfcAirTempTimeSlices', 1)
+fileout.createDimension('nBasalHeatFluxTimeSlices', 1)
# Copy over all of the required grid variables to the new file
@@ -102,14 +107,15 @@
newvar = fileout.createVariable('temperature', datatype, ('Time', 'nCells', 'nVertLevels'))
#newvar = fileout.createVariable('temperature', 'd', ('Time', 'nCells', 'nVertLevelsPlus2'))
newvar[:] = numpy.zeros(newvar.shape)
+
# These boundary conditions are currently part of mesh, and are time independent. If they change, make sure to adjust the dimensions here and in Registry.
-newvar = fileout.createVariable('beta', datatype, ( 'nCells',))
+newvar = fileout.createVariable('betaTimeSeries', datatype, ( 'nCells', 'nBetaTimeSlices', ))
newvar[:] = numpy.zeros(newvar.shape)
-newvar = fileout.createVariable('sfcMassBal', datatype, ('nCells',))
+newvar = fileout.createVariable('sfcMassBalTimeSeries', datatype, ( 'nCells', 'nSfcMassBalTimeSlices', ))
newvar[:] = numpy.zeros(newvar.shape)
-newvar = fileout.createVariable('sfcAirTemp', datatype, ('nCells',))
+newvar = fileout.createVariable('sfcAirTempTimeSeries', datatype, ( 'nCells', 'nSfcAirTempTimeSlices', ))
newvar[:] = numpy.zeros(newvar.shape)
-newvar = fileout.createVariable('basalHeatFlux', datatype, ('nCells',))
+newvar = fileout.createVariable('basalHeatFluxTimeSeries', datatype, ( 'nCells', 'nBasalHeatFluxTimeSlices',))
newvar[:] = numpy.zeros(newvar.shape)
# Assign the global attributes
Modified: branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py
===================================================================
--- branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py 2012-06-07 17:25:03 UTC (rev 1973)
+++ branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py 2012-06-08 19:00:37 UTC (rev 1974)
@@ -193,10 +193,11 @@
sf = 1.0
if netCDF_module == 'netCDF4':
sf = 1.0 # netCDF4 should autoscale!
- beta = outfile.variables['beta']
+ beta = outfile.variables['betaTimeSeries']
print '</font>
<font color="gray">input beta min/max', inbeta[:].min()*sf, inbeta[:].max()*sf
#beta[timelevout,:] = BilinearInterp(x0, y0, inbeta[timelev,:,:]*sf, xCell, yCell)
- beta[:] = BilinearInterp(x0, y0, inbeta[timelev,:,:]*sf, xCell, yCell)
+ beta[:,0] = BilinearInterp(x0, y0, inbeta[timelev,:,:]*sf, xCell, yCell)
+ print 'interim beta min/max', beta[:].min(), beta[:].max()
# Make all beta be positive values
b = beta[:]
b[b<=0.0] = 1.0
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Makefile
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-06-07 17:25:03 UTC (rev 1973)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Makefile 2012-06-08 19:00:37 UTC (rev 1974)
@@ -20,6 +20,7 @@
mpas_land_ice_lifev.o \
mpas_land_ice_sia.o \
mpas_land_ice_global_diagnostics.o \
+ mpas_land_ice_annual_forcing.o \
mpas_ocn_tracer_advection.o \
mpas_ocn_tracer_advection_std.o \
mpas_ocn_tracer_advection_std_hadv.o \
@@ -46,13 +47,15 @@
mpas_land_ice_tendency.o: mpas_land_ice_mask.o mpas_ocn_tracer_advection.o
+mpas_land_ice_annual_forcing.o:
+
mpas_land_ice_lifev.o:
mpas_land_ice_sia.o:
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_ocn_advection.o
+mpas_land_ice_mpas_core.o: mpas_land_ice_global_diagnostics.o mpas_land_ice_time_integration.o mpas_ocn_advection.o mpas_land_ice_annual_forcing.o
# modules needed for FCT tracer advection from ocean core (eventually to be in operators?)
Modified: branches/land_ice_projects/implement_core/src/core_land_ice/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-06-07 17:25:03 UTC (rev 1973)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/Registry 2012-06-08 19:00:37 UTC (rev 1974)
@@ -27,6 +27,8 @@
namelist real land_ice_model config_sealevel 0.0
namelist real land_ice_model config_dynamic_thickness 100.0
% Define the ice thickness below which dynamics are not calculated.
+namelist character land_ice_model config_forcing_frequency Static
+% Valid options are Static, Annual
% The following options are used by the ocean core and may be useful in the future
%namelist logical land_ice_model config_h_ScaleWithMesh false
@@ -80,9 +82,12 @@
% nvertLevels+2 is only needed if attaching b.c. to a 3-d field. We are unlikely to use this, but I am retaining it for now. MJH 04/19/12
% dim nVertLevelsPlus2 nVertLevels+2
dim nVertLevelsPlus1 nVertLevels+1
+dim nSfcMassBalTimeSlices nSfcMassBalTimeSlices
+dim nBetaTimeSlices nBetaTimeSlices
+dim nSfcAirTempTimeSlices nSfcAirTempTimeSlices
+dim nBasalHeatFluxTimeSlices nBasalHeatFluxTimeSlices
-
% VARIABLES =====================
%
% var persistence type name_in_file ( dims ) time_levs iro- name_in_code struct super-array array_class
@@ -181,12 +186,18 @@
% Boundary Conditions (defined as mesh variables because they are read in once and held constant for the simulation - at least for now)
%
% Note: beta should be moved to diagnostic fields when a process model comes online
-var persistent real beta ( nCells ) 0 iro beta mesh - -
-var persistent real basalHeatFlux ( nCells ) 0 iro basalHeatFlux mesh - -
-var persistent real sfcAirTemp ( nCells ) 0 iro sfcAirTemp mesh - -
-var persistent real sfcMassBal ( nCells ) 0 iro sfcMassBal mesh - -
+var persistent real beta ( nCells ) 0 - beta mesh - -
+var persistent real basalHeatFlux ( nCells ) 0 - basalHeatFlux mesh - -
+var persistent real sfcAirTemp ( nCells ) 0 - sfcAirTemp mesh - -
+var persistent real sfcMassBal ( nCells ) 0 - sfcMassBal mesh - -
+% Time-varying forcing for B.C.
+var persistent real betaTimeSeries ( nBetaTimeSlices nCells ) 0 i betaTimeSeries mesh - -
+var persistent real basalHeatFluxTimeSeries ( nBasalHeatFluxTimeSlices nCells ) 0 i basalHeatFluxTimeSeries mesh - -
+var persistent real sfcAirTempTimeSeries ( nSfcAirTempTimeSlices nCells ) 0 i sfcAirTempTimeSeries mesh - -
+var persistent real sfcMassBalTimeSeries ( nSfcMassBalTimeSlices nCells ) 0 i sfcMassBalTimeSeries mesh - -
+
% STATE VARIABLES =====================
% Prognostic variables: read from input, saved in restart, and written to output
Added: branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_annual_forcing.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_annual_forcing.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_annual_forcing.F 2012-06-08 19:00:37 UTC (rev 1974)
@@ -0,0 +1,206 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! ocn_monthly_forcing
+!
+!> \brief MPAS land ice annual forcing
+!> \author Matt Hoffman
+!> \date 06/07/12
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for building the forcing arrays,
+!> if annual forcing is used. Modified from ocn_monthly_forcing.
+!
+!-----------------------------------------------------------------------
+
+module land_ice_annual_forcing
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_timekeeping
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: land_ice_assign_annual_forcing
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine land_ice_assign_annual_forcing
+!
+!> \brief Determines the forcing array used for the annual forcing.
+!> \author Matt Hoffman
+!> \date 06/07/12
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the forcing arrays used later in MPAS.
+!
+!-----------------------------------------------------------------------
+
+ subroutine land_ice_assign_annual_forcing(timeStamp, grid, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type(MPAS_Time_type), intent(in) :: timeStamp
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(inout) :: &
+ grid !< Input: grid information
+
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: Error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+ real (kind=RKIND), dimension(:,:), pointer :: sfcMassBalTimeSeries, sfcAirTempTimeSeries, betaTimeSeries, basalHeatFluxTimeSeries
+ real (kind=RKIND), dimension(:), pointer :: sfcMassBal, sfcAirTemp, beta, basalHeatFlux
+ integer :: iYear, ierr, nCells
+
+
+ err = 0
+
+ sfcMassBalTimeSeries => grid % sfcMassBalTimeSeries % array
+ sfcAirTempTimeSeries => grid % sfcAirTempTimeSeries % array
+ betaTimeSeries => grid % betaTimeSeries % array
+ basalHeatFluxTimeSeries => grid % basalHeatFluxTimeSeries % array
+
+ sfcMassBal => grid % sfcMassBal % array
+ sfcAirTemp => grid % sfcAirTemp % array
+ beta => grid % beta % array
+ basalHeatFlux => grid % basalHeatFlux % array
+
+ nCells = grid % nCells
+
+ if (config_forcing_frequency .eq. 'Static') then
+ iYear = 1 ! Static forcing uses the first entry in the time series at all time steps of the simulation
+ else
+ ! Get the index to use for the time series. This assumes that the starting year
+ ! is year 1 and therefore can be used as the index into the time series.
+ call mpas_get_time(timeStamp, YYYY = iYear, ierr = ierr)
+ err = ierr
+ endif
+
+ ! Call this separately for each field - there will be cases where one field is time varying and others are constant, and this routine will deal with that properly.
+ call assign_annual_forcing_to_1D_array(sfcMassBal, sfcMassBalTimeSeries, iYear, nCells)
+ call assign_annual_forcing_to_1D_array(sfcAirTemp, sfcAirTempTimeSeries, iYear, nCells)
+ call assign_annual_forcing_to_1D_array(beta, betaTimeSeries, iYear, nCells)
+ call assign_annual_forcing_to_1D_array(basalHeatFlux, basalHeatFluxTimeSeries, iYear, nCells)
+
+ end subroutine land_ice_assign_annual_forcing!}}}
+
+ !--------------------------------------------------------------------
+
+
+! Private subroutines=======================
+
+
+!***********************************************************************
+!
+! routine assign_annual_forcing_to_1D_array
+!
+!> \brief Assigns the appropriate annual forcing to the array used by MPAS.
+!> \author Matt Hoffman
+!> \date 06/07/12
+!> \version SVN:$Id$
+!> \details
+!> This routine assigns the forcing for a given field to be used later in MPAS.
+!
+!-----------------------------------------------------------------------
+
+ subroutine assign_annual_forcing_to_1D_array(forcing, forcingAnnual, iYear, nCells)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:) :: forcingAnnual
+ integer, intent(in) :: iYear, nCells
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:) :: forcing
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+ integer :: timeIndex
+
+ ! Make sure that only valid time levels are used.
+ if (iYear .lt. 1) then
+ timeIndex = 1
+ elseif (iYear .gt. size(forcingAnnual, 1) ) then
+ timeIndex = size(forcingAnnual, 1) ! Hold the final value if you run out of time levels - useful if only 1 time level is supplied.
+ else
+ timeIndex = iYear
+ endif
+
+ ! Assign the field
+ forcing(1:nCells) = forcingAnnual( timeIndex, 1:nCells )
+
+ end subroutine assign_annual_forcing_to_1D_array
+
+
+
+!***********************************************************************
+
+end module land_ice_annual_forcing
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
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-06-07 17:25:03 UTC (rev 1973)
+++ branches/land_ice_projects/implement_core/src/core_land_ice/mpas_land_ice_mpas_core.F 2012-06-08 19:00:37 UTC (rev 1974)
@@ -157,12 +157,14 @@
use land_ice_tendency, only: land_ice_diagnostic_solve
use mpas_ocn_tracer_advection
use ocn_advection
+ use land_ice_annual_forcing
implicit none
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
+ type (MPAS_Time_Type) :: currTime
type (dm_info) :: dminfo
integer :: err, err_tmp
@@ -175,6 +177,10 @@
! Call init routines ====
call land_ice_setup_vertical_coords(mesh)
+ currTime = mpas_get_clock_time(clock, MPAS_NOW, err_tmp)
+ err = ior(err, err_tmp)
+ call land_ice_assign_annual_forcing(currTime, mesh, err_tmp)
+ err = ior(err, err_tmp)
! Init for FCT tracer advection
mesh % maxLevelCell % array = mesh % nVertLevels ! Needed for FCT tracer advection
@@ -221,6 +227,7 @@
use mpas_grid_types
use mpas_io_output
use mpas_timer
+ use land_ice_annual_forcing
implicit none
@@ -241,7 +248,9 @@
write(0,*) 'Initial timestep ', timeStamp
write(6,*) 'Initial timestep ', timeStamp
+ call mpas_timer_start("write output frame")
call write_output_frame(output_obj, output_frame, domain)
+ call mpas_timer_stop("write output frame")
! 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(...)
@@ -256,6 +265,14 @@
write(0,*) 'Doing timestep ', timeStamp
write(6,*) 'Doing timestep ', timeStamp
+ call mpas_timer_start("assign forcing fields")
+ block_ptr => domain % blocklist
+ do while(associated(block_ptr))
+ call land_ice_assign_annual_forcing(currTime, block_ptr % mesh, ierr)
+ block_ptr => block_ptr % next
+ end do
+ call mpas_timer_stop("assign forcing fields")
+
call mpas_timer_start("time integration")
call mpas_timestep(domain, itimestep, dt, timeStamp)
call mpas_timer_stop("time integration")
Modified: branches/land_ice_projects/test_cases/dome/setup_dome_initial_conditions.py
===================================================================
--- branches/land_ice_projects/test_cases/dome/setup_dome_initial_conditions.py 2012-06-07 17:25:03 UTC (rev 1973)
+++ branches/land_ice_projects/test_cases/dome/setup_dome_initial_conditions.py 2012-06-08 19:00:37 UTC (rev 1974)
@@ -46,10 +46,11 @@
normalVelocity = gridfile.variables['normalVelocity'][:]
layerThicknessFractions = gridfile.variables['layerThicknessFractions'][:]
temperature = gridfile.variables['temperature'][:]
- beta = gridfile.variables['beta'][:]
- SMB = gridfile.variables['sfcMassBal'][:]
- Tsfc = gridfile.variables['sfcAirTemp'][:]
- G = gridfile.variables['basalHeatFlux'][:]
+ # Get b.c. variables
+ beta = gridfile.variables['betaTimeSeries'][:]
+ SMB = gridfile.variables['sfcMassBalTimeSeries'][:]
+ Tsfc = gridfile.variables['sfcAirTempTimeSeries'][:]
+ G = gridfile.variables['basalHeatFluxTimeSeries'][:]
except:
sys.exit('Error: The grid file specified is either missing or lacking needed dimensions/variables.')
@@ -78,21 +79,21 @@
# boundary conditions
# Sample values to use, or comment these out for them to be 0.
beta[:] = 10000.
-SMB[:] = 2.0/1000.0 * (thickness[:] + bedTopography[:]) - 1.0 # units: m weq/yr, lapse rate of 1 m/yr with 0 at 500 m
+#SMB[:] = 2.0/1000.0 * (thickness[:] + bedTopography[:]) - 1.0 # units: m/yr, lapse rate of 1 m/yr with 0 at 500 m
Tsfc[:] = -5.0/1000.0 * (thickness[:] + bedTopography[:]) # lapse rate of 5 deg / km
G = 0.01
-# Reassign the modified numpy array values back into netCDF variable objects
+# Reassign the modified numpy array values back into netCDF variable objects
gridfile.variables['thickness'][:] = thickness
gridfile.variables['normalVelocity'][:] = normalVelocity
gridfile.variables['bedTopography'][:] = bedTopography
gridfile.variables['temperature'][:] = temperature
gridfile.variables['layerThicknessFractions'][:] = layerThicknessFractions
-gridfile.variables['beta'][:] = beta
-gridfile.variables['sfcMassBal'][:] = SMB
-gridfile.variables['sfcAirTemp'][:] = Tsfc
-gridfile.variables['basalHeatFlux'][:] = G
+gridfile.variables['betaTimeSeries'][:] = beta
+gridfile.variables['sfcMassBalTimeSeries'][:] = SMB
+gridfile.variables['sfcAirTempTimeSeries'][:] = Tsfc
+gridfile.variables['basalHeatFluxTimeSeries'][:] = G
gridfile.close()
print 'Successfully added dome initial conditions to ', gridfilename
</font>
</pre>