;
;---------------------------------------------------
; This script makes plots of biological simulations
; (c) Martin Vancoppenolle, UCL-ASTR, june 2007
; reviewed may 2008
; works for BARROW-ISPOL
;---------------------------------------------------
;

;
;------------------------------------------------------------------------------
; Read input data in exp_id.dat
;------------------------------------------------------------------------------
;

exp_id = 'E001'
indir  = ''
outdir = ''
dt     = 86400
c_bio_model = 'BFMSI'
site        = 'TOURNAI'

OPENR, 11, 'exp_id.dat'
READF, 11, exp_id & READF, 11, indir & READF, 11, outdir & READF, 11, dt & READF, 11, c_bio_model & READF, 11, site & CLOSE, 11

file_name = exp_id+'_contours.ps'

numplot_x = 6       ; number of horizontal plots
numplot_y = 3       ; number of vertical plots

device = 'PS'       ; 'PS' or 'X'

;
;------------------------------------------------------------------------------
; DEFINE OUTPUT PLOT
;------------------------------------------------------------------------------
;
figuresize_x = numplot_x * 4.0 ; figure size on x direction
figuresize_y = numplot_y * 5.0 ; figure size on y direction
set_plot, device
IF ( device EQ 'PS' ) THEN BEGIN
   device, /COLOR, /LANDSCAPE, filename=outdir+file_name, $
           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

!X.MARGIN=[6,3]
!Y.MARGIN=[6,3]

;
;------------------------------------------------------------------------------
; Read files
;------------------------------------------------------------------------------
;

 open_fields_v7, indir, exp_id, c_bio_model, $
                 numt, doy, ts_d, ts_m,                        $
                 h_i, h_s, z_ip, z_ib, s_i, t_i, e_i, PAR, Ra, $
                 dhib, dhisu, dhisi,                           $
                 AoCb, NO3b, PO4b, DSib, chla, eoCb,           $
                 Argb, Argbub, $
                 Oxyb, Oxybub, $
                 DICb, Alkb, CO2b, CO2bub,                     $
                 CO2aq, HCO3m, CO32m, pH, pCO2, Ikab,          $
                 dFeb, aFeb, eFeb,                             $
                 AoCt, NO3t, PO4t, DSit, eoCt, Argt,           $
                 DICt, Alkt, CO2t, Ikat,                       $
                 dFet, aFet, eFet,                             $
                 FDAbd, FDAbpos, FDAsi,                        $
                 FNO3, FNO3bpos, FNO3si,                       $
                 FCO2_atm, FCO2_bub,                           $
                 syn, rsp, lys, rem,                           $
                 lim_lig, lim_no3, lim_po4, lim_dsi, lim_tem, lim_sal
;DEVICE, /CLOSE
;SET_PLOT, "X"
;!P.MULTI=[0,2,2]
;STOP
                
; conversion for carbon units if KRILL is used
IF ( c_bio_model EQ 'KRILL' ) THEN BEGIN

zunit = 12.
AoCb = AoCb * zunit
eoCb = eoCb * zunit
AoCt = AoCt * zunit
eoCt = eoCt * zunit
syn    = syn       * zunit * 86400.
lys    = lys       * zunit * 86400.
rem    = rem       * zunit * 86400.

ENDIF

;
;------------------------------------------------------------------------------
; Extract size of arrays & time axis in months
;------------------------------------------------------------------------------
;
zdummy = SIZE(z_ip)
nlay = zdummy(1)

zdummy = SIZE(z_ib)
nlay_bio = zdummy(1)

zdummy = SIZE(numt)
PRINT, size(numt)
nts = zdummy(1)
numy = LONG ( (nts+364) / 365 )     ; number of years

PRINT, ' nlay     : ', nlay
PRINT, ' nlay_bio : ', nlay_bio
PRINT, ' nts      : ', nts

IF ( site EQ 'BARROW' ) THEN BEGIN
   mon = numt * dt / 86400. / 30.4166 + 1 ; time axis (Barrow)
ENDIF

IF ( site EQ 'ISPOL' ) THEN BEGIN
   mon = numt / 30.41666666 + 1
ENDIF

IF ( site EQ 'YROSIAE' ) THEN BEGIN
   mon = numt / 30.41666666 + 1
