[MITgcm-devel] Merging seaice_growth_if into the main branch
Ian Fenty
ifenty at MIT.EDU
Wed Feb 2 12:17:08 EST 2011
MITgcm devel,
I have gone through the rewrite of seaice_growth and identified a few things that could be added and changed to reconcile it with seaice_growth_if.F. As I see it, the rewritten seaice_growth.F is a significant improvement on the old "legacy" code because: 1) it is now cleaner and has more consistent variable names, 2) there is a clearer and more logical separation of steps, 3) the storage of variables for the adjoint is improved, and 4) some of the parameterizations used in seaice_growth_if are incorporated.
Even though some of the parameterization from seaice_growth_if.F are incorporated, I was able to identify a few that are sufficiently different as to lead to very different results in both forward and adjoint modes. Therefore, I propose the addition of a few well-defined blocks of code, wrapped in #ifdef statements, which permit a closer reproduction of seaice_growth_if.F. In addition, I added a number of diagnostic variables that I found useful for debugging and analysis. The updated seaice_growth.F and seaice_diagnostics_init.F can be found in the contrib: /MITgcm_contrib/ifenty/Fenty_seaice_thermo_code_merge/code_update
A brief overview of the proposed code changes makes up the bulk of this email. I also conducted a more in-depth analysis of the proposed changes which can be downloaded from the MITgcm_contrb directory:
ifenty/Fenty_seaice_thermo_code_merge/documentation/Seaice_Growth_Forward_and_Adjoint_comparisons.pdf
or from a server here at JPL:
http://ecco2.jpl.nasa.gov/data1/meeting/2011/02_02/Seaice_Growth_Forward_and_Adjoint_comparisons.pdf
The aforementioned analysis employs a new verification experiment, a 1-D ocean-sea ice column (MITgcm_contrib/ifenty/Fenty_seaice_thermo_code_merge/verification/1D_ocean_ice_column) and a modified version of the classic lab_sea setup (/MITgcm_contrib/ifenty/Fenty_seaice_thermo_code_merge/verification/lab_sea). It is my hope that these verification experiments be incorporated into the main branch.
It is my belief that if the changes I propose are incorporated into the main branch, all remaining *_if.F routines can be permanently removed from the repository. Anyone using the *_if routines can safely transition to using the new seaice_growth.F and simply define a few new compile-time flags in SEAICE_OPTIONS.h.
In modifying seaice_growth.F, I tried to avoid making any changes that would break the existing verification experiments. Provided that the updated routines don't break anything, who should be the one who brings the changes into the main branch, Dimitris and I or someone at MIT?
Of course, I am open to any comments or suggestions concerning the proposed changes (especially in the naming of the #ifdef flags).
- Ian
** Description of Code Changes **
Because seaice_growth.F is somewhat complicated, changing a parameterization often triggers changes to other parameterizations. Consequently, although the parameterization variants below are relatively self-contained, and one can easily imagine other variants, they are interdependent (unless described as OPTIONAL) and will not work in isolation.
====================================
I. Turbulent ocean-ice heat flux parameterization
#define MCPHEE_OCEAN_ICE_HEAT_FLUX
Thesis: 2.2.3, Appendix A.3.3
====================================
McPhee gives an empirically based closure for unresolved turbulent ocean-ice heat fluxes
(e.g., Air-ice-ocean interaction: turbulent ocean boundary layer exchange processes). In it, the far-field ocean temperature, T, the friction velocity beneath ice, u*,
an observationally-inferred constant the Stanton number, St, are used to parameterize <w'T'>:
<w'T'>= S_t * u_star * (T_o - T_f)
Besides having more empirical support than the old GAMMA timescales that we have been using, parameterizing our fluxes using the language and equations of McPhee will make the MITgcm more comprehensible to the wider sea ice modelling community.
*** New code snippet:
--------------------------------------------
#ifdef MCPHEE_OCEAN_ICE_HEAT_FLUX
c The turbulent ocean-ice heat flux
a_QbyOCN(I,J) = -STANTON_NUMBER * USTAR_BASE * rhoConst *
& HeatCapacity_Cp *(surf_theta - TBC)*
& MixedLayerTurbulenceFactor*hFacC(i,j,kSurface,bi,bj)
c The turbulent ocean-ice heat flux converted to meters
c of potential ice melt
a_QbyOCN(I,J) = a_QbyOCN(I,J) * convertQ2HI
#endif /* MCPHEE_OCEAN_ICE_HEAT_FLUX */
--------------------------------------------
====================================
II. New ice growth criteria in ice-free fraction of grid cell and open water air-sea heat fluxes
#define FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
Thesis: Appendix A.6, A.3.1
====================================
The conditions triggering the expansion of new ice cover in this parameterization are that air-sea heat flux divergence exceeds the potential ocean-ice heat flux convergence in the open-water fraction of the grid cell. Consequently, sea ice forms in grid cells whose average temperature slightly exceeds the seawater freezing point - as observed. Importantly, maintaining the uppermost ocean grid cell slightly above the freezing point ensures that the turbulent ocean-ice heat flux parameterization is always operational, an important requirement for a well-behaved adjoint.
A slight difference from Gael's formulation is that only shortwave radiative fluxes impinging on the ice-free portion of the grid cell are used to when calculating changes to ice thickness due to open-water air-sea heat fluxes. In other words, SW radiative flux convergence beneath the ice warms the ocean and does not melt ice directly.
*** New code snippet:
--------------------------------------------
#ifdef FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
c The sum of ice growth tendencies from open water fluxes
c + any residual ocean-ice fluxes (weighted by the open
c water fraction).
c Using the McPhee parameterization, r_QbyOCN .LE ZERO
c Therefore, initial ice growth is triggered by open water
c heat flux divergence overcoming the melting tendency of
c the ocean-ice heat fluxes.
tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
& (1.0 _d 0 - AREApreTH(i,J))
c A fraction (SWFRACB) of the open-water heat flux includes
c shortwave radiative convergence at depths below the upper
c ocean grid cell. These must be subtracted out.
tmpscal2=SWFRACB * a_QSWbyATM_open(I,J)
tmpscal3 = MAX(0.0 _d 0, tmpscal1-tmpscal2)
d_HEFFbyATMonOCN(I,J)=tmpscal
#endif /* FENTY_DELTA_HEFF_OPEN_WATER_FLUXES */
--------------------------------------------
====================================
III. Changing sea ice concentration
#define FENTY_AREA_EXPANSION_CONTRACTION
Thesis: Appendix A.5
====================================
The three categories of heat flux that change ice thickness (air-sea heat flux, air-ice heat flux, and ocean-ice flux) are logically separated and used to change ice concentration separately, an important requirement for the adjoint.
*** New code snippet:
--------------------------------------------
#ifdef FENTY_AREA_EXPANSION_CONTRACTION
C apply tendency
IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
& (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
AREA(I,J,bi,bj)=max(0. _d 0 , min( 1. _d 0,
& AREA(I,J,bi,bj)+d_AREAbyATM(I,J)
#ifdef FENTY_AREA_EXPANSION_CONTRACTION
& + d_AREAbyOCN(I,J) + d_AREAbyICE(I,J)
#endif
--------------------------------------------
The change of ice thickness by air-ice heat fluxes is given its own name (d_HEFFbyATMonICE) instead of being lumped together with d_HEFFbyATMonOCN
*** New code snippet:
--------------------------------------------
#ifdef FENTY_AREA_EXPANSION_CONTRACTION
tmpscal2 = d_HEFFbyATMOnICE(I,J)
c Ice concentrations are reduced whenever existing ice thins from surface
c heat flux convergence.
if (tmpscal2 .LE. ZERO) then
d_AREAbyICE(I,J) =
& HALF * tmpscal2 * AREApreTH(I,J)/heff_star
endif
#endif /* FENTY_AREA_EXPANSION_CONTRACTION */
--------------------------------------------
Only net air-sea heat flux divergence in the ice-free fraction of the grid cell can increase ice concentration. Using the McPhee parameterization, there is no possible increase of ice thickness when seawater falls below the freezing point. Importantly, using the same parameterization seawater temperature never falls below the freezing point due to the exchange of heat with the atmosphere.
*** New code snippet:
--------------------------------------------
#ifdef FENTY_AREA_EXPANSION_CONTRACTION
tmpscal3 = d_HEFFbyOCNonICE(I,J)
c Sensible heat transfer from the ocean to the sea ice thins it and
c reduces concentrations
if (tmpscal3 .LE.ZERO) then
d_AREAbyOCN(I,J) =
& HALF * tmpscal3 * AREApreTH(I,J)/heff_star
endif
#endif /* FENTY_AREA_EXPANSION_CONTRACTION */
--------------------------------------------
====================================
IV. (OPTIONAL) Smoother variation of the friction velocity (u*) in McPhee parameterization as a function of concentration
#define GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR
used with
#define MCPHEE_OCEAN_ICE_HEAT_FLUX
Thesis: Appendix A.6.1 equation A.34
====================================
In my thesis I used two values of u* depending on the presence of ice in a grid cell. This led to a very low values of ocean-ice heat flux even if only a tiny fraction of the grid cell was ice covered. One alternative is to make u* the area-weighted average of the "typical" values of near-surface u* for the ice-covered and ice-free ocean.
*** New code snippet:
--------------------------------------------
#ifdef GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR
MixedLayerTurbulenceFactor = 12.5 _d 0 -
& 11.5 _d 0 * AREApreTH(I,J)
#else
c If ice is present, MixedLayerTurbulenceFactor = 1.0, else 12.50
IF (AREApreTH(I,J) .GT. 0. _d 0) THEN
MixedLayerTurbulenceFactor = ONE
ELSE
MixedLayerTurbulenceFactor = 12.5 _d 0
ENDIF
#endif /* GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR */
--------------------------------------------
====================================
V. (OPTIONAL) Open water air-sea heat flux convergence melts ice or warms the ocean
#define FENTY_OPEN_WATER_FLUXES_MELT_ICE
used with
#define FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
Thesis: Appendix A.3.1
====================================
In some ice models open water air-sea heat flux convergence is used to melt ice directly and in others it is used to warm the ocean. I provide the option for both.
*** New Code snippet
--------------------------------------------
#ifdef FENTY_OPEN_WATER_FLUXES_MELT_ICE
c Convergent open water heat fluxes melt ice directly,
c bypassing the ocean
tmpscal3 = tmpscal1-tmpscal2
#else
c Convergent open water heat fluxes warm the ocean
c leading to an indirect ice melt.
tmpscal3 = MAX(0.0 _d 0, tmpscal1-tmpscal2)
#endif /* FENTY_OPEN_WATER_FLUXES_MELT_ICE */
d_HEFFbyATMonOCN(I,J)=MAX(tmpscal3, -HEFF(I,J,bi,bj))
d_HEFFbyATMonOCN_open(I,J)=d_HEFFbyATMonOCN(I,J)
r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-d_HEFFbyATMonOCN(I,J)
HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyATMonOCN(I,J)
--------------------------------------------
====================================
VI. New diagnostic variables
====================================
diagName = 'SIdHbATO' : change of ice thickness due to open water air-sea fluxes
diagName = 'SIdHbATC' : change of ice thickness due to air-ice fluxes
diagName = 'SIdHbOCN' : change of ice thickness due to ocean-ice fluxes
diagName = 'SIdA' : total ice concentration change
diagName = 'SIdAbATO' : change of ice area due to open water air-sea fluxes
diagName = 'SIdAbATC' : change of ice area due to air-ice fluxes
diagName = 'SIdAbOCN' : change of ice area due to ocean-ice fluxes
diagName = 'SIaQbATO' : potential ice thickness change due to open water air-sea fluxes
diagName = 'SIaQbATC' : potential ice thickness change due to air-ice fluxes
diagName = 'SIaQbOCN' : potential ice thickness change due ocean-ice fluxes
====================================
VI. Summary
====================================
Three compile-time flags must be defined to ensure that the new parameterizations for 1) ocean-ice heat fluxes, 2) ice growth from open water heat flux divergence, and 3) ice concentration expansion and contraction:
#define MCPHEE_OCEAN_ICE_HEAT_FLUX
#define FENTY_DELTA_HEFF_OPEN_WATER_FLUXES
#define FENTY_AREA_EXPANSION_CONTRACTION)
There are additional parameterization variants that one may use when the above flags are defined. The first is a formulation of u* which gives a slightly smoother solution in forward and adjoint mode (#define GRADIENT_MIXED_LAYER_TURBULENCE_FACTOR). The second allows one to choose whether open water heat flux convergence melts ice directly or indirectly by first warming the ocean (#define FENTY_OPEN_WATER_FLUXES_MELT_ICE).
I suggest the addition of a new variable to keep track of the change of ice thickness due to heat flux convergence at the snow/ice surface: d_HEFFbyATMonICE. I also suggest that the three heat flux terms which change ice thickness (open water air-sea, ocean-ice, and air-ice) be given separate for the purposes of better diagnostics (d_AREAbyATM, d_AREAbyOCN, d_AREAbyICE).
New diagnostics are added to keep track of the various contributions of these fluxes on ice thickness and area.
More information about the MITgcm-devel
mailing list