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 * = 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 (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.