ENDIF

;
;------------------------------------------------------------------------------
; OBS
;------------------------------------------------------------------------------
;
;  OPEN FILE

IF ( site EQ 'BARROW' ) THEN BEGIN

   data_dir = '/Users/ioulianikolskaia/Boulot/SCIENCE/PLOT_SCRIPTS/LIM1D_BIO/IDL/DATA/'
   file_name_2 = 'obs_barrow_seb.txt'
   OPENR, 10, data_dir+file_name_2 ; 'ouverture de l'accès au fichier'
   num_lines = 589 ; nombre de lignes
   jd_obs    = FLTARR(num_lines)
   annee_obs    = FLTARR(num_lines)
   jour_entier = FLTARR(num_lines)
   heure_obs    = FLTARR(num_lines)
   hi_obs    = FLTARR(num_lines)
   hs_obs    = FLTARR(num_lines)
   
   FOR i = 0, num_lines - 1 DO BEGIN
      READF, 10, var1, var2, var3, var4, var5, var6
      jd_obs   (i) = var1 & annee_obs   (i) = var2 & jour_entier(i) = var3
      heure_obs   (i) = var4 & hs_obs   (i) = var5 & hi_obs   (i) = var6
   ENDFOR
   mon_obs    = jd_obs / 30.467 + 1
   
   CLOSE, 10

ENDIF

IF ( site EQ 'ISPOL' ) THEN BEGIN

   data_dir = '/Users/ioulianikolskaia/Boulot/SCIENCE/PLOT_SCRIPTS/LIM1D_BIO/IDL/DATA/'
   file_no = NCDF_OPEN(data_dir+'ISPOL.nc', /NOWRITE)
   ;NCDF_VARGET, file_no, h_i_no, h_i
   day_obs_no  = NCDF_VARID( file_no, 'nday'  )
   hi_ulb_no   = NCDF_VARID( file_no, 'hi_ULB'  )
   hs_ulb_no   = NCDF_VARID( file_no, 'hs_ULB'  )
   hi_awi_no   = NCDF_VARID( file_no, 'hi_AWI'  )
   hs_awi_no   = NCDF_VARID( file_no, 'hs_AWI'  )
   dsi_obs_no  = NCDF_VARID( file_no, 'dsi' )
   s_i_obs_no  = NCDF_VARID( file_no, 'sal_phy' )
   
   NCDF_VARGET, file_no, day_obs_no, day_obs
   NCDF_VARGET, file_no, hi_ulb_no, hi_ulb
   NCDF_VARGET, file_no, hs_ulb_no, hs_ulb
   NCDF_VARGET, file_no, hi_awi_no, hi_awi
   NCDF_VARGET, file_no, hs_awi_no, hs_awi
   NCDF_VARGET, file_no, dsi_obs_no, dsi_obs
   NCDF_VARGET, file_no, s_i_obs_no, s_i_obs
   
   NCDF_CLOSE, file_no
   
   hi_obs = ( hi_ulb + hi_awi ) / 2.0
   hs_obs = ( hs_ulb + hs_awi ) / 2.0

   mon_obs    = day_obs / 30.467 + 1

ENDIF


IF ( site EQ 'YROSIAE' ) THEN BEGIN
   hi_obs =     [ 1.472, 1.609, 1.670, 1.690, 1.708, 1.714, 1.711, 1.677, 1.693 ]
   hi_obs_std = [ 0.023, 0.014, 0.009, 0.011, 0.011, 0.009, 0.016, 0.024, 0.014 ]
   hs_obs = REPLICATE(0.01, 9)
   day_obs =[ 263, 279, 292, 299, 306, 312, 318, 327, 355 ]
   mon_obs    = day_obs / 30.467 + 1
ENDIF

;
;------------------------------------------------------------------------------
; PLOTS
;------------------------------------------------------------------------------
;
!P.MULTI=[0,numplot_x, numplot_y]

; initialize plot positions
i_plot = -1
horpos = 0.02
vertpos = 0.015
dv = 1./FLOAT(numplot_y)
dh = 1./FLOAT(numplot_x)
leveling = 1
cb_length = 0.8*dh


IF ( site EQ 'ISPOL' ) THEN BEGIN
   zmin = -1.30
   zmax = 0.30
