source: branches/2017/dev_v3.20_2017_transport_max/GRAPHICS/YROSIAE_profiles.pro

Last change on this file was 6, checked in by vancop, 8 years ago

initial import of v3.20 /Users/ioulianikolskaia/Boulot/CODES/LIM1D/ARCHIVE/TMP/LIM1D_v3.20/

File size: 37.7 KB
Line 
1;
2;------------------------------------------------------------------------------
3;
4; This script makes plots of biological simulations
5;
6; (c) Martin Vancoppenolle, UCL-ASTR, june 2007
7; reviewed may 2008
8; reviewed feb 2010, for biological model
9; reviewed aug 2014, helsinki
10; reviewed nov 2015, helsinki
11;
12;--------------------
13;--- YROSIAE site ---
14;--------------------
15;
16;------------------------------------------------------------------------------
17;
18miss_val = -9999.
19i_auto = 1             ; 1=use exp_id.dat; 2=use manual mode with sensitivity runs
20plot_type = 'DEFAULT'  ;
21;
22;==============================================================================
23; General Options
24;==============================================================================
25;
26
27IF ( i_auto EQ 1 ) THEN BEGIN ; automatic mode
28
29   nruns = 1
30   control_run = 'E001' & indir = '' & outdir = '' & dt = 86400 & c_bio_model = 'BFMSI' & site = 'TOURNAI'
31   
32   OPENR, 11, 'exp_id.dat'
33   READF, 11, control_run & READF, 11, indir & READF, 11, outdir
34   READF, 11, dt & READF, 11, c_bio_model & READF, 11, site
35   CLOSE, 11
36   
37   PRINT, control_run
38   PRINT, indir
39   PRINT, outdir
40   PRINT, dt
41   PRINT, c_bio_model
42   PRINT, site
43
44ENDIF
45
46;-------------
47; Manual mode
48;-------------
49IF ( i_auto EQ 0 ) THEN BEGIN
50
51   indir_base = '/Users/ioulianikolskaia/Boulot/CODES/BFMLIM_2013/BFM-LIM-2506_2014/LIM/RUN/'
52   outdir = '/Users/ioulianikolskaia/Boulot/SCIENCE/PLOT_SCRIPTS/LIM1D_BIO/IDL/plots/'
53   dt   = 3600.        ; time step
54
55   control_run = 'LIM_KR_ML'
56   sensiti_run = [ ]
57;  sensiti_run =  [ 'LIM_KR_SL', 'LIM_KR_BA' ]
58   dummy = SIZE(sensiti_run)
59   nruns = dummy(1) + 1 ; sensitivity runs + 1 control run
60   print, 'nruns ', nruns
61
62   indir = STRARR(nruns)
63   indir(0) = indir_base+control_run+'/'
64   FOR i = 1, nruns-1 DO BEGIN
65      indir(i) = indir_base+sensiti_run(i-1)+'/'
66   ENDFOR
67
68ENDIF
69
70;
71;==============================================================================
72; Graphic Options
73;==============================================================================
74;
75
76;--------------------------------
77;--- plot_type dependent options
78;--------------------------------
79
80IF ( plot_type EQ 'DEFAULT' ) THEN BEGIN
81
82   colors      = [  0, 200 ,   100   ,   240   ,   240   ,   120   ,   120    ]     
83   linestyle   = [  0,0    ,    0    ,    1    ,     2   ,     2   ,     1    ]
84   thick       = [  5, 5.  ,    5.0  ,    2.0  ,   2.0   ,   2.0   ,    2.0   ]
85   stamp_out   = control_run+'_profiles'
86   cs          = 1.8 ; charsize
87   colorkey    = 'rd'
88   ct          = 0
89   !X.MARGIN   = [8,3]
90   !Y.MARGIN   = [3,3]
91   x_size      = 4.
92   y_size      = 6.
93   add_obs     = 'YES'
94
95ENDIF
96
97file_name = stamp_out+'.ps'
98
99numplot_x = 5       ; number of horizontal plots
100numplot_y = 2       ; number of vertical plots
101
102device = 'PS'       ; 'PS' or 'X'
103
104;--------------------------------
105;--- define output plots
106;--------------------------------
107figuresize_x = numplot_x * x_size ; figure size on x direction
108figuresize_y = numplot_y * y_size ; figure size on y direction
109set_plot, device
110IF ( device EQ 'PS' ) THEN BEGIN
111   device, /COLOR, /LANDSCAPE, filename=outdir+file_name, $
112           XSIZE=figuresize_x,YSIZE=figuresize_y,FONT_SIZE=9.0
113   loadct, ct
114ENDIF
115
116IF ( device EQ 'X' ) THEN BEGIN
117   xsize = 1200
118   ysize = 800
119   init_graphics_x, xsize, ysize, colorkey
120ENDIF
121
122;
123;==============================================================================
124; EXTRACT DATA
125;==============================================================================
126;
127;open_fields_v6, indir(0), control_run, c_bio_model, $
128;                numt, doy, ts_d, ts_m,                        $
129;                h_i, h_s, z_ip, z_ib, s_i, t_i, e_i, PAR, Ra, $
130;                dhib, dhisu, dhisi,                           $
131;                DAFb, NO3b, PO4b, DSib, chla,                 $
132;                DAFt, NO3t, PO4t, DSit,                       $
133;                FDAbd, FDAbpos, FDAsi,                        $
134;                FNO3, FNO3bpos, FNO3si,                       $
135;                syn, lys, rem,                                $
136;                lim_lig, lim_no3, lim_po4, lim_dsi, lim_tem, lim_sal
137               
138 open_fields_v7, indir(0), control_run, c_bio_model, $
139                 numt, doy, ts_d, ts_m,                        $
140                 h_i, h_s, z_ip, z_ib, s_i, t_i, e_i, PAR, Ra, $
141                 dhib, dhisu, dhisi,                           $
142                 DAFb, NO3b, PO4b, DSib, chla, eoCb,           $
143                 Argb, Argbub, $
144                 Oxyb, Oxybub, $
145                 DICb, Alkb, CO2b, CO2bub,                     $
146                 CO2aq, HCO3m, CO32m, pH, pCO2, Ikab,          $
147                 dFeb, aFeb, eFeb,                             $
148                 DAFt, NO3t, PO4t, DSit, eoCt, Argt,           $
149                 DICt, Alkt, CO2t, Ikat,                       $
150                 dFet, aFet, eFet,                             $
151                 FDAbd, FDAbpos, FDAsi,                        $
152                 FNO3, FNO3bpos, FNO3si,                       $
153                 FCO2_atm, FCO2_bub,                           $
154                 syn, lys, rem,                                $
155                 lim_lig, lim_no3, lim_po4, lim_dsi, lim_tem, lim_sal
156
157zsize = SIZE(DAFb)
158nlay = zsize(1)   ; number of layers
159nts  = zsize(2)   ; number of time steps
160
161;
162;------------------------------------------------------------------------------
163; OBS ...
164;------------------------------------------------------------------------------
165;
166
167; observed normalized profiles from Lannuzel et al. (2007)
168
169N_obs = 9 ; station numbers
170N_dpt = 50 ; number of depths
171
172; DOY
173doy_obs =[ 263, 279, 292, 299, 306, 312, 318, 327, 355 ]
174
175; Ice thickness / snow depth
176hi_obs =     [ 1.472, 1.609, 1.670, 1.690, 1.708, 1.714, 1.711, 1.677, 1.693 ]
177hi_obs_std = [ 0.023, 0.014, 0.009, 0.011, 0.011, 0.009, 0.016, 0.024, 0.014 ]
178hs_obs = REPLICATE(0.01, 9)
179
180h_chl     = [ 1.52, 1.60, 1.69, 1.715, 1.75, 1.72, 1.73, 1.70, 1.70  ]
181h_nut     = [ 1.51, 1.60, 1.69, 1.70 , 1.72, 1.735,1.705, 1.70, 1.70 ]
182
183; Depths
184N_chl     = LONARR(N_obs)
185N_chl     =  [ 15, 16, 17, 18, 18, 18, 18, 17, 17 ]
186
187depth_chl = FLTARR(N_obs, N_dpt)
188
189depth_chl(*,*) = miss_val
190depth_chl[0,0:N_chl(0)-1] = [ 5, 15, 25, 35, 45, 55, 65, 76, 87, 97   , 107 , 117  , 127  , 137  , 147 ] / 100.
191depth_chl[1,0:N_chl(1)-1] = [ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95   , 105 , 115  , 125  , 135  , 145  , 155 ] / 100.
192depth_chl[2,0:N_chl(2)-1] = [ 5, 15, 25, 35, 45, 55, 65, 75, 85, 94.5 , 104 , 114  , 124  , 134  , 144  , 154  , 164 ] / 100.
193depth_chl[3,0:N_chl(3)-1] = [ 5, 15, 25, 35, 45, 55, 65, 75, 85, 90.75, 96.5, 106.5, 116.5, 126.5, 136.5, 146.5, 156.5, 166.5 ] / 100.
194depth_chl[4,0:N_chl(4)-1] = [ 5, 15, 25, 35, 45, 55, 65, 75, 85, 92.5 , 100 , 110,   120  , 130  , 140  , 150  , 160  , 170 ] / 100.
195depth_chl[5,0:N_chl(5)-1] = [ 5, 15, 25, 35, 45, 55, 65, 75, 85, 91   , 97  , 107,   117  , 127  , 137  , 147  , 157  , 167 ] / 100.
196depth_chl[6,0:N_chl(6)-1] = [ 5, 15, 25, 35, 45, 55, 65, 65, 75, 85   , 95  , 106.5, 118  , 128  , 138  , 148  , 158  , 168 ] / 100.
197depth_chl[7,0:N_chl(7)-1] = [ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95   , 105.2, 115.2,125.2, 135.2, 145.2, 155.2, 165.2 ] / 100.
198depth_chl[8,0:N_chl(8)-1] = [ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95   , 105  , 115  , 125 , 135  , 145  , 155  , 165   ] / 100.
199
200; Chlorophyll
201rhoi = 917.
202rhow = 1020.
203chla_obs = FLTARR(N_obs, N_dpt)
204chla_obs(*,*) = miss_val
205chla_obs[0,0:N_chl(0)-1] = [ 0.14, 0.36, 0.12, 0.03, 0.14, 0.02, 0.01, 0.03, 0.03, 0.02, 0.19, 1.02, 1.03, 0.45, 13.28]
206chla_obs[1,0:N_chl(1)-1] = [ 0.24, 0.04, 0.06, 0.03, 0.15, 0.06, 0.13, 0.16, 0.04, 0.07, 0.26, 0.62, 0.39, 0.49, 2.85, 45.33]
207chla_obs[2,0:N_chl(2)-1] = [ 0.10, 0.08, 0.32, 0.03, 0.19, 0.06, 0.03, 0.03, 0.04, 0.06, 0.26, 0.26, 0.25, 0.27, 1.25, 2.30, 505.99 ]
208chla_obs[3,0:N_chl(3)-1] = [ 0.33, 0.22, 0.97, 0.00, 0.18, 0.08, 0.06, 0.05, 0.01, 0.42, 0.20, 0.32, 0.30, 0.25, 0.18, 0.26, 1.63, 1440.73 ]
209chla_obs[4,0:N_chl(4)-1] = [ 0.11, 0.12, 0.07, 0.06, 0.12, 0.06, 0.06, 0.04, 0.22, 0.13, 0.49, 0.28, 0.24, 0.36, 0.92, 1.11, 1.98, 501.79 ]
210chla_obs[5,0:N_chl(5)-1] = [ 0.75, 0.20, 0.13, 0.05, 0.03, 0.07, 0.03, 0.05, 0.03, 0.06, 0.08, 0.15, 0.10, 0.05, 0.07, 0.11, 0.81, 2443.38 ]
211chla_obs[6,0:N_chl(6)-1] = [ 0.18, 0.16, 0.06, 0.10, 0.08, 0.04, 0.06, 1.31, 0.03, 0.05, 0.25, 0.16, 0.06, 0.14, 0.05, 0.10, 1.39, 2342.79 ]
212chla_obs[7,0:N_chl(7)-1] = [ 0.35, 0.88, 0.86, 0.53, 0.84, 0.44, 0.96, 0.44, 0.30, 0.34, 1.07, 1.69, 6.05, 0.13, 1.68, 1.18, 1203.62       ]
213chla_obs[8,0:N_chl(8)-1] = [ 0.59, 0.29, 0.33, 0.28, 0.16, 0.17, 0.16, 0.11, 0.13, 0.27, 0.19, 0.22, 0.21, 0.63, 0.89, 2.35, 974.83        ]
214
215zaddr = WHERE(chla_obs NE miss_val)
216chla_obs(zaddr) = chla_obs(zaddr)  * rhoi / rhow
217
218; depths for nuts
219N_nut = REPLICATE(6,9)
220depth_nut = FLTARR(N_obs, 6)
221depth_nut[0, 0:N_nut(0)-1] = [ 12.5, 37.5, 65.05, 103.5, 128.5, 146. ] / 100.
222depth_nut[1, 0:N_nut(1)-1] = [   15, 45  , 75   , 105  , 135  , 155  ] / 100.
223depth_nut[2, 0:N_nut(2)-1] = [   15, 47  , 79.5 , 111  , 143  , 164  ] / 100.
224depth_nut[3, 0:N_nut(3)-1] = [   16, 48  , 80   , 112  , 144  , 165  ] / 100.
225depth_nut[4, 0:N_nut(4)-1] = [   16, 48  , 80.5 , 113.3, 145.8, 167  ] / 100.
226depth_nut[5, 0:N_nut(5)-1] = [16.35, 49.05, 71.75, 114.45, 147.15, 168.5 ] / 100.
227depth_nut[6, 0:N_nut(6)-1] = [16.05, 48.15, 80.25, 112.35, 144.45, 165.5 ] / 100.
228depth_nut[7, 0:N_nut(7)-1] = [   16, 48, 80, 112, 144, 165 ] / 100.
229depth_nut[8, 0:N_nut(8)-1] = [ 17, 49.75, 85, 112.75, 144.25, 165 ] / 100.
230
231; Silicates
232dsi_obs = FLTARR(N_obs, 6)
233dsi_obs[0,0:N_nut(0)-1] = [ 20.10, 18.03, miss_val, 12.78, 12.45, 29.98 ]
234dsi_obs[1,0:N_nut(1)-1] = [ 19.81, 15.01,    14.31, 13.83,  9.68, 16.23 ]
235dsi_obs[2,0:N_nut(2)-1] = [ 17.11, 14.42,    12.32,  6.07,  9.14, 28.56 ]
236dsi_obs[3,0:N_nut(3)-1] = [ 17.23, 14.81,    11.90,  9.92,  4.37, 25.08 ]
237dsi_obs[4,0:N_nut(4)-1] = [ 17.07, 14.81,    12.17, 11.07,  3.49, 23.19 ]
238dsi_obs[5,0:N_nut(5)-1] = [ 17.89, 14.10,    12.72, 10.03,  3.10, 24.25 ]
239dsi_obs[6,0:N_nut(6)-1] = [ 19.54, 14.92,    12.94, 10.36,  2.83, 33.80 ]
240dsi_obs[7,0:N_nut(7)-1] = [ 20.80, 13.99,    10.85,  6.29,  1.29, 17.37 ]
241dsi_obs[8,0:N_nut(8)-1] = [ 16.30, 12.89,     9.86,  2.72,  1.18, 9.83  ]
242
243zaddr = WHERE(dsi_obs NE miss_val)
244dsi_obs(zaddr) = dsi_obs(zaddr)  * rhoi / rhow
245
246; NOX
247nox_obs = FLTARR(N_obs, 6)
248nox_obs[0,0:N_nut(0)-1] = [ 10.00, 6.17, 5.80, 5.98, 3.25, 16.10 ]
249nox_obs[1,0:N_nut(1)-1] = [ 7.14, 3.96, 4.59, 2.28, 0.25, 7.45 ]
250nox_obs[2,0:N_nut(2)-1] = [ 8.64, 6.09, 5.04, 1.05, 1.50, 35.72 ]
251nox_obs[3,0:N_nut(3)-1] = [ 8.13, 5.72, 4.67, 0.78, 0.42, 49.64 ]
252nox_obs[4,0:N_nut(4)-1] = [ 8.15, 7.07, 5.23, 0.83, 0.28, 58.91 ]
253nox_obs[5,0:N_nut(5)-1] = [ 5.80, 5.02, 8.02, 1.19, 0.61, 97.09 ]
254nox_obs[6,0:N_nut(6)-1] = [ 8.68, 5.90, 4.79, 0.84, 0.33, 93.78 ]
255nox_obs[7,0:N_nut(7)-1] = [ 7.39, 4.61, 1.40, 0.50, 0.35, 7.94  ]
256nox_obs[8,0:N_nut(8)-1] = [ 2.74, 1.29, 0.44, 0.21, 0.08, 8.54  ]
257
258zaddr = WHERE(nox_obs NE miss_val)
259nox_obs(zaddr) = nox_obs(zaddr)  * rhoi / rhow
260
261; PO4
262po4_obs = FLTARR(N_obs, 6)
263po4_obs[0,0:N_nut(0)-1] = [ 0.49, 0.38, 0.31, 0.31, 0.29, 0.72 ]
264po4_obs[1,0:N_nut(1)-1] = [ 0.36, 0.24, 0.27, 0.20, 0.15, 1.80 ]
265po4_obs[2,0:N_nut(2)-1] = [ 0.42, 0.31, 0.28, 0.10, 0.23, 8.47 ]
266po4_obs[3,0:N_nut(3)-1] = [ 0.41, 0.28, 0.24, 0.22, 0.14, 16.54 ]
267po4_obs[4,0:N_nut(4)-1] = [ 0.44, 0.33, 0.25, 0.18, 0.11, 22.05 ]
268po4_obs[5,0:N_nut(5)-1] = [ 0.38, 0.33, 0.47, 0.27, 0.16, 37.34 ]
269po4_obs[6,0:N_nut(6)-1] = [ 0.50, 0.34, 0.32, 0.32, 0.19, 35.61 ]
270po4_obs[7,0:N_nut(7)-1] = [ 0.57, 0.38, 0.29, 0.22, 0.16, 9.79  ]
271po4_obs[8,0:N_nut(8)-1] = [ 0.37, 0.21, 0.15, 0.11, 0.11, 10.96 ]
272
273zaddr = WHERE(po4_obs NE miss_val)
274po4_obs(zaddr) = po4_obs(zaddr)  * rhoi / rhow
275
276; physical observations
277
278N_TS = [ 31, 33, 33, 34, 34, 35, 35, 34, 35 ]
279; depth of TS cores
280depth_TS = FLTARR(N_obs, N_dpt)
281depth_TS(*,*) = miss_val
282
283depth_TS[0,0:N_TS(0)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 102.5, 107.5, 112.5, 117.5, 122.5, 127.5, 132.5, 137.5, 142.5, 147.5, 152.5 ] / 100.
284depth_TS[1,0:N_TS(1)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 100  , 105  , 110  , 115  ,   120,   125,   130,   135,   140,   145,   150,   155, 160 ] / 100.
285depth_TS[2,0:N_TS(2)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 99.5, 104.5, 109.5, 114.5, 119.5, 124.5, 129.5, 134.5, 139.9, 144.5, 149.5, 154.5, 159.5, 164.5 ] / 100.
286depth_TS[3,0:N_TS(3)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 102.5, 107.5, 112.5, 117.5, 122.5, 127.5, 132.5, 137.5, 142.5, 147.5, 152.5, 157.5, 162.5, 167.5 ] / 100.
287depth_TS[4,0:N_TS(4)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 102.5, 107.5, 112.5, 117.5, 122.5, 127.5, 132.5, 137.5, 142.5, 147.5, 152.5, 157.5, 162.5, 167.5 ] / 100.
288depth_TS[5,0:N_TS(5)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 94.5, 99.5, 104.5, 109.5, 114.5, 119.5, 124.5, 129.5, 134.5, 139.5, 144.5, 149.5, 154.5, 159.5, 164.5, 169.5 ] / 100.
289depth_TS[6,0:N_TS(6)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 102.5, 107.5, 112.5, 117.5, 122.5, 127.5, 132.5, 137.5, 142.5, 147.5, 152.5, 157.5, 162.5, 167.5, 172.5 ] / 100.
290depth_TS[7,0:N_TS(7)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 102.5, 107.5, 112.5, 117.5, 122.5, 127.5, 132.5, 137.5, 142.5, 147.5, 152.5, 157.5, 162.5, 167.5 ] / 100.
291depth_TS[8,0:N_TS(8)-1] = [ 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 102.5, 107.5, 112.5, 117.5, 122.5, 127.5, 132.5, 137.5, 142.5, 147.5, 152.5, 157.5, 162.5, 167.5, 172.5 ] / 100.
292
293h_TS = FLTARR(N_obs) ; We don't have thicknes of TS cores
294FOR i_obs = 0, N_obs - 1 DO BEGIN
295  zdh         = ( depth_TS(i_obs,N_TS(i_obs)-1) - depth_TS(i_obs,N_TS(i_obs)-2) ) / 2.
296  h_TS(i_obs) = depth_TS(i_obs,N_TS(i_obs)-1) + zdh
297ENDFOR
298
299; T-S
300t_i_obs = FLTARR(N_obs, N_dpt)
301t_i_obs(*,*) = miss_val
302
303t_i_obs[0,0:N_TS(0)-1] = [ -19.8, -18.7, -18.6, -18.2, -17.5, -16.8, -16.4, -16.2, -15.3, -14.7, -14.2, -13.5, -13.1, -12.6, -11.9, -11.5, -12.4, -11.5, -10.8, -10, -9.1, -8.6, -8.1, -7, -5.7, -4.7, -4.4, -3.7, -2.9, -2.4, -2.2 ]
304t_i_obs[1,0:N_TS(1)-1] = [ -12.7, -12.6, -12.8, -12.5, -12.6, -12.2, -12.3, -12.4, -11.9, -11.8, -11.6, -11.1, -10.7, -10.4, -10, -10, -9.9, -9.5, -9.4, -9.3, -8.5, -7.9, -7.3, -6.6, -6.1, -5.5, -5, -4.5, -4.2, -3.6, -2.8, -2.3, -1.9 ]
305t_i_obs[2,0:N_TS(2)-1] = [ -14.6, -14.6, -14.5, -14.1, -13.5, miss_val, -12.6, -11.6, -11.8, -11.5, -11.1, -10.5, -10.2, -9.9, -9.2, -8.9, -8.7, -8.8, -9.1, -6.1, -5.6, -5.2, -4.8, -4.3, -4.2, -4.1, -3.7, -3.6, -3.1, -2.8, -2.4, -1.9, -2 ]
306t_i_obs[3,0:N_TS(3)-1] = [ -12.6, -12.3, -12.2, -11.6, -11.4, -10.9, -10.1, -9.9, -9.6, -9, -8.8, -8.4, -8, -8, -7.9, -7.3, -6.9, -6.4, miss_val, miss_val, -6.2, -5.6, -5, -4.8, -4.5, -3.7, -3.7, -3.2, -2.7, -2.7, -2.6, -1.7, -2.1, -2 ]
307t_i_obs[4,0:N_TS(4)-1] = [ -12.4, -12.9, -12.5, -12.2, -11.8, -11.5, -11.1, -10.5, -10.3, -9.7, -9.2, miss_val, -8.8, -8.2, -8, -7.4, -7.3, -6.7, -6.3, miss_val, -5, -4.6, -4.6, -3.9, -3.6, -3.4, -3.1, -3, -3, -2.5, -2.7, -2.2, -2, -2 ]
308t_i_obs[5,0:N_TS(5)-1] = [ -10.2, -10.1, -10.2, -9.8, -9.4, -9.3, -8.9, miss_val, -8.7, -8.2, -7.8, -7.8, -7.6, -7.2, -6.6, -6.8, -6.4, -6.2, -6.2, miss_val, miss_val, -4.9, -4.6, -4.3, -4, -3.7, miss_val, -3.4, -3.2, -2.8, -2.6, -1.9, -2, -1.8, -2 ]
309t_i_obs[6,0:N_TS(6)-1] = [ -7.6, -8.3, -8.2, -7.7, -7.6, miss_val, -7.1, -6.4, -7, -6.8, -6.8, -6.8, -6.9, -6.8, -6.3, -5.9, miss_val, -5.2, -5, -5, -4.8, -4.5, -4.2, -4, miss_val, -3.7, -3.4, -2.7, -2.7, -2.5, -2, -2, -1.6, -1.6, -1.9 ]
310t_i_obs[7,0:N_TS(7)-1] = [ -4.5, -4.6, -4.7, -4.5, -4.6, miss_val, -4, -3.9, -3.4, -3.4, -3, -3.3, -3.2, -3, -2.7, -2.7, -2.9, -2.4, -2.4, -2.1, -2, -2.3, -1.8, -1.6, -1.7, -1.7, -1.4, -1.1, -1.4, -1.3, -1.4, -0.8, -2, -1.6 ]
311t_i_obs[8,0:N_TS(8)-1] = [ -2.1, -3.3, -3.6, -3.4, -3.1, -3, -3, -2.7, -2.6, miss_val, -2.9, -2.8, -2.2, -1.7, -1.6, -2.3, -1.6, -1.7, -1.7, -1.9, -1.8, -2, -1.8, -1.4, -1.9, -1.4, -1.7, -1.8, miss_val , -1.3, -1, -1, -1.6, miss_val, miss_val ]
312
313; S-physical core
314s_i_obs = FLTARR(N_obs, N_dpt)
315s_i_obs(*,*) = miss_val
316
317s_i_obs[0,0:N_TS(0)-1] = [ 13, 11.5, 11.3, 9.7, 10.1, 9.5, 8.5, 7.6, 7.1, 6.5, 6.8, 6.2, 6.9, 6.7, 7, 7.6, 6.6, 8.2, 5.6, 4.5, 4.8, 4.8, miss_val, 5.1, 5, 4.8, 4.8, 4.7, 4.7, 5.5, 9.2 ]
318s_i_obs[1,0:N_TS(1)-1] = [ 11.4, 10.1, 10.5, miss_val , 10.5, 8.8, 8, 7.1, 6.5, 6.8, 7.3, 7.2, 6.7, 6.1, 5.9, 5.9, 6, 6.2, 5.7, 5.1, 7.3, 5.9, 6, 6, 5, 4.8, 4.8, 5.7, 5.2, 3.7, 6.2, 7.1, 11 ]
319s_i_obs[2,0:N_TS(2)-1] = [ 11.2, 10.5, 10.1, 9.5, 9.6, 7.5, 7.6, 7.4, 6.7, 6.5, 6.5, 6.5, 6.3, 5.8, 5.7, 5.6, 5.9, 5.3, 5, 5.4, 5.3, 4.9, 4.7, 4.2, 4.1, 4.4, 4.8, 4.4, 4.8, 5.3, 5.5, miss_val, 12.5 ]
320s_i_obs[3,0:N_TS(3)-1] = [ 10.8, 10.9, 10.4, 9.4, 9.7, 8.2, 7.6, 7.1, 6.8, 6.7, 7, 7.2, 6.9, 6, 6.7, 5.9, 5.7, 5.3, 4.9, 5.4, 4.9, 5.2, 5, 4.9, 4.3, 4.6, 4.5, 4.7, 4.7, 5, 5.7, 5.4, 5.8, 13.7 ]
321s_i_obs[4,0:N_TS(4)-1] = [ 10.9, 10.8, 10.1, 9, 9.4, 8.2, 8.1, 7.2, 7.3, 6.5, 6.8, 6.7, 6.7, 6.3, miss_val, 6.5, 6.6, 6.9, 6.6, miss_val, 6.1, 5.5, 5.3, 5.2, 6, 5.7, 5.7, 5.5, 5.1, 5, 5.5, 4.8, 5.5, 12.2 ]
322s_i_obs[5,0:N_TS(5)-1] = [ 10.7, 10.6, 10.3, 9.6, 10.1, 8, 7.7, 7, 6.4, 6.7, 6.8, 7, 6.8, 5.4, 5.7, 5.9, 5.7, 5.4, 5, 4.8, 5.1, 5.1, 5, 4.7, 4.2, 4.1, 4.4, 4.3, 4.2, 4.9, 5.6, 5.7, 5, 5.6, 12 ]
323s_i_obs[6,0:N_TS(6)-1] = [ 11.1, 10.1, 10.4, 10.8, 10.3, 9.4, 8.2, 7, 7, 6.9, 7.3, 7.3, 6.7, 6.2, 6.2, 6.2, 6.5, 6.2, 5.5, 5.1, 5.3, 5.8, 5.7, 5.7, 5.8, 5, 4.7, 4.7, 4.7, 4.7, 5.2, 5.4, 4.8, 5.8, 12.2 ]
324s_i_obs[7,0:N_TS(7)-1] = [ 8.4, 9.5, 9.4, 10.1, 8.8, 7.8, 7.3, 7.3, 6.2, 6.7, 6.5, 6.2, 5.7, 5.7, 5.4, 6.1, 5.2, 4.6, 4.2, 4.3, 4.5, 4.3, 4.2, 3.6, 3.7, 3.5, 3.5, 3.3, 3.9, 4.4, 4.3, 4.7, 10.1, miss_val ]
325s_i_obs[8,0:N_TS(8)-1] = [ 5.4, 6.7, 7.2, 7.6, 8.4, 7.2, 6.9, 6.5, 6.1, 6.2, 5.9, 6.1, 5.8, 4.9, 4.8, 4.7, 4.7, 4.4, 3.9, 3.4, 4.1, 3.7, 4.1, 3.8, 3.3, 3.2, 3.6, 3.1, 3.2, 3.4, 4.1, 5.2, 3.9, 4.1, 9 ]
326
327e_i_obs = FLTARR(N_obs, N_dpt)
328e_i_obs(*,*) = miss_val
329zaddr = where( ( s_i_obs NE miss_val ) AND ( t_i_obs NE miss_val ) )
330e_i_obs(zaddr) = - 0.054*s_i_obs(zaddr) / t_i_obs(zaddr)
331
332; brine concentrations
333e_i_nut = FLTARR(N_obs, 6)
334FOR i_obs = 0, N_obs - 1 DO BEGIN
335   zei   = e_i_obs(i_obs,*)
336   zd    = depth_TS(i_obs,*)
337   zaddr = WHERE(zei NE miss_val)
338   zei   = zei(zaddr)
339   zd    = zd(zaddr) / h_TS(i_obs)
340   zei_itp = INTERPOL(zei, zd, depth_nut(i_obs,*)/h_nut(i_obs)); interpolate bvf on nutrient depths
341   e_i_nut(i_obs,*) = zei_itp
342ENDFOR
343
344;DEVICE, /CLOSE
345;SET_PLOT, 'X'
346;!P.MULTI=[0,5,2]
347;FOR i_obs = 0, N_obs -1 DO BEGIN
348;   PLOT, [0., 100.], [-2., 0.], /NODATA, charsize = 2
349;   OPLOT, e_i_obs(i_obs,0:N_TS(i_obs)-1)*100., -depth_TS(i_obs,0:N_TS(i_obs)-1)/h_TS(i_obs), psym = 1
350;   OPLOT, e_i_nut(i_obs,*)*100., -depth_nut(i_obs,*)/h_nut(i_obs)
351;ENDFOR
352;STOP
353
354iaddr = where ( dsi_obs NE miss_val )
355dsi_br_obs = FLTARR(N_obs, 6)
356dsi_br_obs(iaddr) = dsi_obs(iaddr) / e_i_nut(iaddr)
357
358iaddr = where ( po4_obs NE miss_val )
359po4_br_obs = FLTARR(N_obs, 6)
360po4_br_obs(iaddr) = po4_obs(iaddr) / e_i_nut(iaddr)
361
362iaddr = where ( nox_obs NE miss_val )
363nox_br_obs = FLTARR(N_obs, 6)
364nox_br_obs(iaddr) = nox_obs(iaddr) / e_i_nut(iaddr)
365
366;*** Mean observed Profiles
367;--- interpolate T, S and e on standard depths ---
368N_std = 30
369depth_std = FINDGEN(30) / 30.+ 1./60.
370
371t_i_std = FLTARR(N_obs, N_std)
372s_i_std = FLTARR(N_obs, N_std)
373e_i_std = FLTARR(N_obs, N_std)
374chla_std = FLTARR(N_obs, N_std)
375FOR i_obs = 0, N_obs - 1 DO BEGIN
376   zt = t_i_obs(i_obs,0:N_TS(i_obs)-1)
377   zd = depth_TS(i_obs,0:N_TS(i_obs)-1)
378   zaddr = WHERE( zt NE miss_val )
379   zt = zt(zaddr)
380   zd = zd(zaddr) / h_TS(i_obs)
381   ztstd = INTERPOL( zt, zd, depth_std )
382   t_i_std(i_obs,*) = ztstd(*)
383
384   zs = s_i_obs(i_obs,0:N_TS(i_obs)-1)
385   zd = depth_TS(i_obs,0:N_TS(i_obs)-1)
386   zaddr = WHERE( zs NE miss_val )
387   zs = zs(zaddr)
388   zd = zd(zaddr) / h_TS(i_obs)
389   zsstd = INTERPOL( zs, zd, depth_std )
390   s_i_std(i_obs,*) = zsstd(*)
391
392   ze = e_i_obs(i_obs,0:N_TS(i_obs)-1)
393   zd = depth_TS(i_obs,0:N_TS(i_obs)-1)
394   zaddr = WHERE( ze NE miss_val )
395   ze = ze(zaddr)
396   zd = zd(zaddr) / h_TS(i_obs)
397   zestd = INTERPOL( ze, zd, depth_std )
398   e_i_std(i_obs,*) = zestd(*)
399
400   zchl    = chla_obs(i_obs,0:N_chl(i_obs)-1)
401   zd      = depth_chl(i_obs,0:N_chl(i_obs)-1)
402   zaddr   = WHERE( zchl NE miss_val )
403   zchl    = zchl(zaddr)
404   zd      = zd(zaddr) / h_chl(i_obs)
405   zchlstd = INTERPOL( zchl, zd, depth_std )
406   chla_std(i_obs,*) = zchlstd(*)
407
408;  s_i_std = ...
409;  e_i_std = ...
410ENDFOR
411
412;***
413
414t_i_obs_mean = MEAN( t_i_std, dimension = 1 )
415t_i_obs_std = STDDEV( t_i_std, dimension = 1 )
416
417s_i_obs_mean = MEAN( s_i_std, dimension = 1 )
418s_i_obs_std = STDDEV( s_i_std, dimension = 1 )
419
420e_i_obs_mean = MEAN( e_i_std, dimension = 1 )
421e_i_obs_std = STDDEV( e_i_std, dimension = 1 )
422
423;depth_obs_mean = MEAN( depth_obs, DIMENSION = 1 )
424
425chla_obs_mean  = MEAN( chla_std , DIMENSION = 1 )
426chla_obs_std = STDDEV( chla_std, dimension = 1 )
427
428dsi_obs_mean = FLTARR(N_dpt) & dsi_obs_std = FLTARR(N_dpt)
429nox_obs_mean = FLTARR(N_dpt) & nox_obs_std = FLTARR(N_dpt)
430po4_obs_mean = FLTARR(N_dpt) & po4_obs_std = FLTARR(N_dpt)
431depth_nut_mean = FLTARR(N_dpt)
432depth_nut_std = FLTARR(N_obs,6)
433FOR i_obs = 0, N_obs - 1 DO BEGIN
434   depth_nut_std(i_obs,*) = depth_nut(i_obs,*) / h_nut(i_obs)
435ENDFOR
436
437
438FOR i = 0, 5 DO BEGIN
439   i_addr = WHERE( dsi_obs(*,i) NE miss_val )
440   dsi_obs_mean(i) =  MEAN( dsi_obs(i_addr,i) )
441   dsi_obs_std(i)  = STDDEV( dsi_obs(i_addr,i) )
442   nox_obs_mean(i) =  MEAN( nox_obs(i_addr,i) )
443   nox_obs_std(i)  = STDDEV( nox_obs(i_addr,i) )
444   po4_obs_mean(i) =  MEAN( po4_obs(i_addr,i) )
445   po4_obs_std(i)  = STDDEV( po4_obs(i_addr,i) )
446   depth_nut_mean(i) = MEAN( depth_nut_std(i_addr,i) )
447ENDFOR
448
449;*** Corresponding model mean profiles
450t_i_mod_mean = FLTARR(nlay)  & t_i_mod_mean(*) = 0.
451t_i_mod_std  = FLTARR(nlay)  & t_i_mod_std (*) = 0.
452s_i_mod_mean = FLTARR(nlay)  & s_i_mod_mean(*) = 0.
453s_i_mod_std  = FLTARR(nlay)  & s_i_mod_std (*) = 0.
454e_i_mod_mean = FLTARR(nlay)  & e_i_mod_mean(*) = 0.
455e_i_mod_std  = FLTARR(nlay)  & e_i_mod_std (*) = 0.
456chla_mod_mean = FLTARR(nlay) & chla_mod_mean(*) = 0.
457chla_mod_std  = FLTARR(nlay) & chla_mod_std (*) = 0.
458nox_mod_mean = FLTARR(nlay)  & nox_mod_mean(*) = 0.
459nox_mod_std  = FLTARR(nlay)  & nox_mod_std (*) = 0.
460po4_mod_mean = FLTARR(nlay)  & po4_mod_mean(*) = 0.
461po4_mod_std  = FLTARR(nlay)  & po4_mod_std (*) = 0.
462dsi_mod_mean = FLTARR(nlay)  & dsi_mod_mean(*) = 0.
463dsi_mod_std  = FLTARR(nlay)  & dsi_mod_std (*) = 0.
464FOR i = 0, nlay - 1 DO BEGIN
465
466   FOR i_obs = 0, N_obs - 1 DO BEGIN
467      i_mod = WHERE(doy EQ doy_obs(i_obs)-1)
468      t_i_mod_mean(i)  = t_i_mod_mean(i) + t_i(i,i_mod) / FLOAT(N_obs)
469      s_i_mod_mean(i)  = s_i_mod_mean(i) + s_i(i,i_mod) / FLOAT(N_obs)
470      e_i_mod_mean(i)  = e_i_mod_mean(i) + e_i(i,i_mod) / FLOAT(N_obs)
471      chla_mod_mean(i) = chla_mod_mean(i) + chla(i,i_mod) / FLOAT(N_obs)
472      dsi_mod_mean(i)  = dsi_mod_mean(i) + dsib(i,i_mod) / FLOAT(N_obs)
473      nox_mod_mean(i)  = nox_mod_mean(i) + no3b(i,i_mod) / FLOAT(N_obs)
474      po4_mod_mean(i)  = po4_mod_mean(i) + po4b(i,i_mod) / FLOAT(N_obs)
475   ENDFOR
476
477   zsum   = 0
478   zsum_s = 0
479   zsum_e = 0
480   zsum_c = 0
481   zsum_d = 0
482   zsum_n = 0
483   zsum_p = 0
484   FOR i_obs = 0, N_obs - 1 DO BEGIN
485      i_mod = WHERE(doy EQ doy_obs(i_obs)-1)
486      zsum   = zsum + ( t_i(i,i_mod) - t_i_mod_mean(i) ) ^2 / FLOAT(N_obs)
487      zsum_s = zsum_s + ( s_i(i,i_mod) - s_i_mod_mean(i) ) ^2 / FLOAT(N_obs)
488      zsum_e = zsum_e + ( e_i(i,i_mod) - e_i_mod_mean(i) ) ^2 / FLOAT(N_obs)
489      zsum_c = zsum_c + ( chla(i,i_mod) - chla_mod_mean(i) ) ^2 / FLOAT(N_obs)
490      zsum_d = zsum_d + ( dsib(i,i_mod) - dsi_mod_mean(i) ) ^2 / FLOAT(N_obs)
491      zsum_n = zsum_n + ( no3b(i,i_mod) - nox_mod_mean(i) ) ^2 / FLOAT(N_obs)
492      zsum_p = zsum_p + ( po4b(i,i_mod) - po4_mod_mean(i) ) ^2 / FLOAT(N_obs)
493   ENDFOR
494   t_i_mod_std(i)  = SQRT( zsum   )
495   s_i_mod_std(i)  = SQRT( zsum_s )
496   e_i_mod_std(i)  = SQRT( zsum_e )
497   chla_mod_std(i) = SQRT( zsum_c )
498   nox_mod_std(i)  = SQRT( zsum_n )
499   po4_mod_std(i)  = SQRT( zsum_p )
500   dsi_mod_std(i)  = SQRT( zsum_d )
501
502ENDFOR
503
504;-------
505
506
507;
508;==============================================================================
509; DIAGNOSTICS
510;==============================================================================
511;
512phi=findgen(32)*(!PI*2/32.)
513phi = [ phi, phi(0) ]
514usersym, cos(phi), sin(phi), /fill
515
516;
517;==============================================================================
518; PLOTS
519;==============================================================================
520;
521;
522;==============================================================================
523; Time series of profiles
524;==============================================================================
525;
526;---------
527;--- T ---
528;---------
529!P.MULTI=[0,numplot_x, numplot_y]
530FOR i_obs = 0, N_obs - 1 DO BEGIN
531
532   ; prepare plot
533   LOADCT, 0
534   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
535   PLOT, [ -20, 0. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'deg C', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2
536   XYOUTS, -5, -0.1, ztitle
537
538   ; add obs
539   OPLOT, t_i_obs(i_obs,*), - depth_TS(i_obs,*)/h_TS(i_obs), PSYM = 1, THICK = 3, SYMSIZE = 0.8
540
541   ; add model
542   i_mod = WHERE(doy EQ doy_obs(i_obs))
543   zz = FLTARR(nlay)
544   FOR i = 0, nlay - 1 DO BEGIN
545      zz(i) = z_ip(i,i_mod) / h_i(i_mod)
546   ENDFOR
547   LOADCT, 3
548   OPLOT, t_i(*,i_mod), -zz, THICK = 3, COLOR = 150
549
550ENDFOR
551
552
553;---------
554;--- S   
555;---------
556!P.MULTI=[0,numplot_x, numplot_y]
557FOR i_obs = 0, N_obs - 1 DO BEGIN
558
559   ; prepare plot
560   LOADCT, 0
561   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
562   PLOT, [0., 15. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'g/kg', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2
563   XYOUTS, 7., -0.1, ztitle
564
565   ; add obs
566   OPLOT, s_i_obs(i_obs,*), - depth_TS(i_obs,*)/h_TS(i_obs), PSYM = 1, THICK = 3, SYMSIZE = 0.8
567
568   ; add model
569   i_mod = WHERE(doy EQ doy_obs(i_obs))
570   zz = FLTARR(nlay)
571   FOR i = 0, nlay - 1 DO BEGIN
572      zz(i) = z_ip(i,i_mod) / h_i(i_mod)
573   ENDFOR
574   LOADCT, 3
575   OPLOT, s_i(*,i_mod), -zz, THICK = 3, COLOR = 150
576
577ENDFOR
578
579
580;---------
581;--- e   
582;---------
583!P.MULTI=[0,numplot_x, numplot_y]
584FOR i_obs = 0, N_obs - 1 DO BEGIN
585
586   ; prepare plot
587   LOADCT, 0
588   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
589   PLOT, [0., 50. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = '%', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2
590   XYOUTS, 15., -0.1, ztitle
591
592   ; add obs
593   OPLOT, e_i_obs(i_obs,*)*100., - depth_TS(i_obs,*)/h_TS(i_obs), PSYM = 1, THICK = 3, SYMSIZE = 0.8
594
595   ; add model
596   i_mod = WHERE(doy EQ doy_obs(i_obs))
597   zz = FLTARR(nlay)
598   FOR i = 0, nlay - 1 DO BEGIN
599      zz(i) = z_ip(i,i_mod) / h_i(i_mod)
600   ENDFOR
601   LOADCT, 3
602   OPLOT, e_i(*,i_mod), -zz, THICK = 3, COLOR = 150
603
604ENDFOR
605
606
607;--------------
608;--- chla
609;--------------
610!P.MULTI=[0,numplot_x, numplot_y]
611FOR i_obs = 0, N_obs - 1 DO BEGIN
612
613   ; prepare plot
614   LOADCT, 0
615   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
616   PLOT, [ 0.01, 40. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mg chla / m3', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2
617   XYOUTS, 27., -0.1, ztitle
618
619   ; add model
620   i_mod = WHERE(doy EQ doy_obs(i_obs))
621   zz = FLTARR(nlay)
622   FOR i = 0, nlay - 1 DO BEGIN
623      zz(i) = z_ib(i,i_mod) / h_i(i_mod)
624   ENDFOR
625   LOADCT, 9
626   OPLOT, chla(*,i_mod), -zz, THICK = 3, COLOR = 150
627
628   ; add obs
629   zchlaobs = chla_obs(i_obs,0:N_chl(i_obs)-1)
630   zz_obs   = depth_chl(i_obs,0:N_chl(i_obs)-1) / h_chl(i_obs)
631   zchla_obs2mod = INTERPOL( zchlaobs, zz_obs,zz )
632   OPLOT, zchla_obs2mod, -zz, PSYM = 1, THICK = 3, SYMSIZE = 0.8
633;   print, zchla_obs2mod
634
635ENDFOR
636
637
638;------------
639;--- nox 
640;------------
641!P.MULTI=[0,numplot_x, numplot_y]
642
643FOR i_obs = 0, N_obs - 1 DO BEGIN
644
645   ; prepare plot
646   LOADCT, 0
647   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
648   PLOT, [ 0.00, 12. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mmmol N / m3', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2
649   XYOUTS, 2.2, -0.3, ztitle
650
651   ; add obs
652   OPLOT, nox_obs(i_obs,*), - depth_nut(i_obs,*)/h_nut(i_obs), PSYM = 1, THICK = 3, SYMSIZE = 0.8
653
654   ; add model
655   i_mod = WHERE(doy EQ doy_obs(i_obs))
656   zz = FLTARR(nlay)
657   FOR i = 0, nlay - 1 DO BEGIN
658      zz(i) = z_ib(i,i_mod) / h_i(i_mod)
659   ENDFOR
660   LOADCT, 9
661   OPLOT, no3b(*,i_mod), -zz, THICK = 3, COLOR = 150
662
663   PRINT, ' i_mod : ', i_mod
664   PRINT, ' NO3b  : ', no3b(*,i_mod)
665
666ENDFOR
667
668
669;-------------
670;--- po4 
671;-------------
672!P.MULTI=[0,numplot_x, numplot_y]
673
674FOR i_obs = 0, N_obs - 1 DO BEGIN
675
676   ; prepare plot
677   LOADCT, 0
678   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
679   PLOT, [ 0.00, 4.0 ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mmmol P / m3', XTICKFORMAT='(I2)', XMINOR = 2, YMINOR = 2
680   XYOUTS, 0.37, -0.3, ztitle
681
682   ; add obs
683   OPLOT, po4_obs(i_obs,*), - depth_nut(i_obs,*)/h_nut(i_obs) , PSYM = 1, THICK = 3, SYMSIZE = 0.8
684
685   ; add model
686   i_mod = WHERE(doy EQ doy_obs(i_obs))
687   zz = FLTARR(nlay)
688   FOR i = 0, nlay - 1 DO BEGIN
689      zz(i) = z_ib(i,i_mod) / h_i(i_mod)
690   ENDFOR
691   LOADCT, 9
692   OPLOT, po4b(*,i_mod), -zz, THICK = 3, COLOR = 150
693
694ENDFOR
695
696
697
698;--------------
699;--- DSi 
700;--------------
701!P.MULTI=[0,numplot_x, numplot_y]
702
703FOR i_obs = 0, N_obs - 1 DO BEGIN
704
705   ; prepare plot
706   LOADCT, 0
707   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
708   PLOT, [ 0.00, 40. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mmmol Si / m3', XTICKFORMAT='(F4.1)', XMINOR = 2, YMINOR = 2
709   XYOUTS, 3.0, -0.15, ztitle
710
711   ; add obs
712   OPLOT, dsi_obs(i_obs,*), - depth_nut(i_obs,*)/h_nut(i_obs), PSYM = 1, THICK = 3, SYMSIZE = 0.8
713
714   ; add model
715   i_mod = WHERE(doy EQ doy_obs(i_obs))
716   zz = FLTARR(nlay)
717   FOR i = 0, nlay - 1 DO BEGIN
718      zz(i) = z_ib(i,i_mod) / h_i(i_mod)
719   ENDFOR
720   LOADCT, 9
721   OPLOT, dsib(*,i_mod), -zz, THICK = 3, COLOR = 150
722
723ENDFOR
724
725;
726;==============================================================================
727; Time series of profiles, brine concentrations
728;==============================================================================
729;
730
731;--- nox_br
732!P.MULTI=[0,numplot_x, numplot_y]
733
734FOR i_obs = 0, N_obs - 1 DO BEGIN
735
736   ; prepare plot
737   LOADCT, 0
738   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
739   PLOT, [ 0., 300. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mmmol N / m3', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2
740   XYOUTS, 20.0, -0.15, ztitle
741
742   ; add obs
743   OPLOT, nox_br_obs(i_obs,*), - depth_nut(i_obs,*)/h_nut(i_obs), PSYM = 1, THICK = 3, SYMSIZE = 0.8
744
745   ; add model
746   i_mod = WHERE(doy EQ doy_obs(i_obs))
747   zz = FLTARR(nlay)
748   FOR i = 0, nlay - 1 DO BEGIN
749      zz(i) = z_ib(i,i_mod) / h_i(i_mod)
750   ENDFOR
751   LOADCT, 9
752   OPLOT, no3b(*,i_mod)/e_i(*,i_mod)*100., -zz, THICK = 3, COLOR = 150
753   LOADCT, 1
754
755   OPLOT, [ 1.6, 1.6 ], [ -1., 0.], linestyle = 1, color = 150 , thick = 3
756
757ENDFOR
758
759;--- po4_br
760!P.MULTI=[0,numplot_x, numplot_y]
761
762FOR i_obs = 0, N_obs - 1 DO BEGIN
763
764   ; prepare plot
765   LOADCT, 0
766   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
767   PLOT, [ 0.0, 15. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, $
768         SUBTITLE = 'mmmol P / m3', XTICKFORMAT='(I3)', XMINOR = 2, XTICKS = 3 , YMINOR = 2
769   XYOUTS, 10., -0.15, ztitle
770
771   ; add obs
772   OPLOT, po4_br_obs(i_obs,*), - depth_nut(i_obs,*)/h_nut(i_obs), PSYM = 1, THICK = 3, SYMSIZE = 0.8
773
774   ; add model
775   i_mod = WHERE(doy EQ doy_obs(i_obs))
776   zz = FLTARR(nlay)
777   FOR i = 0, nlay - 1 DO BEGIN
778      zz(i) = z_ib(i,i_mod) / h_i(i_mod)
779   ENDFOR
780   LOADCT, 9
781   OPLOT, po4b(*,i_mod)/e_i(*,i_mod)*100., -zz, THICK = 3, COLOR = 150
782   LOADCT, 1
783
784   OPLOT, [ 0.24, 0.24 ], [ -1., 0.], linestyle = 1, color = 150 , thick = 3
785
786ENDFOR
787
788;--- dsi_br
789!P.MULTI=[0,numplot_x, numplot_y]
790
791FOR i_obs = 0, N_obs - 1 DO BEGIN
792
793   ; prepare plot
794   LOADCT, 0
795   ztitle = STRCOMPRESS(STRING(doy_obs(i_obs)), /REMOVE_ALL)
796   PLOT, [ 0., 750. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mmmol Si / m3', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2
797   XYOUTS, 20.0, -0.15, ztitle
798
799   ; add obs
800   OPLOT, dsi_br_obs(i_obs,*), - depth_nut(i_obs,*)/h_nut(i_obs), PSYM = 1, THICK = 3, SYMSIZE = 0.8
801
802   ; add model
803   i_mod = WHERE(doy EQ doy_obs(i_obs))
804   zz = FLTARR(nlay)
805   FOR i = 0, nlay - 1 DO BEGIN
806      zz(i) = z_ib(i,i_mod) / h_i(i_mod)
807   ENDFOR
808   LOADCT, 9
809   OPLOT, dsib(*,i_mod)/e_i(*,i_mod)*100., -zz, THICK = 3, COLOR = 150
810   LOADCT, 1
811
812   OPLOT, [ 3.90, 3.90 ], [ -1., 0.], linestyle = 1, color = 150 , thick = 3
813
814ENDFOR
815
816;
817;==============================================================================
818; Mean profiles
819;==============================================================================
820;
821!P.MULTI=[0,numplot_x, numplot_y]
822
823;--- T ---
824
825; prepare plot
826LOADCT, 0
827ztitle = "Mean T"
828PLOT, [ -20., 0. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'deg C', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2, TITLE = ztitle
829
830; add obs
831OPLOT, t_i_obs_mean, - depth_std, PSYM = 8, THICK = 3, SYMSIZE = 0.8
832FOR i = 0, N_std - 1 DO BEGIN
833   OPLOT, [ t_i_obs_mean(i) - t_i_obs_std(i), t_i_obs_mean(i) + t_i_obs_std(i) ], $
834          [ -depth_std(i), -depth_std(i) ]
835ENDFOR
836LOADCT, 3
837OPLOT, t_i_mod_mean, -zz, THICK = 3, COLOR = 150
838LOADCT, 0
839OPLOT, t_i_mod_mean+t_i_mod_std, -zz, THICK = 1, COLOR = 150
840OPLOT, t_i_mod_mean-t_i_mod_std, -zz, THICK = 1, COLOR = 150
841
842;--- S ---
843
844; prepare plot
845LOADCT, 0
846ztitle = "Mean S"
847PLOT, [0., 12. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'g/kg', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2, TITLE = ztitle
848
849; add obs
850OPLOT, s_i_obs_mean, - depth_std, PSYM = 8, THICK = 3, SYMSIZE = 0.8
851FOR i = 0, N_std - 1 DO BEGIN
852   OPLOT, [ s_i_obs_mean(i) - s_i_obs_std(i), s_i_obs_mean(i) + s_i_obs_std(i) ], $
853          [ -depth_std(i), -depth_std(i) ]
854ENDFOR
855
856; add model
857LOADCT, 3
858OPLOT, s_i_mod_mean, -zz, THICK = 3, COLOR = 150
859LOADCT, 0
860OPLOT, s_i_mod_mean+s_i_mod_std, -zz, THICK = 1, COLOR = 150
861OPLOT, s_i_mod_mean-s_i_mod_std, -zz, THICK = 1, COLOR = 150
862
863;--- e ---
864
865; prepare plot
866LOADCT, 0
867ztitle = "Mean e"
868   PLOT, [0., 50. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = '%', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2, TITLE = ztitle
869
870; add obs
871OPLOT, e_i_obs_mean*100., - depth_std, PSYM = 8, THICK = 3, SYMSIZE = 0.8
872FOR i = 0, N_std - 1 DO BEGIN
873   OPLOT, [ e_i_obs_mean(i) - e_i_obs_std(i), e_i_obs_mean(i) + e_i_obs_std(i) ] * 100., $
874          [ -depth_std(i), -depth_std(i) ]
875ENDFOR
876
877; add model
878LOADCT, 3
879OPLOT, e_i_mod_mean, -zz, THICK = 3, COLOR = 150
880LOADCT, 0
881OPLOT, e_i_mod_mean+e_i_mod_std, -zz, THICK = 1, COLOR = 150
882OPLOT, e_i_mod_mean-e_i_mod_std, -zz, THICK = 1, COLOR = 150
883
884
885
886;--- chl-a ---
887
888; prepare plot
889LOADCT, 0
890ztitle = "Mean chla"
891PLOT, [ 0.01, 40. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mg chla / m3', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2, TITLE = ztitle
892; add obs
893OPLOT, chla_obs_mean, - depth_std, PSYM = 8, THICK = 3, SYMSIZE = 0.8
894FOR i = 0, N_std - 1 DO BEGIN
895   OPLOT, [ chla_obs_mean(i) - chla_obs_std(i), chla_obs_mean(i) + chla_obs_std(i) ], $
896          [ -depth_std(i), -depth_std(i) ]
897ENDFOR
898; plot seawater value
899LOADCT, 1
900OPLOT, [ 0.08 ], [ -1 ], PSYM = 8, THICK = 3, SYMSIZE = 1.6, COLOR = 150
901
902; add model
903LOADCT, 9
904OPLOT, chla_mod_mean, -zz, THICK = 3, COLOR = 150
905LOADCT, 0
906OPLOT, chla_mod_mean+chla_mod_std, -zz, THICK = 1, COLOR = 150
907OPLOT, chla_mod_mean-chla_mod_std, -zz, THICK = 1, COLOR = 150
908
909;--- NOX ---
910
911; prepare plot
912LOADCT, 0
913ztitle = "Mean NOX"
914PLOT, [ 0.00, 35. ], [ -1.0, 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mmmol N / m3', XTICKFORMAT='(I3)', XMINOR = 2, YMINOR = 2, TITLE = ztitle
915; add obs
916OPLOT, nox_obs_mean, - depth_nut_mean, PSYM = 8, THICK = 3, SYMSIZE = 0.8
917FOR i = 0, 5 DO BEGIN
918   OPLOT, [ nox_obs_mean(i) - nox_obs_std(i), nox_obs_mean(i) + nox_obs_std(i) ], $
919          [ -depth_nut_mean(i), -depth_nut_mean(i) ]
920ENDFOR
921
922; plot seawater value
923LOADCT, 1
924OPLOT, [ 31.16 ], [ -1 ], PSYM = 8, THICK = 3, SYMSIZE = 1.6, COLOR = 150
925
926; add model
927LOADCT, 9
928OPLOT, nox_mod_mean, -zz, THICK = 3, COLOR = 150
929LOADCT, 0
930OPLOT, nox_mod_mean+nox_mod_std, -zz, THICK = 1, COLOR = 150
931OPLOT, nox_mod_mean-nox_mod_std, -zz, THICK = 1, COLOR = 150
932
933;--- PO4 ---
934
935; prepare plot
936LOADCT, 0
937ztitle = "Mean PO4"
938PLOT, [ 0.00, 4.0 ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mmmol P / m3', XTICKFORMAT='(I2)', XMINOR = 2, YMINOR = 2, TITLE = ztitle
939; add obs
940OPLOT, po4_obs_mean, - depth_nut_mean, PSYM = 8, THICK = 3, SYMSIZE = 0.8
941FOR i = 0, N_std - 1 DO BEGIN
942   OPLOT, [ po4_obs_mean(i) - po4_obs_std(i), po4_obs_mean(i) + po4_obs_std(i) ], $
943          [ -depth_nut_mean(i), -depth_nut_mean(i) ]
944ENDFOR
945
946; plot seawater value
947LOADCT, 1
948OPLOT, [ 2.05 ], [ -1 ], PSYM = 8, THICK = 3, SYMSIZE = 1.6, COLOR = 150
949
950; add model
951LOADCT, 9
952OPLOT, po4_mod_mean, -zz, THICK = 3, COLOR = 150
953LOADCT, 0
954OPLOT, po4_mod_mean+po4_mod_std, -zz, THICK = 1, COLOR = 150
955OPLOT, po4_mod_mean-po4_mod_std, -zz, THICK = 1, COLOR = 150
956
957;--- DSi ---
958
959; prepare plot
960LOADCT, 0
961ztitle = "Mean DSi"
962PLOT, [ 0.00, 90. ], [ -1., 0. ], CHARSIZE = cs, /NODATA, SUBTITLE = 'mmmol N / m3', XTICKFORMAT='(I2)', XMINOR = 2, YMINOR = 2, TITLE = ztitle
963; add obs
964OPLOT, dsi_obs_mean, - depth_nut_mean, PSYM = 8, THICK = 3, SYMSIZE = 0.8
965FOR i = 0, N_std - 1 DO BEGIN
966   OPLOT, [ dsi_obs_mean(i) - dsi_obs_std(i), dsi_obs_mean(i) + dsi_obs_std(i) ], $
967          [ -depth_nut_mean(i), -depth_nut_mean(i) ]
968ENDFOR
969
970; plot seawater value
971LOADCT, 1
972OPLOT, [ 84.35 ], [ -1 ], PSYM = 8, THICK = 3, SYMSIZE = 1.6, COLOR = 150
973
974; add model
975LOADCT, 9
976OPLOT, dsi_mod_mean, -zz, THICK = 3, COLOR = 150
977LOADCT, 0
978OPLOT, dsi_mod_mean+dsi_mod_std, -zz, THICK = 1, COLOR = 150
979OPLOT, dsi_mod_mean-dsi_mod_std, -zz, THICK = 1, COLOR = 150
980
981;
982;==============================================================================
983; End of the script
984;==============================================================================
985;
986IF ( device EQ 'PS' ) THEN BEGIN
987   DEVICE, /CLOSE
988   SET_PLOT, "X"
989   !P.MULTI=[0,2,2]
990ENDIF
991
992END
993
Note: See TracBrowser for help on using the repository browser.