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	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.