Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/Thermal_properties.pro
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/Thermal_properties.pro	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/Thermal_properties.pro	(revision 10845)
@@ -0,0 +1,269 @@
+;
+;------------------------------------------------------------------------------
+;
+; Plot LIM thermal properties
+; 
+;------------------------------------------------------------------------------
+
+; Model parameters
+; Sample experiments
+N_spl = 3 & Ts_spl = FLTARR(N_spl)
+Ts_spl[*] = [ -15., -10., -5. ] & Tw_spl = 2.
+colors = [ 50., 100., 150. ]
+
+;
+;==============================================================================
+; 1) User parameters
+;==============================================================================
+;
+
+out_dir=''
+file_out   = 'Thermal_properties'
+
+numplot_x = 2       ; number of horizontal plots
+numplot_y = 2       ; number of vertical plots
+ct = 13             ; colortable
+cs = 1.0            ; size of fonts on plots (1=normal)
+th = 2.             ; thickness of curves
+device = 'PS'       ; 'PS' or 'X'
+
+;
+;==============================================================================
+; 2) Initialize graphics
+;==============================================================================
+;
+
+dx = 6 
+dy = 6.
+figuresize_x = numplot_x * dx; figure size on x direction
+figuresize_y = numplot_y * dy; figure size on y direction
+SET_PLOT, DEVICE
+;!P.FONT=1
+TVLCT, [0,255,0,0], [0,0,255,0], [0,0,0,255]
+IF ( device EQ 'PS' ) THEN BEGIN
+;   DEVICE, /COLOR, /LANDSCAPE, filename=out_dir+file_out+'.ps', $
+;           XSIZE=figuresize_x,YSIZE=figuresize_y,FONT_SIZE=12.0, SET_FONT='Helvetica', BITS = 16
+   DEVICE, /COLOR, /LANDSCAPE, filename=out_dir+file_out+'.ps', $
+           XSIZE=figuresize_x,YSIZE=figuresize_y,FONT_SIZE=9.0
+ENDIF
+
+IF ( device EQ 'X' ) THEN BEGIN
+   xsize = 1200
+   ysize = 800
+   colorkey = 'rd'
+   init_graphics_x, xsize, ysize, colorkey
+ENDIF
+
+!P.MULTI=[0,numplot_x,numplot_y]
+!X.MARGIN   = [6,3]
+!Y.MARGIN   = [4,3]
+
+;==============================================================================
+; 3) Computations and plots
+;==============================================================================
+
+; arrays
+N_T = 200
+N_S = 3
+miss_val = 9.99e32
+
+;--- Physical constants
+c0 = 2067.d0
+cw = 3991.d0
+L  = 334000.d0
+mu = 0.054d0
+rhoi = 917.d0
+rhow = 1026.d0
+rhos = 330.d0
+
+Sw = 34.d0
+
+beta1 = 0.09; pringle coefficients
+beta2 = 0.011;
+
+;--- Ice salinity
+;S = ( rhoi - rhos ) / rhoi * Sw
+S = FLTARR(N_S)
+Smin = 5.   
+Smax = 15.
+dS = ( Smax - Smin ) / FLOAT(N_S-1)
+S[0:N_S-1] = Smin + FINDGEN(N_S) * dS
+
+Tf = -mu * S
+
+;--- Temperatures
+T = FLTARR(N_T)
+Tmin = -20.
+Tf   = -mu * Smax
+Tmax = 0. 
+dT = ( Tmax - Tmin ) / FLOAT(N_T-1)
+T[0:N_T-1] = Tmin + FINDGEN(N_T) * dT
+
+;--- Grids
+S_grid = FLTARR(N_S,N_T)
+T_grid = FLTARR(N_S,N_T)
+
+FOR j = 0, N_T - 1 DO BEGIN
+   S_grid(0:N_S-1,j) = S(0:N_S-1)
+ENDFOR
+
+FOR i = 0, N_S - 1 DO BEGIN
+   T_grid(i,0:N_T-1) = T(0:N_T-1)
+ENDFOR
+
+;--- Brine fraction
+phi = -mu * S_grid / T_grid
+
+;--- specific enthalpy
+q_i = c0 * ( T_grid + mu * S_grid ) - L * ( 1 + mu * S_grid / T_grid ) - cw*mu*S_grid
+
+;--- thermal conductivity (pringle)
+kpri = 2.11 + beta1 * S_grid / T_grid - beta2 * T_grid; 
+
+;--- specific heat
+cp = c0 + L*mu*S_grid/(T_grid*T_grid)
+
+;--- PLOTS
+!P.MULTI=[0, numplot_x, numplot_y]
+LOADCT, 13
+
+PLOT, [Tmin, Tmax], [0., 50.], /NODATA, XTITLE='T (C)', YTITLE='phi (%)', CHARSIZE = cs
+FOR i = 0, N_S - 1 DO BEGIN
+   zaddr = where(phi LT 100) 
+   OPLOT, T_grid(i,zaddr), phi(i,zaddr)*100., thick = th, color = i* 50
+ENDFOR
+
+PLOT, [Tmin, Tmax], [-4., 0.], /NODATA, XTITLE='T (C)', YTITLE='q_m (10^5 J/kg)', CHARSIZE = cs
+FOR i = 0, N_S - 1 DO BEGIN
+   zaddr = where(phi LT 100) 
+   OPLOT, T_grid(i,zaddr), q_i(i,zaddr)/1.0e5, thick = th, color = i* 50
+ENDFOR
+
+PLOT, [Tmin, Tmax], [0., 3.], /NODATA, XTITLE='T (C)', YTITLE='k_i (W/m/K)', CHARSIZE = cs
+FOR i = 0, N_S - 1 DO BEGIN
+   zaddr = where(phi LT 100) 
+   OPLOT, T_grid(i,zaddr), kpri(i,zaddr), thick = th, color = i* 50
+ENDFOR
+
+PLOT, [Tmin, Tmax], [0., 20.], /NODATA, XTITLE='T (C)', YTITLE='c_i (10^3 J/kg/K)', CHARSIZE = cs
+FOR i = 0, N_S - 1 DO BEGIN
+   zaddr = where(phi LT 100) 
+   OPLOT, T_grid(i,zaddr), cp(i,zaddr)/1000., thick = th, color = i* 50
+ENDFOR
+
+
+;T[0:N_T-1] = FINDGEN(N_T) * dT + Tmin
+;dT = ( Tmax - Tf) / FLOAT(N/2-1)
+;T[N/2:N-1] = FINDGEN(N/2) * dT + Tf
+
+;;PLOT, T, TITLE = "T", CHARSIZE = cs
+;
+;; Specific heat
+;ci = FLTARR(N) & ci(*) = miss_val
+;cis = FLTARR(N) & cis(*) = miss_val
+;ci(0:N/2) = c0 + L*mu*S/T(0:N/2)/T(0:N/2)
+;cis(0:N/2) = mu*(c0 - L/T(0:N/2))
+;;PLOT, T(0:N/2), ci(0:N/2)/L, CHARSIZE = cs
+;;OPLOT, T(0:N/2), cis(0:N/2)/L, LINESTYLE = 1
+;
+;;--- Brine volume
+;e_i = FLTARR(N) & e_i(*) = 1.
+;i_addr = FINDGEN(N/2+1)
+;e_i(i_addr) = -mu * S / T(i_addr)
+;PLOT, T, e_i, CHARSIZE = cs, TITLE = "e"
+;OPLOT, T, REPLICATE(0.5, N), linestyle = 1
+;T_50 = -mu * S / 0.50 ; temperature at which brine volume is 50%
+;
+;; Notz (2005), Phd thesis, page 44
+;i_addr = WHERE(T LT Tf)
+;Sbr = - 21.4 * T(i_addr) - 0.886 * T(i_addr) * T(i_addr) - 0.0170 * T(i_addr) * T(i_addr) * T(i_addr)
+;zrS = S / Sbr(i_addr)
+;beta_ocs = 0.81
+;
+;zrhoi      = 916.8 - 0.1403 * T(i_addr)              ; pure ice density Pounder(1965)
+;zrhol      = rhow + beta_ocs * ( Sbr(i_addr) - Sw )
+;zrR        = zrhoi / zrhol
+;znum       = 1 - zrS
+;zden       = 1 + zrS * ( zrR - 1 )
+;e_i2 = FLTARR(N) & e_i2(*) = 1.
+;e_i2(i_addr) =  1. - znum / zden
+;;OPLOT, T, e_i2, linestyle = 3
+;
+;;OPLOT, [ T_50, T_50], [ 0., 1.], linestyle = 1
+;
+;;--- Enthalpies
+;Ei = FLTARR(N) & Ei(*) = miss_val
+;Ei2= FLTARR(N) & Ei2(*) = miss_val
+;Es = FLTARR(N) & Es(*) = miss_val
+;Es_wgt = FLTARR(N) & Es(*) = miss_val
+;Ew = FLTARR(N) & Ew(*) = miss_val
+;
+;i_addr = WHERE( T LE 0. )
+;Es(i_addr) = c0 * T(i_addr) - L
+;Es_wgt(i_addr) = Es(i_addr) * rhos / rhoi
+;i_addr = WHERE( T GT 0. )
+;Es(i_addr) = cw * T(i_addr)
+;Es_wgt(i_addr) = Es(i_addr) * ( rhoi - rhos ) / rhoi
+;
+;i_addr = WHERE( T GE Tf )
+;Ew(i_addr) = cw * T(i_addr)
+;i_addr = FINDGEN(N/2+1)
+;
+;Ei(i_addr) = c0 * ( T(i_addr) + mu * S ) - L * ( 1 + mu * S / T(i_addr) ) - cw*mu*S
+;Ei2(i_addr) = - ( 1 - e_i(i_addr) ) * L
+;
+;Emax = MAX(Ew(WHERE(Ew NE miss_val)))
+;
+;PLOT, [Tmin, Tmax], [ MIN(Ei)/1000., Emax/1000. ], /NODATA, TITLE = "Enthalpy", CHARSIZE = cs, XTITLE = "degC", YTITLE = "kJ/kg", YMINOR = 2
+;OPLOT, T, Ei/1000., MAX=0., THICK = 3.;, solid black
+;OPLOT, T, Ei2/1000., MAX=0., COLOR = 150, THICK = 3.; solid grey
+;;OPLOT, T, Es/1000., LINESTYLE = 1, THICK = 3.
+;OPLOT, T, Es_wgt/1000., COLOR = 100, THICK = 3
+;i_addr = WHERE( T GE 0. )
+;OPLOT, T(i_addr), Es_wgt(i_addr)/1000., COLOR = 200, THICK = 3
+;;OPLOT, T, Ew/1000., MAX=Emax/1000., LINESTYLE = 0, THICK = 3.
+;OPLOT, [ T_50, T_50], [ -400., 100.], COLOR = 100
+;OPLOT, [ Tmin, Tmax], [0., 0.], COLOR = 100
+;
+;;--- Sample experiment
+;Es_spl = c0 * Ts_spl - L
+;Ew_spl = cw * Tw_spl
+;
+;LOADCT, ct
+;FOR i = 0, N_spl - 1 DO BEGIN
+;   OPLOT, [ Ts_spl(i) ], [ rhos / rhoi * Es_spl(i) ]/1000., PSYM = 1, SYMSIZE = 2, THICK = 3, COLOR = COLORS(i)
+;ENDFOR
+;LOADCT, 0
+;OPLOT, [ Tw_spl ], [ (rhoi - rhos) / rhoi * Ew_spl ]/1000., PSYM = 1, SYMSIZE = 2, THICK = 3
+;
+;Ei_spl = Es_spl * rhos / rhoi + Ew_spl * ( rhoi - rhos ) / rhoi
+;
+;; Retrieve T
+;zA = ( L *( rhos/rhoi ) - (rhos/rhoi) * c0 * Ts_spl + (rhos-rhoi)/rhoi * cw * Tw_spl ) * rhoi
+;
+;aaa = c0
+;bbb = mu * S * ( c0 - cw ) - L + zA / rhoi
+;ccc = - L * mu * S
+;
+;discr = SQRT ( bbb * bbb - 4 * aaa * ccc )
+;
+;Ti_spl = ( - bbb - discr ) / 2.0 / aaa
+;
+;Ei_test = c0 * ( Ti_spl + mu * S ) - L * ( 1 + mu * S / Ti_spl ) - cw*mu*S
+;PRINT, "Ei_test : ", Ei_test
+;PRINT, "Ei_spl  : ", Ei_spl
+;
+;; Retrieve brine volume
+;e_i_spl = -mu*S/Ti_spl
+;PRINT, "e_i_spl = ", e_i_spl
+;
+;LOADCT, ct
+;FOR i = 0, N_spl - 1 DO BEGIN
+;   OPLOT, [ Ti_spl(i) ], [ Ei_spl(i) ]/1000., PSYM = 1, SYMSIZE = 2, COLOR = COLORS(i), THICK = 3
+;ENDFOR
+;LOADCT, 0
+
+;--------------------------------------------------------------------------------------------------
+DEVICE, /CLOSE
+SET_PLOT, 'X'
+END
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/albedo_cloud_correction.m
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/albedo_cloud_correction.m	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/albedo_cloud_correction.m	(revision 10845)
@@ -0,0 +1,51 @@
+clear all; close all;
+
+% Script to plot the cloud fraction dependency of albedo correction
+
+%------------------
+% Sea Ice
+%------------------
+
+% overcast albedo
+alpha_os = 0:0.001:1.;
+
+% SI3 correction
+f = ( -0.1010 * alpha_os.^2 + 0.1933 * alpha_os - 0.0148 );
+
+% clear sky sea ice albedo
+alpha_cs = alpha_os - f;
+
+% Grenfell and Perovich, Table 3 data
+alpha_GP04 = [ 0.114 0.151 0.141 0.266 0.367 0.309 0.332 0.492 0.678 0.928 ];
+f_GP04     = [ -0.010 -0.022 0.012 0.027 0.020 0.036 0.037 0.056 0.071 0.077 ];
+
+%--------------
+% Ocean
+%--------------
+alpha_oce_cs = 0.05 / ( 1.1 * 0.40^1.4 + 0.15 ) % --- Briegleb and Ramanathan 1972
+alpha_oce_os = 0.06 %--- Payne 1972
+f_oce = alpha_oce_os - alpha_oce_cs
+
+%--------------
+% Plot
+%--------------
+plot(alpha_os,f,'k'); hold on
+zaddr = find( ( alpha_os >= alpha_GP04(3)  )   & ...
+              ( alpha_os <= alpha_GP04(10) ) );
+plot(alpha_os(zaddr),f(zaddr),'k', 'LineWidth', 3)
+plot(alpha_GP04(3:10),f_GP04(3:10),'ksq','MarkerFaceColor', 'k')
+plot(alpha_GP04(5),f_GP04(5),'rsq','MarkerFaceColor','r')
+text(alpha_GP04(5)+0.02,f_GP04(5),'dirty ice')
+plot(alpha_GP04(1:2),f_GP04(1:2),'rsq','MarkerFaceColor','r')
+text(alpha_GP04(1)+0.02,f_GP04(1),'wet tundra')
+text(alpha_GP04(2)+0.02,f_GP04(2),'dry tundra')
+plot(alpha_oce_os,f_oce,'bo','MarkerFaceColor', 'b')
+text(alpha_oce_os+0.02,f_oce,'ocean')
+plot( [ alpha_GP04(3) alpha_GP04(3) ], [ -0.06 0.08 ], 'k--' )
+plot( [ alpha_GP04(10) alpha_GP04(10) ], [ -0.06 0.08 ], 'k--' )
+
+xlabel('\alpha^{os}')
+ylabel('\alpha^{os}-\alpha^{cs}')
+
+set(gca,'fontsize', 16)
+set(gca, 'FontName', 'Helvetica LT Std')
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/albedo_hi.m
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/albedo_hi.m	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/albedo_hi.m	(revision 10845)
@@ -0,0 +1,71 @@
+clear all, close all
+
+% plots SI3 albedo dependence versus ice thickness
+
+z1_href_pnd = 1. / 0.05
+z1_c1 = 1. / ( log(1.5) - log(0.05) );
+z1_c2 = 1. / 0.05;
+z1_c3 = 1. / 0.02
+z1_c4 = 1. / 0.03
+
+rn_alb_idry = 0.60;
+rn_alb_imlt = 0.50;
+rn_alb_sdry = 0.85
+rn_alb_smlt = 0.75
+rn_alb_oce  = 0.06;
+rn_alb_dpnd = 0.27;
+
+h_i = 0:0.001:3;
+
+%--- Bare ice albedo (for hi > 150cm)
+zalb_ice    = h_i;
+zalb_ice(:) = rn_alb_idry;
+
+%--- Bare ice albedo (for hi < 150cm)
+zaddr = find( ( 0.05 < h_i ) & ( h_i <= 1.5 ) ); % 5cm < hi < 150cm
+zalb_ice(zaddr) = zalb_ice(zaddr) + ( 0.18 - zalb_ice(zaddr) ) * z1_c1 .* ( log(1.5) - log(h_i(zaddr)) );
+zaddr = find( h_i <= 0.05 );
+zalb_ice(zaddr) = rn_alb_oce  + ( 0.18 - rn_alb_oce ) * z1_c2 * h_i(zaddr);
+
+figure; 
+
+
+subplot(1,3,1)
+
+plot(h_i,zalb_ice,'k','LineWidth', 3)
+xlabel('h_i(m)')
+ylabel('\alpha')
+ylim([0.1 0.9])
+set(gca,'fontsize', 16)
+set(gca, 'FontName', 'Helvetica LT Std')
+
+
+
+%--- Snow covered albedo
+
+h_s = 0:0.001:0.6;
+zalb_ice    = 0.6
+zalb_snw    = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * exp( -h_s * z1_c3 )
+
+subplot(1,3,2)
+plot(h_s,zalb_snw,'k','LineWidth', 3)
+xlabel('h_s(m)')
+ylabel('\alpha')
+ylim([0.1 0.9])
+set(gca,'fontsize', 16)
+set(gca, 'FontName', 'Helvetica LT Std')
+
+%--- Ponded ice albedo
+zalb_ice    = 0.5
+h_pnd = 0:0.001:0.6;
+zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * exp( -h_pnd * z1_href_pnd )
+
+subplot(1,3,3)
+plot(h_pnd,zalb_pnd,'k','LineWidth', 3)
+xlabel('h_{pnd}(m)')
+ylabel('\alpha')
+ylim([0.1 0.9])
+set(gca,'fontsize', 16)
+set(gca, 'FontName', 'Helvetica LT Std')
+
+
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/ice_cats.pro
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/ice_cats.pro	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/Figures/scripts/ice_cats.pro	(revision 10845)
@@ -0,0 +1,157 @@
+
+;
+;==============================================================================
+; 1) User parameters
+;==============================================================================
+;
+
+out_dir=''
+file_out   = 'ThicknessSpace'
+
+numplot_x = 1       ; number of horizontal plots
+numplot_y = 1       ; number of vertical plots
+ct = 13             ; colortable
+cs = 2.0            ; size of fonts on plots (1=normal)
+th = 2.             ; thickness of curves
+device = 'PS'       ; 'PS' or 'X'
+
+;
+;==============================================================================
+; 2) Initialize graphics
+;==============================================================================
+;
+
+dx = 12.
+dy = 8.
+figuresize_x = numplot_x * dx; figure size on x direction
+figuresize_y = numplot_y * dy; figure size on y direction
+SET_PLOT, DEVICE
+;!P.FONT=1
+TVLCT, [0,255,0,0], [0,0,255,0], [0,0,0,255]
+IF ( device EQ 'PS' ) THEN BEGIN
+;   DEVICE, /COLOR, /LANDSCAPE, filename=out_dir+file_out+'.ps', $
+;           XSIZE=figuresize_x,YSIZE=figuresize_y,FONT_SIZE=12.0, SET_FONT='Helvetica', BITS = 16
+   DEVICE, /COLOR, /LANDSCAPE, filename=out_dir+file_out+'.ps', $
+           XSIZE=figuresize_x,YSIZE=figuresize_y,FONT_SIZE=9.0
+ENDIF
+
+IF ( device EQ 'X' ) THEN BEGIN
+   xsize = 1200
+   ysize = 800
+   colorkey = 'rd'
+   init_graphics_x, xsize, ysize, colorkey
+ENDIF
+
+!P.MULTI=[0,numplot_x,numplot_y]
+!X.MARGIN   = [7,3]
+!Y.MARGIN   = [4,3]
+
+;
+;==============================================================================
+; 2) Initialize graphics
+;==============================================================================
+;
+
+jpl = 5
+jpl1 = 10
+
+;-----------------------
+; original CICE version
+;-----------------------
+
+hi_max_cice_5 = FLTARR(jpl+1)
+hi_max_cice_5(0) = 0
+zc1 =  3. / FLOAT(jpl)
+zc2 = 10. * zc1 
+zc3 =  3.
+
+FOR jl = 1, jpl  DO BEGIN
+   zx1 = FLOAT( jl-1 ) / FLOAT( jpl )
+   hi_max_cice_5(jl) = hi_max_cice_5(jl-1) + zc1 + zc2 * (1. + TANH( zc3 * (zx1 - 1. ) ) ) 
+ENDFOR
+
+hi_max_cice_10 = FLTARR(jpl1+1)
+hi_max_cice_10(0) = 0
+zc1 =  3. / FLOAT(jpl1)
+zc2 = 10. * zc1 
+zc3 =  3.
+
+FOR jl = 1, jpl1  DO BEGIN
+   zx1 = FLOAT( jl-1 ) / FLOAT( jpl1 )
+   hi_max_cice_10(jl) = hi_max_cice_10(jl-1) + zc1 + zc2 * (1. + TANH( zc3 * (zx1 - 1. ) ) ) 
+ENDFOR
+
+;-------------------
+; LIM3 distribution
+;-------------------
+alpha = 1./20.
+
+; 5 categories, Arctic ice
+hmean = 2.5
+hL = 3.*hmean
+
+hi_max_1 = FLTARR(jpl + 1)
+hi_max_1(0) = 0.
+
+FOR jl = 1, jpl DO BEGIN
+   znum = jpl * ( hL+1 )^alpha
+   zden = ( jpl - jl ) * ( hL+1 )^alpha + jl
+   hi_max_1(jl) = ( znum / zden )^(1./alpha) - 1
+ENDFOR
+
+; 10 categories, Arctic ice
+hi_max_2 = FLTARR(jpl1 + 1)
+hi_max_2(0) = 0.
+FOR jl = 1, jpl1 DO BEGIN
+   znum = jpl1 * ( hL+1 )^alpha
+   zden = ( jpl1 - jl ) * ( hL+1 )^alpha + jl
+   hi_max_2(jl) = ( znum / zden )^(1./alpha) - 1
+ENDFOR
+
+; 5 categories, Spitzberg
+hmean = 1.0
+hL = 3.*hmean
+
+hi_max_3 = FLTARR(jpl + 1)
+hi_max_3(0) = 0.
+
+FOR jl = 1, jpl DO BEGIN
+   znum = jpl * ( hL+1 )^alpha
+   zden = ( jpl - jl ) * ( hL+1 )^alpha + jl
+   hi_max_3(jl) = ( znum / zden )^(1./alpha) - 1
+ENDFOR
+
+; 10 categories, Spitzberg
+hmean = 1.0
+hi_max_4 = FLTARR(jpl1 + 1)
+hi_max_4(0) = 0.
+
+FOR jl = 1, jpl1 DO BEGIN
+   znum = jpl1 * ( hL+1 )^alpha
+   zden = ( jpl1 - jl ) * ( hL+1 )^alpha + jl
+   hi_max_4(jl) = ( znum / zden )^(1./alpha) - 1
+ENDFOR
+
+!P.MULTI = [ 0, numplot_x, numplot_y ]
+
+PLOT, [0.1, 10.], [0., 7.], /NODATA, charsize = cs, XSTYLE = 1, YSTYLE = 1
+OPLOT, hi_max_cice_5, FINDGEN(jpl+1)*0.+1, psym = 4, symsize = 1, thick = 2
+OPLOT, hi_max_cice_10, FINDGEN(jpl1+1)*0.+2, psym = 4, symsize = 1, thick = 2, color =150
+OPLOT, hi_max_1, FINDGEN(jpl+1)*0.+3   , psym = 1, symsize = 1, thick = 2, color = 0
+OPLOT, hi_max_2, FINDGEN(jpl1+1)*0.+4  , psym = 1, symsize = 1, thick = 2, color = 150
+OPLOT, hi_max_3, FINDGEN(jpl+1)*0.+5   , psym = 2, symsize = 1, thick = 2, color = 0
+OPLOT, hi_max_4, FINDGEN(jpl1+1)*0.+6   , psym = 2, symsize = 1, thick = 2, color = 150
+
+PLOT, FINDGEN(jpl+1)/jpl, hi_max_cice_5, psym = 4, symsize = 1, thick = 2, charsize = cs
+OPLOT, FINDGEN(jpl1+1)/jpl1, hi_max_cice_10, psym = 4, symsize = 1, thick = 2, color = 150
+OPLOT, FINDGEN(jpl+1)/jpl, hi_max_1,   psym = 1, symsize = 1, thick = 2
+OPLOT, FINDGEN(jpl1+1)/jpl1, hi_max_2, psym = 1, symsize = 1, thick = 2, color = 150
+OPLOT, FINDGEN(jpl+1)/jpl, hi_max_3,   psym = 2, symsize = 1, thick = 2 
+OPLOT, FINDGEN(jpl1+1)/jpl1, hi_max_4, psym = 2, symsize = 1, thick = 2, color = 150
+
+
+DEVICE, /CLOSE
+SET_PLOT, 'X'
+!P.MULTI=[0,1,1]
+END
+
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namalb
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namalb	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namalb	(revision 10845)
@@ -0,0 +1,10 @@
+!------------------------------------------------------------------------------
+&namalb         !   albedo parameters
+!------------------------------------------------------------------------------
+   !                                  !                          !  obs range (cloud-sky)
+   rn_alb_sdry      =   0.85          !  dry snow albedo         :  0.85 -- 0.87
+   rn_alb_smlt      =   0.75          !  melting snow albedo     :  0.72 -- 0.82
+   rn_alb_idry      =   0.60          !  dry ice albedo          :  0.54 -- 0.65
+   rn_alb_imlt      =   0.50          !  bare puddled ice albedo :  0.49 -- 0.58
+   rn_alb_dpnd      =   0.27          !  ponded ice albedo       :  0.10 -- 0.30 
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdia
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdia	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdia	(revision 10845)
@@ -0,0 +1,9 @@
+!------------------------------------------------------------------------------
+&namdia         !   Diagnostics
+!------------------------------------------------------------------------------
+   ln_icediachk     = .false.         !  check online the heat, mass & salt budgets (T) or not (F)
+   ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F)
+   ln_icectl        = .false.         !  ice points output for debug (T or F)
+   iiceprt          =  10             !  i-index for debug
+   jiceprt          =  10             !  j-index for debug
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn	(revision 10845)
@@ -0,0 +1,17 @@
+!------------------------------------------------------------------------------
+&namdyn         !   Ice dynamics
+!------------------------------------------------------------------------------
+   ln_dynFULL       = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection)
+   ln_dynADV        = .false.         !  dyn.: only advection w prescribed vel.(rn_uvice + advection)
+      rn_uice       =   0.00001       !        prescribed ice u-velocity
+      rn_vice       =   0.            !        prescribed ice v-velocity
+   rn_ishlat        =   2.            !  free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2)
+   ln_landfast      = .false.         !  landfast ice parameterization (T or F)                           
+      rn_gamma      =   0.15          !     fraction of ocean depth that ice must reach to initiate landfast
+                                      !        recommended range: [0.1 ; 0.25]
+      rn_icebfr     =  10.            !     maximum bottom stress per unit area of contact [N/m2]                 
+                                      !        a very large value ensures ice velocity=0 even with a small contact area
+                                      !        recommended range: ?? (should be greater than atm-ice stress => >0.1 N/m2)
+      rn_lfrelax    =   1.e-5         !     relaxation time scale to reach static friction [s-1]
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_adv
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_adv	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_adv	(revision 10845)
@@ -0,0 +1,7 @@
+!------------------------------------------------------------------------------
+&namdyn_adv     !   Ice advection
+!------------------------------------------------------------------------------
+   ln_adv_Pra       = .false.         !  Advection scheme (Prather)
+   ln_adv_UMx       = .true.          !  Advection scheme (Ultimate-Macho)
+      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order)
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_rdgrft
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_rdgrft	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_rdgrft	(revision 10845)
@@ -0,0 +1,25 @@
+!------------------------------------------------------------------------------
+&namdyn_rdgrft  !   Ice ridging/rafting
+!------------------------------------------------------------------------------
+          ! -- ice_rdgrft_strength -- !
+   ln_str_H79       = .true.          !  ice strength param.: Hibler_79   => P = pstar*<h>*exp(-c_rhg*A)                      
+      rn_pstar      =   2.0e+04       !     ice strength thickness parameter [N/m2]
+      rn_crhg       =   20.0          !     ice strength conc. parameter (-)
+                   ! -- ice_rdgrft -- !
+   rn_csrdg         =   0.5           !  fraction of shearing energy contributing to ridging
+              ! -- ice_rdgrft_prep -- !
+   ln_partf_lin     = .false.         !  Linear ridging participation function (Thorndike et al, 1975)
+      rn_gstar      =   0.15          !     fractional area of thin ice being ridged 
+   ln_partf_exp     = .true.          !  Exponential ridging participation function (Lipscomb, 2007)
+      rn_astar      =   0.03          !     exponential measure of ridging ice fraction [set to 0.05 if hstar=100]
+   ln_ridging       = .true.          !  ridging activated (T) or not (F)
+      rn_hstar      =  25.0           !     determines the maximum thickness of ridged ice [m] (Hibler, 1980)
+      rn_porordg    =   0.3           !     porosity of newly ridged ice (Lepparanta et al., 1995)
+      rn_fsnwrdg    =   0.5           !     snow volume fraction that survives in ridging
+      rn_fpndrdg    =   1.0           !     pond fraction that survives in ridging (small a priori)
+   ln_rafting       = .true.          !  rafting activated (T) or not (F)
+      rn_hraft      =   0.75          !     threshold thickness for rafting [m]
+      rn_craft      =   5.0           !     squeezing coefficient used in the rafting function
+      rn_fsnwrft    =   0.5           !     snow volume fraction that survives in rafting
+      rn_fpndrft    =   1.0           !     pond fraction that survives in rafting (0.5 a priori)
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_rhg
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_rhg	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namdyn_rhg	(revision 10845)
@@ -0,0 +1,11 @@
+!------------------------------------------------------------------------------
+&namdyn_rhg     !   Ice rheology
+!------------------------------------------------------------------------------
+   ln_rhg_EVP       = .true.          !  EVP rheology
+      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017)
+      rn_creepl     =   2.0e-9        !     creep limit [1/s]
+      rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve          
+      nn_nevp       = 120             !     number of EVP subcycles                             
+      rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
+                                      !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namforcing
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namforcing	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namforcing	(revision 10845)
@@ -0,0 +1,13 @@
+!------------------------------------------------------------------------------
+&namforcing     !   Ice surface forcing
+!------------------------------------------------------------------------------
+   rn_cio           =   5.0e-03       !  ice-ocean drag coefficient (-)
+   rn_blow_s        =   0.66          !  mesure of snow blowing into the leads
+                                      !     = 1 => no snow blowing, < 1 => some snow blowing
+   nn_flxdist       =  -1             !  Redistribute heat flux over ice categories
+                                      !     =-1  Do nothing (needs N(cat) fluxes)
+                                      !     = 0  Average N(cat) fluxes then apply the average over the N(cat) ice
+                                      !     = 1  Average N(cat) fluxes then redistribute over the N(cat) ice using T-ice and albedo sensitivity
+                                      !     = 2  Redistribute a single flux over categories
+   nice_jules       = 0               !  Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE)
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namini
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namini	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namini	(revision 10845)
@@ -0,0 +1,25 @@
+!------------------------------------------------------------------------------
+&namini         !   Ice initialization
+!------------------------------------------------------------------------------
+   ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
+   ln_iceini_file   = .false.         !  netcdf file provided for initialization (T) or not (F)
+   rn_thres_sst     =   2.0           !  max delta temp. above Tfreeze with initial ice = (sst - tfreeze)
+   rn_hts_ini_n     =   0.3           !  initial real snow thickness (m), North
+   rn_hts_ini_s     =   0.3           !        "            "             South
+   rn_hti_ini_n     =   3.0           !  initial real ice thickness  (m), North
+   rn_hti_ini_s     =   1.0           !        "            "             South
+   rn_ati_ini_n     =   0.9           !  initial ice concentration   (-), North
+   rn_ati_ini_s     =   0.9           !        "            "             South
+   rn_smi_ini_n     =   6.3           !  initial ice salinity     (g/kg), North
+   rn_smi_ini_s     =   6.3           !        "            "             South
+   rn_tmi_ini_n     = 270.            !  initial ice/snw temperature (K), North
+   rn_tmi_ini_s     = 270.            !        "            "             South
+
+   sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hts = 'Ice_initialization'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_ati = 'Ice_initialization'    , -12 ,'ati'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_tsu = 'Ice_initialization'    , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_tmi = 'Ice_initialization'    , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_smi = 'Ice_initialization'    , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   cn_dir='./'
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namitd
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namitd	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namitd	(revision 10845)
@@ -0,0 +1,9 @@
+!------------------------------------------------------------------------------
+&namitd         !   Ice discretization
+!------------------------------------------------------------------------------
+   ln_cat_hfn       = .true.          !  ice categories are defined by a function following rn_himean**(-0.05)
+      rn_himean     =   2.0           !  expected domain-average ice thickness (m)
+   ln_cat_usr       = .false.         !  ice categories are defined by rn_catbnd below (m)
+      rn_catbnd     =   0.,0.45,1.1,2.1,3.7,6.0  
+   rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/nampar
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/nampar	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/nampar	(revision 10845)
@@ -0,0 +1,10 @@
+!------------------------------------------------------------------------------
+&nampar         !   Generic parameters
+!------------------------------------------------------------------------------
+   jpl              =   5             !  number of ice  categories
+   nlay_i           =   2             !  number of ice  layers
+   nlay_s           =   1             !  number of snow layers (only 1 is working)
+   nn_virtual_itd   =   0             !  virtual ITD mono-category parameterizations (1-3 => jpl = 1 only) or not (0)
+                                      !     2: activate enhanced thermal conductivity only --- temporary option
+                                      !     3: activate virtual thin ice melting only      ---  temporary option
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd	(revision 10845)
@@ -0,0 +1,8 @@
+!------------------------------------------------------------------------------
+&namthd         !   Ice thermodynamics
+!------------------------------------------------------------------------------
+   ln_icedH         = .true.          !  activate ice thickness change from growing/melting (T) or not (F)
+   ln_icedA         = .true.          !  activate lateral melting param. (T) or not (F)
+   ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F)
+   ln_icedS         = .true.          !  activate brine drainage (T) or not (F)
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_da
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_da	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_da	(revision 10845)
@@ -0,0 +1,11 @@
+!------------------------------------------------------------------------------
+&namthd_da      !   Ice lateral melting
+!------------------------------------------------------------------------------
+   rn_beta          =   1.0           !  coef. beta for lateral melting param. Recommended range=[0.8-1.2]
+                                      !   => decrease = more melt and melt peaks toward higher concentration (A~0.5 for beta=1 ; A~0.8 for beta=0.2)
+                                      !      0.3 = best fit for western Fram Strait and Antarctica
+                                      !      1.4 = best fit for eastern Fram Strait      
+   rn_dmin          =   8.            !  minimum floe diameter for lateral melting param. Recommended range=[6-10]
+                                      !   => 6  vs 8m = +40% melting at the peak (A~0.5)
+                                      !      10 vs 8m = -20% melting
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_do
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_do	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_do	(revision 10845)
@@ -0,0 +1,9 @@
+!------------------------------------------------------------------------------
+&namthd_do      !   Ice growth in open water
+!------------------------------------------------------------------------------
+   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_hnewice
+   ln_frazil        = .false.         !  Frazil ice parameterization (ice collection as a function of wind)
+      rn_maxfraz    =   1.0           !     maximum fraction of frazil ice collecting at the ice base
+      rn_vfraz      =   0.417         !     thresold drift speed for frazil ice collecting at the ice bottom (m/s)
+      rn_Cfraz      =   5.0           !     squeezing coefficient for frazil ice collecting at the ice bottom
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_pnd
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_pnd	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_pnd	(revision 10845)
@@ -0,0 +1,10 @@
+!------------------------------------------------------------------------------
+&namthd_pnd     !   Melt ponds
+!------------------------------------------------------------------------------
+   ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+      ln_pnd_fwb    = .false.         !     melt ponds store freshwater or not
+   ln_pnd_CST       = .false.         !  activate constant melt ponds
+      rn_apnd       =   0.2           !     prescribed pond fraction, at Tsu=0
+      rn_hpnd       =   0.05          !     prescribed pond depth, at Tsu=0
+   ln_pnd_alb       = .false.         !  melt ponds affect albedo or not
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_sal
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_sal	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_sal	(revision 10845)
@@ -0,0 +1,15 @@
+!------------------------------------------------------------------------------
+&namthd_sal     !   Ice salinity
+!------------------------------------------------------------------------------
+   nn_icesal        =   2             !  ice salinity option
+                                      !     1: constant ice salinity (S=rn_icesal)
+                                      !     2: varying salinity parameterization S(z,t)
+                                      !     3: prescribed salinity profile S(z), Schwarzacher, 1959
+   rn_icesal        =   4.            !    (nn_icesal=1) ice salinity (g/kg)
+   rn_sal_gd        =   5.            !  restoring ice salinity, gravity drainage (g/kg)
+   rn_time_gd       =   1.73e+6       !  restoring time scale, gravity drainage  (s)
+   rn_sal_fl        =   2.            !  restoring ice salinity, flushing (g/kg)
+   rn_time_fl       =   8.64e+5       !  restoring time scale, flushing (s)
+   rn_simax         =  20.            !  maximum tolerated ice salinity (g/kg)
+   rn_simin         =   0.1           !  minimum tolerated ice salinity (g/kg)
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_zdf
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_zdf	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/namelists/namthd_zdf	(revision 10845)
@@ -0,0 +1,10 @@
+!------------------------------------------------------------------------------
+&namthd_zdf     !   Ice heat diffusion
+!------------------------------------------------------------------------------
+   ln_zdf_BL99      = .true.          !  Heat diffusion follows Bitz and Lipscomb 1999
+   ln_cndi_U64      = .false.         !  sea ice thermal conductivity: k = k0 + beta.S/T            (Untersteiner, 1964)
+   ln_cndi_P07      = .true.          !  sea ice thermal conductivity: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007)
+   rn_cnd_s         =   0.31          !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
+                                      !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
+   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m]
+/
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_main/SI3_manual.bbl
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_main/SI3_manual.bbl	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_main/SI3_manual.bbl	(revision 10845)
@@ -0,0 +1,252 @@
+\begin{thebibliography}{47}
+\providecommand{\natexlab}[1]{#1}
+\providecommand{\url}[1]{\texttt{#1}}
+\providecommand{\urlprefix}{URL }
+\expandafter\ifx\csname urlstyle\endcsname\relax
+  \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
+  \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
+  \urlstyle{rm}\Url}\fi
+\providecommand{\eprint}[2][]{\url{#2}}
+
+\bibitem[{Assur(1958)}]{Assur58}
+Assur, A., 1958: Composition of sea ice and its tensile strength. \textit{A.
+  Nat. Acad. Sci./Nat. Res. Council, Arctic Sea Ice}, \textbf{598}, 106--138.
+
+\bibitem[{Babko et~al.(2002)Babko, Rothrock, and Maykut}]{Babkoetal02}
+Babko, O., D.~A. Rothrock, and G.~A. Maykut, 2002: Role of rafting in the
+  mechanical redistribution of sea ice thickness. \textit{Journal of
+  Geophysical Research}, \textbf{107}, 3113, doi:10.1029/1999JC000\,190,
+  \doi{doi:10.1029/1999JC000190}.
+
+\bibitem[{Bitz et~al.(2001)Bitz, Holland, Weaver, and Eby}]{Bitzetal01}
+Bitz, C.~M., M.~M. Holland, A.~J. Weaver, and M.~Eby, 2001: Simulating the
+  ice-thickness distribution in a coupled climate model. \textit{Journal of
+  Geophysical Research}, \textbf{106}, 2441--2463.
+
+\bibitem[{Bitz and Lipscomb(1999)}]{BitzLipscomb99}
+Bitz, C.~M. and W.~H. Lipscomb, 1999: An energy-conserving thermodynamic model
+  of sea ice. \textit{Journal of Geophysical Research}, \textbf{104},
+  15,669--15,677.
+
+\bibitem[{Bouillon et~al.(2013)Bouillon, Fichefet, Legat, and
+  Madec}]{Bouillonetal13}
+Bouillon, S., T.~Fichefet, V.~Legat, and G.~Madec, 2013: The
+  elastic-viscous-plastic method revisited. \textit{Ocean Modelling},
+  \textbf{71}, 2--12.
+
+\bibitem[{Brandt et~al.(2005)Brandt, Warren, Worby, and
+  Grenfell}]{Brandtetal05}
+Brandt, R.~E., S.~G. Warren, A.~P. Worby, and T.~C. Grenfell, 2005: {Surface
+  albedo of the Antarctic sea ice zone}. \textit{Journal of Climate},
+  \textbf{18}, 3606--3622.
+
+\bibitem[{Coon et~al.(1974)Coon, Maykut, Pritchard, Rothrock, and
+  Thorndike}]{Coonetal74}
+Coon, M.~D., G.~A. Maykut, R.~S. Pritchard, D.~A. Rothrock, and A.~S.
+  Thorndike, 1974: Modeling the pack ice as an elastic-plastic material.
+  \textit{AIDJEX Bulletin}, \textbf{24}, 1--105.
+
+\bibitem[{Ebert and Curry(1993)}]{EbertCurry93}
+Ebert, E.~E. and J.~A. Curry, 1993: An intermediate one-dimensional
+  thermodynamic sea ice model for investigating ice-atmosphere interactions.
+  \textit{Journal of Geophysical Research}, \textbf{98}, 10,085--10,109.
+
+\bibitem[{Fichefet and Morales~Maqueda(1997)}]{FichefetMaqueda97}
+Fichefet, T. and M.~A. Morales~Maqueda, 1997: Sensitivity of a global sea ice
+  model to the treatment of ice thermodynamics and dynamics. \textit{Journal of
+  Geophysical Research}, \textbf{102}, 12,609--12,646.
+
+\bibitem[{Flato and Hibler(1995)}]{FlatoHibler95}
+Flato, G.~M. and W.~D. Hibler, 1995: Ridging and strength in modeling the
+  thickness distribution of arctic sea ice. \textit{Journal of Geophysical
+  Research}, \textbf{100}, 18\,611--18\,626.
+
+\bibitem[{Girard et~al.(2009)Girard, Weiss, Molines, Barnier, and
+  Bouillon}]{Girardetal09}
+Girard, L., J.~Weiss, J.~Molines, B.~Barnier, and S.~Bouillon, 2009: Evaluation
+  of high-resolution sea ice models on the basis of statistical and scaling
+  properties of arctic sea ice drift and deformation. \textit{Journal of
+  Geophysical Research}, \textbf{114}, \doi{10.1029/2008JC005182}.
+
+\bibitem[{Grenfell and Maykut(1977)}]{GrenfellMaykut77}
+Grenfell, T.~C. and G.~A. Maykut, 1977: The optical properties of ice and snow
+  in the arctic basin. \textit{Journal of Glaciology}, \textbf{18}, 445--463.
+
+\bibitem[{Grenfell and Perovich(2004)}]{GrenfellPerovich04}
+Grenfell, T.~C. and D.~K. Perovich, 2004: Seasonal and spatial evolution of
+  albedo in a snow-ice-land-ocean environment. \textit{J. Geophys. Res.},
+  \textbf{109~(C01001)}, \doi{10.1029/2003JC001866}.
+
+\bibitem[{Haapala(2000)}]{Haapala00}
+Haapala, J., 2000: On the modelling of ice-thickness redistribution.
+  \textit{Journal of Glaciology}, \textbf{46~(154)}, 427--437.
+
+\bibitem[{Hibler(1979)}]{Hibler79}
+Hibler, W.~D., 1979: A dynamic thermodynamic sea ice model. \textit{Journal of
+  Physical Oceanography}, \textbf{9}, 815--846.
+
+\bibitem[{Hibler(1980)}]{Hibler80}
+---{}---{}---, 1980: Modeling a variable thickness sea ice cover.
+  \textit{Monthly Weather Review}, \textbf{108}, 1943--1973.
+
+\bibitem[{H{{\o{}}}yland(2002)}]{Hoyland02}
+H{{\o{}}}yland, 2002: Consolidation of first-year sea ice ridges.
+  \textit{Journal of Geophysical Research}, \textbf{107}, 3062,
+  \doi{doi:10.1029/2000JC000526.}
+
+\bibitem[{{IOC, SCOR and IAPSO}(2010)}]{TEOS-10}
+{IOC, SCOR and IAPSO}, 2010: The international thermodynamic equation of
+  seawater - 2010: Calculation and use of thermodynamic properties.
+  Intergovernmental Oceanographic Commission, Manuals and Guides~56, UNESCO.
+
+\bibitem[{Kovacs(1996)}]{Kovacs96}
+Kovacs, A., 1996: Sea ice, part i. bulk salinity versus ice floe thickness.
+  Tech. Rep. 96-7, Cold Regions Research and Engineering Laboratory, Hanover,
+  New Hampshire.
+
+\bibitem[{Kreyscher et~al.(1999)Kreyscher, Harder, Lemke, and
+  Flato}]{Kreyscheretal00}
+Kreyscher, M., M.~Harder, P.~Lemke, and G.~M. Flato, 1999: Results of the sea
+  ice model intercomparison project: Evaluation of sea ice rheology scheme for
+  use in climate simulations. \textit{Jounal of Geophysical Research},
+  \textbf{105}, 11,299--11,320.
+
+\bibitem[{Lecomte et~al.(2015)Lecomte, Fichefet, Flocco, Schroeder, and
+  Vancoppenolle}]{Lecomteetal15}
+Lecomte, O., T.~Fichefet, D.~Flocco, D.~Schroeder, and M.~Vancoppenolle, 2015:
+  Interactions between wind-blown snow redistribution and melt ponds in a
+  coupled ocean-sea ice model. \textit{Ocean Modelling}, \textbf{87}, 67--80.
+
+\bibitem[{Lecomte et~al.(2011)Lecomte, Fichefet, Vancoppenolle, and
+  Nicolaus}]{Lecomteetal11}
+Lecomte, O., T.~Fichefet, M.~Vancoppenolle, and M.~Nicolaus, 2011: A new snow
+  thermodynamic scheme for large-scale sea-ice models. \textit{Annals of
+  Glaciology}, \textbf{52}, 337--346.
+
+\bibitem[{Lepp\"aranta(2005)}]{Lepparanta05}
+Lepp\"aranta, M., 2005: \textit{The Drift of Sea ice}. Springer, 266 pp. pp.
+
+\bibitem[{Lepp\"aranta et~al.(1995)Lepp\"aranta, Lensu, and
+  Witch}]{Lepparantaetal95}
+Lepp\"aranta, M., M.~Lensu, and P.~K.~B. Witch, 1995: The life story of a
+  first-year sea ice ridge. \textit{Cold Regions Science and Technology},
+  \textbf{23(3)}, 279--290.
+
+\bibitem[{Lipscomb(2001)}]{Lipscomb01}
+Lipscomb, W.~H., 2001: Remapping the thickness distribution in sea ice models.
+  \textit{Journal of Geophysical Research}, \textbf{106}, 13\,989--14\,000.
+
+\bibitem[{Lipscomb et~al.(2007)Lipscomb, Hunke, Maslowski, and
+  Jackaci}]{Lipscombetal07}
+Lipscomb, W.~H., E.~Hunke, W.~Maslowski, and J.~Jackaci, 2007: Ridging,
+  strength, and stability in high-resolution sea ice models. \textit{Journal of
+  Geophysical Research}, \textbf{112}, C03S9, \doi{doi:10.1029/2005JC003355}.
+
+\bibitem[{Massonnet et~al.(2018)Massonnet, Barth\'el\'emy, Worou, Fichefet,
+  Vancoppenolle, and Rousset}]{Massonnetetal18b}
+Massonnet, F., A.~Barth\'el\'emy, K.~Worou, T.~Fichefet, M.~Vancoppenolle, and
+  C.~Rousset, 2018: Insights on the discretization of the ice thickness
+  distribution in large-scale sea ice models. \textit{submitted}.
+
+\bibitem[{Maykut(1986)}]{Maykut86}
+Maykut, G.~A., 1986: \textit{The Geophysics of Sea Ice}, NATO ASI Series.
+  Series B, Physics., Vol. 146, chap. The surface heat and mass balance,
+  395--463. Plenum Press, New York.
+
+\bibitem[{Maykut and McPhee(1995)}]{MaykutMcPhee95}
+Maykut, G.~A. and M.~G. McPhee, 1995: Solar heating of the arctic mixed layer.
+  \textit{Journal of Geophysical Research}, \textbf{100}, 24,691--24,703.
+
+\bibitem[{Maykut and Thorndike(1973)}]{MaykutThorndike73}
+Maykut, G.~A. and A.~S. Thorndike, 1973: An approach to coupling the dynamics
+  and thermodynamics of arctic sea ice. \textit{AIDJEX Bulletin}, \textbf{21},
+  23--29.
+
+\bibitem[{Maykut and Untersteiner(1971)}]{MaykutUntersteiner71}
+Maykut, G.~A. and N.~Untersteiner, 1971: Some results from a time-dependent
+  thermodynamic model of sea ice. \textit{Journal of Geophysical Research},
+  \textbf{76}, 1550--1575.
+
+\bibitem[{McPhee(1992)}]{McPhee92}
+McPhee, M., 1992: Turbulent heat flux in the upper ocean under sea ice.
+  \textit{Journal of Geophysical Research}, \textbf{97(C4)}, 5365--5379.
+
+\bibitem[{Mellor and Kantha(1989)}]{MellorKantha89}
+Mellor, G.~L. and L.~Kantha, 1989: An ice-ocean coupled model. \textit{Journal
+  of Geophysical Research}, \textbf{94}, 10,937--10,954.
+
+\bibitem[{Notz et~al.(2016)Notz, Jahn, Holland, Hunke, Massonnet, Stroeve,
+  Tremblay, and Vancoppenolle}]{Notzetal16}
+Notz, D., A.~Jahn, M.~Holland, E.~Hunke, F.~Massonnet, J.~Stroeve, B.~Tremblay,
+  and M.~Vancoppenolle, 2016: {The CMIP6 Sea-Ice Model Intercomparison Project
+  (SIMIP): understanding sea ice through climate-model simulations}.
+  \textit{Geoscientific Model Development}, \textbf{9 (9)}, 3427--3446.
+
+\bibitem[{Parmerter(1975)}]{Parmerter75}
+Parmerter, R.~R., 1975: A model of simple rafting in sea ice. \textit{Journal
+  of Geophysical Research}, \textbf{80}, 1948--1952.
+
+\bibitem[{Perovich et~al.(2002)Perovich, Grenfell, Light, and
+  Hobbs}]{Perovichetal02alb}
+Perovich, D.~K., T.~C. Grenfell, B.~Light, and P.~V. Hobbs, 2002: Seasonal
+  evolution of the albedo of multiyear arctic sea ice. \textit{Journal of
+  Geophysical Research}, \textbf{107}, 8044, doi:10.1029/2000JC000\,438,
+  \doi{doi:10.1029/2000JC000438}.
+
+\bibitem[{Prather(1986)}]{Prather86}
+Prather, M.~J., 1986: Numerical advection by conservation of second-order
+  moments. \textit{Journal of Geophysical Research}, \textbf{91}, 6671--6681.
+
+\bibitem[{Pringle et~al.(2007)Pringle, Eicken, Trodahl, and
+  Backstrom}]{Pringleetal07}
+Pringle, D.~J., H.~Eicken, H.~J. Trodahl, and L.~Backstrom, 2007: Thermal
+  conductivity of landfast antarctic and arctic sea ice. \textit{Journal of
+  Geophysical Research}, \textbf{112}, C04017, \doi{10.1029/2006JC003641}.
+
+\bibitem[{Schmidt et~al.(2004)Schmidt, Bitz, Mikolajewicz, and
+  Tremblay}]{Schmidtetal04}
+Schmidt, G., C.~M. Bitz, U.~Mikolajewicz, and L.-B. Tremblay, 2004: Ice-ocean
+  boundary conditions for coupled models. \textit{Ocean Modelling}, \textbf{7},
+  59--74, \doi{doi:10.1016/S1463-5003(03)00030-1}.
+
+\bibitem[{Shine and Henderson-Sellers(1985)}]{ShineHenderson85}
+Shine, K.~P. and A.~Henderson-Sellers, 1985: The sensitivity of a thermodynamic
+  sea ice model to changes in surface albedo parameterization. \textit{Journal
+  of Geophysical Research}, \textbf{90}, 2243--2250.
+
+\bibitem[{Thorndike et~al.(1975)Thorndike, Rothrock, Maykut, and
+  Colony}]{Thorndikeetal75}
+Thorndike, A.~S., D.~A. Rothrock, G.~A. Maykut, and R.~Colony, 1975: The
+  thickness distribution of sea ice. \textit{Journal of Geophysical Research},
+  \textbf{80}, 4501--4513.
+
+\bibitem[{Tuhkuri and Lensu(2002)}]{TuhkuriLensu02}
+Tuhkuri, J. and M.~Lensu, 2002: Laboratory tests on ridging and rafting of ice
+  sheets. \textit{Journal of Geophysical Research}, \textbf{107}, 3125,
+  doi:10.1029/2001JC000\,848, \doi{doi:10.1029/2001JC000848}.
+
+\bibitem[{Vancoppenolle et~al.(2009)Vancoppenolle, Fichefet, and
+  Goosse}]{Vancoppenolleetal09b}
+Vancoppenolle, M., T.~Fichefet, and H.~Goosse, 2009: {Simulating the mass
+  balance and salinity of Arctic and Antarctic sea ice. 2. Sensitivity to
+  salinity processes}. \textit{Ocean Modelling}, \textbf{27 (1--2)}, 54--69,
+  \doi{10.1016/j.ocemod.2008.11.003}.
+
+\bibitem[{Weiss(2013)}]{Weiss13}
+Weiss, J., 2013: \textit{Drift, Deformation and Fracture of Sea Ice: A
+  Perspective across Scales}. Springer.
+
+\bibitem[{Worster(1992)}]{Worster92}
+Worster, M.~G., 1992: \textit{Interactive Dynamics of Convection and
+  Solidification}, chap. The dynamics of mushy layers, 113--138. Kluwer.
+
+\bibitem[{Zhang and Rothrock(2001)}]{ZhangRothrock01}
+Zhang, J. and D.~Rothrock, 2001: A thickness and enthalpy distribution sea ice
+  model. \textit{Journal of Physical Oceanography}, \textbf{31}, 2986--3001.
+
+\bibitem[{Zuo and Oerlemans(1996)}]{ZuoOerlemans96}
+Zuo, Z. and J.~Oerlemans, 1996: Modelling albedo and specific balance of the
+  greenland ice sheet: calculations for the sondre stromfjord transect.
+  \textit{Journal of Glaciology}, \textbf{42}, 305--317.
+
+\end{thebibliography}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_main/SI3_manual.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_main/SI3_manual.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_main/SI3_manual.tex	(revision 10845)
@@ -0,0 +1,129 @@
+%% ==============================================================================
+%% NEMO_book.tex: build the NEMO Reference Manual
+%% ==============================================================================
+
+
+%% Run ../build_NEMO_manual.sh to generate the PDF version (with pdflatex)
+
+%% Document layout
+%% ==============================================================================
+
+\documentclass[a4paper, 11pt]{book}
+
+
+%% Preamble
+%% ==============================================================================
+
+\usepackage{../../tex_main/NEMO_manual}
+\usepackage{../../tex_sub/NEMO_minted}
+
+\usepackage{bm}
+
+\makeindex
+
+
+%% 
+%% ==============================================================================
+
+\begin{document}
+
+
+%% Frontmatter
+%% ==============================================================================
+
+%% Frontpage
+
+\title{
+\vspace{1.5cm}
+\rule{400pt}{1.5pt} \\
+\vspace{0.45cm}
+ {\Huge SI$^3$ \\}
+ {\Large -- Sea Ice modelling Integrated Initiative --}\\
+ {\Large The NEMO Sea Ice engine}\\ 
+\date{\today}
+\author{The NEMO Sea Ice Working Group \\
+\\
+\small \textit{Note du P\^ole de mod\'elisation de l'Institut Pierre-Simon Laplace No 31} \\
+\small \textcolor{red}{First draft documentation. Needs completion and rewriting}. \\
+\\
+ISSN No 1288-1619\\
+}
+\rule{400pt}{1.5pt}
+}
+
+\maketitle						% generate the title
+
+\frontmatter
+
+%% ToC i.e. Table of Contents
+
+\dominitoc
+\tableofcontents
+
+
+%% Mainmatter
+%% ==============================================================================
+
+\mainmatter
+
+\subfile{../tex_sub/todolist}
+
+%% Abstract - Foreword
+
+\subfile{../tex_sub/abstract_foreword}
+
+% Introduction
+
+\subfile{../tex_sub/introduction}
+
+% Chapters
+
+\subfile{../tex_sub/chap_model_basics}
+
+\subfile{../tex_sub/chap_domain}
+
+\subfile{../tex_sub/chap_dynamics}
+
+\subfile{../tex_sub/chap_transport}
+
+\subfile{../tex_sub/chap_ridging_rafting}
+
+\subfile{../tex_sub/chap_radiative_transfer}
+
+\subfile{../tex_sub/chap_thermo}
+
+\subfile{../tex_sub/chap_interfaces}
+
+\subfile{../tex_sub/chap_output_diagnostics}
+
+\subfile{../tex_sub/chap_single_category_use}
+
+\subfile{../tex_sub/chap_bdy_agrif}
+
+\subfile{../tex_sub/chap_miscellaneous}
+
+%% Appendix
+
+%\appendix
+%% Backmatter
+%% ==============================================================================
+
+\backmatter
+
+%% Bibliography
+
+\cleardoublepage
+\phantomsection
+\addcontentsline{toc}{chapter}{Bibliography}
+
+\bibliography{../tex_main/SI3_manual.bbl}
+%\bibliography{/Users/ioulianikolskaia/Boulot/LIBRARY/BIBTEX/Martin_biblio.bib} 
+
+%% Index
+
+\cleardoublepage
+\phantomsection
+\addcontentsline{toc}{chapter}{Index}
+\printindex
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/abstract_foreword.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/abstract_foreword.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/abstract_foreword.tex	(revision 10845)
@@ -0,0 +1,44 @@
+
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Abstract
+% ================================================================
+
+\chapter*{Abstract}
+
+SI$^3$ (Sea Ice modelling Integrated Initiative) is the sea ice engine of NEMO (Nucleus for European Modelling of the Ocean). It is adapted to regional and global sea ice and climate problems. It is intended to be a flexible tool for studying sea ice and its interactions with the other components of the Earth System over a wide range of space and time scales. SI$^3$ is based on the Arctic Ice Dynamics Joint EXperiment (AIDJEX) framework \citep{Coonetal74}, combining the ice thickness distribution framework, the conservation of horizontal momentum, an elastic-viscous plastic rheology, and energy-conserving halo-thermodynamics. Prognostic variables are the two-dimensional horizontal velocity field, ice volume, area, enthalpy, salt content, snow volume and enthalpy. In the horizontal direction, the model uses a curvilinear orthogonal grid. In the vertical direction, the model uses equally-spaced layers. In thickness space, the model uses thickness categories with prescribed boundaries. Various physical and numerical choices are available to describe sea ice physics. SI$^3$ is interfaced with the NEMO ocean engine, and, via the OASIS coupler, with several atmospheric general circulation models. It also supports two-way grid embedding via the AGRIF software.
+
+%\vspace{-40pt}
+
+{\small
+} 
+
+% ================================================================
+% Disclaimer
+% ================================================================
+\chapter*{Disclaimer}
+
+Like all components of NEMO, the sea ice component is developed under the \href{http://www.cecill.info/}{CECILL license}, 
+which is a French adaptation of the GNU GPL (General Public License). Anyone may use it 
+freely for research purposes, and is encouraged to communicate back to the NEMO team 
+its own developments and improvements. The model and the present document have been 
+made available as a service to the community. We cannot certify that the code and its manual 
+are free of errors. Bugs are inevitable and some have undoubtedly survived the testing phase. 
+Users are encouraged to bring them to our attention. The authors assume no responsibility 
+for problems, errors, or incorrect usage of NEMO.
+
+\vspace{1cm}
+SI3 reference in papers and other publications is as follows:
+\vspace{0.5cm}
+
+The NEMO Sea Ice Working Group, 2018: SI$^3$ -- Sea Ice modelling Integrated Initiative -- The NEMO Sea Ice Engine, \textit{Note du P\^ole de mod\'{e}lisation}, Institut Pierre-Simon Laplace (IPSL), France, 
+No XX, ISSN No 1288-1619.\\
+
+\vspace{0.5cm}
+Additional information can be found on \href{http://www.nemo-ocean.eu/}{www.nemo-ocean.eu}.
+\vspace{0.5cm}
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/ametsoc.bst
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/ametsoc.bst	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/ametsoc.bst	(revision 10845)
@@ -0,0 +1,1662 @@
+%%
+%% This is file `ametsoc7.bst',
+%% generated with the docstrip utility.
+%%
+%% The original source files were:
+%%
+%% merlin.mbs  (with options: `ay,nat,nm-rev1,ed-rev,jnrlst,nmdash,nmd-3,dt-beg,yr-com,yrp-col,thtit-a,vol-2bf,vnum-sp,volp-com,pgsep-s,num-xser,numser,ser-vol,bkpg-x,pg-bk,pre-edn,agu-doi,doi,edparc,bkedcap,edby,edbyy,blk-tit,in-x,pp,ed,abr,ednx,ord,jabr,eprint,url,url-nt,em-it,nfss,')
+%% ----------------------------------------
+%% *** Bibliography style file for ALL AMS Journals...version 1.0  ***
+%% *** Brian Papa - American Meteorological Society ***
+%% 
+%% Copyright 1994-2004 Patrick W Daly
+ % ===============================================================
+ % IMPORTANT NOTICE:
+ % This bibliographic style (bst) file has been generated from one or
+ % more master bibliographic style (mbs) files, listed above.
+ %
+ % This generated file can be redistributed and/or modified under the terms
+ % of the LaTeX Project Public License Distributed from CTAN
+ % archives in directory macros/latex/base/lppl.txt; either
+ % version 1 of the License, or any later version.
+ % ===============================================================
+ % Name and version information of the main mbs file:
+ % \ProvidesFile{merlin.mbs}[2004/02/09 4.13 (PWD, AO, DPC)]
+ %   For use with BibTeX version 0.99a or later
+ %-------------------------------------------------------------------
+ % This bibliography style file is intended for texts in ENGLISH
+ % This is an author-year citation style bibliography. As such, it is
+ % non-standard LaTeX, and requires a special package file to function properly.
+ % Such a package is    natbib.sty   by Patrick W. Daly
+ % The form of the \bibitem entries is
+ %   \bibitem[Jones et al.(1990)]{key}...
+ %   \bibitem[Jones et al.(1990)Jones, Baker, and Smith]{key}...
+ % The essential feature is that the label (the part in brackets) consists
+ % of the author names, as they should appear in the citation, with the year
+ % in parentheses following. There must be no space before the opening
+ % parenthesis!
+ % With natbib v5.3, a full list of authors may also follow the year.
+ % In natbib.sty, it is possible to define the type of enclosures that is
+ % really wanted (brackets or parentheses), but in either case, there must
+ % be parentheses in the label.
+ % The \cite command functions as follows:
+ %   \citet{key} ==>>                Jones et al. (1990)
+ %   \citet*{key} ==>>               Jones, Baker, and Smith (1990)
+ %   \citep{key} ==>>                (Jones et al., 1990)
+ %   \citep*{key} ==>>               (Jones, Baker, and Smith, 1990)
+ %   \citep[chap. 2]{key} ==>>       (Jones et al., 1990, chap. 2)
+ %   \citep[e.g.][]{key} ==>>        (e.g. Jones et al., 1990)
+ %   \citep[e.g.][p. 32]{key} ==>>   (e.g. Jones et al., p. 32)
+ %   \citeauthor{key} ==>>           Jones et al.
+ %   \citeauthor*{key} ==>>          Jones, Baker, and Smith
+ %   \citeyear{key} ==>>             1990
+ %---------------------------------------------------------------------
+
+ENTRY
+  { address
+    archive
+    author
+    booktitle
+    chapter
+    doi
+    edition
+    editor
+    eid
+    eprint
+    howpublished
+    department
+    institution
+    journal
+    key
+    month
+    note
+    number
+    organization
+    pages
+    publisher
+    school
+    series
+    title
+    type
+    url
+    volume
+    year
+  }
+  {}
+  { label extra.label sort.label short.list }
+INTEGERS { output.state before.all mid.sentence after.sentence after.block }
+FUNCTION {init.state.consts}
+{ #0 'before.all :=
+  #1 'mid.sentence :=
+  #2 'after.sentence :=
+  #3 'after.block :=
+}
+STRINGS { s t}
+FUNCTION {output.nonnull}
+{ 's :=
+  output.state mid.sentence =
+    { ", " * write$ }
+    { output.state after.block =
+        { add.period$ write$
+          newline$
+          "\newblock " write$
+        }
+        { output.state before.all =
+            'write$
+            { add.period$ " " * write$ }
+          if$
+        }
+      if$
+      mid.sentence 'output.state :=
+    }
+  if$
+  s
+}
+FUNCTION {output}
+{ duplicate$ empty$
+    'pop$
+    'output.nonnull
+  if$
+}
+FUNCTION {output.check}
+{ 't :=
+  duplicate$ empty$
+    { pop$ "empty " t * " in " * cite$ * warning$ }
+    'output.nonnull
+  if$
+}
+FUNCTION {fin.entry}
+{ add.period$
+  write$
+  newline$
+}
+
+FUNCTION {new.block}
+{ output.state before.all =
+    'skip$
+    { after.block 'output.state := }
+  if$
+}
+FUNCTION {new.sentence}
+{ output.state after.block =
+    'skip$
+    { output.state before.all =
+        'skip$
+        { after.sentence 'output.state := }
+      if$
+    }
+  if$
+}
+FUNCTION {add.blank}
+{  " " * before.all 'output.state :=
+}
+
+FUNCTION {date.block}
+{
+  ":" *
+  add.blank
+}
+
+FUNCTION {not}
+{   { #0 }
+    { #1 }
+  if$
+}
+FUNCTION {and}
+{   'skip$
+    { pop$ #0 }
+  if$
+}
+FUNCTION {or}
+{   { pop$ #1 }
+    'skip$
+  if$
+}
+FUNCTION {new.block.checkb}
+{ empty$
+  swap$ empty$
+  and
+    'skip$
+    'new.block
+  if$
+}
+FUNCTION {field.or.null}
+{ duplicate$ empty$
+    { pop$ "" }
+    'skip$
+  if$
+}
+FUNCTION {emphasize}
+{ duplicate$ empty$
+    { pop$ "" }
+    { "\textit{" swap$ * "}" * }
+  if$
+}
+FUNCTION {bolden}
+{ duplicate$ empty$
+    { pop$ "" }
+    { "\textbf{" swap$ * "}" * }
+  if$
+}
+FUNCTION {tie.or.space.prefix}
+{ duplicate$ text.length$ #3 <
+    { "~" }
+    { " " }
+  if$
+  swap$
+}
+
+FUNCTION {capitalize}
+{ "u" change.case$ "t" change.case$ }
+
+FUNCTION {space.word}
+{ " " swap$ * " " * }
+ % Here are the language-specific definitions for explicit words.
+ % Each function has a name bbl.xxx where xxx is the English word.
+ % The language selected here is ENGLISH
+FUNCTION {bbl.and}
+{ "and"}
+
+FUNCTION {bbl.etal}
+{ "et~al." }
+
+FUNCTION {bbl.editors}
+{ "Eds." }
+
+FUNCTION {bbl.editor}
+{ "Ed." }
+
+FUNCTION {bbl.edby}
+{ "edited by" }
+
+FUNCTION {bbl.edition}
+{ "ed." }
+
+FUNCTION {bbl.volume}
+{ "Vol." }
+
+FUNCTION {bbl.of}
+{ "of" }
+
+FUNCTION {bbl.number}
+{ "No." }
+
+FUNCTION {bbl.nr}
+{ "No." }
+
+FUNCTION {bbl.in}
+{ "in" }
+
+FUNCTION {bbl.pages}
+{ "pp." }
+
+FUNCTION {bbl.page}
+{ "p." }
+
+FUNCTION {bbl.chapter}
+{ "chap." }
+
+FUNCTION {bbl.techrep}
+{ "Tech. Rep." }
+
+FUNCTION {bbl.mthesis}
+{ "Master's thesis" }
+
+FUNCTION {bbl.phdthesis}
+{ "Ph.D. thesis" }
+
+FUNCTION {bbl.first}
+{ "1st" }
+
+FUNCTION {bbl.second}
+{ "2d" }
+
+FUNCTION {bbl.third}
+{ "3d" }
+
+FUNCTION {bbl.fourth}
+{ "4th" }
+
+FUNCTION {bbl.fifth}
+{ "5th" }
+
+FUNCTION {bbl.st}
+{ "st" }
+
+FUNCTION {bbl.nd}
+{ "d" }
+
+FUNCTION {bbl.rd}
+{ "d" }
+
+FUNCTION {bbl.th}
+{ "th" }
+
+MACRO {jan} {"Jan."}
+
+MACRO {feb} {"Feb."}
+
+MACRO {mar} {"Mar."}
+
+MACRO {apr} {"Apr."}
+
+MACRO {may} {"May"}
+
+MACRO {jun} {"Jun."}
+
+MACRO {jul} {"Jul."}
+
+MACRO {aug} {"Aug."}
+
+MACRO {sep} {"Sep."}
+
+MACRO {oct} {"Oct."}
+
+MACRO {nov} {"Nov."}
+
+MACRO {dec} {"Dec."}
+
+FUNCTION {eng.ord}
+{ duplicate$ "1" swap$ *
+  #-2 #1 substring$ "1" =
+     { bbl.th * }
+     { duplicate$ #-1 #1 substring$
+       duplicate$ "1" =
+         { pop$ bbl.st * }
+         { duplicate$ "2" =
+             { pop$ bbl.nd * }
+             { "3" =
+                 { bbl.rd * }
+                 { bbl.th * }
+               if$
+             }
+           if$
+          }
+       if$
+     }
+   if$
+}
+
+MACRO {acmcs} {"ACM Comput. Surv."}
+
+MACRO {acta} {"Acta Inf."}
+
+MACRO {cacm} {"Commun. ACM"}
+
+MACRO {ibmjrd} {"IBM J. Res. Dev."}
+
+MACRO {ibmsj} {"IBM Syst.~J."}
+
+MACRO {ieeese} {"IEEE Trans. Software Eng."}
+
+MACRO {ieeetc} {"IEEE Trans. Comput."}
+
+MACRO {ieeetcad}
+ {"IEEE Trans. Comput. Aid. Des."}
+
+MACRO {ipl} {"Inf. Process. Lett."}
+
+MACRO {jacm} {"J.~ACM"}
+
+MACRO {jcss} {"J.~Comput. Syst. Sci."}
+
+MACRO {scp} {"Sci. Comput. Program."}
+
+MACRO {sicomp} {"SIAM J. Comput."}
+
+MACRO {tocs} {"ACM Trans. Comput. Syst."}
+
+MACRO {tods} {"ACM Trans. Database Syst."}
+
+MACRO {tog} {"ACM Trans. Graphic."}
+
+MACRO {toms} {"ACM Trans. Math. Software"}
+
+MACRO {toois} {"ACM Trans. Office Inf. Syst."}
+
+MACRO {toplas} {"ACM Trans. Progr. Lang. Syst."}
+
+MACRO {tcs} {"Theor. Comput. Sci."}
+
+FUNCTION {bibinfo.check}
+{ swap$
+  duplicate$ missing$
+    {
+      pop$ pop$
+      ""
+    }
+    { duplicate$ empty$
+        {
+          swap$ pop$
+        }
+        { swap$
+          pop$
+        }
+      if$
+    }
+  if$
+}
+FUNCTION {bibinfo.warn}
+{ swap$
+  duplicate$ missing$
+    {
+      swap$ "missing " swap$ * " in " * cite$ * warning$ pop$
+      ""
+    }
+    { duplicate$ empty$
+        {
+          swap$ "empty " swap$ * " in " * cite$ * warning$
+        }
+        { swap$
+          pop$
+        }
+      if$
+    }
+  if$
+}
+FUNCTION {format.eprint}
+{ eprint duplicate$ empty$
+    'skip$
+    { "\eprint"
+      archive empty$
+        'skip$
+        { "[" * archive * "]" * }
+      if$
+      "{" * swap$ * "}" *
+    }
+  if$
+}
+STRINGS  { bibinfo}
+INTEGERS { nameptr namesleft numnames }
+
+FUNCTION {format.names}
+{ 'bibinfo :=
+  duplicate$ empty$ 'skip$ {
+  's :=
+  "" 't :=
+  #1 'nameptr :=
+  s num.names$ 'numnames :=
+  numnames 'namesleft :=
+    { namesleft #0 > }
+    { s nameptr
+      duplicate$ #1 >
+        { "{f.~}{vv~}{ll}{, jj}" }
+        { "{vv~}{ll}{, f.}{, jj}" }
+      if$
+      format.name$
+      bibinfo bibinfo.check
+      't :=
+      nameptr #1 >
+        {
+          namesleft #1 >
+            { ", " * t * }
+            {
+              numnames #2 >
+                { "," * }
+                'skip$
+              if$
+              s nameptr "{ll}" format.name$ duplicate$ "others" =
+                { 't := }
+                { pop$ }
+              if$
+              t "others" =
+                {
+                  " " * bbl.etal *
+                }
+                {
+                  bbl.and
+                  space.word * t *
+                }
+              if$
+            }
+          if$
+        }
+        't
+      if$
+      nameptr #1 + 'nameptr :=
+      namesleft #1 - 'namesleft :=
+    }
+  while$
+  } if$
+}
+FUNCTION {format.names.ed}
+{
+  format.names
+}
+FUNCTION {format.key}
+{ empty$
+    { key field.or.null }
+    { "" }
+  if$
+}
+
+FUNCTION {format.authors}
+{ author "author" format.names
+}
+FUNCTION {get.bbl.editor}
+{ editor num.names$ #1 > 'bbl.editors 'bbl.editor if$ }
+
+FUNCTION {format.editors}
+{ editor "editor" format.names duplicate$ empty$ 'skip$
+    {
+      "," *
+      " " *
+      get.bbl.editor
+      capitalize
+   "(" swap$ * ")" *
+      *
+    }
+  if$
+}
+FUNCTION {format.book.pages}
+{ pages "pages" bibinfo.check
+  duplicate$ empty$ 'skip$
+    { " " * bbl.pages * }
+  if$
+}
+FUNCTION {format.doi}
+{ doi "doi" bibinfo.check
+  duplicate$ empty$ 'skip$
+    {
+      "\doi{" swap$ * "}" *
+    }
+  if$
+}
+FUNCTION {format.note}
+{
+  url empty$
+    'skip$
+    { "\urlprefix\url{" url * "}" * output }
+  if$
+ note empty$
+    { "" }
+    { note #1 #1 substring$
+      duplicate$ "{" =
+        'skip$
+        { output.state mid.sentence =
+          { "l" }
+          { "u" }
+        if$
+        change.case$
+        }
+      if$
+      note #2 global.max$ substring$ * "note" bibinfo.check
+    }
+  if$
+}
+
+FUNCTION {format.title}
+{ title
+  duplicate$ empty$ 'skip$
+    { "t" change.case$ }
+  if$
+  "title" bibinfo.check
+}
+FUNCTION {format.full.names}
+{'s :=
+ "" 't :=
+  #1 'nameptr :=
+  s num.names$ 'numnames :=
+  numnames 'namesleft :=
+    { namesleft #0 > }
+    { s nameptr
+      "{vv~}{ll}" format.name$
+      't :=
+      nameptr #1 >
+        {
+          namesleft #1 >
+            { ", " * t * }
+            {
+              s nameptr "{ll}" format.name$ duplicate$ "others" =
+                { 't := }
+                { pop$ }
+              if$
+              t "others" =
+                {
+                  " " * bbl.etal *
+                }
+                {
+                  numnames #2 >
+                    { "," * }
+                    'skip$
+                  if$
+                  bbl.and
+                  space.word * t *
+                }
+              if$
+            }
+          if$
+        }
+        't
+      if$
+      nameptr #1 + 'nameptr :=
+      namesleft #1 - 'namesleft :=
+    }
+  while$
+}
+
+FUNCTION {author.editor.key.full}
+{ author empty$
+    { editor empty$
+        { key empty$
+            { cite$ #1 #3 substring$ }
+            'key
+          if$
+        }
+        { editor format.full.names }
+      if$
+    }
+    { author format.full.names }
+  if$
+}
+
+FUNCTION {author.key.full}
+{ author empty$
+    { key empty$
+         { cite$ #1 #3 substring$ }
+          'key
+      if$
+    }
+    { author format.full.names }
+  if$
+}
+
+FUNCTION {editor.key.full}
+{ editor empty$
+    { key empty$
+         { cite$ #1 #3 substring$ }
+          'key
+      if$
+    }
+    { editor format.full.names }
+  if$
+}
+
+FUNCTION {make.full.names}
+{ type$ "book" =
+  type$ "inbook" =
+  or
+    'author.editor.key.full
+    { type$ "proceedings" =
+        'editor.key.full
+        'author.key.full
+      if$
+    }
+  if$
+}
+
+FUNCTION {output.bibitem}
+{ newline$
+  "\bibitem[{" write$
+  label write$
+  ")" make.full.names duplicate$ short.list =
+     { pop$ }
+     { * }
+   if$
+  "}]{" * write$
+  cite$ write$
+  "}" write$
+  newline$
+  ""
+  before.all 'output.state :=
+}
+
+FUNCTION {if.digit}
+{ duplicate$ "0" =
+  swap$ duplicate$ "1" =
+  swap$ duplicate$ "2" =
+  swap$ duplicate$ "3" =
+  swap$ duplicate$ "4" =
+  swap$ duplicate$ "5" =
+  swap$ duplicate$ "6" =
+  swap$ duplicate$ "7" =
+  swap$ duplicate$ "8" =
+  swap$ "9" = or or or or or or or or or
+}
+FUNCTION {n.separate}
+{ 't :=
+  ""
+  #0 'numnames :=
+  { t empty$ not }
+  { t #-1 #1 substring$ if.digit
+      { numnames #1 + 'numnames := }
+      { #0 'numnames := }
+    if$
+    t #-1 #1 substring$ swap$ *
+    t #-2 global.max$ substring$ 't :=
+    numnames #5 =
+      { duplicate$ #1 #2 substring$ swap$
+        #3 global.max$ substring$
+        "\," swap$ * *
+      }
+      'skip$
+    if$
+  }
+  while$
+}
+FUNCTION {n.dashify}
+{
+  n.separate
+  't :=
+  ""
+    { t empty$ not }
+    { t #1 #1 substring$ "-" =
+        { t #1 #2 substring$ "--" = not
+            { "--" *
+              t #2 global.max$ substring$ 't :=
+            }
+            {   { t #1 #1 substring$ "-" = }
+                { "-" *
+                  t #2 global.max$ substring$ 't :=
+                }
+              while$
+            }
+          if$
+        }
+        { t #1 #1 substring$ *
+          t #2 global.max$ substring$ 't :=
+        }
+      if$
+    }
+  while$
+}
+
+FUNCTION {word.in}
+{ "" }
+
+FUNCTION {format.date}
+{ year "year" bibinfo.check duplicate$ empty$
+    {
+      "empty year in " cite$ * "; set to ????" * warning$
+       pop$ "????"
+    }
+    'skip$
+  if$
+  extra.label *
+  before.all 'output.state :=
+  ", " swap$ *
+}
+FUNCTION {format.btitle}
+{ title "title" bibinfo.check
+  duplicate$ empty$ 'skip$
+    {
+      emphasize
+    }
+  if$
+}
+FUNCTION {either.or.check}
+{ empty$
+    'pop$
+    { "can't use both " swap$ * " fields in " * cite$ * warning$ }
+  if$
+}
+FUNCTION {format.bvolume}
+{ volume empty$
+    { "" }
+    { bbl.volume volume tie.or.space.prefix
+      "volume" bibinfo.check * *
+      series "series" bibinfo.check
+      duplicate$ empty$ 'pop$
+%%% NO ITALIC SERIES        { emphasize ", " * swap$ * }
+      {", " * swap$ * }
+      if$
+      "volume and number" number either.or.check
+    }
+  if$
+}
+FUNCTION {format.number.series}
+{ volume empty$
+    { number empty$
+        { series field.or.null }
+        { series empty$
+            { number "number" bibinfo.check }
+            { output.state mid.sentence =
+                { bbl.number }
+                { bbl.number capitalize }
+              if$
+              number tie.or.space.prefix "number" bibinfo.check * *
+              bbl.in space.word *
+              series "series" bibinfo.check *
+            }
+          if$
+        }
+      if$
+    }
+    { "" }
+  if$
+}
+FUNCTION {is.num}
+{ chr.to.int$
+  duplicate$ "0" chr.to.int$ < not
+  swap$ "9" chr.to.int$ > not and
+}
+
+FUNCTION {extract.num}
+{ duplicate$ 't :=
+  "" 's :=
+  { t empty$ not }
+  { t #1 #1 substring$
+    t #2 global.max$ substring$ 't :=
+    duplicate$ is.num
+      { s swap$ * 's := }
+      { pop$ "" 't := }
+    if$
+  }
+  while$
+  s empty$
+    'skip$
+    { pop$ s }
+  if$
+}
+
+FUNCTION {convert.edition}
+{ extract.num "l" change.case$ 's :=
+  s "first" = s "1" = or
+    { bbl.first 't := }
+    { s "second" = s "2" = or
+        { bbl.second 't := }
+        { s "third" = s "3" = or
+            { bbl.third 't := }
+            { s "fourth" = s "4" = or
+                { bbl.fourth 't := }
+                { s "fifth" = s "5" = or
+                    { bbl.fifth 't := }
+                    { s #1 #1 substring$ is.num
+                        { s eng.ord 't := }
+                        { edition 't := }
+                      if$
+                    }
+                  if$
+                }
+              if$
+            }
+          if$
+        }
+      if$
+    }
+  if$
+  t
+}
+
+FUNCTION {format.edition}
+{ edition duplicate$ empty$ 'skip$
+    {
+      convert.edition
+      output.state mid.sentence =
+        { "l" }
+        { "t" }
+      if$ change.case$
+      "edition" bibinfo.check
+      " " * bbl.edition *
+    }
+  if$
+}
+INTEGERS { multiresult }
+FUNCTION {multi.page.check}
+{ 't :=
+  #0 'multiresult :=
+    { multiresult not
+      t empty$ not
+      and
+    }
+    { t #1 #1 substring$
+      duplicate$ "-" =
+      swap$ duplicate$ "," =
+      swap$ "+" =
+      or or
+        { #1 'multiresult := }
+        { t #2 global.max$ substring$ 't := }
+      if$
+    }
+  while$
+  multiresult
+}
+FUNCTION {format.pages}
+{ pages duplicate$ empty$ 'skip$
+    { duplicate$ multi.page.check
+        {
+          n.dashify
+        }
+        {
+        }
+      if$
+      "pages" bibinfo.check
+    }
+  if$
+}
+FUNCTION {format.journal.pages}
+{ pages duplicate$ empty$ 'pop$
+    { swap$ duplicate$ empty$
+        { pop$ pop$ format.pages }
+        {
+          ", " *
+          swap$
+          n.dashify
+          "pages" bibinfo.check
+          *
+        }
+      if$
+    }
+  if$
+}
+FUNCTION {format.journal.eid}
+{ eid "eid" bibinfo.check
+  duplicate$ empty$ 'pop$
+    { swap$ duplicate$ empty$ 'skip$
+      {
+          ", " *
+      }
+      if$
+      swap$ *
+    }
+  if$
+}
+FUNCTION {format.vol.num.pages}
+{ volume field.or.null
+  duplicate$ empty$ 'skip$
+    {
+      "volume" bibinfo.check
+    }
+  if$
+  number "number" bibinfo.check duplicate$ empty$ 'skip$
+    {
+      swap$ duplicate$ empty$
+        { "there's a number but no volume in " cite$ * warning$ }
+        'skip$
+      if$
+      swap$
+      "~(" swap$ * ")" *
+    }
+  if$ *
+  bolden
+  eid empty$
+    { format.journal.pages }
+    { format.journal.eid }
+  if$
+}
+
+FUNCTION {format.chapter.pages}
+{ chapter empty$
+    'format.pages
+    { type empty$
+        { bbl.chapter }
+        { type "l" change.case$
+          "type" bibinfo.check
+        }
+      if$
+      chapter tie.or.space.prefix
+      "chapter" bibinfo.check
+      * *
+      pages empty$
+        'skip$
+        { ", " * format.pages * }
+      if$
+    }
+  if$
+}
+
+FUNCTION {format.booktitle}
+{
+  booktitle "booktitle" bibinfo.check
+  emphasize
+}
+FUNCTION {format.in.ed.booktitle}
+{ format.booktitle duplicate$ empty$ 'skip$
+    {
+      editor "editor" format.names.ed duplicate$ empty$ 'pop$
+        {
+          get.bbl.editor
+          swap$ "," *
+          " " * swap$ *
+          swap$
+          "," *
+          " " * swap$
+          * }
+      if$
+      word.in swap$ *
+    }
+  if$
+}
+FUNCTION {format.thesis.type}
+{ type duplicate$ empty$
+    'pop$
+    { swap$ pop$
+      "t" change.case$ "type" bibinfo.check
+    }
+  if$
+}
+FUNCTION {format.tr.number}
+{ number "number" bibinfo.check
+  type duplicate$ empty$
+    { pop$ bbl.techrep }
+    'skip$
+  if$
+  "type" bibinfo.check
+  swap$ duplicate$ empty$
+    { pop$ "t" change.case$ }
+    { tie.or.space.prefix * * }
+  if$
+}
+FUNCTION {format.article.crossref}
+{
+  word.in
+  " \cite{" * crossref * "}" *
+}
+FUNCTION {format.book.crossref}
+{ volume duplicate$ empty$
+    { "empty volume in " cite$ * "'s crossref of " * crossref * warning$
+      pop$ word.in
+    }
+    { bbl.volume
+      capitalize
+      swap$ tie.or.space.prefix "volume" bibinfo.check * * bbl.of space.word *
+    }
+  if$
+  " \cite{" * crossref * "}" *
+}
+FUNCTION {format.incoll.inproc.crossref}
+{
+  word.in
+  " \cite{" * crossref * "}" *
+}
+FUNCTION {format.org.or.pub}
+{ 't :=
+  ""
+  address empty$ t empty$ and
+    'skip$
+    {
+      t empty$
+        { address "address" bibinfo.check *
+        }
+        { t *
+          address empty$
+            'skip$
+            { ", " * address "address" bibinfo.check * }
+          if$
+        }
+      if$
+    }
+  if$
+}
+FUNCTION {format.publisher.address}
+{ publisher "publisher" bibinfo.warn format.org.or.pub
+}
+
+FUNCTION {format.organization.address}
+{ organization "organization" bibinfo.check format.org.or.pub
+}
+
+STRINGS {oldname}
+
+FUNCTION {name.or.dash}
+{ 's :=
+   oldname empty$
+     { s 'oldname := s }
+     { s oldname =
+         { "---{}---{}---" }
+         { s 'oldname := s }
+       if$
+     }
+   if$
+}
+
+FUNCTION {article}
+{ output.bibitem
+  format.authors "author" output.check
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title "title" output.check
+  new.sentence
+  crossref missing$
+    {
+      journal
+      "journal" bibinfo.check
+      emphasize
+      "journal" output.check
+      format.vol.num.pages output
+  format.doi output
+    }
+    { format.article.crossref output.nonnull
+      format.pages output
+    }
+  if$
+  format.note output
+  format.eprint output
+  fin.entry
+}
+FUNCTION {book}
+{ output.bibitem
+  author empty$
+    { format.editors "author and editor" output.check
+      editor format.key output
+      add.blank
+      name.or.dash
+    }
+    { format.authors output.nonnull
+      name.or.dash
+      crossref missing$
+        { "author and editor" editor either.or.check }
+        'skip$
+      if$
+    }
+  if$
+  format.date "year" output.check
+  date.block
+  format.btitle "title" output.check
+  crossref missing$
+    { format.bvolume output
+  new.sentence
+      format.edition output
+      format.number.series output
+      format.publisher.address output
+    }
+    {
+  new.sentence
+      format.book.crossref output.nonnull
+    }
+  if$
+  format.book.pages output
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+FUNCTION {booklet}
+{ output.bibitem
+  format.authors output
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title "title" output.check
+  new.sentence
+  howpublished "howpublished" bibinfo.check output
+  address "address" bibinfo.check output
+  format.book.pages output
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+
+FUNCTION {inbook}
+{ output.bibitem
+  author empty$
+    { format.editors "author and editor" output.check
+      editor format.key output
+      name.or.dash
+    }
+    { format.authors output.nonnull
+      name.or.dash
+      crossref missing$
+        { "author and editor" editor either.or.check }
+        'skip$
+      if$
+    }
+  if$
+  format.date "year" output.check
+  date.block
+  format.btitle "title" output.check
+  crossref missing$
+    {
+      format.bvolume output
+      format.chapter.pages "chapter and pages" output.check
+  new.sentence
+      format.edition output
+      format.number.series output
+      format.publisher.address output
+    }
+    {
+      format.chapter.pages "chapter and pages" output.check
+  new.sentence
+      format.book.crossref output.nonnull
+    }
+  if$
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+
+FUNCTION {incollection}
+{ output.bibitem
+  format.authors "author" output.check
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title "title" output.check
+  new.sentence
+  crossref missing$
+    { format.in.ed.booktitle "booktitle" output.check
+      format.publisher.address output      
+      format.bvolume output
+      format.number.series output
+      format.edition output
+      format.chapter.pages output
+    }
+    { format.incoll.inproc.crossref output.nonnull
+      format.chapter.pages output
+    }
+  if$
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+FUNCTION {inproceedings}
+{ output.bibitem
+  format.authors "author" output.check
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title "title" output.check
+  new.sentence
+  crossref missing$
+    { format.in.ed.booktitle "booktitle" output.check
+      publisher empty$
+        { format.organization.address output }
+        { organization "organization" bibinfo.check output
+          format.publisher.address output
+        }
+      if$
+      format.bvolume output
+      format.pages output
+      format.number.series output
+      howpublished "howpublished" bibinfo.check output
+    }
+    { format.incoll.inproc.crossref output.nonnull
+      format.pages output
+    }
+  if$
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+FUNCTION {conference} { inproceedings }
+FUNCTION {manual}
+{ output.bibitem
+  format.authors output
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.btitle "title" output.check
+  new.sentence
+  organization "organization" bibinfo.check output
+  address "address" bibinfo.check output
+  format.edition output
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+
+FUNCTION {mastersthesis}
+{ output.bibitem
+  format.authors "author" output.check
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title
+  "title" output.check
+  new.sentence
+  bbl.mthesis format.thesis.type output.nonnull
+  department "department" bibinfo.warn output
+  school "school" bibinfo.warn output
+  address "address" bibinfo.check output
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+
+FUNCTION {misc}
+{ output.bibitem
+  format.authors output
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title output
+  new.sentence
+  howpublished "howpublished" bibinfo.check output
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+FUNCTION {phdthesis}
+{ output.bibitem
+  format.authors "author" output.check
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title
+  "title" output.check
+  new.sentence
+  bbl.phdthesis format.thesis.type output.nonnull
+  school "school" bibinfo.warn output
+  address "address" bibinfo.check output
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+
+FUNCTION {proceedings}
+{ output.bibitem
+  format.editors output
+  editor format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.btitle "title" output.check
+  format.bvolume output
+  format.number.series output
+  publisher empty$
+    { format.organization.address output }
+    { organization "organization" bibinfo.check output
+      format.publisher.address output
+    }
+  if$
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+
+FUNCTION {techreport}
+{ output.bibitem
+  format.authors "author" output.check
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title
+  "title" output.check
+  new.sentence
+  format.tr.number output.nonnull
+  institution "institution" bibinfo.warn output
+  address "address" bibinfo.check output
+  format.doi output
+  format.note output
+  format.eprint output
+  fin.entry
+}
+
+FUNCTION {unpublished}
+{ output.bibitem
+  format.authors "author" output.check
+  author format.key output
+  name.or.dash
+  format.date "year" output.check
+  date.block
+  format.title "title" output.check
+  format.doi output
+  format.note "note" output.check
+  format.eprint output
+  fin.entry
+}
+
+FUNCTION {default.type} { misc }
+READ
+FUNCTION {sortify}
+{ purify$
+  "l" change.case$
+}
+INTEGERS { len }
+FUNCTION {chop.word}
+{ 's :=
+  'len :=
+  s #1 len substring$ =
+    { s len #1 + global.max$ substring$ }
+    's
+  if$
+}
+FUNCTION {format.lab.names}
+{ 's :=
+  "" 't :=
+  s #1 "{vv~}{ll}" format.name$
+  s num.names$ duplicate$
+  #2 >
+    { pop$
+      " " * bbl.etal *
+    }
+    { #2 <
+        'skip$
+        { s #2 "{ff }{vv }{ll}{ jj}" format.name$ "others" =
+            {
+              " " * bbl.etal *
+            }
+            { bbl.and space.word * s #2 "{vv~}{ll}" format.name$
+              * }
+          if$
+        }
+      if$
+    }
+  if$
+}
+
+FUNCTION {author.key.label}
+{ author empty$
+    { key empty$
+        { cite$ #1 #3 substring$ }
+        'key
+      if$
+    }
+    { author format.lab.names }
+  if$
+}
+
+FUNCTION {author.editor.key.label}
+{ author empty$
+    { editor empty$
+        { key empty$
+            { cite$ #1 #3 substring$ }
+            'key
+          if$
+        }
+        { editor format.lab.names }
+      if$
+    }
+    { author format.lab.names }
+  if$
+}
+
+FUNCTION {editor.key.label}
+{ editor empty$
+    { key empty$
+        { cite$ #1 #3 substring$ }
+        'key
+      if$
+    }
+    { editor format.lab.names }
+  if$
+}
+
+FUNCTION {calc.short.authors}
+{ type$ "book" =
+  type$ "inbook" =
+  or
+    'author.editor.key.label
+    { type$ "proceedings" =
+        'editor.key.label
+        'author.key.label
+      if$
+    }
+  if$
+  'short.list :=
+}
+
+FUNCTION {calc.label}
+{ calc.short.authors
+  short.list
+  "("
+  *
+  year duplicate$ empty$
+     { pop$ "????" }
+     { purify$ #-1 #4 substring$ }
+  if$
+  *
+  'label :=
+}
+
+FUNCTION {sort.format.names}
+{ 's :=
+  #1 'nameptr :=
+  ""
+  s num.names$ 'numnames :=
+  numnames 'namesleft :=
+    { namesleft #0 > }
+    { s nameptr
+      "{vv{ } }{ll{ }}{  f{ }}{  jj{ }}"
+      format.name$ 't :=
+      nameptr #1 >
+        {
+          "   "  *
+          namesleft #1 = t "others" = and
+            { "zzzzz" * }
+            { t sortify * }
+          if$
+        }
+        { t sortify * }
+      if$
+      nameptr #1 + 'nameptr :=
+      namesleft #1 - 'namesleft :=
+    }
+  while$
+}
+
+FUNCTION {sort.format.title}
+{ 't :=
+  "A " #2
+    "An " #3
+      "The " #4 t chop.word
+    chop.word
+  chop.word
+  sortify
+  #1 global.max$ substring$
+}
+FUNCTION {author.sort}
+{ author empty$
+    { key empty$
+        { "to sort, need author or key in " cite$ * warning$
+          ""
+        }
+        { key sortify }
+      if$
+    }
+    { author sort.format.names }
+  if$
+}
+FUNCTION {author.editor.sort}
+{ author empty$
+    { editor empty$
+        { key empty$
+            { "to sort, need author, editor, or key in " cite$ * warning$
+              ""
+            }
+            { key sortify }
+          if$
+        }
+        { editor sort.format.names }
+      if$
+    }
+    { author sort.format.names }
+  if$
+}
+FUNCTION {editor.sort}
+{ editor empty$
+    { key empty$
+        { "to sort, need editor or key in " cite$ * warning$
+          ""
+        }
+        { key sortify }
+      if$
+    }
+    { editor sort.format.names }
+  if$
+}
+FUNCTION {presort}
+{ calc.label
+  label sortify
+  "    "
+  *
+  type$ "book" =
+  type$ "inbook" =
+  or
+    'author.editor.sort
+    { type$ "proceedings" =
+        'editor.sort
+        'author.sort
+      if$
+    }
+  if$
+  #1 entry.max$ substring$
+  'sort.label :=
+  sort.label
+  *
+  "    "
+  *
+  title field.or.null
+  sort.format.title
+  *
+  #1 entry.max$ substring$
+  'sort.key$ :=
+}
+
+ITERATE {presort}
+SORT
+STRINGS { last.label next.extra }
+INTEGERS { last.extra.num number.label }
+FUNCTION {initialize.extra.label.stuff}
+{ #0 int.to.chr$ 'last.label :=
+  "" 'next.extra :=
+  #0 'last.extra.num :=
+  #0 'number.label :=
+}
+FUNCTION {forward.pass}
+{ last.label label =
+    { last.extra.num #1 + 'last.extra.num :=
+      last.extra.num int.to.chr$ 'extra.label :=
+    }
+    { "a" chr.to.int$ 'last.extra.num :=
+      "" 'extra.label :=
+      label 'last.label :=
+    }
+  if$
+  number.label #1 + 'number.label :=
+}
+FUNCTION {reverse.pass}
+{ next.extra "b" =
+    { "a" 'extra.label := }
+    'skip$
+  if$
+  extra.label 'next.extra :=
+  extra.label
+  duplicate$ empty$
+    'skip$
+    { "{\natexlab{" swap$ * "}}" * }
+  if$
+  'extra.label :=
+  label extra.label * 'label :=
+}
+EXECUTE {initialize.extra.label.stuff}
+ITERATE {forward.pass}
+REVERSE {reverse.pass}
+FUNCTION {bib.sort.order}
+{ sort.label
+  "    "
+  *
+  year field.or.null sortify
+  *
+  "    "
+  *
+  title field.or.null
+  sort.format.title
+  *
+  #1 entry.max$ substring$
+  'sort.key$ :=
+}
+ITERATE {bib.sort.order}
+SORT
+FUNCTION {begin.bib}
+{ preamble$ empty$
+    'skip$
+    { preamble$ write$ newline$ }
+  if$
+  "\begin{thebibliography}{" number.label int.to.str$ * "}" *
+  write$ newline$
+  "\providecommand{\natexlab}[1]{#1}"
+  write$ newline$
+  "\providecommand{\url}[1]{\texttt{#1}}"
+  write$ newline$
+  "\providecommand{\urlprefix}{URL }"
+  write$ newline$
+  "\expandafter\ifx\csname urlstyle\endcsname\relax"
+  write$ newline$
+  "  \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else"
+  write$ newline$
+  "  \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup \urlstyle{rm}\Url}\fi"
+  write$ newline$
+  "\providecommand{\eprint}[2][]{\url{#2}}"
+  write$ newline$
+}
+EXECUTE {begin.bib}
+EXECUTE {init.state.consts}
+ITERATE {call.type$}
+FUNCTION {end.bib}
+{ newline$
+  "\end{thebibliography}" write$ newline$
+}
+EXECUTE {end.bib}
+%% End of customized bst file
+%%
+%% End of file `ametsoc7.bst'.
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_bdy_agrif.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_bdy_agrif.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_bdy_agrif.tex	(revision 10845)
@@ -0,0 +1,16 @@
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 11 - BDY and AGRIF with SI$^3$
+% ================================================================
+
+\chapter{BDY and AGRIF with SI$^3$}
+\label{chap:REG}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_domain.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_domain.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_domain.tex	(revision 10845)
@@ -0,0 +1,106 @@
+
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 2 Ñ Domain
+% ================================================================
+
+\chapter{Time, space and thickness space domain}
+\label{chap:DOM}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+Having defined the model equations in previous Chapter, we need now to choose the numerical discretization.  In the present chapter, we provide a general description of the SI$^3$ discretization strategy, in terms of time, space and thickness, which is considered as an extra independent variable. 
+
+Sea ice state variables are typically expressed as:
+\begin{equation}
+X(ji,jj,\textcolor{gray}{jk},jl).
+\end{equation}
+$ji$ and $jj$ are x-y spatial indices, as in the ocean. $jk=1, ..., nlay\_i$ corresponds to the vertical coordinate system in sea ice (ice layers), and only applies to vertically-resolved quantities (ice enthalpy and salinity). $jl=1, ..., jpl$ corresponds to the ice categories, discretizing thickness space.
+
+\section{Time domain}
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : Time Stepping
+%
+\begin{figure}[ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=6cm,angle=-00]{../Figures/time_stepping.png}
+\caption{Schematic representation of time stepping in SI$^3$, assuming $nn\_fsbc=5$.}
+\label{ice_scheme}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+The sea ice time stepping is synchronized with that of the ocean. Because of the potentially large numerical cost of sea ice physics, in particular rheology, SI$^3$ can be called every nn\_fsbc time steps (namsbc in \textit{namelist\_ref}). The sea ice time step is therefore $rdt\_ice = rdt * nn\_fsbc$. In terms of quality, the best value for \textit{nn\_fsbc} is 1, providing full consistency between sea ice and oceanic fields. Larger values (typically 2 to 5) can be used but numerical instabilities can appear because of the progressive decoupling between the state of sea ice and that of the ocean, hence changing $nn\_fsbc$ must be done carefully.
+
+Ice dynamics (rheology, advection, ridging/rafting) and thermodynamics are called successively. To avoid pathological situations, thermodynamics were chosen to be applied on fields that have been updated by dynamics, in a somehow semi-implicit procedure.
+
+There are a few iterative / subcycling procedures throughout the code, notably for rheology, advection, ridging/ rafting and the diffusion of heat. In some cases, the arrays at the beginning of the sea ice time step are required. Those are referred to as $X\_b$.
+
+\section{Spatial domain}
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIGx : Vertical grid
+%
+%
+\begin{figure}[!ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=10cm,angle=-00]{../Figures/thermogrid.eps}
+\caption{\footnotesize{Vertical grid of the model, used to resolve vertical temperature and salinity profiles}}\label{fig_dom_icelayers}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+The horizontal indices $ji$ and $jj$ are handled as for the ocean in NEMO, assuming C-grid discretization and in most cases a finite difference expression for scale factors.
+
+The vertical index $jk=1, ..., nlay\_i$ is used for enthalpy (temperature) and salinity. In each ice category, the temperature and salinity profiles are vertically resolved over $nlay\_i$ equally-spaced layers. The number of snow layers can currently only be set to $nlay\_s=1$ (Fig. \ref{fig_dom_icelayers}).
+
+To increase numerical efficiency of the code, the two horizontal dimensions of an array $X(ji,jj,jk,jl)$ are collapsed into one (array $X\_1d(ji,jk,jl)$) for thermodynamic computations, and re-expanded afterwards.
+
+\forfile{../namelists/nampar}
+
+\section{Thickness space domain}
+
+\forfile{../namelists/namitd}
+
+Thickness space is discretized using $jl=1, ..., jpl$ thickness categories, with prescribed boundaries $hi\_max(jl-1),hi\_max(jl)$. Following \cite{Lipscomb01}, ice thickness can freely evolve between these boundaries. The number of ice categories $jpl$ can be adjusted from the namelist ($nampar$).
+
+There are two means to specify the position of the thickness boundaries of ice categories. The first option (ln\_cat\_hfn) is to use a fitting function that places the category boundaries between 0 and 3$\overline h$, with $\overline h$ the expected mean ice thickness over the domain (namelist parameter rn\_himean), and with a greater resolution for thin ice (Fig. \ref{fig_dom_icecats}). More specifically, the upper limits for ice in category $jl=1, ..., jpl-1$ are:
+\begin{eqnarray} 
+hi\_max(jl) = \biggr ( \frac{jl \cdot (3\overline h + 1 )^{\alpha}}{ (jpl-jl)(3 \overline h + 1)^{\alpha} + jl }\biggr )^{\alpha^{-1}} - 1,
+\end{eqnarray}
+with $hi\_max(0)$=0 m and $\alpha = 0.05$. The last category has no upper boundary, so that it can contain arbitrarily thick ice.
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIGx : Ice categories
+%
+%
+\begin{figure}[!ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=6cm,angle=-00]{../Figures/ice_cats.eps}
+\caption{\footnotesize{Boundaries of the model ice thickness categories (m) for varying number of categories and prescribed mean thickness ($\overline h$). The formerly used $tanh$ formulation is also depicted.}}\label{fig_dom_icecats}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+The other option (ln\_cat\_usr) is to specify category boundaries by hand using rn\_catbnd. The first category must always be thickner than rn\_himin (0.1 m by default).
+
+The choice of ice categories is important, because it constraints the ability of the model to resolve the ice thickness distribution. The latest study \citep{Massonnetetal18b} recommends to use at least 5 categories, which should include one thick ice with lower bounds at $\sim$4 m and $\sim$2 m for the Arctic and Antarctic, respectively, for allowing the storage of deformed ice.
+
+With a fixed number of cores, the cost of the model linearly increases with the number of ice categories. Using $jpl=1$ single ice category is also much cheaper than with 5 categories, but seriously deteriorates the ability of the model to grow and melt ice. Indeed, thin ice thicknes faster than thick ice, and shrinks more rapidly as well. When nn\_virtual\_itd=1 ($jpl$ = 1 only), two parameterizations are activated to compensate for these shortcomings. Heat conduction and areal decay of melting ice are adjusted to closely approach the 5 categories case.
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_dynamics.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_dynamics.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_dynamics.tex	(revision 10845)
@@ -0,0 +1,35 @@
+
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 3 Ñ Ice dynamics
+% ================================================================
+
+\chapter{Ice dynamics}
+\label{chap:DYN}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+In this chapter, we describe how the momentum equation is solved.
+
+I guess we could be brief. Gurvan and Clem should take the lead to write this.
+
+1) aEVP approach, revisited by Bouillon et al (2013).
+
+2) Short description of the numerical method. NEMO scale factors.
+
+3) Landfast ice option
+
+4) Boundary conditions (free-slip, no-slip, ...)
+
+4) Practical use of namelist parameters (nevp, telast, ...).
+
+\end{document}
+
+
+
+
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_interfaces.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_interfaces.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_interfaces.tex	(revision 10845)
@@ -0,0 +1,30 @@
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 8 Ñ Interfaces
+% ================================================================
+
+\chapter{Ice-atmosphere and ice-ocean interfaces}
+\label{chap:INT}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+In this chapter, we give information on the representation of interfaces.
+
+\section{Ice-ocean interface}
+
+\section{Ice-atmosphere interface}
+
+Drags?
+
+Blowing snow parameter
+
+Flux redistributor
+
+Jules coupling
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_miscellaneous.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_miscellaneous.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_miscellaneous.tex	(revision 10845)
@@ -0,0 +1,18 @@
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 12 -Miscellaneous
+% ================================================================
+
+\chapter{Miscellaneous topics}
+\label{chap:MIS}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+getting started, namelists, configs, ...
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_model_basics.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_model_basics.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_model_basics.tex	(revision 10845)
@@ -0,0 +1,457 @@
+
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+% ================================================================
+% Chapter 1 Ñ Model Basics
+% ================================================================
+
+\chapter{Model Basics}
+\label{chap:MB}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+%-------------------------------------------------------------------------------------------------------------------------
+%
+% 2.1 Assumptions
+\section{Rationale and assumptions}
+%
+%-------------------------------------------------------------------------------------------------------------------------
+\begin{itemize}
+\item Separation of horizontal dynamics and vertical thermodynamics;
+\item Dynamics: sea ice is a viscous-plastic continuum;
+\item Thermodynamics: sea ice is a mushy layer covered by snow
+\item Subgrid-scale physics
+\end{itemize}
+
+\subsection{Scales, thermodynamics and dynamics}
+Because sea ice is much wider -- $\mathcal{O}$(100-1000 km) -- than thick -- $\mathcal{O}$(1 m) -- ice drift can be considered as purely horizontal: vertical motions around the hydrostatic equilibrium position are negligible. The same scaling argument justifies the assumption that heat exchanges are purely vertical\footnote{The latter assumption is probably less valid, because the horizontal scales of temperature variations are $\mathcal{O}$(10-100 m)}. It is on this basis that thermodynamics and dynamics are separated and rely upon different frameworks and sets of hypotheses: thermodynamics use the ice thickness distribution \citep{Thorndikeetal75} and the mushy-layer \citep{Worster92} frameworks, whereas dynamics assume continuum mechanics \citep[e.g.,][]{Lepparanta05}. Thermodynamics and dynamics interact by two means: first, advection impacts state variables; second, the horizontal momentum equation depends, among other things, on the ice state.
+
+\subsection{Subgrid scale variations}
+
+Sea ice properties -- in particular ice thickness -- feature important changes at horizontal scales $\mathcal{O}$(1m) \citep{Thorndikeetal75}. An explicit representation of these variations is not and will not be -- at least in the next twenty years or so -- accessible to large-scale sea ice models. Yet important features, such as energy exchanges through the ice, quite non-linearly depend on ice thickness \citep{Maykut86}; whereas ice motion depends on the presence of open water, thin and thick ice at the very least, suggesting that subgrid-scale variations in ice properties must be accounted for, at least in a statistical fashion \citep{MaykutThorndike73}. 
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : Ice scheme
+%
+\begin{figure}[ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=10cm,angle=-00]{../Figures/ice_scheme.png}
+\caption{Representation of the ice pack, using multiple categories with specific ice concentration ($a_l, l=1, 2, ..., L$), thickness ($h^i_l$), snow depth ($h^s_l$), vertical temperature and salinity profiles ($T^i_{kl}$, $S^{*}_{kl}$) and a single ice velocity vector ($\bm{u}$).}
+\label{ice_scheme}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+The \textit{multi-category} framework \citep{MaykutThorndike73} addresses this issue by treating the ice thickness as an independent variable next to spatial coordinates and time, and introducing a thickness distribution\footnote{$g(h)$, termed the \textit{ice thickness distribution} is the density of probability of ice thickness \citep{Thorndikeetal75}.} $g(h)$ as the main prognostic model field. In the discrete world, the thickness distribution is converted into $L$ thickness categories. Ice thickness categories occupy a fraction of each grid cell, termed ice concentration ($a_l, l=1, 2, ..., L$), with specific thickness and properties. 
+
+The \textit{single-category} framework \citep{Hibler79} tackles the subgrid-scale issue by drastically simplifying the ice thickness distribution. The grid cell is divided into open water and sea ice characterized by a single ice concentration $A$ and mean thickness $H$. Single-category models (in particular LIM2) typically add parameterizations to represent the effects of unresolved ice thickness distribution on ice growth and melt \citep[see, e.g.][]{MellorKantha89,FichefetMaqueda97}.
+
+SI$^3$ provides the choice between either a multi- or a single-category framework. The default mode is multi-category. The single-category mode can be activated by setting the number of categories ($jpl=1$) and by activating the virtual thickness distribution parameterizations (nn\_monocat=1).
+
+\subsection{Thermodynamic formulation}
+
+%
+%--------------------------------------------------------------------------------------------------------------------
+%
+% TABLE x Physical constants
+%
+%
+\begin{table}[ht]
+\begin{footnotesize}
+\caption{Thermodynamic constants of the model.}
+\begin{center}
+\begin{tabular}{@{}llccl}\hline
+ & Description & Value & Units & Ref \\ \hline
+$c_i$ (cpic) & Pure ice specific heat & 2067 & J/kg/K & ? \\
+$c_w$ (rcp) & Seawater specific heat & 3991 & J/kg/K & \cite{TEOS-10} \\
+$L$ (lfus) & Latent heat of fusion (0$^\circ$C) & 334000 & J/kg/K & \cite{BitzLipscomb99} \\
+$\rho_i$ (rhoic) & Sea ice density & 917 & kg/m$^3$ & \cite{BitzLipscomb99} \\
+$\rho_s$ (rhosn) & Snow density & 330 & kg/m$^3$ & \cite{MaykutUntersteiner71} \\
+$\mu$ (tmut) & Linear liquidus coefficient & 0.054 & $^\circ$C/(g/kg) & \cite{Assur58} \\
+\hline
+\end{tabular}
+\end{center}
+\label{PhysicalConstants_table}
+\end{footnotesize}
+\end{table}
+%
+%--------------------------------------------------------------------------------------------------------------------
+%
+
+Ice thermodynamics are formulated assuming that sea ice is covered by snow. Within each thickness category, both snow and sea ice are horizontally uniform, hence each thickness category has a specific ice thickness ($h_i^l$) and snow depth ($h_s^l$). Snow is assumed to be fresh, with constant density and thermal conductivity. Sea ice is assumed to be a \textit{mushy layer}\footnote{Mushy layers are two-phase, two-component porous media.} \citep{Worster92} of constant density, made of pure ice and brine in thermal equilibrium, related by a linear liquidus relationship \citep{BitzLipscomb99}. A vertically-averaged bulk salinity $S_l$ uniquely characterizes brine fraction for each thickness category, and changes through  time from a simple parametrization of brine drainage. The linear vertical salinity profile ($S^*_{kl}$) is reconstructed from the vertical mean \citep{Vancoppenolleetal09b}. The diffusion of heat affects the vertical temperature profile, discretized on a unique layer of snow and multiple ice layers (typically 2-5) for each category, whereas thermal properties depend on local brine fraction. Growth and melt rates are computed, also for each ice category. The choice of the main thermodynamic constants  is described in Tab. \ref{PhysicalConstants_table}.
+
+
+
+%
+%--------------------------------------------------------------------------------------------------------------------
+%
+% TABLE x Global variables
+%
+%
+\begin{table}[ht]
+\begin{footnotesize}
+\caption{LIM global variables.}
+\begin{center}
+\begin{tabular}{@{}llcc}\hline
+Symbol   & Description & Units & Code name \\ \hline
+$\mathbf{u}$ & Sea ice velocity & [$m.s^{-1}$] & \textit{u\_ice, v\_ice (ji,jj)} \\
+$\mathbf{\sigma}$ & Stress tensor & [$Pa.m$] & \textit{stress1\_i, stress2\_i} \\
+ & & & \textit{stress12\_i (ji,jj)} \\
+$a_l$ & Concentration of sea ice in category $l$ & [-] & \textit{a\_i(ji,jj,jl)} \\
+$v^i_l$ & Volume of sea ice per unit area in category $l$ & [$m$] & \textit{v\_i(ji,jj,jl)} \\
+$v^s_l$ & Volume of snow per unit area in category $l$ & [$m$] & \textit{v\_s(ji,jj,jl)} \\
+$e^i_{kl}$ & Sea ice enthalpy per unit area in layer $k$ and category $l$ & [$J.m^{-2}$] & \textit{e\_i(ji,jj,jk,jl)} \\
+$e^s_{l}$ & Snow enthalpy per unit area in category $l$ & [$J.m^{-2}$] & \textit{e\_s(ji,jj,jl)} \\
+$M^s_l$ & Sea ice salt content in category $l$ & [$g/kg.m$] & \textit{smv\_i(ji,jj,jl)}\\
+%$O_l$ & Sea ice age content in category $l$ & $[days.m]$ \\
+\hline
+\end{tabular}
+\end{center}
+\label{GVariables_table}
+\end{footnotesize}
+\end{table}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+ Temperature, salinity, ice thickness, and snow depth are not extensive variables and therefore not conservative. Hence, conservative, extensive variables, must be introduced to ensure mass, salt and energy conservation. There are several back-and-forth conversions from extensive (conservative) state variables (see Table \ref{GVariables_table}) to intensive state variables of practical use (Table \ref{IVariables_table}).
+ 
+ %--------------------------------------------------------------------------------------------------------------------
+%
+% TABLE x Equivalent variables
+%
+%
+\begin{table}[ht]
+\begin{footnotesize}
+\caption{Intensive variables of practical use.}
+\begin{center}
+\begin{tabular}{@{}llcc}\hline
+Symbol   & Description & Units \\ \hline
+$h^i_l=v^i_l/g^i_l$ & Ice thickness & [$m$] \\
+$h^s_l=v^s_l/g^i_l$ & Snow depth & [$m$] \\
+$q^i_{m,kl}=e^i_{l,k}/(h^i_l/N)$ & Ice specific enthalpy & [$J.kg^{-1}$]\\
+$q^s_{m,l}=e^s_{l}/h^s_l$ & Snow specific enthalpy & [$J.kg^{-1}$] \\
+$T^i_{kl}=T(q^i_{kl})$ & Ice temperature & [$K$] \\
+$T^s_{l} =T(q^s_{l}) $ & Snow temperature & [$K$] \\
+$\overline S^i_{l}=M^s_l / v^i_l$ & Vertically-averaged bulk ice salinity & [g/kg] \\
+$S^*_{kl}$ & Depth-dependent ice salinity & [g/kg] \\
+$\phi_{kl}$ & Brine fraction & [-] \\
+%$o^i_{l}=O_l / g^i_l $ & Ice age & [$days$] \\
+\hline
+\end{tabular}
+\end{center}
+\label{IVariables_table}
+\end{footnotesize}
+\end{table}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+\subsection{Dynamic formulation}
+
+The formulation of ice dynamics is based on the continuum approach. The latter holds provided the drift ice particles are much larger than single ice floes, and much smaller than typical gradient scales. This compromise is rarely achieved in practice \citep{Lepparanta05}. Yet the continuum approach generates a convenient momentum equation for the horizontal ice velocity vector $\mathbf{u}=(u,v)$, which can be solved with classical numerical methods (here, finite differences on the NEMO C-grid). The most important term in the momentum equation is internal stress. We follow the viscous-plastic (VP) rheological framework \citep{Hibler79}, assuming that sea ice has no tensile strength but responds to compressive and shear deformations in a plastic way. In practice, the elastic-viscous-plastic (EVP) technique of  \citep{Bouillonetal13} is used, more convient numerically than VP.  It is well accepted that the VP rheology and its relatives are the minimum complexity to get reasonable ice drift patterns \citep{Kreyscheretal00}, but fail at generating the observed deformation patterns \citep{Girardetal09}. This is a long-lasting problem: what is the ideal rheological model for sea ice and how it should be applied are still being debated \citep[see, e.g.][]{Weiss13}. 
+
+%-------------------------------------------------------------------------------------------------------------------------
+%
+\section{Thickness distribution framework}
+%
+%-------------------------------------------------------------------------------------------------------------------------
+
+We first present the essentials of the thickness distribution framework \citep{Thorndikeetal75}. Consider a given region of area $R$ centered at spatial coordinates $(\bm{x})$ at a given time $t$. $R$ could be e.g. a model grid cell. The ice thickness distribution $g(\mathbf{x},t, h)$ is introduced as follows:
+\begin{linenomath}
+\begin{align}
+g(h) &= \lim_{dh \to 0} \frac{dA(h,h+dh)}{dh},
+\end{align}
+\end{linenomath}
+where $dA(h,h+dh)$ is the surface fraction of all parts of $R$ with ice thickness between $h$ and $h+dh$. Using this definition, the spatial structure of ice thickness is lost (see Fig. \ref{fig_g_h}), and $h$ becomes an extra independent variable, next to spatial coordinates and time, that can be thought as random. $g$ is by definition normalized to 1. The conservation of area, expressed in terms of $g(h)$, is given by \citep{Thorndikeetal75}: 
+\begin{linenomath}
+\begin{align}
+\frac{\partial g}{\partial t} = - \nabla .(g \mathbf{u}) - \frac{\partial}{\partial h} (f g) + \psi,
+\label{eq:gt}
+\end{align}
+\end{linenomath}
+where the terms on the right hand side refer to horizontal transport, thermodynamic transport in thickness space ($f$, m/s is the growth/melt rate), and mechanical redistribution rate, e.g. by ridging and rafting, where $\psi$ must conserve ice area and volume by construction.
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : Ice thickness distribution
+%
+\begin{figure}[ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=6cm,angle=-00]{../Figures/g_h.png}
+\caption{Representation of the relation between real thickness profiles and the ice thickness distribution function $g(h)$}
+\label{fig_g_h}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+In numerical implementations, the thickness distribution is discretized into several thickness categories, with specific ice concentration $a_l$ and ice volume per area $v_l^i$:
+\begin{linenomath}
+\begin{align}
+a_l &= \int _{H^*_{l-1}}^{H^*_l} dh \cdot g(h) , \\
+v_l^i &= \int _{H^*_{l-1}}^{H^*_l} dh \cdot h \cdot g(h).
+\end{align}
+\end{linenomath}
+Ice volume per area is the extensive counterpart for ice thickness, connected with volume through $h_l^i = v_l^i / a_l$. Evolution equations for extensive variables can be readily derived from equation \ref{eq:gt} by integration between thickness boundaries of the $l^{th}$ category \citep{Bitzetal01}. This applies to all model extensive variables (see Table \ref{GVariables_table}). For ice area, this reads:
+\begin{linenomath}
+\begin{align}
+\frac{\partial a_l}{\partial t} = - \bm{\nabla} \cdot (a_l \mathbf{u}) + \Theta^a_l + \int_{H^*_{l-1}}^{H^*_l} dh \psi.
+\label{eq:gt}
+\end{align}
+\end{linenomath}
+wher $\Theta^a_l$ refers to the effect of thermodynamics. Enthalpy is a particular case because it also has a vertical depth dependence $z$, which corresponds to $K$ vertical layers of equal thickness. The solution adopted here, following from \cite{ZhangRothrock01}, is that enthalpy from the individual layers are conserved separately. This is a practical solution, for lack of better.
+
+One of the major actions of LIM is to resolve conservation equations for all extensive variables that characterize the ice state.  Let us now connect this detailed information with classical sea ice fields. The ice concentration $A$ and the ice volume per area\footnote{Ice volume per area is equivalent to the grid-cell averaged ice thickness.} $V_i$ (m) directly derive from $g$: 
+\begin{linenomath}
+\begin{align}
+ A(\bm{x},t) &=\int_{0^+}^{\infty} dh \cdot g(h,\bm{x},t) \sim A_{ij} = \sum_{l=1}^L a_{ijl}, & \\
+ V_i(\bm{x},t)&=\int_{0}^{\infty} dh \cdot g(h,\bm{x},t) \cdot h \sim V^i_{ij} = \sum_{l=1}^L v^i_{ijl}. & \\
+\end{align}
+\end{linenomath}
+where the $0^+$ boundary implies that the means exclude open water. The mean ice thicknesses $H_i$ (m) is:
+\begin{linenomath}
+\begin{align}
+H_i = V_i/A,
+\end{align}
+\end{linenomath}
+whereas the open water fraction is simply $1-A$.
+
+\section{Governing equations}
+
+Let us now readily present the set of the LIM "governing" equations in the framework of the assumptions developed above. The conservation of horizontal momentum reads:
+\begin{linenomath}
+\begin{align}
+m \frac{\partial \bm{u}} {\partial t} & = \bm{\nabla}\cdot\bm{\sigma} +A \left(\bm{\tau}_{a}+\bm{\tau}_{w}\right) - m f \bm{k} \times \bm{u} - m g \bm{\nabla}{\eta},
+\label{a}
+\end{align}
+\end{linenomath}
+where $m=\rho_i V_i + \rho_s V_s $ is the ice and snow mass per unit area, $\bm{u}$ is the ice velocity, $\bm{\sigma}$ is the internal stress tensor, $\bm{\tau}_a$ and $\bm{\tau}_w$ are the air and ocean stresses, respectively, $f$ is the Coriolis parameter, $\bm{k}$ is a unit vector pointing upwards, $g$ is the gravity acceleration and $\eta$ is the ocean surface elevation. The EVP approach used in LIM \citep{Bouillonetal13} gives the stress tensor as a function of the strain rate tensor $\dot{\bm{\epsilon}}$ and some of the sea ice state variables:
+\begin{linenomath}
+\begin{align}
+\bm{\sigma} & = \bm{\sigma} (\dot{ \bm{\epsilon}}, \text{ice state}).
+\end{align}
+\end{linenomath}
+To the exception of velocity and internal stress, all extensive variables in Table \ref{GVariables_table} follow a conservation equation of the form:
+\begin{linenomath}
+\begin{align}
+\frac{\partial X}{\partial t} & = - \nabla .(\mathbf{u}X) + \Theta^X + \Psi^X,
+\label{eq:Xt}
+\end{align}
+\end{linenomath}
+including the effets of transport, thermodynamics ($\Theta^X$) and mechanical redistribution ($\Psi^X$). Solving these $jpl.(4+2.jpk)$ equations gives the temporal evolution of $\bm{u}$, $\bm{\sigma}$ and the rest of the global (extensive) variables listed in Table \ref{GVariables_table}.
+
+\section{Ice Dynamics}
+
+Dynamical processes include the conservation of momentum, rheology, transport and mechanical redistribution.
+To resolve the momentum equation, atmospheric stress is taken either as forcing or from an atmospheric model, oceanic stress and sea surface elevation from the ocean model, the Coriolis term is trivial. The last term, the divergence of the internal stress tensor $\sigma$, is the most critical term in the momentum equation and requires a rheological formulation. The EVP approach used in LIM gives the stress tensor components as \citep{Bouillonetal13}:
+\begin{linenomath}
+\begin{align}
+\sigma_{ij} & = \frac{P}{2(\Delta + \Delta_{min}) }
+\biggr [ ( \dot{\epsilon}_{kk} - \Delta) \delta_{ij} + \frac{1}{e^2} ( 2 \dot{\epsilon}_{ij} - \dot{\epsilon}_{kk} \delta_{ij} ) \biggr ],
+\end{align}
+\end{linenomath}
+where $\Delta$ is a particular measure of the deformation rate, $\Delta_{min}$ a parameter determining a smooth transition from pure viscous vlow ($\Delta<<\Delta_{min}$) to pure plastic flow ($\Delta >> \Delta_{min}$), and $e$ is a parameter giving the ratio between the maximum compressive stress and twice the maximum shear stress. In the pure plastic regime, the stress principal components should lie on the edge of an elliptical yield curve (Fig. \ref{fig_yield}). In the viscous regime, they are within the ellipse. The ice strength $P$ determines the plastic failure criterion and connects the momentum equation with the state of the sea ice. $P$ is not well constrained and must be parameterized. The heuristic option of \cite{Hibler79} was here adopted as a reference formulation:
+\begin{linenomath}
+\begin{align}
+P & = P^* V_i e^{-C(1-A)},
+\end{align}
+\end{linenomath}
+where $P^*$ and $C$ are empirical constants (see Table \ref{Modelparameters_table} for the values of the main model parameters).
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : Ice thickness distribution
+%
+\begin{figure}[ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=6cm,angle=-00]{../Figures/yield_curve.png}
+\caption{Elliptical yield curve used in the VP rheologies, drawn in the space of the principal components of the stress tensor ($\sigma_1$ and $\sigma_2$).}
+\label{fig_yield}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+%
+%--------------------------------------------------------------------------------------------------------------------
+%
+% TABLE x Model parameters %
+%
+\begin{table}[ht]
+\begin{footnotesize}
+\caption{Main model parameters.}
+\begin{center}
+\begin{tabular}{@{}llccl}\hline
+ & Description & Value & Units & Ref \\ \hline
+ 
+% Parameters
+$P*$ (rn\_pstar) & ice strength thickness param. & 20000 & N/m2 & - \\
+$C$ (rn\_crhg) & ice strength concentration param. & 20 & - & \citep{Hibler79} \\
+$H^*$ (rn\_hstar) & maximum ridged ice thickness param. & 25 & m & \citep{Lipscombetal07} \\
+$p$ (rn\_por\_rdg) & porosity of new ridges & 0.3 & - & \citep{Lepparantaetal95} \\
+$amax$ (rn\_amax) & maximum ice concentration & 0.999 & - & -\\
+$h_0$ (rn\_hnewice) & thickness of newly formed ice & 0.1 & m & - \\
+ 
+\hline
+\end{tabular}
+\end{center}
+\label{Modelparameters_table}
+\end{footnotesize}
+\end{table}
+%
+%--------------------------------------------------------------------------------------------------------------------
+%
+
+
+Transport connects the horizontal velocity fields and the rest of the ice properties. LIM assumes that the ice properties in the different thickness categories are transported at the same velocity. The scheme of \cite{Prather86}, based on the conservation of 0, 1$^{st}$ and 2$^{nd}$ order moments in $x-$ and $y-$directions,  is used, with some numerical diffusion if desired. Whereas this scheme is accurate, nearly conservative, it is also quite expensive since, for each advected field, five moments need to be advected, which proves CPU consuming, in particular when multiple categories are used. Other solutions are currently explored.
+
+The dissipation of energy associated with plastic failure under convergence and shear is accomplished by rafting (overriding of two ice plates) and ridging (breaking of an ice plate and subsequent piling of the broken ice blocks into pressure ridges). Thin ice preferentially rafts whereas thick ice preferentially ridges \citep{TuhkuriLensu02}. Because observations of these processes are limited, their representation in LIM is rather heuristic. The amount of ice that rafts/ridges depends on the strain rate tensor invariants (shear and divergence) as in \citep{FlatoHibler95}, while the ice categories involved are determined by a participation function favouring thin ice \citep{Lipscombetal07}. The thickness of ice being deformed ($h'$) determines whether ice rafts ($h'<$ 0.75 m) or ridges ($h'>$ 0.75 m), following \cite{Haapala00}. The deformed ice thickness is $2h'$ after rafting, and is distributed between $2h'$ and $2 \sqrt{H^*h'}$ after ridging, where $H^* = 25$ m \citep{Lipscombetal07}. Newly ridged ice is highly porous, effectively trapping seawater. To represent this, a prescribed volume fraction (30\%) of newly ridged ice \citep{Lepparantaetal95} incorporates mass, salt and heat are extracted from the ocean. Hence, in contrast with other models, the net thermodynamic ice production during convergence is not zero in LIM, since mass is added to sea ice during ridging. Consequently, simulated new ridges have high temperature and salinity as observed \citep{Hoyland02}. A fraction of snow (50 \%) falls into the ocean during deformation.
+
+\section{Ice thermodynamics}
+In this section, we develop the underlying principles of the thermodynamic formulation, summarized in the term $\Theta^X$, where $X$ refers to all extensive state variables. $\Theta^X$ includes the contributions of transport in thickness space and thermodynamic source and sink terms.
+
+\subsection{Transport in thickness space}
+Transport in thickness space describes how vertical growth and melt moves ice state variables among the different thicknesses at a velocity $f(h)$, the net ice growth/melt rate, which needs to be first computed.  In discretized form, this term moves ice properties between neighbouring categories. The linear remapping scheme of \cite{Lipscomb01} is used. This scheme is semi-lagrangian, second-order, is less diffusive and converges faster than other options.
+
+\subsection{Thermodynamic source and sink terms}
+
+Since heat, salt and mass are strongly inter-dependent for sea ice, the thermodynamic source and sink terms are treated together. They include the changes in extensive sea ice state associated with thermodynamic processes. The latter are separated in two main parts: (i) open water fraction processes, where atmosphere and ocean are in direct interaction; and (ii) vertical ice thermodynamic processes, driven by surface snow/ice-atmosphere and basal ice-ocean exchanges, for each thickness category. For each part, first, the energy available or lost is specified. Then the impact on mass exchanges is evaluated. The latter part requires to specify how sea ice and snow responds to energy supply or loss, which is achieved through the enthalpy formulation.
+
+%\subsubsection{Atmosphere-ice-ocean heat and mass exchanges}
+%
+%The atmosphere-ice heat exchanges are either specified in external files (forcing) or computed from an atmospheric model. The following heat and mass fluxes should be specified:
+%\begin{itemize}
+%\item Net solar flux, for open water ($Q^{sol}$) and ice categories ($Q^{sol,ice}_{l}$);
+%\item Net non-solar flux, for open water ($Q^{nsol}$) and ice categories ($Q^{nsol,ice}_{l}$);
+%\item Solid and liquid precipitation (kg/m$^2$/s).
+%\end{itemize}
+%
+%In terms of mass (freshwater), the ocean gets:
+%\begin{itemize}
+%\item All liquid precipitation;
+%\item The part of solid precipitation that is not intercepted by the ice, accounting for the part that is blown off the sea ice by winds;
+%\item The net melting minus freezing mass flux, as diagnosed from thermodynamic sea ice calculations.
+%\end{itemize}
+%
+%In terms of energy, the ocean gets:
+%\begin{itemize}
+%\item The net atmospheric solar flux, minus what is intercepted by the ice; plus what is transmitted below the ice; 
+%\item The net atmospheric non-solar heat flux, minus what is used by the sea ice; 
+%\item The sensible heat exchange associated with ice-ocean and atmosphere-ocean mass fluxes. 
+%\end{itemize}
+%
+%In terms of salt, the ocean gets the net salt release minus uptake associated with sea ice growth, melt and brine drainage. There is no atmosphere-ocean salt flux.
+
+\subsubsection{Enthalpy formulation}
+
+% Add snow and water
+
+A first overarching aspect of the thermodynamic calculations is the specification of the response of sea ice to energy supply. This is achieved by defining the internal energy (or enthalpy\footnote{Wording it internal energy or enthalpy is equivalent since pressure effects are not considered.}). This ultimately relies on the response of the phase composition to salinity and temperature changes. The enthalpy formulation used in LIM is based on the following assumptions:
+\begin{itemize}
+\item Sea ice is gas-free, composed solely of pure ice and saline brine, characterized by brine fraction $\phi$;
+\item brine and pure ice are in thermodynamic equilibrium;
+\item the salinity-dependence of the freezing point is linear (linear liquidus);
+\item the density of the sea ice (ice+brine) medium is constant ($\rho_i$).
+\end{itemize}
+%In this context, the bulk salinity of a small control sea ice volume is related to the salinity of brine ($S_{br}$):
+%\begin{linenomath}
+%\begin{align}
+%S = \phi S_{br}
+%\label{eq_S_Sbr}
+%\end{align}
+%\end{linenomath}
+%Because of thermal equilibrium and linear liquidus, brine salinity directly depends on local temperature: $T=-\mu S_{br}$. Hence, brine fraction is given by
+%\begin{linenomath}
+%\begin{align}
+%\phi = -\frac{\mu S}{T}.
+%\end{align}
+%\end{linenomath}
+%This equation relates brine fraction to bulk salinity and temperature, which therefore can be seen as the two main state variables in this linear mushy layer formulation for sea ice thermodynamics. 
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : Thermal properties
+%
+\begin{figure}[ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=8cm,angle=-00]{../Figures/Thermal_properties.png}
+\caption{Thermal properties of sea ice vs temperature for different bulk salinities: brine fraction, specific enthalpy, thermal conductivity, and effective specific heat.}
+\label{fig_thermal_properties}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+Based on these, brine fraction reduces to $\phi = -\mu S/T$ (see Fig. \ref{fig_thermal_properties}), where $\mu$ relates the freezing point of brine to salinity, and one can derive the specific enthalpy $q_m(S,T)$, defined as the energy required to warm and melt a unit control volume of sea ice at temperature $T$ (in Celsius) and salinity $S$ until 0$^\circ$ C, taken as a reference zero-energy level \citep{BitzLipscomb99,Schmidtetal04}:
+\begin{linenomath}
+\begin{align}
+q_m(S,T) = \biggr [ c_i ( T + \mu S) - L \biggr ( 1 + \frac{\mu S}{T} \biggr ) - c_w \mu S \biggr ]
+\end{align}
+\end{linenomath}
+where $c_i$ is pure ice specific heat, $L$ is latent heat of fusion at 0$^\circ$C, and $c_w$ is water specific heat. The first term expresses the warming of solid ice. The second term expresses internal change in brine fraction, which is often the largest because the Stefan number  ($c_i T/L$) is generally small. The last term gives the warming of the remaining water from $T_{fr} = -\mu S$ until 0$^\circ$C. Similar, but simpler and linear expressions for snow and water can be derived.
+
+The second overarching aspect is that all growth and melt processes must be calculated consistently with the enthalpy formulation. Energetics of phase transitions are handled using the formalism of \cite{Schmidtetal04}. For each phase transition, initial and final states (temperature and salinity) are defined, and the ice-to-ocean mass flux to the ice $F_m$ (kg/s) relates to the energy gain or loss $\Delta Q$ through:
+\begin{linenomath}
+\begin{align}
+\Delta Q / \Delta q_m = F_m,
+\label{eq_phasechange}
+\end{align}
+\end{linenomath}
+where $\Delta q_m$ is the change in specific enthalpy involved in the considered phase transition, from initial to final state.
+
+
+\subsubsection{Open water processes}
+
+As part of the sea ice thermodynamic calculations, a heat budget estimate for the uppermost ocean level ($B^{opw}$) must be included, to compute the rate of new ice formation or the contribution of sensible heat to bottom melting. $B^{opw}$ includes:
+\begin{itemize}
+\item the absorption of a fraction $f_1^{qsr}$ of solar radiation (given by radiative transfer component of the ocean model);
+\item the non-solar heat flux absorbed at the surface;
+\item the sensible heat content of precipitation
+\item the sensible heat flux from the ocean to the sea ice ($A.F_w$)
+\end{itemize}
+Other contributions are not assumed not to contribute. The ocean-to-ice sensible heat flux is formulated the bulk formula of \citep{McPhee92}.
+
+% Ice growth
+If $B^{opw}$ is such that the SST would decrease below the freezing point, the remainder of the heat is used to form new ice. The heat loss is converted into a volume of new ice $v_0$. The thickness $h_0$ of the new ice grown during a sea ice time step depends on unresolved small-scale currents and waves and is prescribed. The fraction $a_0=v_0/h_0$ is computed accordingly. The salinity of this new ice $S_0$ is given by the salinity-thickness empirical relationship of \cite{Kovacs96}. The temperature assumed for this new ice is the local freezing point. If by contrast $B^{opw}$ is positive and there still is ice in the grid cell, then $B^{opw}$ is directly redirected to bottom melting. This argument follows from \cite{MaykutMcPhee95}, who found that most of solar heat absorbed in the surface waters is converted into melting. In practise, this prevents the SST to be above freezing as long ice is present.
+
+$B^{opw}$ can be seen as a predictor of the heat budget of the first ocean level. As such, it only helps to compute new ice formation and the extra bottom melt in summer, but is not part of the conservation of heat in the model. To ensure heat conservation, the heat effectively contributing to changing sea ice is removed from the non-solar flux sent to the ocean. This includes: (i) the heat loss used for ice formation, (ii) the heat gain used to melt ice, and (iii) the sensible heat given by the ocean to the ice. Finally, because ice dynamics are not able to maintain the small amount of open water that is observed, a maximum ice fraction ($amax, <1$) is prescribed.
+
+\subsubsection{Vertical ice thermodynamic processes}
+
+The second part of the computations regard the computation of purerly vertical processes in the ice-covered part of the grid cell, similarly for each ice category. 
+
+\textbf{Surface melt, basal growth and melt and diffusion of heat}. The surface melt rate, as well as the basal growth / melt rate depend on the energy budget at the upper and lower interfaces, respectively, between the external fluxes either from the atmosphere or the ocean, and the internal conduction fluxes. The internal conduction fluxes depend on the internal temperature profile, which is determined by solving the enthalpy equation:
+\begin{linenomath}
+\begin{align}
+\rho \frac{\partial q_m}{\partial t} = - \frac{\partial}{\partial z} (F_c + F_r ).
+\end{align}
+\end{linenomath}
+which state that the local change in enthalpy is given by the divergence of the vertical conduction ($F_c=-k(S,T) \partial T / \partial z$) and radiation ($F_r$) fluxes. $\rho$ is the density of ice or snow. Re-expressed as a function of temperature, this becomes the heat diffusion equation.  This equation is non-linear in $T$, because of $q$ and $k$, and its main specificity is that internal melting requires large amounts of energy near the freezing point. The thermal conductivity is formulated following \cite{Pringleetal07}, empirically accounting for the reduction of thermal conductivity at large brine fractions. 
+
+At the ice base, we assume that the temperature is at the local freezing point.  Ice grows or melt if the heat balance between the oceanic sensible heat flux ($F_w$) and internal conduction is negative or positive.
+
+At the ice surface, the boundary condition on the heat diffusion equation is:
+\begin{linenomath}
+\begin{align}
+Q^{sr}  + Q^{ns}(T_{su})= F_c + Q^{sum}.
+\end{align}
+\end{linenomath}
+where $Q^{sr}$ and $Q^{ns}$ are the net downwelling atmospheric solar and non-solar flux components. If the solution of this equation without melting gives a surface temperature ($T_{su}$) below 0$^\circ$ C, then there is no melting and the heat available for surface melting $Q_{sum}=0$.  Otherwise $T_{su}$ is capped at 0$^\circ$ C and $Q_{sum}$ is calculated as a residual.
+
+\textbf{Radiation}. Radiation contributes to the surface and internal heat budget. The radiative transfer scheme is currently basic, composed of surface albedo, transmission through the ice interior and attenuation with vertical depth. The albedo is computed empirically as a function of ice thickness, snow depth and surface temperature, using a reformulation of the parameterization of \cite{ShineHenderson85}. When snow is present, all the absorbed radiation is transformed into sensible heat available for conduction or melting. Over snow-free ice, a fraction of solar radiation is transmitted below the surface and attenuates exponentially with depth, until it reaches the base of the ice.
+
+\textbf{Growth and melt processes}.  Snow grows from precipitation and loses mass from melting and snow-ice conversion once the snow base is below sea level. Sea ice grows and melts by various means. Ice forms by congelation or melt at the base, can melt at the surface and form from snow-to-ice conversion at the snow-ice interface if the latter is below sea level. Some new ice is also added to the system when seawater is trapped into newly formed pressure ridges.
+
+\textbf{Salt dynamics}. Bulk salinity is empirically parameterized, as a function of salt uptake during growth, gravity drainage and flushing. The shape of the vertical profile depends on the bulk salinity \citep{Vancoppenolleetal09b}.
+
+\textbf{Single-category parameterizations}. If the single-category representation is adopted, then two parameterizations can be activated, following \citep{FichefetMaqueda97}. First, the thermal conductivity of both ice and snow is multiplied by a factor $>1$ accounting for the unresolved thin ice, effectively increasing the ice growth rate. Second, to account for the loss of thin ice in summer, the ice concentration is reduced in proportion to the loss of ice thickness. Both parameterizations have been tuned to match the results in multi-category mode.
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_output_diagnostics.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_output_diagnostics.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_output_diagnostics.tex	(revision 10845)
@@ -0,0 +1,65 @@
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 9 - Outputs and diagnostics
+% ================================================================
+
+\chapter{Output and diagnostics}
+\label{chap:DIA}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+\section{SIMIP diagnostics}
+
+The SIMIP protocol \citep{Notzetal16} was designed for CMIP6, to standardize sea ice model outputs in climate simulations. We tried to follow the data request as closely as possible. Outputs are in most cases directly managed with XIOS2 in \textbf{limwri.F90}, but not always. In the code, output fields keep their native LIM reference name. 
+
+A corresponding entry exists in \textbf{field\_def\_nemo-lim.xml}, where fields are given their SIMIP specifications (standard name, long name, units). At the end of the file the fields are gathered in the field groups \textbf{SIday\_fields}, \textbf{SImon\_fields} and \textbf{SImon\_scalar} for separation of the daily (SIday) and monthly (SImon) requests. 
+
+In \textbf{file\_def\_nemo-lim.xml}, the daily, monthly and scalar output files are created.
+
+In the reference xml files, the largest possible SIMIP-based diagnostics with LIM are distributed among the field groups. If some fields are to be discarded, the best way to do so is to remove them from the field groups in  field\_def\_nemo-lim.xml.
+
+\subsection{Missing SIMIP fields}
+
+About 90\% of the SIMIP fields can be output. Below is the list of the missing fields and why they are missing. \\
+
+\textbf{1. Fields that are not part of the sea ice representation in LIM3.6}
+\begin{itemize}
+\item sisnconc (snow area fraction), siitdsnconc (snow area fractions in thickness categories);
+\item simpconc (meltpond area fraction), simpmass (melt pond mass per area), simprefrozen (thickness of refrozen ice on ponds);
+\item sirdgconc (ridged ice area fraction), sirdgmass (ridged ice thickness);
+\item sidmasslat (lateral sea ice melt rate);
+\item sndmasswindrif (snow mass change through wind drift of snow);
+\end{itemize}
+
+\textbf{2. Fields which value is trivial}
+\begin{itemize}
+\item sipr (rainfall over sea ice): all rain falls in open water;
+\item sidragtop (atmospheric drag over sea ice): namelist parameter;
+\item sidragbot (oceanic drag over sea ice): namelist parameter
+\end{itemize}
+
+\textbf{3. Fields that belong to the atmospheric component}
+\begin{itemize}
+\item siflswdtop, siflswutop, siflswdbot, sifllwdtop, sifllwutop, siflsenstop, sifllatstop (surface energy budget components)
+\end{itemize}
+
+Ice thickness and snow depth were masked below 5\% ice concentration, because below this value, they become meaninglessly large in LIM. This is notably because of the Prather advection scheme. We hope to fix these issues for our next release. For similar reasons, the ice age is masked below 15\% concentration.
+
+Fluxes through straits and passages were not directly implemented. Instead, ice mass, snow mass, and ice area transports were implemented as 2D arrays, for x- and y- directions. A python script is available to derive the fluxes through straits and passages from full 2D arrays for ORCA2 and eORCA1 grids.
+
+\subsection{Links}
+
+\begin{itemize}
+\item \textcolor{blue}{\href{http://www.geosci-model-dev.net/9/3427/2016/}{Paper of Notz et al}};
+\item \textcolor{blue}{\href{http://clipc-services.ceda.ac.uk/dreq/u/SIMIP.html}{SIMIP CMIP6 data request page}};
+\item \textcolor{blue}{\href{http://www.climate-cryosphere.org/activities/targeted/simip}{SIMIP description on CliC website}}.
+\end{itemize}
+
+\section{Conservation checks}
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_radiative_transfer.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_radiative_transfer.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_radiative_transfer.tex	(revision 10845)
@@ -0,0 +1,210 @@
+
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+
+
+
+% ================================================================
+% Chapter 6 Ñ Radiative transfer
+% ================================================================
+
+\chapter{Radiative transfer}
+\label{chap:RAD}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+Radiative transfer in SI$^3$ currently reduces to the parameterization of solar radiation partitionning through the snow/ice/open water system, treated using a single wavelength band. This will likely be improved in future versions of the code. In this chapter, we first explain how solar radiation is partionned in the snow-ice system, then describe how, solar radiation-wise, the snow-ice system is framed in the context of the atmosphere-ice-ocean boundary.
+
+\section{Solar radiation partitionning in the snow-ice system}
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : Radiation balance %
+\begin{figure}[!ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=6cm,angle=-00]{../Figures/radiative_transfer.png}
+\caption{Partitionning of solar radiation in the snow-ice system, as represented in SI$^3$.}
+\label{fig_radiative_transfer}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+Solar radiation in the snow-ice system is represented following the principles of \cite{MaykutUntersteiner71}, see Fig.\ref{fig_radiative_transfer}, using a unique band of solar radiation. Incident solar radiation (W/m$^2$, counted per unit ice area - not per grid cell area) is specified in the SBC routines and is a priori category dependent, because multiple atmosphere-surface reflexions are frequent in polar regions imply that incident radiation depends on the surface albedo and therefore surface state. 
+
+Net solar radiation qsr\_ice(i,j,l) is obtained by substracting the reflected part of the incident radiation using the surface albedo $\alpha(i,j,l)$, parameterized as a function of environmental conditions. 
+
+The subsequent attenuation of solar radiation through the snow-ice system is represented assuming the presence of a highly diffusive surface scattering layer, absorbing a fraction $i_o$ of net solar radiation, which is transformed into sensible heat, contributing to the surface energy balance. 
+
+The remainder of solar radiation, qtr\_ice\_top(i,j,l), is transmitted below the surface and attenuates following Beer-Lambert law. The part of solar radiation that is absorbed on its path to the base of the ice is given as sensible heat to the snow/ice system, via a source term in the heat diffusion equation. The rest of solar radiation that reaches the ice base, qtr\_ice\_bot(i,j,l), is transmitted to the ocean. 
+
+In the rest of this section, we describe how the albedo, the surface transmission parameter $i_o$ and the attenuation of solar radiation are parameterized.
+
+\subsection{Surface albedo}
+
+The surface albedo determines the amount of solar radiation that is reflected by the ice surface, hence also net solar radiation. The philosophy of the parameterization of surface albedo is the following: each ice category has its own albedo value $\alpha(i,j,l)$, determined as a function of cloud fraction, ice thickness, snow depth, melt pond fraction and depth, using observation-based empirical fits. 
+
+The original \cite{ShineHenderson85} parameterization had a few inconsistencies and flaws that the revisited parameterization described hereafter fixes. In particular, the dependencies of albedos on ice thickness, snow depth and cloud fraction have been revised in the light of recent observational constraints \citep{Brandtetal05,GrenfellPerovich04}. In addition, the asymptotic properties of albedo are better specified and now fully consistent with oceanic values. Finally, the effect of melt ponds has been included \citep{Lecomteetal15}.
+
+The user has control on 5 reference namelist values, which describe the asymptotic values of albedo of snow and ice for dry and wet conditions, as well as the deep ponded-ice albedo. Observational surveys, in particular during SHEBA in the Arctic \citep{Perovichetal02alb} and further additional experiments \citep{GrenfellPerovich04}, as well as by \cite{Brandtetal05} in the Antarctic, have provided relatively strong constraints on the surface albedo. In this context, the albedo can hardly be used as the main model tuning parameter, at least outside of these observation-based bounds (see namalb for reference values).
+
+\forfile{../namelists/namalb}
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : cloud correction
+%
+\begin{figure}[ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=10cm,angle=-00]{../Figures/albedo_cloud_correction.png}
+\caption{Albedo correction $\Delta \alpha$ as a function of overcast sky (diffuse light) albedo $\alpha_os$, from field observations \cite[][their Table 3]{GrenfellPerovich04} (squares) and 2nd-order fit (Eq. \ref{eq_albedo_cloud_correction}). Red squares represent the irrelevant data points excluded from the fit. For indication, the amplitude of the correction used in the ocean component is also depicted (blue circle).}
+% ocean uses 0.06 for overcast sky (Payne 74) and Briegleb and Ramanathan parameterization
+\label{fig_albedo_cloud_correction}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+Because the albedo is not an intrinsic optical property, it depends on the type of light (diffuse of direct), which is practically handled by weighting the clear (cs) and overcast (os) skies values by cloud fraction $c(i,j)$ \citep{FichefetMaqueda97}:
+\begin{equation}
+\alpha(i,j,l) = [ 1 - c(i,j) ] \cdot \alpha_{cs} (i,j,l) + c (i,j) \cdot \alpha_{os}(i,j,l).
+\end{equation}
+For concision, we drop the spatial and category indices hereafter. \cite{GrenfellPerovich04} observations at Point Barrow, on the Alaskan Coast, suggest that clear and overcast sky albedos are directly related through
+\begin{equation}
+\alpha_{cs} = \alpha_{os} - \Delta \alpha(\alpha_{os}).
+\end{equation}
+The relation between $\Delta \alpha$ and $ \alpha_{os}$ can well be handled using a 2$^{nd}$-order polynomial fit (Fig. \ref{fig_albedo_cloud_correction}):
+\begin{equation}
+\Delta \alpha =  ( -0.1010 \cdot \alpha_{os}^2 + 0.1933 \cdot \alpha_{os} - 0.0148 ).
+\label{eq_albedo_cloud_correction}
+\end{equation}
+Overcast sky surface albedo is used as a reference, from which the clear-sky value is derived.
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : thickness, snow depth, pond depth dependencies
+%
+\begin{figure}[ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=4cm,angle=-00]{../Figures/albedo_dependencies.png}
+\caption{Example albedo dependencies on ice thickness, snow depth and pond depth, as parameterized in SI$^3$.}
+\label{fig_albedo_dependencies}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+The second important parameter that controls surface albedo is surface type. In each category, we assume that three types of surfaces can coexist (bare, snow-covered and ponded ice), with respective fractions $f_{ice}$, $f_{snw}$ and $f_{pnd}$ summing to 1. Then the overcast albedo is expressed as
+\begin{equation}
+\alpha_{os}(i,j,l) = f_{ice} \cdot \alpha_{ice} + f_{snw} \cdot \alpha_{snw} + f_{pnd}\cdot \alpha_{pnd}
+\end{equation}
+with a specific albedo value for each surface type. 
+
+The surface fractions $f_{ice}$, $f_{snw}$ and $f_{pnd}$ are currently crudely parameterized: if snow is present ($h_s>0$), then $f_{snw}=1$ and $f_{ice}=f_{pnd}=0$. In the absence of snow, $f_{pnd}$ is either specified or calculated (depending on melt pond options in nampnd), and $f_{ice}=1.-f_{pnd}$. Admittedly, more refined parameterizations of $f_{snw}$ could improve the realism of the model. Note finally that the dependence of surface albedo on the presence of melt ponds can be included or not (namelist parameter ln\_pnd\_alb). If the latter is set to false, $f_{pnd}$ is always assumed zero in the albedo computations.
+
+Works by \cite{Brandtetal05} and references therein, indicate that the dependence of the albedo of bare ice on ice thickness depends is linear/logarithmic/constant from thin to thick ice. Hence, the following expressions capture  the essence of their works:
+\begin{eqnarray}
+\alpha_{ice} = 
+\begin{cases}
+\alpha_{ice}^{\infty} & \text{ if } h_i > 1.5, \\
+\alpha_{ice}^{\infty} + ( 0.18 - \alpha_{ice}^{\infty} ) \cdot \frac{ln(1.5) - ln(h_i)}{ln(1.5) - ln(0.05)} & \text{ if } 0.05 < h_i ,<= 1.5 \\
+\alpha_{oce}  + ( 0.18 - \alpha_{oce} ) h_i /0.05 & \text{ if } h_i < 0.05.
+\end{cases}
+\end{eqnarray}
+The thick-ice constant albedo value depends on whether the surface is dry or melting:
+\begin{eqnarray}
+\alpha_{ice}^{\infty} =
+\begin{cases}
+\alpha_{i,dry} & \text{ if } T_{su} < T_{fr} \\
+\alpha_{i,mlt} & \text{ if } T_{su} = T_{fr},
+\end{cases}
+\end{eqnarray}
+values that are to be specified from the namelist.
+
+\cite{GrenfellPerovich04} suggest that the dependence of surface albedo on snow depth is exponential,
+\begin{eqnarray}
+\alpha_{snw} = \alpha_{snw}^{\infty} - ( \alpha_{snw}^{\infty} - \alpha_{ice} ) * exp( -h_s / h_s^{ref} ),
+\end{eqnarray}
+where $h_s^{ref} = 0.02$ $(0.03)$ m for dry (wet) snow. As for bare ice, the deep-snow asymptotic albedo also depends on whether the surface is dry or melting:
+\begin{eqnarray}
+\alpha_{snw}^{\infty} =
+\begin{cases}
+\alpha_{s,dry} & \text{ if } T_{su} < T_{fr} \\
+\alpha_{s,mlt} & \text{ if } T_{su} = T_{fr},
+\end{cases}
+\end{eqnarray}
+values that are to be specified from the namelist.
+
+Based on ideas developed from melt ponds on continental ice \citep{ZuoOerlemans96}, the albedo of ponded ice was proposed to follow \citep{Lecomteetal11}:
+\begin{eqnarray}
+\alpha_{pnd} = \alpha_{dpnd} - ( \alpha_{dpnd} -\alpha_{ice} ) \cdot exp( -h_{pnd} / 0.05 )
+\end{eqnarray}
+$\alpha_{dpnd}$ is a namelist parameter. \cite{EbertCurry93} also use such dependency for their multi-spectral albedo.
+
+The dependencies of surface albedo on ice thickness, snow depth and pond depth are illustrated in Fig. \ref{fig_albedo_dependencies}.
+
+\subsection{Transmission below the snow/ice surface}
+
+The transmitted solar radiation below the surface is represented following \cite{FichefetMaqueda97} and \cite{MaykutUntersteiner71}:
+\begin{eqnarray}
+qtr\_ice\_top(i,j,l) = i_o(i,j) qsr\_ice(i,j,l),
+\end{eqnarray}
+where $i_o=0$ in presence of snow, and depends on cloud fraction otherwise, based on works of \cite{GrenfellMaykut77}. This parameterization needs to be re-evaluated and likely updated.
+
+\subsection{Attenuation and transmission below the ice/ocean interface}
+
+Attenuation of solar radiation through the ice follows Beer-Lambert law. In practise, we assume that irradiance below layer $k$ is given by
+
+\begin{eqnarray}
+radtr\_i(i,j,k,l) = qtr\_ice\_top(i,j,l) \cdot exp(-\kappa_i z),
+\end{eqnarray}
+where $\kappa_i = 1$ m$^{-1}$ is the exponential attenuation coefficient (namelist parameter rn\_kappa\_i). Hence, at the ice base, remains below the $l^{th}$ category a transmitted flux: 
+\begin{eqnarray}
+qtr\_ice\_bot(i,j,l) = qtr\_ice\_top(i,j,l) \cdot exp(-\kappa_i h_i).
+\end{eqnarray}
+
+\section{Solar radiation: framing sea ice at the ocean-atmosphere boundary}
+
+How solar radiation transfer through sea ice is framed into the atmosphere-ice-ocean is nearly identical but not exactly the same in forced and coupled mode (see Fig. \ref{fig_radiative_transfer}.
+
+The basic principle of the computation is that the irradiant flux given to the ocean model (qsr) is computed as the average flux per grid cell area (qsr\_tot) minus what is given to the sea ice ($\sum a(l) qsr\_ice(l)$), plus what is transmitted below sea ice $\sum$ qtr\_ice\_bot(jl) (see at the base of Fig. \ref{fig_radiative_transfer}). Such formulation ensures heat conservation by construction.
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : Radiation balance %
+\begin{figure}[!ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=8cm,angle=-00]{../Figures/radiation_atm_ice_oce.png}
+\caption{Framing solar radiation transfer through sea ice into the atmosphere-ice-ocean context.}
+% ocean uses 0.06 for overcast sky (Payne 74) and Briegleb and Ramanathan parameterization
+\label{fig_radiative_transfer}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+\subsection{Forced mode}
+
+In forced-atmosphere mode, it is the incoming solar irradiance fluxes above the ocean and sea ice (categories) that are specified (from files) or computed (from bulk formulae), and constitute the basis of solar radiation transfer computations. Then the net solar fluxes above open water (qsr\_oce) and ice categories (qsr\_ice) are obtained by multiplication by $1-\alpha$. qsr\_tot is then diagnosed as a weighted sum of qsr\_oce and the qsr\_ice(jl)'s.
+
+\subsection{Coupled mode}
+
+In coupled-atmosphere mode, qsr\_tot and qsr\_ice have to be provided by the atmospheric model, whereas qsr\_oce is diagnosed from qsr\_ice and qsr\_tot.
+
+Some atmospheric models enable \textit{tiling} and can provide solar fluxes over individual ice categories. For such atmospheric models, net solar radiation fluxes are directly useable by SI$^3$ (nn\_flxdist =  -1). Other models cannot do tiling, being only able to provide a net solar flux above all ice categories, seen as a single surface type. For such models a first option is to give the net solar flux above sea ice identically to all sea ice categories (nn\_flxdist = 0). Yet a better option is to redistribute the mean solar flux above sea ice $< qsr_ice >$ above categories (nn\_flxdist = 2) using the following scaling, conserving heat by construction:
+
+\begin{eqnarray}
+qsr\_ice(jl) = < qsr\_ice > \frac{1 - \alpha(jl)}{1 - < \alpha >}
+\end{eqnarray}
+
+where $ < \alpha >$ is the albedo averaged over the ice categories. Note that for testing, the flux redistributor can be emulated in forced mode (nn\_flxdist = 1).
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_ridging_rafting.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_ridging_rafting.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_ridging_rafting.tex	(revision 10845)
@@ -0,0 +1,156 @@
+
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 5 Ñ Ridging and rafting
+% ================================================================
+
+\chapter{Ridging and rafting}
+\label{chap:RDG}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+This chapter focuses on how LIM solves the red part of the general equation:
+\begin{equation}
+\frac{\partial X}{\partial t} = - \nabla .(\mathbf{u} X)  \color{black} \Theta^X  + \color{red} \Psi^X,
+\label{eq:Xt}
+\end{equation}
+where $X$ refers to any global sea ice state variable.
+
+Divergence and shear open the ice pack and create ice of zero thickness. Convergence and shear consumes thin ice and create thicker ice by mechanical deformation. The redistribution functions $\Psi^X$ describe how opening and mechanical deformation redistribute the global ice state variables into the various ice thickness categories. 
+
+The fundamental redistribution function is $\Psi^g$, which accounts for area redistribution. The other redistribution functions $\Psi^X$ associated with other state variables will derive naturally. The redistribution function $\Psi^g$ should first ensure area conservation. By integrating the evolution equation for $g(h)$ over all thicknesses, recalling that $\int_0^\infty g(h) = 1$, and that the total areal change due to thermodynamics must be zero, e.g. $\int_0^\infty \partial (fg) / \partial h=0$, then the area conservation reads: 
+\begin{equation}
+\int_0^\infty h \Psi^g dh = \nabla \cdot \mathbf{u}.
+\end{equation}
+
+Second, we must say something about volume conservation, and this will be done more specifically later. Following \cite{Thorndikeetal75}, we separate the $\Psi^X$'s into \textbf{(i)} \textit{dynamical inputs}, \textbf{(ii)} \textit{participation functions}, i.e., how much area of ice with a given thickness participates to mechanical deformation \textbf{(iii)} \textit{transfer functions}, i.e., where in thickness space the ice is transferred after deformation.
+
+\section{Dynamical inputs}
+
+A general expression of $\Psi^g$, the mechanical redistribution function associated to the ice concentration, was proposed by \cite{Thorndikeetal75}:
+\begin{equation}
+\Psi^g = \vert \dot{\epsilon} \vert [ \alpha_o(\theta) \delta(h) + \alpha_d(\theta) w_d(h,g) ],
+\label{eq:psi}
+\end{equation}
+which is convenient to separate the dependence in $\mathbf{u}$ from those in $g$ and $h$. The first and second terms on the right-hand side correspond to opening and deformation, respectively. $\vert \dot{\epsilon} \vert = (\dot{\epsilon}_{I}^2+ \dot{\epsilon}_{II}^2)^{1/2}$, where $\dot{\epsilon}_{I}=\nabla \cdot \mathbf{u}$ and $\dot{ \epsilon}_{II}$ are the strain rate tensor invariants; $\theta=\textrm{atan}(\dot{\epsilon_{II}} / \dot{\epsilon_{I}} )$. $w_d(h,g)$, the deformation mode will be discussed in the next section. $\vert \dot{\epsilon} \vert \alpha_o$ and $\vert \dot{\epsilon} \vert \alpha_d$ are called the lead opening and closing rates, respectively.
+
+The \textbf{dynamical} inputs of the mechanical redistribution in LIM are:
+\begin{itemize}
+\item $\vert \dot{\epsilon} \vert \alpha_o$, the opening rate,
+\item $\vert \dot{\epsilon} \vert \alpha_d$, the net closing rate.
+\end{itemize}
+Following \cite{Thorndikeetal75}, we choose $\int_0^\infty w_d(h,g) = - 1$. In order to satisfy area conservation, the relation $\vert \dot{\epsilon} \vert \alpha_o - \vert \dot{\epsilon} \vert \alpha_d = \nabla \cdot \mathbf{u}$ must be verified. In the model, there are two ways to compute the divergence of the velocity field. A first way is to use the velocity components ($\dot{\epsilon}_I = \nabla \cdot \mathbf{u}\vert^{rhg}$) as computed after the rheology (superscript $rhg$). Another way is to derive it from the horizontal transport of ice concentration and open water fraction. In principle, the equality $A^o + \sum_{l=1}^L g^i_L = 1$ should always be verified. However, after ice transport (superscript $trp$), this is not the case, and one can diagnose a velocity divergence using the departure from this equality: $\nabla \cdot \mathbf{u}\vert^{trp} = ( 1 - A^o - \sum_{l=1}^L g^i_L )/ \Delta t$. In general, we will use $\dot{\epsilon}_I$ unless otherwise stated.
+
+The \textbf{net closing rate} is written as a sum of two terms representing the energy dissipation by shear and convergence \citep{FlatoHibler95}:
+\begin{equation}
+\vert \dot{\epsilon} \vert \alpha_d (\theta) = C_s \frac{1}{2} ( \Delta - \vert \dot{\epsilon}_I \vert) - \textrm{min} (\dot{\epsilon}_I,0),
+\end{equation}
+where $\Delta$ is a measure of deformation (defined in the rheology section). The factor $C_s=0.5$ (\textit{Cs} in \textit{namelist\_ice}) is added to allow for energy sinks other than ridge building (e.g., sliding friction) during shear. In case of convergence, the closing rate must be large enough to satisfy area conservation after ridging, so we take:
+\begin{equation}
+\vert \dot{\epsilon} \vert \alpha_d (\theta) = \text{max}(\vert \dot{\epsilon} \vert \alpha_d (\theta), - \nabla \cdot \mathbf{u} \vert^{trp}) \quad \text{if } \nabla \cdot \mathbf{u} < 0.
+\end{equation}
+
+The \textbf{opening rate} is obtained by taking the difference:
+\begin{equation}
+\vert \dot{\epsilon} \vert \alpha_o = \vert \dot{\epsilon} \vert \alpha_d = \nabla \cdot \mathbf{u}\vert^{trp}
+\end{equation}
+
+\section{The two deformation modes: ridging and rafting}
+
+The deformation mode is separated into ridging $w^{ri}$ and rafting $w^{ra}$ modes:
+\begin{equation}
+w^d(h,g) = w^{ri}(g,h) + w^{ra}(g,h).
+\end{equation}
+
+\textbf{Rafting} is the piling of two ice sheets on top of each other. Rafting doubles the participating ice thickness and is a volume-conserving process. \cite{Babkoetal02} concluded that rafting plays a significant role during initial ice growth in fall, therefore we included it into the model. 
+
+\textbf{Ridging} is the piling of a series of broken ice blocks into pressure ridges. Ridging redistributes participating ice on a various range of thicknesses. Ridging does not conserve ice volume, as pressure ridges are porous. Therefore, the volume of ridged ice is larger than the volume of new ice being ridged. In the model, newly ridged is has a prescribed porosity $p=30\%$ (\textit{ridge\_por} in \textit{namelist\_ice}), following observations \citep{Lepparantaetal95,Hoyland02}. The importance of ridging is now since the early works of \citep{Thorndikeetal75}.
+
+The deformation modes are formulated using \textbf{participation} and \textbf{transfer} functions with specific contributions from ridging and rafting:
+\begin{equation}
+w_d(h,g) = - [ b^{ra}(h) + b^{ri}(h) ] g(h) + n^{ra}(h)+n^{ri}(h).
+\label{eq:wd}
+\end{equation}
+$b^{ra}(h)$ and $b^{ri}(h)$ are the rafting and ridging participation functions. They determine which regions of the ice thickness space participate in the redistribution. $n^{ra}(h)$ and $n^{ri}(h)$, called transfer functions, specify how thin, deformation ice is redistributed onto thick, deformed ice. Participation and transfer functions are normalized in order to conserve area.
+
+\section{Participation functions}
+
+We assume that the participation of ice in redistribution does not depend upon whether the deformation process is rafting or ridging. Therefore, the participation functions can be written as follows:
+\begin{equation}
+b^{ra}(h)= \beta(h) b(h),
+\label{eq:bra}
+\end{equation}
+\begin{equation}
+b^{ri}(h) = [1-\beta(h)] b(h),
+\label{eq:bri}
+\end{equation}
+where $b(h)$ is an exponential weighting function with an e-folding scale $a*$ \citep{Lipscombetal07} (\textit{astar} in \textit{namelist\_ice}) which preferentially apportions the thinnest available ice to ice deformation:
+\begin{equation}
+b(h) = \frac{\text{exp} [ - G(h) / a^\star ]}{a^\star[1-\text{exp}(-1/a^\star)]} ,
+\end{equation}
+It is numerically more stable than the original version of \cite{Thorndikeetal75}. This scheme is still present in the code and can be activated using \textit{partfun\_swi} from \textit{namelist\_ice}, with the associated parameter \textit{Gstar}.
+
+$\beta(h)$ partitions deformation ice between rafted and ridged ice. $\beta(h)$ is formulated following \cite{Haapala00}, using the \cite{Parmerter75} law, which states that, under a critical participating ice thickness $h_P$, ice rafts, otherwise it ridges:
+\begin{equation}
+\beta(h)= \frac{ \textrm{tanh} [ -C_{ra} (h-h_P) ] + 1 } {2},
+\end{equation}
+where $C_{ra}=5$ m$^{-1}$ (\textit{Craft} in \textit{namelist\_ice})  and $h_P=0.75$ m (\textit{hparmeter} in \textit{namelist\_ice}) \citep{Haapala00,Babkoetal02}. The $tanh$ function is used to smooth the transition between ridging and rafting.  If namelist parameter \textit{raftswi} is set to 0, ice only ridges and does not raft. 
+
+\section{Transfer functions}
+
+The rafting transfer function assumes a doubling of ice thickness :
+\begin{equation}
+n^{ra}(h)=\frac{1}{2}\int_0^\infty \delta(h-2h') b(h') g(h') dh,
+\end{equation}
+where $\delta$ is the Dirac delta function. 
+
+The ridging transfer function is :
+\begin{equation}
+n^{ri}(h)=\int_0^\infty \gamma(h',h) (1+p) b(h') g(h') dh.
+\label{eq:nri}
+\end{equation}
+The redistributor $\gamma(h',h)$ specifies how area of thickness $h'$ is redistributed on area of thickness $h$. We follow \citep{Hibler80} who constructed a rule, based on observations, that forces all ice participating in ridging with thickness $h'$ to be linearly distributed between ice that is between $2h'$ and $2\sqrt{H^*h'}$ thick, where $H^\star=100$ m (\textit{Hstar} in \textit{namelist\_ice}). This in turn determines how to construct the ice volume redistribution function $\Psi^v$. Volumes equal to participating area times thickness are removed from thin ice. They are redistributed following Hibler's rule. The factor $(1+p)$ accounts for initial ridge porosity $p$ (\textit{ridge\_por} in \textit{namelist\_ice}, defined as the fractional volume of seawater initially included into ridges. In many previous models, the initial ridge porosity has been assumed to be 0, which is not the case in reality since newly formed ridges are porous, as indicated by in-situ observations \citep{Lepparantaetal95,Hoyland02}. In other words, LIM3 creates a higher volume of ridged ice with the same participating ice.
+
+For the numerical computation of the integrals, we have to compute several temporary values:
+\begin{itemize}
+\item The thickness of rafted ice $h^{ra}_l = 2 h^i_l$
+\item The mean thickness of ridged ice $h^{ri,mean}_l=\text{max} ( \sqrt{H^\star h^i_l} , h^i_l \cdot 1.1 )$
+\item The minimum thickness of ridged ice $h^{ri,min}_l = \text{min} [ 2*h^i_l, 0.5\cdot(h^{ri,mean}_l+h^i_l) ]$
+\item The maximum thickness of ridged ice $h^{ri,min}_l = 2 h^{ri,mean}_l - h^{ri,min}_l$
+\item The mean rate of thickening of ridged ice $k^{ri}_l = h^{ri,mean}_l / h^{i}_l$
+\end{itemize}
+
+\section{Ridging shift}
+
+%\textit{Subroutine}: \textbf{lim\_itd\_ridgeshift}
+
+The numerical computation of the impact of mechanical redistribution on ice concentration involves:
+\begin{itemize}
+\item A normalization factor that ensures volume conservation (\textit{aksum})
+\item The removal of ice participating in deformation (including the closing of open water)
+\item The addition of deformed ice
+\end{itemize}
+
+For ice concentrations, the numerical procedure reads:
+\begin{equation}
+\Delta g^i_l = C^{net} \Delta t \biggr [ - ( b^{ri}_l + b^{ra}_l) + \sum_{l_2=1}^L \biggr ( f^{ra}_{l,l_2}  \frac{b^{ra}_{l_2}}{k^{ra}} + f^{ri}_{l,l_2} \frac{ b^{ri}_{l_2} }{k^{ri}_{l_2}} \biggr ) \biggr ]
+\end{equation}
+\begin{itemize}
+\item $C^{net}$ is the normalized closing rate ($\vert \dot{\epsilon} \vert \alpha^d/\textit{aksum}$)
+\item $b^{ri}_l$ and $b^{ra}_l$ are the area participating into redistribution for category $l$
+\item $f^{ra}_{l,l_2}$ and $f^{ri}_{l,l_2}$ are the fractions of are of category $l$ being redistributed into category $l_2$
+\item $k^{ra}$ is the rate of thickening of rafted ice (=2)
+\end{itemize}
+
+Because of the nonlinearities involved in the integrals, the ridging procedure has to be iterated until $A^\star = A^{ow} + \sum_{l=1}^L g^i_l = 1$.
+
+\section{Mechanical redistribution for other global ice variables}
+
+The other global ice state variables redistribution functions $\Psi^X$ are computed based on $\Psi^g$ for the ice age content and on $\Psi^{v^i}$ for the remainder (ice enthalpy and salt content, snow volume and enthalpy). The general principles behind this derivation are described in Appendix A of \cite{Bitzetal01}. A fraction $f_s=0.5$ (\textit{fsnowrdg} and \textit{fsnowrft} in \textit{namelist\_ice}) of the snow volume and enthalpy is assumed to be lost during ridging and rafting and transferred to the ocean. The contribution of the seawater trapped into the porous ridges is included in the computation of the redistribution of ice enthalpy and salt content (i.e., $\Psi^{e^i}$ and $\Psi^{M^s}$). During this computation, seawater is supposed to be in thermal equilibrium with the surrounding ice blocks. Ridged ice desalination induces an implicit decrease in internal brine volume, and heat supply to the ocean, which accounts for ridge consolidation as described by \cite{Hoyland02}. The inclusion of seawater in ridges does not imply any net change in ocean salinity. The energy used to cool down the seawater trapped in porous ridges until the seawater freezing point is rejected into the ocean.
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_single_category_use.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_single_category_use.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_single_category_use.tex	(revision 10845)
@@ -0,0 +1,22 @@
+
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 10 Ñ Single-category use
+% ================================================================
+
+\chapter{Using SI$^3$ with a single category}
+
+\label{chap:INT}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+\section{Enhanced conduction}
+\section{Virtual thin ice melting (fake lateral melting)}
+\section{Ridging and rafting with a single ice category}
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_thermo.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_thermo.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_thermo.tex	(revision 10845)
@@ -0,0 +1,88 @@
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 7 Ñ Ice thermodynamics
+% ================================================================
+
+\chapter{Ice thermodynamics}
+\label{chap:THD}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+
+% ================================================================
+% Thermodynamics % ================================================================
+
+Mass, energy and salt are intertwined for sea ice. Referred to as thermodynamics, or more exactly halo-thermodynamics. Mushy-layer theory.
+
+\section{Open water and new ice formation}
+
+% There is everything in here, but there are too much details.
+% Move things around into chapters and summarize.
+
+As part of the sea ice computations, a heat budget of the uppermost oceanic level is estimated. This heat budget is used if negative to compute the production of new ice or, if positive, for bottom melting. 
+
+%--------------------------------------------------------------------------------------------------------------------
+%
+% FIG x : Open water energy budget
+%
+\begin{figure}[ht]
+\begin{center}
+\vspace{0cm}
+\includegraphics[height=10cm,angle=-00]{../Figures/Openwater_eb.png}
+\caption{Scheme of the estimate of the heat budget of the first ocean level.}
+\label{fig_yield}
+\end{center}
+\end{figure}
+%
+%--------------------------------------------------------------------------------------------------------------------
+
+We estimate the heat budget of the first ocean level ($B^{opw}$) assuming four contributions, namely:
+\begin{itemize}
+\item The absorption of a fraction $f_1^{qsr}$ of solar radiation (given by radiative transfer component of the ocean model);
+\item The non-solar heat flux absorbed at the surface;
+\item The sensible heat content of precipitation
+\item The sensible heat loss to the sea ice
+\end{itemize}
+The estimated heat budget thus reads:
+\begin{linenomath}
+\begin{align}
+B^{opw} = Q^{sr} ( 1 - A ) f_1^{qsr} + Q^{ns} ( 1 - A ) + Q^{emp} - A F_w,
+\end{align}
+\end{linenomath}
+Hence there is no consideration of the entrainment of heat at the base of the first ocean level, or of solar radiation transmitted below the ice. The ocean-to-ice turbulent sensible heat flux is formulated following \citep{McPhee92}
+\begin{linenomath}
+\begin{align}
+F_w = \rho_0 c_w C_h u^* ( SST - T_b )
+\end{align}
+\end{linenomath}
+where $\rho_0$ is the reference ocean density, $c_w$ is the seawater specific heat, $C_h = 7.5 \times 10^{-3}$ is a heat transfer coefficient, and $u^*=\sqrt{\tau_{iw}/\rho_0}$. There are two additional conditions, the oceanic heat flux cannot be negative. Second, $F_w$ cannot exceed the heat content of the first ocean level.
+
+% Ice growth
+If the $B^{opw}$ is such that the SST would decrease below the freezing point, the remainder of the heat is used to form new ice. The heat loss is converted into mass through [\ref{eq_phasechange}], giving by multiplication by density a volume of new ice $v_0$. The thickness $h_0$ of the new ice grown during a sea ice time step depends on unresolved small currents and waves and is prescribed. The fraction $a_0=v_0/h_0$ is computed accordingly. The salinity of this new ice $S_0$ is given by the S-h empirical relationship of \cite{Kovacs96}. The temperature assumed for this new ice is the local freezing point.
+
+% Ice melt
+If there is ice in the grid cell and that the $B^{opw}$ is positive, it is directly given attributed to the heat available for bottom melting. This argument follows from \cite{MaykutMcPhee95}, who found that most of solar heat absorbed in the surface waters is converted into melting. In practise, this means that the SST cannot go above freezing as long ice is present.
+
+The heat loss used for ice formation, heat gain used to melt ice and the sensible heat given by the ocean to the ice, are all removed from the non-solar heat flux transmitted to the ocean.
+
+Because ice dynamics are not able to maintain the small amount of open water that is observed, a maximum ice fraction ($amax, <1$) is prescribed.
+
+\section{Diffusion of heat}
+
+\section{Vertical growth and melt}
+
+\section{Desalination}
+
+\section{Remapping}
+
+\section{Transport in thickness space}
+
+\section{True lateral melting}
+
+\end{document}
+
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_transport.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_transport.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/chap_transport.tex	(revision 10845)
@@ -0,0 +1,27 @@
+
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Chapter 4 Ñ Ice transport
+% ================================================================
+
+\chapter{Ice transport}
+\label{chap:TRP}
+\minitoc
+
+\newpage
+$\ $\newline    % force a new line
+
+As soon as ice dynamics are activated (\textit{ln\_dyn\_xxx}), all extensive state variables are to be advected following the horizontal velocity field.
+
+\section{Second order moments conserving (Prather 1986) scheme (\textit{ln\_adv\_Pra})}
+
+The scheme of  \cite{Prather86} explicitly computes the conservation of second-order moments of the spatial distribution of global sea ice state variables. This scheme preserves positivity of the transported variables and is practically non-diffusive. It is also computationally expensive, however it allows to localize the ice edge quite accurately. As the scheme is conditionally stable, the time step is split into two parts if the ice drift is too fast, based on the CFL criterion. 
+
+State variables per unit grid cell area are first multiplied by grid cell area. Then, for each state variable, the 0$^{th}$ (mean), 1$^{st}$ (x, y) and 2$^{nd}$ (xx, xy, yy) order moments of the spatial distribution are transported. At 1st time step, all moments are zero (if prescribed initial state); or read from a restart file, and then evolve through the course of the run. Therefore, for each global variable, 5 additional tracers have to be kept into memory and written in the restart file, which significantly increases the required memory. Advection following x and y are computed independently. The succession order of x- and y- advection is reversed every day.
+
+\section{5$^{th}$ order flux-corrected transport scheme (UM5)}
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/introduction.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/introduction.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/introduction.tex	(revision 10845)
@@ -0,0 +1,59 @@
+
+% ================================================================
+% INTRODUCTION
+% ================================================================
+
+\chapter*{Introduction}
+
+\textcolor{red}{[ \textit{July 2018} ]} \\
+
+%Near the poles of the Earth, the seas and oceans freeze when seawater at the freezing point loses heat. The resulting forms of saline ice are collectively called \textit{sea ice} \citep{WMO70}, reaching up to a few meters in thickness, where as sea ice coverage is about 5\% of the global ocean, about 30 million square kilometers. All sea ice characteristics vary over a wide range of spatio-temporal scales, reflecting changes in heat, mass and momentum exchanges with the atmosphere and the ocean; the clearest temporal signal being an ample seasonal cycle. Sea ice formation and melting affects water mass formation in the ocean \citep{GoosseFichefet99}. It not only impacts, but also reflects the state of the climate system \citep{Budyko69,NotzStroeve16}.  Sea ice also affects marine life, water chemistry and human activities in polar regions. Local populations use sea ice for travelling and hunting, whereas navigation and resource exploitation are dependent on sea ice conditions. For such reasons, ocean modelling systems, including NEMO, must include a sea ice component.
+
+The sea Ice Modelling Integrated Initiative (SI$^3$) is the sea ice engine of the Nucleus for European Modelling of the Ocean (NEMO). It is intended to be a flexible tool for studying sea ice and its interactions with the other components of the Earth System over a wide range of space and time scales. SI$^3$ is a curvilinear grid, finite-difference implementation of the classical AIDJEX\footnote{AIDJEX=\textbf{A}rctic \textbf{I}ce \textbf{D}ynamics \textbf{J}oint \textbf{EX}periment} model \citep{Coonetal74}, combining the conservation of momentum for viscous-plastic continuum, energy and salt-conserving halo-thermodynamics, an explicit representation of subgrid-scale ice thickness variations, snow and melt ponds. An option to switch back to the \textit{single-category} (or \textit{2-level}) framework of \cite{Hibler79} provides a cheap sea ice modelling solution.
+
+SI$^3$ is the result of the recommendation of the Sea Ice Working Group (SIWG) to reduce duplication and better use development resources. SI$^3$ merges the capabilities of the 3 formerly used NEMO sea ice models (CICE, GELATO and LIM). The \textbf{3} in SI$^3$ refers to the three formerly used sea ice models. It also refers to linkages between 3 different media (ocean, ice, snow). The model can be spelt 'SI3' in situations where the superscript could be problematic (i.e., within code and svn repository etc.) The model name would be pronounced as 'si-cube' for short (or 'sea ice cubed' for slightly longer). 
+
+% Limitations & scope
+%There are limitations to the applicability of models such as SI$^3$. The continuum approach is not invalid for grid cell size above at least 1 km, below which sea ice particles may include just a few floes, which is not sufficient \citep{Lepparanta05}. Second, one must remember that our current knowledge of sea ice is not as complete as for the ocean: there are no fundamental equations such as Navier Stokes equations for sea ice. Besides, important features and processes span widely different scales, such as brine inclusions (1 $\mu$m-1 mm) \citep{PerovichGow96}, horizontal thickness variations (1 m-100 km) \citep{Percivaletal08}, deformation and fracturing (10 m-1000 km) \citep{Marsanetal04}. These impose complicated and often subjective subgrid-scale treatments. All in all, there is more empirism in sea ice models than in ocean models. 
+
+In order to handle all the subsequent required subjective choices, we applied the following guidelines or principles:
+\begin{itemize}
+\item Sea ice is frozen seawater that is in tight interaction with the underlying ocean. This close connexion suggests that the sea ice and ocean model components must be as consistent as possible. In practice, this is materialized by the close match between LIM and NEMO, in terms of numerical choices, regarding the grid (Arakawa C-type) and the numerical discretization (finite differences with NEMO scale factors).
+\item It is useful to be able to either prescribe the atmospheric state or to use an atmospheric model. For consistency and simplicity of the code, we choose to use formulations as close as possible in both cases.
+\item Different resolutions and time steps can be used. There are parameters that depend on such choices. We thrieved to achieve a resolution and time-step independent code, by imposing a priori scaling on the resolution / time step dependence of such parameters.
+\item Energy, mass and salt must be conserved as much as possible.
+\end{itemize}
+
+This manual is organised as follows. \\
+
+List of chapters... \\
+
+There are no more CPP keys in the code. \\
+
+Namelists and output management follow NEMO guidelines. \\
+
+Changes between releases. \\
+
+The list of people that should be acknowledged is too long, but a great number of people or more exactly a number of great people contributed to the code and should be gratefully acknowledged. As for today,
+the SIWG members are \textit{(July 2018)}.
+
+\begin{itemize}
+\begin{footnotesize}
+\item Yevgeny Aksenov (NOCS, Southampton, UK)
+\item Ed Blockley (Met Office, Exeter, UK, co-chair)
+\item Matthieu Chevallier (CNRM-GAME, M\'et\'eo France, Toulouse)
+\item Danny Feltham (CPOM, Reading, UK)
+\item Thierry Fichefet (UCL, Louvain-la-Neuve, Belgium)
+\item Gilles Garric (Mercator-Oc\'ean, Toulouse, France)
+\item Paul Holland (BAS, Cambridge, UK)
+\item Dorotea Iovino (CMCC, Bologna, Italy)
+\item Gurvan Madec (LOCEAN, CNRS, Paris, France)
+\item Fran\c cois Massonnet (UCL, Louain-la-Neuve, Belgium)
+\item Jeff Ridley (Met Office, Exeter, UK)
+\item Cl\'ement Rousset (LOCEAN, CNRS, Paris, France)
+\item David Salas (CNRM-GAME, M\'et\'eo France, Toulouse)
+\item David Schroeder (CPOM, Reading, UK)
+\item Steffen Tietsche (ECMWF, Reading, UK)
+\item Martin Vancoppenolle (LOCEAN, CNRS, Paris, France, co-chair) 
+\end{footnotesize}
+\end{itemize}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/todolist.tex
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/todolist.tex	(revision 10845)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/doc/si3_doc/tex_sub/todolist.tex	(revision 10845)
@@ -0,0 +1,56 @@
+\documentclass[../../tex_main/NEMO_manual]{subfiles}
+
+\begin{document}
+
+% ================================================================
+% Dummy Chapter - To do list
+% ================================================================
+
+\chapter*{To-do-list for documentation}
+\label{chap:TDL}
+
+\begin{itemize}
+
+\item \textcolor{gray}{Documentation infrastructure ready by Jul 19, 2018.}
+
+\item \textcolor{gray}{Abstract: v1 written.}
+
+\item \textcolor{gray}{Disclaimer: v1 written.}
+
+\item Introduction: list of chapters and change between realeases to be written (1h)
+
+\item \textcolor{gray}{1. Model basics: v1 written.}
+
+\item \textcolor{gray}{Add namelist infrastructure as in NEMO doc (done)}
+
+\item \textcolor{gray}{2. Time, space and thickness space domain}
+
+\item 3. Ice dynamics: to be written (Clem \& Martin, 3h)
+
+\item 4. Ice transport: UM5 description to be written (Clem, 3h)
+
+\item 5. Ridging and rafting: V1 from previous doc. Should be refreshed. A bit too long. All namelist parameters are not explained (namdyn\_rdgrft).
+
+\item \textcolor{gray}{6. Radiative transfer: Full rewriting required (Martin \& Clem, 3-6h). Available elements from previous documentation, but mostly obsolete.}
+
+\item 7. Thermodynamics: Lots of rewriting required. Available elements from previous documentation, but mostly obsolete (Martin, 6h)
+
+\item 8. Interfaces: Full writing required 
+
+\item 9. Output and diagnostics: Full writing required
+
+\item 10. Single category (Martin, 3h). Full writing required.
+
+\item 11. BDY-AGRIF (Clem, 3h). Full writing required.
+
+\item 12. Miscellaneous. Full writing required.
+
+\item Melt ponds are not currently considered in the plan.
+
+\item Make sure all namelist parameters are properly described.
+
+\item Peer-review by SIWG members.
+
+\end{itemize}
+
+\end{document}
Index: /NEMO/branches/UKMO/dev_r10037_GPU/src/OCE/ZDF/zdftke.F90
===================================================================
--- /NEMO/branches/UKMO/dev_r10037_GPU/src/OCE/ZDF/zdftke.F90	(revision 10844)
+++ /NEMO/branches/UKMO/dev_r10037_GPU/src/OCE/ZDF/zdftke.F90	(revision 10845)
@@ -707,5 +707,5 @@
       IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
       IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
-      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
+      IF( ( nn_htau < 0   .OR.  nn_htau > 1 ) .AND. nn_htau .NE. 4 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 4 ' )
       IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
       !
@@ -725,4 +725,14 @@
          CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
             htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
+         CASE( 4 )                                 ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( gphit(ji,jj) <= 0._wp ) THEN
+                     htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
+                  ELSE
+                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
+                  ENDIF
+               END DO
+            END DO
          END SELECT
       ENDIF
