ifenty/Fenty_Thermo_Code_Updates/README.txt

I went through the seaice_*_if.F codes in the repository and made some updates to seaice_budget_ice_if.F and seaice_growth_if.F, brought the code into spec as per Jean-Michel's suggestions (remove unnecessary overlaps, send myIter and myThid to seaice_budget_ice, removed seaice_budget_ocean_if.F, etc.), and clarified most of the comments/in line documentation.  

In addition, I configured two verification-type setups.  The first is a modification of the existing verification/lab sea setup and the other is basically a 1-D column.  In both one can run forward and adjoint calculations to compare 

If you have the time, please look at the code updates or try them in one of your setups.  I would like to know if you find any problems.

The verification experiments and code updates can be found in the /MITgcm_contrib/ifenty/Fenty_Thermo_Code_Updates

Code updates (which include changes to SEAICE_PARAM.h, seaice_readparams.F, and seaice_summary.F), are found in the directory code_updates.

===== Verification Experiments =====

The lab sea verification setup is found in code_20x16_verification_experiment and input_20x16_verification_experiment while the the 1d colum is in code_1D_verification_experiment and input_1D_verification_experiment.  In addition, there are code_XXX_variation_1 in which the default thermodynamic codes are used instead.

The 1D setup uses a vertical T,S profile from the central Arctic that An made and atmospheric forcing which is similar to the Labrador Sea.  Using the default thermo code, the adjoint NaN's in a run over 11,000 time steps.

Since the handling of salt in the salt plume is somewhat different between the default and _if codes, useSalt_PLUME is set to false by default in data.pkg to faciliate comparison between the _if and default codes.
Note: regardless of which thermodynamic ice codes are used, the adjoint run loops enlessly unless SEAICE_ALLOW_DYNAMICS is defined, for reasons which are not clear to me.


===== Major code changes/thoughts =====

1) Ocean-ice heat fluxes
2) salt plume
3) ice salinity
4) ice age

=== Ocean-ice heat fluxes ===

Now the  model can use two alternative schemes to 
calculate turbulent ocean-ice heat fluxes.  
Each is described in turn

= METHOD 1 =

The first uses an empircal relationship between the
'far-field' ocean temperature and ocean to sea ice turbulent heat fluxes
given by McPhee:

Flux  = rho_sw * cp_sw * <w'T'>
<w'T'>= S_t *  u_star * (T_o - T_f)

Where,
rho_sw : seawater density (kg m^-3)
cp_sw  : seawater heat capcity (J kg^-1 K^-1)
S_t    : Stanton Number ~ 0.0056 (dimensionless)
u_star : friction velocity beneath ice ~ 0.015 m s^-1
T_o, T_f : The 'far-field' ocean temperature and ice
freezing point (deg C)

Using typical values, ice advected over waters warmer by 1 degree
will be subjected to a basal heat flux of ~ 345 W m^-2.

In the model, the criteria for ice growth is that the air-sea
heat loss must exceed the potential ocean-ice heat flux (i.e., more
potential ice production over one time step than melt).
Using the above parameterization, ice will form in waters which
are much warmer than its salinity-determined freezing point
using typical wintertime conditions in, say, the North Atlantic.

Therefore, when no ice is present, an additional factor is applied
to the above <w'T'> (mixedLayerTurbulenceFactor) to represent the idea that
turbulent mixing at and near the surface of an ice-free ocean
is much greater than mixing beneath an ice-covered one.

A factor of 12.5 is chosen for mixedLayerTurbulenceFactor.  Consequently,
for each 0.1 degree above freezing of the upper ocean grid cell,
the air-sea heat losses must exceed ~ 515 W m^-2 before ice forms.

Once ice gains a foothold in a grid cell, the McPhee parameterization
is immediately invoked.

= METHOD 2 =

The second scheme (the original version)
subsumes (S_t * u_star) into a single variable (SEAICE_gamma_t)
in the following way:

Flux  = rho_sw * cp_sw * GAMMA
GAMMA = dRf(1)/SEAICE_gamma_t * (T_o - T_f)

Where,
dRf(1)  : depth of surface level grid cell (originally
assumed to be 10 m)
SEAICE_gamma_t : an empirical parameter (~ 3 days in seconds)

Using typical values, ice advected into a grid cell
that is warmer than the sea ice freezing point by 1 degree
will be subjected to a basal heat flux of ~ 160 W m^-2.   Therefore,
as written, one should choose a SEAICE_gamma_t = 1.4 days in seconds
to have flux magnitudes match those of the McPhee parameterization.  
Of course, even with this shorter SEAICE_gamma_t, ice will form
in grid cells which are far above freezing (> 1 degree) when
subjected to air-sea heat fluxes of ~ 345 W m^-2.  (To get around
this I choose SEAICE_gamma_t to be ~ 10^4 seconds for my thesis, 
but I really should have just added a mixedLayerTurbulenceFactor 
instead.)

