Index: /utils/figures/SI3/Figures/scripts/Thermal_properties.pro
===================================================================
--- /utils/figures/SI3/Figures/scripts/Thermal_properties.pro	(revision 11010)
+++ /utils/figures/SI3/Figures/scripts/Thermal_properties.pro	(revision 11010)
@@ -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: /utils/figures/SI3/Figures/scripts/albedo_cloud_correction.m
===================================================================
--- /utils/figures/SI3/Figures/scripts/albedo_cloud_correction.m	(revision 11010)
+++ /utils/figures/SI3/Figures/scripts/albedo_cloud_correction.m	(revision 11010)
@@ -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: /utils/figures/SI3/Figures/scripts/albedo_hi.m
===================================================================
--- /utils/figures/SI3/Figures/scripts/albedo_hi.m	(revision 11010)
+++ /utils/figures/SI3/Figures/scripts/albedo_hi.m	(revision 11010)
@@ -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: /utils/figures/SI3/Figures/scripts/ice_cats.pro
===================================================================
--- /utils/figures/SI3/Figures/scripts/ice_cats.pro	(revision 11010)
+++ /utils/figures/SI3/Figures/scripts/ice_cats.pro	(revision 11010)
@@ -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
+