ENDIF

IF ( site EQ 'BARROW' ) THEN BEGIN
   zmin = -1.5
   zmax = 0.5
ENDIF

IF ( site EQ 'YROSIAE' ) THEN BEGIN
   zmin = -2.0
   zmax = 0.0
ENDIF

;-------------------
; Salinity contours
;-------------------
i_plot   = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 13
levels = [ 0,1,2,3,4,5,6,7,8,9,10,11,12]
colors = [ 241,221,201,181,161,141,121,101,81,61,41,21,1]

ct = 1 ; colortable
tmin = MIN(mon)
tmax = MAX(mon)

title = ' Salinity '

icecontour, s_i, h_i, h_s, nlevels, colors, levels, ct, z_ip, mon, nlay, nts, $
            tmin, tmax, zmin, zmax, title

LOADCT, 13
oplot, mon_obs, -hi_obs, PSYM=1, SYMSIZE=1.0, THICK = 3, COLOR =255
oplot, mon_obs, hs_obs, PSYM=1, SYMSIZE=1.0, THICK = 3, COLOR = 255
LOADCT, ct

; colorbar
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
format='(f5.0)'
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;--------------
; Brine volume
;--------------
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 11
levels = [   0.,   2.,  4.,   5.,   8.,  10.,  15.,20., 25., 30., 50.]
colors = [   20,  60,  100,  150,  160,  170, 180, 190, 200, 210, 250 ]
;colors = [  150, 140,  130,  100,   90,   80,  70,  50,  40,  30,  20 ]

ct = 33 ; colortable
title = ' Brine volume '

icecontour, e_i, h_i, h_s, nlevels, colors, levels, ct, z_ip, mon, nlay, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
format='(i3)'
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;-----------------
; Rayleigh number
;-----------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 11
levels = [ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10. ]
colors = [ 20,  40,   60,   80,  100,  160, 180, 200, 220, 240, 250 ]

ct = 11 ; colortable
title = ' Rayleigh number '

; contours
icecontour, Ra, h_i, h_s, nlevels, colors, levels, ct, z_ip, mon, nlay, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(i2)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;--------------
; PAR         
;--------------
;
!P.MULTI=[0,numplot_x, numplot_y]
i_plot = 0
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
levels = [  0.0,0.015, 0.1,  1.0,  10.,  15.,  20.,  30., 50., 80.,120., 150. ]
levels = [  0.0,0.015, 0.1,  1.0,  10.,  15.,  20.,  30., 40., 50., 60.,  70. ]
levels = [  0.0,0.05, 0.1,  2.5,  5., 7.5,    10.,  12.5, 15., 17.5, 20. ]
colors = [   20,  40,   60,   80,  140,  160, 180, 190, 200, 210,  220 ]

ct = 3 ; colortable
title = ' PAR (muE/m2/s)'

icecontour, PAR , h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

LOADCT, 1
oplot, mon_obs, -hi_obs, PSYM=1, SYMSIZE=1.0, THICK = 3, COLOR = 150
oplot, mon_obs, hs_obs, PSYM=1, SYMSIZE=1.0, THICK = 3, COLOR = 150
LOADCT, ct

; colorbar
format='(f5.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;----------------------
; Temperature contours
;----------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
levels = [ -25., -18, -16., -14., -12., -10., -8., -6., -4., -2., 0. ]
;levels = [ -10., -9.,  -8.,  -7.,  -6.,  -5., -4., -3., -2., -1., 0. ]
colors = [   20,  40,   60,   80,  100,  120, 140, 160, 180, 200, 220 ]

ct = 3 ; colortable
title = ' Temperature '

icecontour, t_i, h_i, h_s, nlevels, colors, levels, ct, z_ip, mon, nlay, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(f5.0)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;-------------------------
; Brine Salinity contours
;-------------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

sigma = s_i * 100. / e_i
sigma = - t_i / 0.053

nlevels = 10
fmin = 0. & fmax = 300.
cmin = 75 & cmax = 250
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

levels(nlevels) = 400. ; trick trick to avoid meaningless huge values on ze plot

ct = 1 ; colortable
title = ' Brine Salinity '