The original scheme is accessed by defining  #SEAICE_USE_ORIGINAL_HEAT FLUX


=== Seaice salinity ===

At present, melting and growing ice is assumed to have a 
constant salinity specified by SEAICE_salinity_fixed.  I plan to 
put in Dimitris' variable salt fraction and HSALT to store 
the total mass of salt stored within each grid cell's ice, but
it would take some more thinking on my part to be sure that I 
conserve salt correctly.

=== Salt Plumes ===

I assume that ice production in leads can generate
salt plumes.  The fraction of the salt sent to the plume package
from ice produced in leads (defined by the open water
fraction) is a nonlinear function of ice concentration, AREA.

Specifically, the function is a logistic curve (sigmoid) with a range and
domain {0,1}.  The function, f(AREA), has a single free parameter,
SEAICE_plumeInflectionPoint, the inflection point of the curve.

By construction, the function has the following properties:
f(1) \approx 1.0
f(SEAICE_plumeInflectionPoint) = 0.5
f(0) \approx 0.0 (when SEAICE_plumeInflectionPoint \geq 0.5)
f(0) > 0.0 (when SEAICE_plumeInflectionPoint < 0.5)

As AREA --> 1, the open water fraction is confined to
narrow leads implying that new ice production is spatially non-uniform,
therby violating the assumption of spatially uniform forcing over
the grid cell required by KPP.

To treat the vertical overturning of the high-salinity brine
in a more physically realistic way, the salt produced
in the leads is sent to depth via the plume package.  To assure
only narrow leads generate plumes, choose a relatively high
SEAICE_plumeInflectionPoint.  For testing purposes, I chose 0.8
but I don't know what the "right" value should be.

=== Ice Age ===

