s_phys_pkgs/text/thsice.tex

% \documentclass[12pt]{article}
% \usepackage{amssymb}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  \usepackage{graphics}


% \oddsidemargin -4mm \evensidemargin 0mm
% \textwidth 165mm
% \textheight 230mm
% \topmargin -2mm \headsep -2mm
% \renewcommand{\baselinestretch}{1.5}
% \begin{document}


\def\deg{$^o$}
%%%--------------------------------------%%%
\subsection{THSICE: The Thermodynamic Sea Ice Package}
\label{sec:pkg:thsice}
\begin{rawhtml}
<!-- CMIREDIR:package_thsice: -->
\end{rawhtml}

{\bf Important note:}
This document has been written by Stephanie Dutkiewicz
and describes an earlier implementation of the sea-ice package.
This needs to be updated to reflect the recent changes (JMC).

\noindent
This thermodynamic ice model is based on the 3-layer model by Winton (2000).
and the energy-conserving LANL CICE model (Bitz and Lipscomb, 1999).
The model considers two equally thick ice layers; the upper layer has
a variable specific heat resulting from brine pockets,
the lower layer has a fixed heat capacity. A zero heat capacity snow
layer lies above the ice. Heat fluxes at the top and bottom
surfaces are used to calculate the change in ice and snow layer
thickness. Grid cells of the ocean model are 
either fully covered in ice or are open water. There is
a provision to parametrize ice fraction (and leads) in this package.
Modifications are discussed in small font following the
subroutine descriptions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{1cm}

\noindent
The ice model is called from {\it thermodynamics.F}, subroutine
{\it ice\_forcing.F} is called in place of {\it external\_forcing\_surf.F}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{1cm}
\noindent
{\bf \underline{subroutine ICE\_FORCING}}

\noindent
In {\it ice\_forcing.F}, we calculate the freezing potential of the
ocean model surface layer of water:
\[
  {\bf frzmlt} = (T_f - SST) \frac{c_{sw} \rho_{sw} \Delta z}{\Delta t}
\]
where $c_{sw}$ is seawater heat capacity, 
$\rho_{sw}$ is the seawater density, $\Delta z$
is the ocean model upper layer thickness and $\Delta t$ is the model (tracer)
timestep. The freezing temperature, $T_f=\mu S$ is a function of the
salinity.


1) Provided there is no ice present and {\bf frzmlt} is less than 0,
   the surface tendencies of wind, heat and freshwater are calculated
   as usual (ie. as in {\it external\_forcing\_surf.F}).

2) If there is ice present in the grid cell
   we call the main ice model routine {\it ice\_therm.F} (see below).
   Output from this routine gives net heat and freshwater flux 
   affecting the top of the ocean.

Subroutine {\it ice\_forcing.F} uses these values to find the 
sea surface tendencies
in grid cells. When there is ice present,  
the surface stress tendencies are
set to zero; the ice model is purely thermodynamic and the
effect of ice motion on the sea-surface is not examined.

Relaxation of surface $T$ and $S$ is only allowed equatorward
of {\bf relaxlat} (see {\bf DATA.ICE below}), and no relaxation is
allowed under the ice at any latitude.

\noindent
{\tiny (Note that there is provision for allowing grid cells to have both
open water and seaice; if {\bf compact} is between  0 and 1)}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{1cm}
\noindent
{\bf {\underline{ subroutine ICE\_FREEZE}}}

This routine is called from {\it thermodynamics.F}
after the new temperature calculation, {\it calc\_gt.F}, 
but before {\it calc\_gs.F}.
In {\it ice\_freeze.F}, any ocean upper layer grid cell
with no ice cover, but with temperature below freezing,
$T_f=\mu S$ has ice initialized.
We calculate {\bf frzmlt} from all the grid cells in
the water column that have a temperature less than
freezing. In this routine, any water below the surface
that is below freezing is set to $T_f$.
A call to
{\it ice\_start.F} is made if {\bf frzmlt} $>0$, 
and salinity tendancy is updated for brine release.

\noindent
{\tiny (There is a provision for fractional ice:
In the case where the grid cell has less ice coverage than
{\bf icemaskmax} we allow {\it ice\_start.F} to be called).}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{1cm}
\noindent
{\bf {\underline{ subroutine ICE\_START}}}

