3.1. Implementation

Icepack is written in FORTRAN90 and runs on platforms using UNIX, LINUX, and other operating systems. The code is not parallelized. (CHANGE IF OPENMP IS IMPLEMENTED)

Icepack consists of the sea ice column physics code, contained in the columnphysics/ directory, and a configuration/ directory that includes a driver for testing the column physics and a set of scripts for configuring the tests. Icepack is designed such that the column physics code may be used by a host sea ice model without direct reference to the driver or scripts, although these may be consulted for guidance when coupling the column physics code to the host sea ice model (CICE may also be useful for this.) Information about the interface between the column physics and the driver or host sea ice model is located in the Initialization and Forcing section.

3.1.1. Directory structure

The present code distribution includes source code for the column physics, source code for the driver, and the scripts. Forcing data is available from the ftp site. The directory structure of Icepack is as follows. All columnphysics filename have a prefix of icepack_ and all driver files are prefixed with icedrv_*.

COPYRIGHT.pdf

copyright notice

LICENSE.pdf

license for using and sharing the code

DistributionPolicy.pdf

policy for using and sharing the code

README.md

basic information and pointers

columnphysics/

columnphysics source code, see Icepack Column Physics

configuration/scripts/

support scripts, see Scripts Implementation

configuration/driver/

icepack driver code, see Driver Implementation

doc/

documentation

icepack.setup

main icepack script for creating cases

A case (compile) directory is created upon initial execution of the script icepack.setup at the user-specified location provided after the -c flag. Executing the command ./icepack.setup -h provides helpful information for this tool.

3.1.2. Grid and boundary conditions

The driver configures a collection of grid cells on which the column physics code will be run. This “horizontal” grid is a vector of length nx, with a minimum length of 4. The grid vector is initialized with different sea ice conditions, such as open water, a uniform slab of ice, a multi-year ice thickness distribution with snow, and land. For simplicity, the same forcing values are applied to all grid cells.

Icepack includes two vertical grids. The basic vertical grid contains nilyr equally spaced grid cells. History variables available for column output are ice and snow temperature, Tinz and Tsnz. These variables also include thickness category as a fourth dimension.

In addition, there is a bio-grid that can be more finely resolved and includes additional nodes for boundary conditions. It is used for solving the brine height variable \(h_b\) and for discretizing the vertical transport equations of biogeochemical tracers. The bio-grid is a non-dimensional vertical grid which takes the value zero at \(h_b\) and one at the ice–ocean interface. The number of grid levels is specified during compilation by setting the variable NBGCLYR equal to an integer (\(n_b\)) .

Ice tracers and microstructural properties defined on the bio-grid are referenced in two ways: as bgrid \(=n_b+2\) points and as igrid\(=n_b+1\) points. For both bgrid and igrid, the first and last points reference \(h_b\) and the ice–ocean interface, respectively, and so take the values \(0\) and \(1\), respectively. For bgrid, the interior points \([2, n_b+1]\) are spaced at \(1/n_b\) intervals beginning with bgrid(2) = \(1/(2n_b)\). The igrid interior points \([2, n_b]\) are also equidistant with the same spacing, but physically coincide with points midway between those of bgrid.

3.1.3. Initialization and Forcing

Icepack’s parameters and variables are initialized in several steps. Many constants and physical parameters are set in icepack_parameters.F90. In the current driver implementation, a namelist file is read to setup the model. Namelist values are given default values in the code, which may then be changed when the input file icepack_in is read. Other physical constants, numerical parameters, and variables are first set in initialization routines for each ice model component or module. Then, if the ice model is being restarted from a previous run, core variables are read and reinitialized in restartfile, while tracer variables needed for specific configurations are read in separate restart routines associated with each tracer or specialized parameterization. Finally, albedo and other quantities dependent on the initial ice state are set. Some of these parameters will be described in more detail in the Tables of Namelist Options.

Two namelist variables control model initialization, ice_ic and restart. Setting ice_ic = ‘default’ causes the model to run using initial values set in the code and the namelist. To start from a file filename, set restart = .true. and ice_ic = filename. When restarting using the Icepack driver, for simplicity the tracers are assumed to be set the same way (on/off) as in the run that created the restart file; i.e. that the restart file contains exactly the information needed for the new run. CICE is more flexible in this regard.