icecontour, sigma, h_i, h_s, nlevels, colors, levels, ct, z_ip, mon, nlay, nts, $
            tmin, tmax, zmin, zmax, title

levels(nlevels) = 300. ; trick trick

; colorbar
format='(i3)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;---------------
; Dissolved NO3 
;---------------
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
colors = [   20,  40,   60,   80,  100,  120, 140, 160, 180, 200, 220 ]
IF ( site EQ 'BARROW' ) THEN BEGIN
   levels = [   0.,   0.2,  0.4,   0.6,   0.8,  1.0,  1.2, 1.4, 1.6, 1.8, 2.0]
ENDIF
IF ( site EQ 'ISPOL' ) THEN BEGIN
   levels = [   0.,     2.,   4.,    6.,    8.,  10.,   12.,14., 16., 18., 20.]
   nlevels = 7.
   levels = [ 0., 1.6, 4., 6., 8., 10., 12., 20.]
   colors = [ 20, 60, 100, 140, 180, 220, 240 ]
ENDIF
IF ( site EQ 'YROSIAE' ) THEN BEGIN
   nlevels = 7.
   levels = [ 0., 1.6, 4., 6., 8., 10., 12., 20.]
   colors = [ 20, 60, 100, 140, 180, 220, 240 ]
ENDIF

ct =  5 ; colortable
title = ' NO3 (mmol/m3)'

icecontour, no3b, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I2)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;---------------
; Dissolved PO4 
;---------------
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 11
IF ( site EQ 'BARROW' ) THEN BEGIN
   levels = [   0.,   0.1,  0.2,   0.3,   0.4,  0.5,  0.6, 0.7, 0.8, 0.9, 1.0]
   colors = [   20,  40,   60,   80,  100,  120, 140, 160, 180, 200, 220 ]
ENDIF
IF ( site EQ 'ISPOL' ) THEN BEGIN
   levels = [   0.,   0.2,  0.4,   0.6,   0.8,  1.0,  1.2, 1.4, 1.6, 1.8, 2.0]
   nlevels = 7
   levels = [ 0., 0.24, 0.4, 0.6, 0.8, 1., 2. ]
   colors = [ 20, 60, 100, 140, 180, 220, 240 ]
ENDIF
IF ( site EQ 'YROSIAE' ) THEN BEGIN
   nlevels = 7
   levels = [ 0., 0.24, 0.4, 0.6, 0.8, 1., 2. ]
   colors = [ 20, 60, 100, 140, 180, 220, 240 ]
ENDIF

ct = 5 ; colortable
title = ' PO4 (mmol/m3)'