\noindent
The energy available from freezing
the sea surface is brought into this routine as {\bf esurp}.
The enthalpy of the 2 layers of any new ice is calculated as:
\begin{eqnarray}
q_1 & = & -c_{i}*T_f + L_i \nonumber \\
q_2 & = & -c_{f}T_{mlt}+ c_{i}(T_{mlt}-T{f}) + L_i(1-\frac{T_{mlt}}{T_f} 
\nonumber \\
\end{eqnarray}
where  $c_f$ is specific heat of liquid fresh water, $c_i$ is the
specific heat of fresh ice, $L_i$ is latent heat of freezing, 
$\rho_i$ is density of ice and
$T_{mlt}$ is melting temperature of ice with salinity of 1.
The height of a new layer of ice is
\[
  h_{i new} = \frac{{\bf esurp} \Delta t}{qi_{0av}}
\]
where $qi_{0av}=-\frac{\rho_i}{2} (q_1+q_2)$.

The surface skin temperature $T_s$ and ice temperatures
$T_1$, $T_2$ and the sea surface temperature are set at $T_f$.

\noindent
{\tiny ( There is provision for fractional ice:
new ice is formed over open water; the first freezing in the cell
must have a height of {\bf himin0}; this determines the ice
fraction {\bf compact}. If there is already ice in the grid cell,
the new ice must have the same height and the new ice fraction
is 
\[
i_f=(1-\hat{i_f}) \frac{h_{i new}}{h_i}
\]
where $\hat{i_f}$ is ice fraction from previous timestep
and $h_i$ is current ice height. Snow is redistributed 
over the new ice fraction. The ice fraction is
not allowed to become larger than {\bf iceMaskmax} and
if the ice height is above {\bf hihig} then freezing energy
comes from the full grid cell,  ice growth does not occur
under orginal ice due to freezing water.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{1cm}
\noindent
{\bf {\underline{subroutine ICE\_THERM}}}

\noindent
The main subroutine of this package is {\it ice\_therm.F} where the
ice temperatures are calculated and the changes in ice and snow
thicknesses are determined. Output provides the net heat and fresh 
water fluxes that force the top layer of the ocean model.

If the current ice height is less than {\bf himin} then
the ice layer is set to zero and the ocean model upper layer temperature
is allowed to drop lower than its freezing temperature; and atmospheric
fluxes are allowed to effect the grid cell.
If the ice height is greater than  {\bf himin} we proceed with
the ice model calculation.

We follow the procedure
of Winton (1999) -- see equations 3 to 21 -- to calculate
the surface and internal ice temperatures. 
The surface temperature is found from the balance of the
flux at the surface $F_s$, the shortwave heat flux absorbed by the ice, 
{\bf fswint}, and
the upward conduction of heat through the snow and/or ice, $F_u$.
We linearize $F_s$ about the surface temperature, $\hat{T_s}$, 
at the previous timestep (where \mbox{}$\hat{ }$ indicates the value at
the  previous timestep):
\[
F_s (T_s) = F_s(\hat{T_s}) + \frac{\partial F_s(\hat{T_s)}}{\partial T_s}
(T_s-\hat{T_s})
\]
where, 
\[
F_s  =  F_{sensible}+F_{latent}+F_{longwave}^{down}+F_{longwave}^{up}+ (1-
\alpha) F_{shortwave}
\]
and
\[
 \frac{d F_s}{dT} = \frac{d F_{sensible}}{dT} + \frac{d F_{latent}}{dT}
+\frac{d F_{longwave}^{up}}{dT}.
\]
$F_s$ and $\frac{d F_s}{dT}$ are currently calculated from the {\bf BULKF} 
package described separately, but could also be provided by an atmospheric
model. The surface albedo is calculated from the ice height and/or 
surface temperature (see below, {\it srf\_albedo.F}) and the 
shortwave flux absorbed in the ice is
\[
{\bf fswint} = (1-e^{\kappa_i h_i})(1-\alpha) F_{shortwave}
\]
where $\kappa_i$ is bulk extinction coefficient.

The conductive flux to the surface is
\[
F_u=K_{1/2}(T_1-T_s)
\]
where $K_{1/2}$ is the effective conductive coupling of the snow-ice
layer between the surface and the mid-point of the upper layer of ice
$
K_{1/2}=\frac{4 K_i K_s}{K_s h_i + 4 K_i h_s}
$.
$K_i$ and $K_s$ are constant thermal conductivities of seaice and snow.

From the above equations we can develop a system of equations to
find the skin surface temperature, $T_s$ and the two ice layer
temperatures (see Winton, 1999, for details). We solve these
equations iteratively until the change in $T_s$ is small.
When the surface temperature is greater then
the melting temperature of the surface, the temperatures are
recalculated setting $T_s$ to 0.  The enthalpy
of the ice layers are calculated in order to keep track of the energy in the
ice model. Enthalpy is defined, here, as the energy required to melt a
unit mass of seaice with temperature $T$.
For the upper layer (1) with brine pockets  and
the lower fresh layer (2):
\begin{eqnarray}
q_1 & = & - c_f T_f + c_i (T_f-T)+ L_{i}(1-\frac{T_f}{T})
\nonumber \\
q_2 & = & -c_i T+L_i \nonumber
\end{eqnarray}
where $c_f$ is specific heat of liquid fresh water, $c_i$ is the
specific heat of fresh ice, and $L_i$ is latent heat of melting fresh ice.


From the new ice temperatures, we can calculate
the energy flux at the surface available for melting (if $T_s$=0)
and the energy at the ocean-ice interface for either melting or freezing.
\begin{eqnarray}
E_{top} &  =  & (F_s- K_{1/2}(T_s-T_1) ) \Delta t
\nonumber \\
E_{bot} &= & (\frac{4K_i(T_2-T_f)}{h_i}-F_b) \Delta t
\nonumber
\end{eqnarray}
where $F_b$ is the heat flux at the ice bottom due to the sea surface
temperature variations from freezing.
If $T_{sst}$ is above freezing, $F_b=c_{sw} \rho_{sw} 
\gamma (T_{sst}-T_f)u^{*}$, $\gamma$ is the heat transfer coefficient
and $u^{*}=QQ$ is frictional velocity between ice 
and water. If $T_{sst}$ is below freezing,
$F_b=(T_f - T_{sst})c_f \rho_f \Delta z /\Delta t$ and set $T_{sst}$
to $T_f$. We also
include the energy from lower layers that drop below freezing,
and set those layers to $T_f$.

If $E_{top}>0$ we melt snow from the surface, if all the snow is melted
and there is energy left, we melt the ice. If the ice is all gone
and there is still energy left, we apply the left over energy to 
heating the ocean model upper layer (See Winton, 1999, equations 27-29).
Similarly if $E_{bot}>0$ we melt ice from the bottom. If all the ice
is melted, the snow is melted (with energy from the ocean model upper layer
if necessary). If $E_{bot}<0$ we grow ice at the bottom
\[
\Delta h_i = \frac{-E_{bot}}{(q_{bot} \rho_i)}
\]
where $q_{bot}=-c_{i} T_f + L_i$ is the enthalpy of the new ice,
The enthalpy of the second ice layer, $q_2$ needs to be modified:
\[
q_2 = \frac{ \hat{h_i}/2 \hat{q_2} + \Delta h_i q_{bot} }
        {\hat{h_i}/{2}+\Delta h_i}
\]

If there is a ice layer and the overlying air temperature is
below 0$^o$C then any precipitation, $P$ joins the snow layer:
\[
\Delta h_s  = -P \frac{\rho_f}{\rho_s} \Delta t, 
\]
$\rho_f$ and $\rho_s$ are the fresh water and snow densities.
Any evaporation, similarly, removes snow or ice from the surface.
We also calculate the snow age here, in case it is needed for
the surface albedo calculation (see {\it srf\_albedo.F} below).

For practical reasons we limit the ice growth to {\bf hilim}
and snow is limited to {\bf hslim}. We converts any
ice and/or snow above these limits back to water, maintaining the salt
balance. Note however, that heat is not conserved in this
conversion; sea surface temperatures below the ice are not
recalculated.

If the snow/ice interface is below the waterline, snow is converted
to ice (see Winton, 1999, equations 35 and 36). The subroutine
{\it new\_layers\_winton.F}, described below, repartitions the ice into
equal thickness layers while conserving energy.

The subroutine {\it ice\_therm.F} now calculates the heat and fresh
water fluxes affecting the ocean model surface layer. The heat flux:
\[
q_{net}= {\bf fswocn} - F_{b} - \frac{{\bf esurp}}{\Delta t}
\]
is composed of the shortwave flux that has passed through the
ice layer and is absorbed by the water, {\bf fswocn}$=QQ$,
the ocean flux to the ice $F_b$,
and the surplus energy left over from the melting, {\bf esurp}.
The fresh water flux is determined from the amount of
fresh water and salt in the ice/snow system before and after the
timestep.

\noindent
{\tiny (There is a provision for fractional ice:
If ice height is above {\bf hihig} then all energy from freezing at
sea surface is used only in the open water aparts of the cell (ie.
$F_b$ will only have the conduction term).
The melt energy is partitioned by {\bf frac\_energy} between melting
ice height and ice extent. However, once ice height drops below
{\bf himon0} then all energy melts ice extent.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{1cm}

\noindent
{\bf {\underline{subroutine SFC\_ALBEDO} } }

\noindent
The routine {\it ice\_therm.F} calls this routine to determine
the surface albedo. There are two calculations provided here:

\noindent
{\bf 1)} from LANL CICE model
\[ \alpha = f_s \alpha_s + (1-f_s) (\alpha_{i_{min}}
         + (\alpha_{i_{max}}- \alpha_{i_{min}}) (1-e^{-h_i/h_{\alpha}}))
\]
where $f_s$ is 1 if there is snow, 0 if not; the snow albedo, 
$\alpha_s$ has two values
depending on whether $T_s<0$ or not; $\alpha_{i_{min}}$ and 
$\alpha_{i_{max}}$ are ice albedos for thin melting ice, and
thick bare ice respectively, and $h_{\alpha}$ is a scale
height.

\noindent
{\bf 2)} From GISS model (Hansen et al 1983)
\[
 \alpha = \alpha_i e^{-h_s/h_a} + \alpha_s (1-e^{-h_s/h_a})
\]
where $\alpha_i$ is a constant albedo for bare ice, $h_a$
is a scale height and $\alpha_s$ is a variable snow albedo.
\[
\alpha_s = \alpha_1 + \alpha_2 e^{-\lambda_a a_s}
\]
where $\alpha_1$ is a constant, $\alpha_2$ depends on $T_s$,
$a_s$ is the snow age, and $\lambda_a$ is a scale frequency.
The snow age is calculated in {\it ice\_therm.F} and is given
in equation 41 in Hansen et al (1983).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{1cm}

\noindent
{\bf {\underline{subroutine NEW\_LAYERS\_WINTON}}}

\noindent
The subroutine
{\it new\_layers\_winton.F} repartitions the ice into
equal thickness layers while conserving energy. We pass
to this subroutine, the ice layer enthalpies after
melting/growth and the new height of the ice layers.
The ending layer height should be half the sum of the
new ice heights from {\it ice\_therm.F}. The enthalpies
of the ice layers are adjusted accordingly to maintain
total energy in the ice model. If layer 2 height is
greater than layer 1 height then layer 2 gives ice to
layer 1 and:
\[
q_1=f_1 \hat{q_1} + (1-f1) \hat{q_2}
\]
where $f_1$ is the fraction of the new to old upper layer heights.
$T_1$ will therefore also have changed.
Similarly for when ice layer height 2 is less than
layer 1 height, except here we need to to be careful
that the new $T_2$ does not fall below the melting temperature.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{1cm}

\noindent
{\bf {\underline{Initializing subroutines}}}

\noindent
{\it ice\_init.F}:
Set ice variables to zero, or reads in pickup information
from {\bf pickup.ic} (which was written out in {\it checkpoint.F})

\noindent
{\it ice\_readparms.F}:
Reads {\bf data.ice}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{1cm}

\noindent
{\bf {\underline{Diagnostic subroutines}}}

\noindent
{\it ice\_ave.F}:
Keeps track of means of the ice variables

\noindent
{\it ice\_diags.F}:
Finds averages and writes out diagnostics

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{1cm}

\noindent
{\bf {\underline{Common Blocks}}}

\noindent
{\it ICE.h}: Ice Varibles, also 
{\bf relaxlat} and {\bf startIceModel}

\noindent
{\it ICE\_DIAGS.h}: matrices for diagnostics: averages of fields
from {\it ice\_diags.F}

\noindent
{\it BULKF\_ICE\_CONSTANTS.h} (in {\bf BULKF} package): 
all the parameters need by the ice model

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{1cm}

\noindent
{\bf {\underline{Input file DATA.ICE}}}

\noindent
Here we need to set {\bf StartIceModel}: which is 1 if the
model starts from no ice; and 0 if there is a pickup file
with the ice matrices ({\bf pickup.ic}) which is read
in {\it ice\_init.F} and written out in {\it checkpoint.F}.
The parameter {\bf relaxlat} defines the latitude poleward
of which there is no relaxing of surface $T$ or $S$ to
observations. This avoids the relaxation forcing the ice
model at these high latitudes.

\noindent
({\tiny Note: {\bf hicemin} is set to 0 here. If the
provision for allowing grid cells to have both
open water and seaice is ever implemented, this would
be greater than 0})

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{1cm}

\noindent
{\bf {\underline{Important Notes}}}

\noindent
{\bf 1)} heat fluxes have different signs in the ocean and ice
models.

\noindent
{\bf 2)} {\bf StartIceModel} must be changed in {\bf data.ice}:
1 (if starting from no ice), 0 (if using pickup.ic file).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\vspace{1cm}

\noindent 
{\bf {\underline{References}}}

\noindent
Bitz, C.M. and W.H. Lipscombe, 1999: An Energy-Conserving
Thermodynamic Model of Sea Ice.
{\it Journal of Geophysical Research}, 104, 15,669 -- 15,677.

\vspace{.2cm}

\noindent
Hansen, J., G. Russell, D. Rind, P. Stone, A. Lacis, S. Lebedeff,
R. Ruedy and L.Travis, 1983: Efficient Three-Dimensional
Global Models for Climate Studies: Models I and II.
{\it Monthly Weather Review}, 111, 609 -- 662.

\vspace{.2cm}

\noindent
Hunke, E.C and W.H. Lipscomb, circa 2001: CICE: the Los Alamos
Sea Ice Model Documentation and Software User's Manual.
LACC-98-16v.2.\\
(note: this documentation is no longer available as CICE has progressed
to a very different version 3)


\vspace{.2cm}

\noindent
Winton, M, 2000: A reformulated Three-layer Sea Ice Model.
{\it Journal of Atmospheric and Ocean Technology}, 17, 525 -- 531.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \end{document}
1	jmc	1.1	% \documentclass[12pt]{article}
2			% \usepackage{amssymb}
3
4			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5			%% \usepackage{graphics}
6
7
8			% \oddsidemargin -4mm \evensidemargin 0mm
9			% \textwidth 165mm
10			% \textheight 230mm
11			% \topmargin -2mm \headsep -2mm
12			% \renewcommand{\baselinestretch}{1.5}
13			% \begin{document}
14
15
16			\def\deg{$^o$}
17			%%%--------------------------------------%%%
18	molod	1.5	\subsection{THSICE: The Thermodynamic Sea Ice Package}
19	edhill	1.3	\label{sec:pkg:thsice}
20			\begin{rawhtml}
21			<!-- CMIREDIR:package_thsice: -->
22			\end{rawhtml}
23	jmc	1.1
24			{\bf Important note:}
25			This document has been written by Stephanie Dutkiewicz
26			and describes an earlier implementation of the sea-ice package.
27			This needs to be updated to reflect the recent changes (JMC).
28
29			\noindent
30			This thermodynamic ice model is based on the 3-layer model by Winton (2000).
31			and the energy-conserving LANL CICE model (Bitz and Lipscomb, 1999).
32			The model considers two equally thick ice layers; the upper layer has
33			a variable specific heat resulting from brine pockets,
34			the lower layer has a fixed heat capacity. A zero heat capacity snow
35			layer lies above the ice. Heat fluxes at the top and bottom
36			surfaces are used to calculate the change in ice and snow layer
37			thickness. Grid cells of the ocean model are
38			either fully covered in ice or are open water. There is
39			a provision to parametrize ice fraction (and leads) in this package.
40			Modifications are discussed in small font following the
41			subroutine descriptions.
42
43			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44
45			\vspace{1cm}
46
47			\noindent
48			The ice model is called from {\it thermodynamics.F}, subroutine
49			{\it ice\_forcing.F} is called in place of {\it external\_forcing\_surf.F}.
50
51			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52
53			\vspace{1cm}
54			\noindent
55			{\bf \underline{subroutine ICE\_FORCING}}
56
57			\noindent
58			In {\it ice\_forcing.F}, we calculate the freezing potential of the
59			ocean model surface layer of water:
60			\[
61			{\bf frzmlt} = (T_f - SST) \frac{c_{sw} \rho_{sw} \Delta z}{\Delta t}
62			\]
63			where $c_{sw}$ is seawater heat capacity,
64			$\rho_{sw}$ is the seawater density, $\Delta z$
65			is the ocean model upper layer thickness and $\Delta t$ is the model (tracer)
66			timestep. The freezing temperature, $T_f=\mu S$ is a function of the
67			salinity.
68
69
70			1) Provided there is no ice present and {\bf frzmlt} is less than 0,
71			the surface tendencies of wind, heat and freshwater are calculated
72			as usual (ie. as in {\it external\_forcing\_surf.F}).
73
74			2) If there is ice present in the grid cell
75			we call the main ice model routine {\it ice\_therm.F} (see below).
76			Output from this routine gives net heat and freshwater flux
77			affecting the top of the ocean.
78
79			Subroutine {\it ice\_forcing.F} uses these values to find the
80			sea surface tendencies
81			in grid cells. When there is ice present,
82			the surface stress tendencies are
83			set to zero; the ice model is purely thermodynamic and the
84			effect of ice motion on the sea-surface is not examined.
85
86			Relaxation of surface $T$ and $S$ is only allowed equatorward
87			of {\bf relaxlat} (see {\bf DATA.ICE below}), and no relaxation is
88			allowed under the ice at any latitude.
89
90			\noindent
91			{\tiny (Note that there is provision for allowing grid cells to have both
92			open water and seaice; if {\bf compact} is between 0 and 1)}
93
94			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95			\vspace{1cm}
96			\noindent
97			{\bf {\underline{ subroutine ICE\_FREEZE}}}
98
99			This routine is called from {\it thermodynamics.F}
100			after the new temperature calculation, {\it calc\_gt.F},
101			but before {\it calc\_gs.F}.
102			In {\it ice\_freeze.F}, any ocean upper layer grid cell
103			with no ice cover, but with temperature below freezing,
104			$T_f=\mu S$ has ice initialized.
105			We calculate {\bf frzmlt} from all the grid cells in
106			the water column that have a temperature less than
107			freezing. In this routine, any water below the surface
108			that is below freezing is set to $T_f$.
109			A call to
110			{\it ice\_start.F} is made if {\bf frzmlt} $>0$,
111			and salinity tendancy is updated for brine release.
112
113			\noindent
114			{\tiny (There is a provision for fractional ice:
115			In the case where the grid cell has less ice coverage than
116			{\bf icemaskmax} we allow {\it ice\_start.F} to be called).}
117
118			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119
120			\vspace{1cm}
121			\noindent
122			{\bf {\underline{ subroutine ICE\_START}}}
123
124			\noindent
125			The energy available from freezing
126			the sea surface is brought into this routine as {\bf esurp}.
127			The enthalpy of the 2 layers of any new ice is calculated as:
128			\begin{eqnarray}
129			q_1 & = & -c_{i}*T_f + L_i \nonumber \\
130			q_2 & = & -c_{f}T_{mlt}+ c_{i}(T_{mlt}-T{f}) + L_i(1-\frac{T_{mlt}}{T_f}
131			\nonumber \\
132			\end{eqnarray}
133			where $c_f$ is specific heat of liquid fresh water, $c_i$ is the
134			specific heat of fresh ice, $L_i$ is latent heat of freezing,
135			$\rho_i$ is density of ice and
136			$T_{mlt}$ is melting temperature of ice with salinity of 1.
137			The height of a new layer of ice is
138			\[
139			h_{i new} = \frac{{\bf esurp} \Delta t}{qi_{0av}}
140			\]
141			where $qi_{0av}=-\frac{\rho_i}{2} (q_1+q_2)$.
142
143			The surface skin temperature $T_s$ and ice temperatures
144			$T_1$, $T_2$ and the sea surface temperature are set at $T_f$.
145
146			\noindent
147			{\tiny ( There is provision for fractional ice:
148			new ice is formed over open water; the first freezing in the cell
149			must have a height of {\bf himin0}; this determines the ice
150			fraction {\bf compact}. If there is already ice in the grid cell,
151			the new ice must have the same height and the new ice fraction
152			is
153			\[
154			i_f=(1-\hat{i_f}) \frac{h_{i new}}{h_i}
155			\]
156			where $\hat{i_f}$ is ice fraction from previous timestep
157			and $h_i$ is current ice height. Snow is redistributed
158			over the new ice fraction. The ice fraction is
159			not allowed to become larger than {\bf iceMaskmax} and
160			if the ice height is above {\bf hihig} then freezing energy
161			comes from the full grid cell, ice growth does not occur
162			under orginal ice due to freezing water.
163			}
164			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
165
166			\vspace{1cm}
167			\noindent
168			{\bf {\underline{subroutine ICE\_THERM}}}
169
170			\noindent
171			The main subroutine of this package is {\it ice\_therm.F} where the
172			ice temperatures are calculated and the changes in ice and snow
173			thicknesses are determined. Output provides the net heat and fresh
174			water fluxes that force the top layer of the ocean model.
175
176			If the current ice height is less than {\bf himin} then
177			the ice layer is set to zero and the ocean model upper layer temperature
178			is allowed to drop lower than its freezing temperature; and atmospheric
179			fluxes are allowed to effect the grid cell.
180			If the ice height is greater than {\bf himin} we proceed with
181			the ice model calculation.
182
183			We follow the procedure
184			of Winton (1999) -- see equations 3 to 21 -- to calculate
185			the surface and internal ice temperatures.
186			The surface temperature is found from the balance of the
187			flux at the surface $F_s$, the shortwave heat flux absorbed by the ice,
188			{\bf fswint}, and
189			the upward conduction of heat through the snow and/or ice, $F_u$.
190			We linearize $F_s$ about the surface temperature, $\hat{T_s}$,
191			at the previous timestep (where \mbox{}$\hat{ }$ indicates the value at
192			the previous timestep):
193			\[
194			F_s (T_s) = F_s(\hat{T_s}) + \frac{\partial F_s(\hat{T_s)}}{\partial T_s}
195			(T_s-\hat{T_s})
196			\]
197			where,
198			\[
199			F_s = F_{sensible}+F_{latent}+F_{longwave}^{down}+F_{longwave}^{up}+ (1-
200			\alpha) F_{shortwave}
201			\]
202			and
203			\[
204			\frac{d F_s}{dT} = \frac{d F_{sensible}}{dT} + \frac{d F_{latent}}{dT}
205			+\frac{d F_{longwave}^{up}}{dT}.
206			\]
207			$F_s$ and $\frac{d F_s}{dT}$ are currently calculated from the {\bf BULKF}
208			package described separately, but could also be provided by an atmospheric
209			model. The surface albedo is calculated from the ice height and/or
210			surface temperature (see below, {\it srf\_albedo.F}) and the
211			shortwave flux absorbed in the ice is
212			\[
213			{\bf fswint} = (1-e^{\kappa_i h_i})(1-\alpha) F_{shortwave}
214			\]
215			where $\kappa_i$ is bulk extinction coefficient.
216
217			The conductive flux to the surface is
218			\[
219			F_u=K_{1/2}(T_1-T_s)
220			\]
221			where $K_{1/2}$ is the effective conductive coupling of the snow-ice
222			layer between the surface and the mid-point of the upper layer of ice
223			$
224			K_{1/2}=\frac{4 K_i K_s}{K_s h_i + 4 K_i h_s}
225			$.
226			$K_i$ and $K_s$ are constant thermal conductivities of seaice and snow.
227
228			From the above equations we can develop a system of equations to
229			find the skin surface temperature, $T_s$ and the two ice layer
230			temperatures (see Winton, 1999, for details). We solve these
231			equations iteratively until the change in $T_s$ is small.
232			When the surface temperature is greater then
233			the melting temperature of the surface, the temperatures are
234			recalculated setting $T_s$ to 0. The enthalpy
235			of the ice layers are calculated in order to keep track of the energy in the
236			ice model. Enthalpy is defined, here, as the energy required to melt a
237			unit mass of seaice with temperature $T$.
238			For the upper layer (1) with brine pockets and
239			the lower fresh layer (2):
240			\begin{eqnarray}
241			q_1 & = & - c_f T_f + c_i (T_f-T)+ L_{i}(1-\frac{T_f}{T})
242			\nonumber \\
243			q_2 & = & -c_i T+L_i \nonumber
244			\end{eqnarray}
245			where $c_f$ is specific heat of liquid fresh water, $c_i$ is the
246			specific heat of fresh ice, and $L_i$ is latent heat of melting fresh ice.
247
248
249
250			From the new ice temperatures, we can calculate
251			the energy flux at the surface available for melting (if $T_s$=0)
252			and the energy at the ocean-ice interface for either melting or freezing.
253			\begin{eqnarray}
254			E_{top} & = & (F_s- K_{1/2}(T_s-T_1) ) \Delta t
255			\nonumber \\
256			E_{bot} &= & (\frac{4K_i(T_2-T_f)}{h_i}-F_b) \Delta t
257			\nonumber
258			\end{eqnarray}
259			where $F_b$ is the heat flux at the ice bottom due to the sea surface
260			temperature variations from freezing.
261			If $T_{sst}$ is above freezing, $F_b=c_{sw} \rho_{sw}
262			\gamma (T_{sst}-T_f)u^{*}$, $\gamma$ is the heat transfer coefficient
263			and $u^{*}=QQ$ is frictional velocity between ice
264			and water. If $T_{sst}$ is below freezing,
265			$F_b=(T_f - T_{sst})c_f \rho_f \Delta z /\Delta t$ and set $T_{sst}$
266			to $T_f$. We also
267			include the energy from lower layers that drop below freezing,
268			and set those layers to $T_f$.
269
270			If $E_{top}>0$ we melt snow from the surface, if all the snow is melted
271			and there is energy left, we melt the ice. If the ice is all gone
272			and there is still energy left, we apply the left over energy to
273			heating the ocean model upper layer (See Winton, 1999, equations 27-29).
274			Similarly if $E_{bot}>0$ we melt ice from the bottom. If all the ice
275			is melted, the snow is melted (with energy from the ocean model upper layer
276			if necessary). If $E_{bot}<0$ we grow ice at the bottom
277			\[
278			\Delta h_i = \frac{-E_{bot}}{(q_{bot} \rho_i)}
279			\]
280			where $q_{bot}=-c_{i} T_f + L_i$ is the enthalpy of the new ice,
281			The enthalpy of the second ice layer, $q_2$ needs to be modified:
282			\[
283			q_2 = \frac{ \hat{h_i}/2 \hat{q_2} + \Delta h_i q_{bot} }
284			{\hat{h_i}/{2}+\Delta h_i}
285			\]
286
287			If there is a ice layer and the overlying air temperature is
288			below 0$^o$C then any precipitation, $P$ joins the snow layer:
289			\[
290			\Delta h_s = -P \frac{\rho_f}{\rho_s} \Delta t,
291			\]
292			$\rho_f$ and $\rho_s$ are the fresh water and snow densities.
293			Any evaporation, similarly, removes snow or ice from the surface.
294			We also calculate the snow age here, in case it is needed for
295			the surface albedo calculation (see {\it srf\_albedo.F} below).
296
297			For practical reasons we limit the ice growth to {\bf hilim}
298			and snow is limited to {\bf hslim}. We converts any
299			ice and/or snow above these limits back to water, maintaining the salt
300			balance. Note however, that heat is not conserved in this
301			conversion; sea surface temperatures below the ice are not
302			recalculated.
303
304			If the snow/ice interface is below the waterline, snow is converted
305			to ice (see Winton, 1999, equations 35 and 36). The subroutine
306			{\it new\_layers\_winton.F}, described below, repartitions the ice into
307			equal thickness layers while conserving energy.
308
309			The subroutine {\it ice\_therm.F} now calculates the heat and fresh
310			water fluxes affecting the ocean model surface layer. The heat flux:
311			\[
312			q_{net}= {\bf fswocn} - F_{b} - \frac{{\bf esurp}}{\Delta t}
313			\]
314			is composed of the shortwave flux that has passed through the
315			ice layer and is absorbed by the water, {\bf fswocn}$=QQ$,
316			the ocean flux to the ice $F_b$,
317			and the surplus energy left over from the melting, {\bf esurp}.
318			The fresh water flux is determined from the amount of
319			fresh water and salt in the ice/snow system before and after the
320			timestep.
321
322			\noindent
323			{\tiny (There is a provision for fractional ice:
324			If ice height is above {\bf hihig} then all energy from freezing at
325			sea surface is used only in the open water aparts of the cell (ie.
326			$F_b$ will only have the conduction term).
327			The melt energy is partitioned by {\bf frac\_energy} between melting
328			ice height and ice extent. However, once ice height drops below
329			{\bf himon0} then all energy melts ice extent.}
330
331			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
332			\vspace{1cm}
333
334			\noindent
335			{\bf {\underline{subroutine SFC\_ALBEDO} } }
336
337			\noindent
338			The routine {\it ice\_therm.F} calls this routine to determine
339			the surface albedo. There are two calculations provided here:
340
341			\noindent
342			{\bf 1)} from LANL CICE model
343			\[ \alpha = f_s \alpha_s + (1-f_s) (\alpha_{i_{min}}
344			+ (\alpha_{i_{max}}- \alpha_{i_{min}}) (1-e^{-h_i/h_{\alpha}}))
345			\]
346			where $f_s$ is 1 if there is snow, 0 if not; the snow albedo,
347			$\alpha_s$ has two values
348			depending on whether $T_s<0$ or not; $\alpha_{i_{min}}$ and
349			$\alpha_{i_{max}}$ are ice albedos for thin melting ice, and
350			thick bare ice respectively, and $h_{\alpha}$ is a scale
351			height.
352
353			\noindent
354			{\bf 2)} From GISS model (Hansen et al 1983)
355			\[
356			\alpha = \alpha_i e^{-h_s/h_a} + \alpha_s (1-e^{-h_s/h_a})
357			\]
358			where $\alpha_i$ is a constant albedo for bare ice, $h_a$
359			is a scale height and $\alpha_s$ is a variable snow albedo.
360			\[
361			\alpha_s = \alpha_1 + \alpha_2 e^{-\lambda_a a_s}
362			\]
363			where $\alpha_1$ is a constant, $\alpha_2$ depends on $T_s$,
364			$a_s$ is the snow age, and $\lambda_a$ is a scale frequency.
365			The snow age is calculated in {\it ice\_therm.F} and is given
366			in equation 41 in Hansen et al (1983).
367
368			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
369
370			\vspace{1cm}
371
372			\noindent
373			{\bf {\underline{subroutine NEW\_LAYERS\_WINTON}}}
374
375			\noindent
376			The subroutine
377			{\it new\_layers\_winton.F} repartitions the ice into
378			equal thickness layers while conserving energy. We pass
379			to this subroutine, the ice layer enthalpies after
380			melting/growth and the new height of the ice layers.
381			The ending layer height should be half the sum of the
382			new ice heights from {\it ice\_therm.F}. The enthalpies
383			of the ice layers are adjusted accordingly to maintain
384			total energy in the ice model. If layer 2 height is
385			greater than layer 1 height then layer 2 gives ice to
386			layer 1 and:
387			\[
388			q_1=f_1 \hat{q_1} + (1-f1) \hat{q_2}
389			\]
390			where $f_1$ is the fraction of the new to old upper layer heights.
391			$T_1$ will therefore also have changed.
392			Similarly for when ice layer height 2 is less than
393			layer 1 height, except here we need to to be careful
394			that the new $T_2$ does not fall below the melting temperature.
395
396			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
397
398			\vspace{1cm}
399
400			\noindent
401			{\bf {\underline{Initializing subroutines}}}
402
403			\noindent
404			{\it ice\_init.F}:
405			Set ice variables to zero, or reads in pickup information
406			from {\bf pickup.ic} (which was written out in {\it checkpoint.F})
407
408			\noindent
409			{\it ice\_readparms.F}:
410			Reads {\bf data.ice}
411
412			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
413
414			\vspace{1cm}
415
416			\noindent
417			{\bf {\underline{Diagnostic subroutines}}}
418
419			\noindent
420			{\it ice\_ave.F}:
421			Keeps track of means of the ice variables
422
423			\noindent
424			{\it ice\_diags.F}:
425			Finds averages and writes out diagnostics
426
427			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
428			\vspace{1cm}
429
430			\noindent
431			{\bf {\underline{Common Blocks}}}
432
433			\noindent
434			{\it ICE.h}: Ice Varibles, also
435			{\bf relaxlat} and {\bf startIceModel}
436
437			\noindent
438			{\it ICE\_DIAGS.h}: matrices for diagnostics: averages of fields
439			from {\it ice\_diags.F}
440
441			\noindent
442			{\it BULKF\_ICE\_CONSTANTS.h} (in {\bf BULKF} package):
443			all the parameters need by the ice model
444
445			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
446			\vspace{1cm}
447
448			\noindent
449			{\bf {\underline{Input file DATA.ICE}}}
450
451			\noindent
452			Here we need to set {\bf StartIceModel}: which is 1 if the
453			model starts from no ice; and 0 if there is a pickup file
454			with the ice matrices ({\bf pickup.ic}) which is read
455			in {\it ice\_init.F} and written out in {\it checkpoint.F}.
456			The parameter {\bf relaxlat} defines the latitude poleward
457			of which there is no relaxing of surface $T$ or $S$ to
458			observations. This avoids the relaxation forcing the ice
459			model at these high latitudes.
460
461			\noindent
462			({\tiny Note: {\bf hicemin} is set to 0 here. If the
463			provision for allowing grid cells to have both
464			open water and seaice is ever implemented, this would
465			be greater than 0})
466
467			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
468			\vspace{1cm}
469
470			\noindent
471			{\bf {\underline{Important Notes}}}
472
473			\noindent
474			{\bf 1)} heat fluxes have different signs in the ocean and ice
475			models.
476
477			\noindent
478			{\bf 2)} {\bf StartIceModel} must be changed in {\bf data.ice}:
479			1 (if starting from no ice), 0 (if using pickup.ic file).
480
481			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
482
483			\vspace{1cm}
484
485			\noindent
486			{\bf {\underline{References}}}
487
488			\noindent
489			Bitz, C.M. and W.H. Lipscombe, 1999: An Energy-Conserving
490			Thermodynamic Model of Sea Ice.
491			{\it Journal of Geophysical Research}, 104, 15,669 -- 15,677.
492
493			\vspace{.2cm}
494
495			\noindent
496			Hansen, J., G. Russell, D. Rind, P. Stone, A. Lacis, S. Lebedeff,
497			R. Ruedy and L.Travis, 1983: Efficient Three-Dimensional
498			Global Models for Climate Studies: Models I and II.
499			{\it Monthly Weather Review}, 111, 609 -- 662.
500
501			\vspace{.2cm}
502
503			\noindent
504			Hunke, E.C and W.H. Lipscomb, circa 2001: CICE: the Los Alamos
505			Sea Ice Model Documentation and Software User's Manual.
506			LACC-98-16v.2.\\
507			(note: this documentation is no longer available as CICE has progressed
508			to a very different version 3)
509
510
511			\vspace{.2cm}
512
513			\noindent
514			Winton, M, 2000: A reformulated Three-layer Sea Ice Model.
515			{\it Journal of Atmospheric and Ocean Technology}, 17, 525 -- 531.
516
517
518
519			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
520			% \end{document}