When the model is not running from a restart file (i.e., ice_ic = ‘default’ and restart = .false.), there are namelist options that control the initial snow depth , ice thicknesses and mixed layer temperature (sst_init). For the slab-initialized grid cell (nx = 2), the run starts with a single ice thickness category having 100% ice cover . hsno_init_slab and hi_init_slab define the initial snow depth and ice thickness for that ice thickness category. For the itd-initialized grid cell (nx = 3), the ice thickness in each category is set to the midpoint of that category’s ice thickness range (excluding the last category, which is set to 1 m thicker than the lower bound). The area fraction of each category is set according to a normalized, downward-facing parabolic function of ice thickness, where the maximum of the parabola is hbar_init_itd and the area fraction of open water is zero. All thickness categories are initialized with a snow depth of hsno_init_itd.

For stand-alone runs, routines in icedrv_forcing.F90 read and interpolate data from files. The namelist variables precalc_forc, atm_data_type, atm_data_format, ocn_data_type, and ocn_data_format control how the forcing data is handled. If precalc_forc = .false. and the atm/ocn_data_type = ‘bin’, when init_forcing is called, a subroutine for the specific dataset (e.g., atm_CFS) stores the forcing data in the *_data arrays in essentially the same format that the raw data is present in, without timestamp information. Then, at each timestep the get_forcing subroutine has a code block for each forcing dataset that contains the forcing’s time basis and interpolates the forcing data to the given timestep. The forcing data that is available in ‘bin’ format are intended merely for code testing, not scientific results.

If precalc_forc = .true., atm/ocn_data_type = ‘MDF’, and atm/ocn_data_format = ‘nc’, then init_forcing reads data from netCDF files formatted according to the Merged Data File (MDF) conventions, which includes timestamp information. During initialization, the forcing data is averaged/interpolated to the Icepack timestep and stored the the *_data arrays. The get_forcing subroutine then simply queries the *_data arrays at each timestep. MDF forcing data are expected to come from observations, and hence may contain missing data and may be present at a shorter or longer sampling interval than the Icepack timestep. To handle variable frequencies and missing data, forcing data are first temporally-averaged for each timestep and then interpolated. For a give data variable (var_data), The MDF_average subroutine takes the average of all forcing datapoints within each timestep +- 0.5 dt excluding missing values and stores the results in var_data. If there is no valid data within 0.5 dt of a given timestep (e.g., most timesteps if dt is much smaller than the sampling interval) then a missing value is placed in var_data(timestep). Then, the MDF_interpolate subroutine linearly interpolates missing values in var_data. The MDF conventions were developed by the Year of Polar Prediction supersite Model Intercomparison Project (Uttal et al., 2024) and a python toolbox is available to build MDF files from raw data. The ocn_MDF subroutine currently assumes that the oceanic heat flux convergence (qdp) is equal to the turbulent heat flux over the thermocline.

If no ocean forcing is provided, namelist variables provide constant values of the ocean mixed layer salinity (sss_fixed), thickness (hmix_fixed), and oceanic heat flux convergence (qdp_fixed). If forcing data is provided then these variables are ignored.

3.1.4. Choosing an appropriate time step

Transport in thickness space imposes a restraint on the time step, given by the ice growth/melt rate and the smallest range of thickness among the categories, \(\Delta t<\min(\Delta H)/2\max(f)\), where \(\Delta H\) is the distance between category boundaries and \(f\) is the thermodynamic growth rate. For the 5-category ice thickness distribution used as the default in this distribution, this is not a stringent limitation: \(\Delta t < 19.4\) hr, assuming \(\max(f) = 40\) cm/day.

3.1.5. Model output

The Icepack model provides diagnostic output files, binary or netCDF restart files, and a primitive netCDF history file capability. The sea ice model CICE provides more extensive options for model output, including many derived output variables.

3.1.5.1. Diagnostic files

Icepack writes diagnostic information for each grid cell as a separate file, ice_diag.*, identified by the initial ice state of the grid cell (no ice, slab, land, etc).

3.1.5.2. Restart files