icecontour, PO4b, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(f4.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;--------------
; Dissolved Si
;--------------
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 11
;levels = [   0.,   2.,  4.,   6.,   8.,  10.,  12.,14., 16., 18., 20.]
;levels = [   0.,     4.,   8.,   12.,   16.,  20.,   24.,28., 32., 36., 40.]
levels = [   0.,   1.,  2.,   3.,   4.,  5.,  6., 7., 8., 9., 10.]
IF ( site EQ 'BARROW' ) THEN BEGIN
   levels = [   0., 0.5,   1.,  1.5,  2.0, 2.5, 3.0,3.5, 4., 4.5, 5.]
   colors = [   20,  40,   60,   80,  100,  120, 140, 160, 180, 200, 220 ]
ENDIF
IF ( site EQ 'ISPOL' ) THEN BEGIN
   levels = [   0.,   2.,  4.,   6.,   8.,  10.,  12.,14., 16., 18., 20.]
   nlevels = 7.
   levels = [ 0., 3.9, 8., 12., 16., 20., 24., 40.]
   colors = [ 20, 60, 100, 140, 180, 220, 240 ]
ENDIF
IF ( site EQ 'YROSIAE' ) THEN BEGIN
   nlevels = 7.
   levels = [ 0., 3.9, 8., 12., 16., 20., 24., 40.]
   colors = [ 20, 60, 100, 140, 180, 220, 240 ]
ENDIF

ct = 5 ; colortable
title = ' DSi (mmol/m3)'

icecontour, dsib, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(f4.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;------------------
; light limitation
;------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 1.001
cmin = 0 & cmax = 250
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3  ; colortable
title = ' Light limitation '

; contours
icecontour, lim_lig, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(f3.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;------------------
; Temp limitation
;------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 1.0001
cmin = 0 & cmax = 250
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3  ; colortable
title = ' Temp limitation '

; contours
icecontour, lim_tem, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f3.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format
;
;------------------
; Salt limitation
;------------------
;

i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 1.0001
cmin = 0 & cmax = 250
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3  ; colortable
title = ' Salt limitation '

; contours
icecontour, lim_sal, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f3.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;------------------
; NO3 limitation
;------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 1.0001
cmin = 0 & cmax = 250
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3 ; colortable
title = ' NO3 limitation '

; contours
icecontour, lim_no3, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f3.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;------------------
; PO4 limitation
;------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 1.0001
cmin = 0 & cmax = 250
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3 ; colortable
title = ' PO4 limitation '

; contours
icecontour, lim_po4, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f3.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;------------------
; DSI limitation
;------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 1.0001
cmin = 0 & cmax = 250
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3 ; colortable
title = ' DSi limitation '

; contours
icecontour, lim_dsi, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f3.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;------------------
; 3-nut limitation
;------------------
;

zlim_nut = FLTARR(nlay, nts)
FOR i = 0, nts - 1 DO BEGIN
   FOR layer = 0, nlay - 1 DO BEGIN
      zlim_nut(layer,i) = MIN( [ lim_dsi(layer,i), lim_no3(layer,i), lim_po4(layer,i) ] ) 
   ENDFOR
ENDFOR

i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 1.0001
cmin = 0 & cmax = 250
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3 ; colortable
title = ' Nutrient limitation '

; contours
icecontour, zlim_nut, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f3.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;----------------------
; Algal organic carbon
;----------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 8
fmin = 0.0 & fmax = MAX(AoCb)
fmin = 0.0 & fmax = 1000.
cmin = 255. & cmax = 80.
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 8 ; colortable
title = ' AoCb (mgC/m3)'

; contours
icecontour, AoCb, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f5.0)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.0155 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;-------------------------
; Detrital organic carbon
;-------------------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 8
fmin = 0.0 & fmax = MAX(eOCb)
fmin = 0.0 & fmax = 1000.
cmin = 255. & cmax = 80.
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 8 ; colortable
title = ' eoCb (mgC/m3)'

; contours
icecontour, eOCb, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f5.0)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.0155 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format


;
;-----------------
; chla
;-----------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1
; contour parameters
nlevels = 10
fmin = 0.0 & fmax = MAX(chla)
fmin = 0.0 & fmax = 80.
cmin = 255. & cmax = 80.
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

IF ( site EQ 'ISPOL' ) THEN BEGIN
;   zmax = 0.
;    levels = [   0.,   5., 10.,  15.,  20.,  25.,  30., 35., 40., 45., 50.] ; CHL_LONG_001
;   levels = [   0.,  2.5,  5.,  7.5,  10., 12.5,  15.,17.5, 20.,22.5, 25.] ; CTL_LONG_01
;   levels = [   0.,  35., 70., 105., 140., 175., 210., 245., 280., 315., 350.] ; PCS_LON_008
;   levels = [ 0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100. ] ; PCS_LONG_009
;   levels = [ 0., 30., 60., 90., 120., 150., 180., 210., 240., 270., 300. ]
ENDIF

ct = 8 ; colortable
title = ' chla (mg/m3) '

; contours
icecontour, chla, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f5.0)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;-----------------
; Synthesis
;-----------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1
; contour parameters
nlevels = 8 
fmin = 0.00 & fmax = MAX(syn) ; LONG(MAX(syn)) + 1
cmin = 255 & cmax = 80
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 8 ; colortable
title = ' syn (mgC/m3/d) '

; contours
icecontour, syn, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f5.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;-----------------
; Lysis
;-----------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1
; contour parameters
nlevels = 8
fmin = 0.0 & fmax = MAX(lys) ; LONG(MAX(lys)) + 1
cmin = 255 & cmax = 80
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors
;IF ( zmax LE 0.0 ) THEN levels = [ 0, 1, 2, 3 ,4,5,6]

ct = 8 ; colortable
title = ' lys (mgC/m3/d) '