Adding ice age means just adding a short snippet at the end
of seaice_growth.F.  I propose that adding to ice age be done
in a separate subroutine which takes as arguments, AREA, AREANm1, 
HEFF, HEFFNm1, and whatever other variables from seaice_growth 
it needs.  Interfacing with a separate subroutine will make
for a cleaner seaice growth code in the long run.
1	ifenty	1.1	I went through the seaice_*_if.F codes in the repository and made some updates to seaice_budget_ice_if.F and seaice_growth_if.F, brought the code into spec as per Jean-Michel's suggestions (remove unnecessary overlaps, send myIter and myThid to seaice_budget_ice, removed seaice_budget_ocean_if.F, etc.), and clarified most of the comments/in line documentation.
2
3			In addition, I configured two verification-type setups. The first is a modification of the existing verification/lab sea setup and the other is basically a 1-D column. In both one can run forward and adjoint calculations to compare
4
5			If you have the time, please look at the code updates or try them in one of your setups. I would like to know if you find any problems.
6
7			The verification experiments and code updates can be found in the /MITgcm_contrib/ifenty/Fenty_Thermo_Code_Updates
8
9			Code updates (which include changes to SEAICE_PARAM.h, seaice_readparams.F, and seaice_summary.F), are found in the directory code_updates.
10
11			===== Verification Experiments =====
12
13			The lab sea verification setup is found in code_20x16_verification_experiment and input_20x16_verification_experiment while the the 1d colum is in code_1D_verification_experiment and input_1D_verification_experiment. In addition, there are code_XXX_variation_1 in which the default thermodynamic codes are used instead.
14
15			The 1D setup uses a vertical T,S profile from the central Arctic that An made and atmospheric forcing which is similar to the Labrador Sea. Using the default thermo code, the adjoint NaN's in a run over 11,000 time steps.
16
17			Since the handling of salt in the salt plume is somewhat different between the default and _if codes, useSalt_PLUME is set to false by default in data.pkg to faciliate comparison between the _if and default codes.
18			Note: regardless of which thermodynamic ice codes are used, the adjoint run loops enlessly unless SEAICE_ALLOW_DYNAMICS is defined, for reasons which are not clear to me.
19
20
21			===== Major code changes/thoughts =====
22
23			1) Ocean-ice heat fluxes
24			2) salt plume
25			3) ice salinity
26			4) ice age
27
28			=== Ocean-ice heat fluxes ===
29
30			Now the model can use two alternative schemes to
31			calculate turbulent ocean-ice heat fluxes.
32			Each is described in turn
33
34			= METHOD 1 =
35
36			The first uses an empircal relationship between the
37			'far-field' ocean temperature and ocean to sea ice turbulent heat fluxes
38			given by McPhee:
39
40			Flux = rho_sw * cp_sw * <w'T'>
41			<w'T'>= S_t * u_star * (T_o - T_f)
42
43			Where,
44			rho_sw : seawater density (kg m^-3)
45			cp_sw : seawater heat capcity (J kg^-1 K^-1)
46			S_t : Stanton Number ~ 0.0056 (dimensionless)
47			u_star : friction velocity beneath ice ~ 0.015 m s^-1
48			T_o, T_f : The 'far-field' ocean temperature and ice
49			freezing point (deg C)
50
51			Using typical values, ice advected over waters warmer by 1 degree
52			will be subjected to a basal heat flux of ~ 345 W m^-2.
53
54			In the model, the criteria for ice growth is that the air-sea
55			heat loss must exceed the potential ocean-ice heat flux (i.e., more
56			potential ice production over one time step than melt).
57			Using the above parameterization, ice will form in waters which
58			are much warmer than its salinity-determined freezing point
59			using typical wintertime conditions in, say, the North Atlantic.
60
61			Therefore, when no ice is present, an additional factor is applied
62			to the above <w'T'> (mixedLayerTurbulenceFactor) to represent the idea that
63			turbulent mixing at and near the surface of an ice-free ocean
64			is much greater than mixing beneath an ice-covered one.
65
66			A factor of 12.5 is chosen for mixedLayerTurbulenceFactor. Consequently,
67			for each 0.1 degree above freezing of the upper ocean grid cell,
68			the air-sea heat losses must exceed ~ 515 W m^-2 before ice forms.
69
70			Once ice gains a foothold in a grid cell, the McPhee parameterization
71			is immediately invoked.
72
73			= METHOD 2 =
74
75			The second scheme (the original version)
76			subsumes (S_t * u_star) into a single variable (SEAICE_gamma_t)
77			in the following way:
78
79			Flux = rho_sw * cp_sw * GAMMA
80			GAMMA = dRf(1)/SEAICE_gamma_t * (T_o - T_f)
81
82			Where,
83			dRf(1) : depth of surface level grid cell (originally
84			assumed to be 10 m)
85			SEAICE_gamma_t : an empirical parameter (~ 3 days in seconds)
86
87			Using typical values, ice advected into a grid cell
88			that is warmer than the sea ice freezing point by 1 degree
89			will be subjected to a basal heat flux of ~ 160 W m^-2. Therefore,
90			as written, one should choose a SEAICE_gamma_t = 1.4 days in seconds
91			to have flux magnitudes match those of the McPhee parameterization.
92			Of course, even with this shorter SEAICE_gamma_t, ice will form
93			in grid cells which are far above freezing (> 1 degree) when
94			subjected to air-sea heat fluxes of ~ 345 W m^-2. (To get around
95			this I choose SEAICE_gamma_t to be ~ 10^4 seconds for my thesis,
96			but I really should have just added a mixedLayerTurbulenceFactor
97			instead.)
98
99			The original scheme is accessed by defining #SEAICE_USE_ORIGINAL_HEAT FLUX
100
101
102			=== Seaice salinity ===
103
104			At present, melting and growing ice is assumed to have a
105			constant salinity specified by SEAICE_salinity_fixed. I plan to
106			put in Dimitris' variable salt fraction and HSALT to store
107			the total mass of salt stored within each grid cell's ice, but
108			it would take some more thinking on my part to be sure that I
109			conserve salt correctly.
110
111			=== Salt Plumes ===
112
113			I assume that ice production in leads can generate
114			salt plumes. The fraction of the salt sent to the plume package
115			from ice produced in leads (defined by the open water
116			fraction) is a nonlinear function of ice concentration, AREA.
117
118			Specifically, the function is a logistic curve (sigmoid) with a range and
119			domain {0,1}. The function, f(AREA), has a single free parameter,
120			SEAICE_plumeInflectionPoint, the inflection point of the curve.
121
122			By construction, the function has the following properties:
123			f(1) \approx 1.0
124			f(SEAICE_plumeInflectionPoint) = 0.5
125			f(0) \approx 0.0 (when SEAICE_plumeInflectionPoint \geq 0.5)
126			f(0) > 0.0 (when SEAICE_plumeInflectionPoint < 0.5)
127
128			As AREA --> 1, the open water fraction is confined to
129			narrow leads implying that new ice production is spatially non-uniform,
130			therby violating the assumption of spatially uniform forcing over
131			the grid cell required by KPP.
132
133			To treat the vertical overturning of the high-salinity brine
134			in a more physically realistic way, the salt produced
135			in the leads is sent to depth via the plume package. To assure
136			only narrow leads generate plumes, choose a relatively high
137			SEAICE_plumeInflectionPoint. For testing purposes, I chose 0.8
138			but I don't know what the "right" value should be.
139
140			=== Ice Age ===
141
142			Adding ice age means just adding a short snippet at the end
143			of seaice_growth.F. I propose that adding to ice age be done
144			in a separate subroutine which takes as arguments, AREA, AREANm1,
145			HEFF, HEFFNm1, and whatever other variables from seaice_growth
146			it needs. Interfacing with a separate subroutine will make
147			for a cleaner seaice growth code in the long run.