Icepack provides restart data in binary unformatted format or netCDF. The restart files created by the Icepack driver contain all of the variables needed for a full, exact restart. The filename begins with the character string ‘iced.’ and is placed in the directory specified by the namelist variable restart_dir. The restart dump frequency is given by the namelist variable dumpfreq. The namelist variable ice_ic contains the pointer to the filename from which the restart data is to be read and the namelist option restart must be set to .true. to use the file. dump_last namelist can also be set to true to trigger restarts automatically at the end of runs. The default restart file format is binary, set in namelist with restart_format = ‘bin’. For netCDF, set restart_format = ‘nc’ or use icepack.setup -s restcdf.

The default configuration of Icepack does not support netCDF. If netCDF restart files are desired, the USE_NETCDF C preprocessor directive must be set during compilation. This is done by setting ICE_IOTYPE to netcdf in icepack.settings or using the icepack.setup -s option ionetcdf. If netCDF is used on a particular machine, the machine env and Macros file must support compilation with netCDF.

3.1.5.3. History files

Icepack has a primitive netCDF history capability that is turned on with the history_format namelist. When history_format is set to ‘nc’, history files are created for each run with a naming convention of icepack.h.yyyymmdd.nc in the run directory history directory. The yyyymmdd is the start date for each run. Use icepack.setup -s histcdf to turn on netCDF history files automatically.

When Icepack history files are turned on, data for a set of fixed fields is written to the history file for each column at every timestep without ability to control fields, frequencies, or temporal averaging. All output fields are hardwired into the implementation in configuration/driver/icedrv_history.F90 file. The netCDF file does NOT meet netCDF CF conventions and is provided as an amenity in the standalone Icepack model. Users are free to modify the output fields or extend the implementation and are encouraged to share any updates with the Consortium.

The default configuration of Icepack does not support netCDF. If netCDF history files are desired, the USE_NETCDF C preprocessor directive must be set during compilation. This is done by setting ICE_IOTYPE to netcdf in icepack.settings or using the icepack.setup -s option ionetcdf. If netCDF is used on a particular machine, the machine env and Macros file must support compilation with netCDF.

Note that some of the ponds history fields are not yet implemented for the topo ponds option.

3.1.5.4. Biogeochemistry History Fields

History output is not provided with Icepack. This documentation indicates what is available for output and is implemented in CICE.

Table Biogeochemical History variables lists the biogeochemical tracer history flags along with a short description and the variable or variables saved. Not listed are flags appended with _ai, i.e. f_fbio_ai. These fields are identical to their counterpart. i.e. f_fbio, except they are averaged by ice area.

Biogeochemical History variables