; contours
icecontour, lys, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title
; colorbar
format='(f5.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;--------------
; Argon bulk
;--------------
;
!P.MULTI=[0,numplot_x, numplot_y]
i_plot = 0
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 10. ;MAX(Argb)
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' Arg (mmol/m3)'

icecontour, Argb , h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I2)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;---------------
; Argon bubbles
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 10. ; MAX(Argbub)
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' Arg bubbles (mmol/m3) '

icecontour, Argbub, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I2)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;--------------
; Oxygen bulk
;--------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 200. ;MAX(Argb)
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' Oxy (mmol/m3)'

icecontour, Oxyb , h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I3)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;----------------
; Oxygen bubbles
;----------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 200.; MAX(Argbub)
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' Oxy bubbles (mmol/m3) '

icecontour, Oxybub, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I2)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format


;
;---------------
; DIC
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 600. ; MAX(DICb)
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' DIC (mmol/m3) '

icecontour, DICb, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I4)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;---------------
; CO2aq
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 20. ; MAX(CO2aq)
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' CO2aq (mmol/m3) '

icecontour, CO2aq, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I4)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;---------------
; HCO3m
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 600.
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' HCO3m (mmol/m3) '

icecontour, HCO3m, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I4)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;---------------
; CO32m
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax =  40.
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' CO32m (mmol/m3) '

icecontour, CO32m, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I4)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;---------------
; pCO2 
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 5
fmin = 0.0 & fmax =12000.
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' pCO2 in brine (muatm) '

icecontour, pCO2, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(E7.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;---------------
; pH   
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 7.5 & fmax = 9.5
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' pH '

icecontour, pH, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(F3.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format


;
;---------------
; Alk
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax = 600. ; MAX(Alkb)
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' Alk (mmol/m3) '

icecontour, Alkb, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I4)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;---------------
; Ika
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin =  0.0 & fmax = 120. ; MAX(Alkb)
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' Ika (mmol/m3) '

icecontour, Ikab, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(I4)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format


; DEVICE, /CLOSE
; SET_PLOT, "X"
; !P.MULTI=[0,2,2]
; STOP

;
;---------------
; CO2bub
;---------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
fmin = 0.0 & fmax =  60. ;  MAX(CO2bub) ;1.0
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 1 ; colortable
title = ' CO2bub (mmol/m3) '

icecontour, CO2bub, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(F3.1)'
format='(I2)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;--------------
; dFe         
;--------------
;
!P.MULTI=[0,numplot_x, numplot_y]
i_plot = 0
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
nlevels = 10
fmin = 0.0 & fmax = MAX(dFeb) ;
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3 ; colortable
title = ' dFe (mumol/m3)'

icecontour, dFeb, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(f5.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;--------------
; aFe         
;--------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
nlevels = 10
fmin = 0.0 & fmax = MAX(aFeb) ;
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3 ; colortable
title = ' aFe (mumol/m3)'

icecontour, aFeb, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(f5.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;--------------
; eFe         
;--------------
;
i_plot = i_plot + 1
i_plot_h = i_plot MOD numplot_x
i_plot_v = i_plot / numplot_x + 1

; contour parameters
nlevels = 10
nlevels = 10
fmin = 0.0 & fmax = MAX(eFeb) ;
cmin = 000 & cmax = 255
do_levels, fmin, fmax, cmin, cmax, nlevels, levels, colors

ct = 3 ; colortable
title = ' eFe (mumol/m3)'

icecontour, eFeb, h_i, h_s, nlevels, colors, levels, ct, z_ib, mon, nlay_bio, nts, $
            tmin, tmax, zmin, zmax, title

; colorbar
format='(f5.1)'
horpos   = 0.02 + i_plot_h * dh
vertpos  = 0.015 + ( numplot_y - i_plot_v ) * dv ; position of colorbar.
colorbarn, horpos, vertpos, cb_length, 0.02, MIN(levels), MAX(levels), leveling, nlevels+1, levels, colors, format

;
;------------------------------------------------------------------------------
; End of the script
;------------------------------------------------------------------------------
;
IF ( device EQ 'PS' ) THEN BEGIN
   DEVICE, /CLOSE
   SET_PLOT, "X"
   !P.MULTI=[0,2,2]
ENDIF
END