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$}
%%%--------------------------------------%%%
\section{Thermodynamic Sea Ice Package: ``thsice''}

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