History Flag	Definition	Variable(s)	Units
f_fiso_atm	atmospheric water isotope deposition flux	fiso_atm	kg m\(^{-2}\) s\(^{-1}\)
f_fiso_ocn	water isotope flux from ice to ocean	fiso_ocn	kg m\(^{-2}\) s\(^{-1}\)
f_iso	isotope mass (snow and ice)	isosno, isoice	kg/kg
f_faero_atm	atmospheric aerosol deposition flux	faero_atm	kg m\(^{-2}\) s\(^{-1}\)
f_faero_ocn	aerosol flux from ice to ocean	faero_ocn	kg m\(^{-2}\) s\(^{-1}\)
f_aero	aerosol mass (snow and ice ssl and int)	aerosnossl, aerosnoint,aeroicessl, aeroiceint	kg/kg
f_fbio	biological ice to ocean flux	fN, fDOC, fNit, fAm,fDON,fFep\(^a\), fFed\(^a\), fSil,fhum, fPON, fDMSPd,fDMS, fDMSPp, fzaero	mmol m\(^{-2}\) s\(^{-1}\)
f_zaero	bulk z-aerosol mass fraction	zaero	kg/kg
f_bgc_S	DEPRECATED	bgc_S	ppt
f_bgc_N	bulk algal N concentration	bgc_N	mmol m\(^{-3}\)
f_bgc_C	bulk algal C concentration	bgc_C	mmol m\(^{-3}\)
f_bgc_DOC	bulk DOC concentration	bgc_DOC	mmol m\(^{-3}\)
f_bgc_DON	bulk DON concentration	bgc_DON	mmol m\(^{-3}\)
f_bgc_DIC	bulk DIC concentration	bgc_DIC	mmol m\(^{-3}\)
f_bgc_chl	bulk algal chlorophyll concentration	bgc_chl	mg chl m\(^{-3}\)
f_bgc_Nit	bulk nitrate concentration	bgc_Nit	mmol m\(^{-3}\)
f_bgc_Am	bulk ammonium concentration	bgc_Am	mmol m\(^{-3}\)
f_bgc_Sil	bulk silicate concentration	bgc_Sil	mmol m\(^{-3}\)
f_bgc_DMSPp	bulk particulate DMSP concentration	bgc_DMSPp	mmol m\(^{-3}\)
f_bgc_DMSPd	bulk dissolved DMSP concentration	bgc_DMSPd	mmol m\(^{-3}\)
f_bgc_DMS	bulk DMS concentration	bgc_DMS	mmol m\(^{-3}\)
f_bgc_Fe	bulk dissolved and particulate iron conc.	bgc_Fed, bgc_Fep	\(\mu\,\)mol m\(^{-3}\)
f_bgc_hum	bulk humic matter concentration	bgc_hum	mmol m\(^{-3}\)
f_bgc_PON	bulk passive mobile tracer conc.	bgc_PON	mmol m\(^{-3}\)
f_upNO	Total algal \({\mbox{NO$_3$}}\) uptake rate	upNO	mmol m\(^{-2}\) d\(^{-1}\)
f_upNH	Total algal \({\mbox{NH$_4$}}\) uptake rate	upNH	mmol m\(^{-2}\) d\(^{-1}\)
f_bgc_ml	upper ocean tracer concentrations	ml_N, ml_DOC, ml_Nit,ml_Am, ml_DON, ml_Fep\(^b\),ml_Fed\(^b\), ml_Sil, ml_hum, ml_PON,ml_DMS, ml_DMSPd, ml_DMSPp	mmol m\(^{-3}\)
f_bTin	ice temperature on the bio grid	bTizn	\(^o\)C
f_bphi	ice porosity on the bio grid	bphizn	%
f_iDin	brine eddy diffusivity on the interface bio grid	iDin	m\(^{2}\) d\(^{-1}\)
f_iki	ice permeability on the interface bio grid	ikin	mm\(^{2}\)
f_fbri	ratio of brine tracer height to ice thickness	fbrine	1
f_hbri	brine tracer height	hbrine	m
f_zfswin	internal ice PAR on the interface bio grid	zfswin	W m\(^{-2}\)
f_bionet	brine height integrated tracer concentration	algalN_net, algalC_net, chl_net, pFe\(^c\)_net, dFe\(^c\)_net, Sil_net, Nit_net, Am_net, hum_net, PON_net, DMS_net, DMSPd_net, DMSPp_net, DOC_net, zaero_net, DON_net	mmol m\(^{-2}\)
f_biosnow	snow integrated tracer concentration”	algalN_snow, algalC_snow,chl_snow, pFe\(^c\)_snow, dFe\(^c\)_snow,Sil_snow, Nit_snow, Am_snow, hum_snow, PON_snow, DMS_snow, DMSPd_snow, DMSPp_snow, DOC_snow, zaero_snow, DON_snow	mmol m\(^{-2}\)
f_grownet	Net specific algal growth rate	grow_net	m d\(^{-1}\)
f_PPnet	Net primary production	PP_net	mgC m\(^{-2}\) d\(^{-1}\)
f_algalpeak	interface bio grid level of peak chla	peak_loc	1
f_zbgc_frac	mobile fraction of tracer	algalN_frac, chl_frac, pFe_frac,dFe_frac, Sil_frac, Nit_frac,Am_frac, hum_frac, PON_frac,DMS_frac, DMSPd_frac, DMSPp_frac,DOC_frac, zaero_frac, DON_frac	1

\(^a\) units are \(\mu\)mol m\(^{-2}\) s\(^{-1}\)

\(^b\) units are \(\mu\)mol m\(^{-3}\)

\(^c\) units are \(\mu\)mol m\(^{-2}\)