New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
std_plot_all.pro in branches/2014/dev_r4650_UKMO14.11_SETTE_OBSASM/NEMOGCM/CONFIG/ORCA2_LIM_OBS/IDL_scripts – NEMO

source: branches/2014/dev_r4650_UKMO14.11_SETTE_OBSASM/NEMOGCM/CONFIG/ORCA2_LIM_OBS/IDL_scripts/std_plot_all.pro @ 4751

Last change on this file since 4751 was 4751, checked in by djlea, 10 years ago

Changes to include an OBS test in SETTE. At the moment this uses an example profile observation.
File size: 30.9 KB
Line 
1pro std_plot_all, doplot = doplot, _extra = ex
2
3  compile_opt idl2, strictarrsubs
4
5@common
6@std_common
7                                ; scripts for nemo v3_2 and v3_3
8
9  PRINT, ''
10  PRINT, '  ############################################'
11  PRINT, ''
12  PRINT, '                    LAUNCH of std_plots'
13  PRINT, ''
14  PRINT, '  ############################################'
15  PRINT, ''
16;
17  std_iodir_data    = isadirectory(getenv('DIR_DATA'),     title = 'path of data in NetCdf format')
18  std_iodir_climato = isadirectory(getenv('DIR_CLIMATO'),  title = 'path of climatological data')
19  std_iodir_mask    = isadirectory(getenv('DIR_MASK'),     title = 'path of mask files (ex: subbasins)')
20; meshmask
21  std_file_mesh = isafile(getenv('FILE_MESH_MASK'),        title = 'mesh_mask', iodir = std_iodir_mask)
22  std_file_msksub = isafile(getenv('FILE_MASK_SUBDOMAIN'), title = 'sub-bassin masks', iodir = std_iodir_mask)
23
24; climatologies
25  std_file_Levitus_T =  isafile(getenv('FILE_TEMP_3D'),    title = 'Levitus_T', iodir = std_iodir_climato)
26  std_file_Levitus_S =  isafile(getenv('FILE_SAL_3D'),     title = 'Levitus_S', iodir = std_iodir_climato)
27  std_file_reynolds  =  isafile(getenv('FILE_SST'),        title = 'Reynolds', iodir = std_iodir_climato)
28  std_file_oaflux    =  isafile(getenv('FILE_FLUX'),       title = 'oaflux', iodir = std_iodir_climato)
29  std_file_mld       =  isafile(getenv('FILE_MLD'),        title = 'Mixed layer depth', iodir = std_iodir_climato)
30  std_file_ice       =  isafile(getenv('FILE_ICE'),        title = 'ICE', iodir = std_iodir_climato)
31  std_file_snow_arc  =  isafile(getenv('FILE_SNOW_ARC'),   title = 'SNOW_ARC', iodir = std_iodir_climato)
32  std_file_snow_ant  =  isafile(getenv('FILE_SNOW_ANT'),   title = 'SNOW_ANT', iodir = std_iodir_climato)
33
34  IF strlowcase(getenv('FILE_GEOHEAT')) EQ 'no' THEN std_file_geoheat = 'no' $
35  ELSE std_file_geoheat =  isafile(getenv('FILE_GEOHEAT'), title = 'Geothermal heating', iodir = std_iodir_climato)
36;
37  allrec =  1 - keyword_set(long(getenv('READ_ONLY_FIRST_RECORD')))
38 
39; Output run experience1
40  std_file1_T     = isafile(getenv('FILE1_T'), title = 'exp1 grid T input file', iodir = std_iodir_data)
41  std_file1_U     = isafile(getenv('FILE1_U'), title = 'exp1 grid U input file', iodir = std_iodir_data)
42  std_file1_V     = isafile(getenv('FILE1_V'), title = 'exp1 grid V input file', iodir = std_iodir_data)
43  std_file1_I     = isafile(getenv('FILE1_I'), title = 'exp1 ice    input file', iodir = std_iodir_data)
44 
45; Output run experience2
46  std_file2_T     = isafile(getenv('FILE2_T'), title = 'exp2 grid T input file', iodir = std_iodir_data)
47  std_file2_U     = isafile(getenv('FILE2_U'), title = 'exp2 grid U input file', iodir = std_iodir_data)
48  std_file2_V     = isafile(getenv('FILE2_V'), title = 'exp2 grid V input file', iodir = std_iodir_data)
49  std_file2_I     = isafile(getenv('FILE2_I'), title = 'exp2 ice    input file', iodir = std_iodir_data)
50
51  PRINT, ''
52  PRINT, '  std_iodir_data : ' + std_iodir_data
53  PRINT, '  std_file1T : ' + std_file1_T
54  PRINT, '  std_file1U : ' + std_file1_U
55  PRINT, '  std_file1V : ' + std_file1_V
56;  PRINT, ' std_file1W : ' + std_file1_W
57  PRINT, '  std_file2I : ' + std_file1_I
58  PRINT, '  std_file2T : ' + std_file2_T
59  PRINT, '  std_file2U : ' + std_file2_U
60  PRINT, '  std_file2V : ' + std_file2_V
61;  PRINT, ' std_file2W : ' + std_file2_W
62  PRINT, '  std_file2I : ' + std_file2_I
63  PRINT, ''
64
65;#########################################################################
66;##########################  Load Grids   ################################
67;#########################################################################
68; load the grid
69  load_orca, std_file_mesh 
70; reading variables
71  masknp = read_ncdf('tmaskutil', file = std_file_mesh, /nostruct, /cont_nofill)
72;#########################################################################
73;############################  Read Data  ################################
74;#########################################################################
75;
76  allrec =  1; - keyword_set(long(getenv('READ_ONLY_FIRST_RECORD')))
77;
78;;; 3D ;;;
79; temperature
80  T1 = read_ncdf(getenv('VAR1_T'), allrecords = allrec, direc = 't', filename = std_file1_T )
81  IF std_file2_T NE std_file1_T THEN BEGIN
82    T2 = read_ncdf(getenv('VAR2_T'), allrecords = allrec, direc = 't', filename = std_file2_T )
83  ENDIF ELSE T2 = {arr:-1}
84  TLev = read_ncdf(getenv('VAR_TEMP_3D'), filename = std_file_Levitus_T )
85  TRey = read_ncdf(getenv('VAR_SST'), filename = std_file_reynolds )
86
87; salinity
88  S1 = read_ncdf(getenv('VAR1_S'), allrecords = allrec, direc = 't', filename = std_file1_T )
89  IF std_file2_T NE std_file1_T THEN BEGIN
90    S2 = read_ncdf(getenv('VAR2_S'), allrecords = allrec, direc = 't', filename = std_file2_T )
91  ENDIF ELSE S2 = {arr:-1}
92  SLev = read_ncdf(getenv('VAR_SAL_3D'), filename = std_file_Levitus_S )
93
94;;; 2D ;;;
95; Net Downward heat flux
96  Q1 = read_ncdf(getenv('VAR1_QNET'), allrecords = allrec, direc = 't', filename = std_file1_T )
97  IF std_file2_T NE std_file1_T THEN BEGIN
98    Q2 = read_ncdf(getenv('VAR2_QNET'), allrecords = allrec, direc = 't', filename = std_file2_T )
99  ENDIF ELSE Q2 = {arr:-1}
100; Geothermal heating
101  IF std_file_geoheat EQ 'no' THEN geo = {arr:float(getenv('VAR_GEOHEAT'))} $
102  ELSE geo = read_ncdf(getenv('VAR_GEOHEAT'), filename =  std_file_geoheat )
103  geo = geo.arr*1.e-3          ; convert into W/m2
104;climatology
105  QNET = read_ncdf(getenv('VAR_FLUX'), filename = std_file_oaflux )
106
107; erp (evaporation damping)
108  ERP1 = read_ncdf(getenv('VAR1_ERP'), allrecords = allrec, direc = 't', filename = std_file1_T )
109  ERP1 = {arr:ERP1.arr * 86400., unit:'mm/day', grid:'T'}
110  IF std_file2_T NE std_file1_T THEN BEGIN
111     ERP2 = read_ncdf(getenv('VAR2_ERP'), allrecords = allrec, direc = 't', filename = std_file2_T )
112     ERP2 = {arr:ERP2.arr * 86400., unit:'mm/day', grid:'T'}
113  ENDIF ELSE ERP2 = {arr:-1}
114 
115; emp (evaporation minus precipitation)
116  EMP1 = read_ncdf(getenv('VAR1_EMP'), allrecords = allrec, direc = 't', filename = std_file1_T )
117  EMP1 = {arr:EMP1.arr * 86400., unit:'mm/day', grid:'T'}
118   IF std_file2_T NE std_file1_T THEN BEGIN
119     EMP2 = read_ncdf(getenv('VAR2_EMP'), allrecords = allrec, direc = 't', filename = std_file2_T )
120     EMP2 = {arr:EMP2.arr * 86400., unit:'mm/day', grid:'T'}
121   ENDIF ELSE EMP2 = {arr:-1}
122   
123 ;mixed layer depth
124   MLD1 = read_ncdf(getenv('VAR1_MLD'), allrecords = allrec, direc = 't', filename = std_file1_T ) ; 10 m
125   IF std_file2_T NE std_file1_T THEN BEGIN
126     MLD2 = read_ncdf(getenv('VAR2_MLD'), allrecords = allrec, direc = 't', filename = std_file2_T ) ; 10 m
127   ENDIF ELSE MLD2 = {arr:-1}
128 ;climatology
129   MLD = read_ncdf(getenv('VAR_MLD'), filename = std_file_mld )
130
131 ; velocities
132   U1 = read_ncdf(getenv('VAR1_U'), allrecords = allrec, direc = 't', filename = std_file1_U )
133   ; old formulation: we tested variable name
134   ; IF strlowcase(getenv('VAR1_U')) EQ 'uocetr_eff' OR strlowcase(getenv('VAR1_U')) EQ 'vozoeftr' THEN BEGIN
135   IF strlowcase(U1.unit) EQ 'm3/s' THEN BEGIN
136     ;IF it is a transport it is transofrmed in velocity
137     U1.arr = U1.arr / e3u_3d(/e2) * umask()
138     U1.unit = 'm/s'
139   ENDIF
140   IF std_file2_U NE std_file1_U THEN BEGIN
141     U2 = read_ncdf(getenv('VAR2_U'), allrecords = allrec, direc = 't', filename = std_file2_U )
142     ; old formulation: we tested variable name
143     ; IF strlowcase(getenv('VAR2_U')) EQ 'uocetr_eff' OR strlowcase(getenv('VAR2_U')) EQ 'vozoeftr' THEN BEGIN
144     IF strlowcase(U2.unit) EQ 'm3/s' THEN BEGIN
145       U2.arr = U2.arr / e3u_3d(/e2) * umask()
146       U2.unit = 'm/s'
147     ENDIF
148   ENDIF ELSE U2 = {arr:-1}
149;
150   V1 = read_ncdf(getenv('VAR1_V'), allrecords = allrec, direc = 't', filename = std_file1_V )
151   ; old formulation: we tested variable name
152   ; IF strlowcase(getenv('VAR1_V')) EQ 'vocetr_eff' OR strlowcase(getenv('VAR1_V')) EQ 'vomeeftr' THEN BEGIN
153   IF strlowcase(V1.unit) EQ 'm3/s' THEN BEGIN ; test on unit to understand if it is a transport or velocity
154     ;IF it is a transport it is transofrmed in velocity
155     V1.arr = V1.arr / e3v_3d(/e1) * vmask()
156     V1.unit = 'm/s'
157   ENDIF
158   IF std_file2_V NE std_file1_V THEN BEGIN
159     V2 = read_ncdf(getenv('VAR2_V'), allrecords = allrec, direc = 't', filename = std_file2_V )
160     ; old formulation
161     ; IF strlowcase(getenv('VAR2_V')) EQ 'vocetr_eff' OR strlowcase(getenv('VAR2_V')) EQ 'vozoeftr' THEN BEGIN
162     IF strlowcase(V2.unit) EQ 'm3/s' THEN BEGIN
163       V2.arr = V2.arr / e3v_3d(/e1) * vmask()
164       V2.unit = 'm/s'
165     ENDIF
166   ENDIF ELSE V2 = {arr:-1}
167
168; ice
169  Ithi_1 = read_ncdf(getenv('VAR1_Ithick'), allrecords = allrec, filename = std_file1_I )
170  caldat, time, mm
171  april = where(mm EQ 4, cnt)
172  Ithi_april_1 = {arr:1./float(cnt) * total(reform(Ithi_1.arr[*, *, temporary(april)],nxt,nyt,cnt), 3), unit:Ithi_1.unit}
173  jan = where(mm EQ 1, cnt)
174  sept = where(mm EQ 9, cnt)
175  Ithi_jan_1 = {arr:1./float(cnt) * total(reform(Ithi_1.arr[*, *, temporary(jan)],nxt,nyt,cnt), 3), unit:Ithi_1.unit}
176  Ithi_sept_1 = {arr:1./float(cnt) * total(reform(Ithi_1.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Ithi_1.unit}
177  undefine, Ithi_1
178;
179  IF std_file2_I NE std_file1_I THEN BEGIN
180     Ithi_2 = read_ncdf(getenv('VAR2_Ithick'), allrecords = allrec, filename = std_file2_I )
181     caldat, time, mm
182     april = where(mm EQ 4, cnt)
183     Ithi_april_2 = {arr:1./float(cnt) * total(reform(Ithi_2.arr[*, *, temporary(april)],nxt,nyt,cnt), 3), unit:Ithi_2.unit}
184     jan = where(mm EQ 1, cnt)
185     sept = where(mm EQ 9, cnt)
186     Ithi_sept_2 = {arr:1./float(cnt) * total(reform(Ithi_2.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Ithi_2.unit}
187     Ithi_jan_2 = {arr:1./float(cnt) * total(reform(Ithi_2.arr[*, *, temporary(jan)],nxt,nyt,cnt), 3), unit:Ithi_2.unit}
188     undefine, Ithi_2
189  ENDIF ELSE BEGIN
190     Ithi_april_2 = {arr:-1}
191     Ithi_sept_2 = {arr:-1}
192     Ithi_jan_2 = {arr:-1}
193  ENDELSE
194;
195  Ifra_1 = read_ncdf(getenv('VAR1_Ifrac'), allrecords = allrec, filename = std_file1_I )
196  caldat, time, mm
197  jan = where(mm EQ 1, cnt)
198  Ifra_jan_1 = {arr:1./float(cnt) * total(reform(Ifra_1.arr[*, *, temporary(jan)],nxt,nyt,cnt), 3), unit:Ifra_1.unit}
199  febr = where(mm EQ 2, cnt)
200  Ifra_febr_1 = {arr:1./float(cnt) * total(reform(Ifra_1.arr[*, *, temporary(febr)],nxt,nyt,cnt), 3), unit:Ifra_1.unit}
201  march = where(mm EQ 3, cnt)
202  Ifra_march_1 = {arr:1./float(cnt) * total(reform(Ifra_1.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Ifra_1.unit}
203  april = where(mm EQ 4, cnt)
204  Ifra_april_1 = {arr:1./float(cnt) * total(reform(Ifra_1.arr[*, *, temporary(april)],nxt,nyt,cnt), 3), unit:Ifra_1.unit}
205  sept = where(mm EQ 9, cnt)
206  Ifra_sept_1 = {arr:1./float(cnt) * total(reform(Ifra_1.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Ifra_1.unit}
207  undefine, Ifra_1
208;
209  IF std_file2_I NE std_file1_I THEN BEGIN
210    Ifra_2 = read_ncdf(getenv('VAR2_Ifrac'), allrecords = allrec, filename = std_file2_I )
211    caldat, time, mm
212    jan = where(mm EQ 1, cnt)
213    Ifra_jan_2 = {arr:1./float(cnt) * total(reform(Ifra_2.arr[*, *, temporary(jan)],nxt,nyt,cnt), 3), unit:Ifra_2.unit}
214    febr = where(mm EQ 2, cnt)
215    Ifra_febr_2 = {arr:1./float(cnt) * total(reform(Ifra_2.arr[*, *, temporary(febr)],nxt,nyt,cnt), 3), unit:Ifra_2.unit}
216    march = where(mm EQ 3, cnt)
217    Ifra_march_2 = {arr:1./float(cnt) * total(reform(Ifra_2.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Ifra_2.unit}
218    april = where(mm EQ 4, cnt)
219    Ifra_april_2 = {arr:1./float(cnt) * total(reform(Ifra_2.arr[*, *, temporary(april)],nxt,nyt,cnt), 3), unit:Ifra_2.unit}
220    sept = where(mm EQ 9, cnt)
221    Ifra_sept_2 = {arr:1./float(cnt) * total(reform(Ifra_2.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Ifra_2.unit}
222    undefine, Ifra_2
223  ENDIF ELSE BEGIN
224    Ifra_jan_2 = {arr:-1}
225    Ifra_febr_2 = {arr:-1}
226    Ifra_march_2 = {arr:-1}
227    Ifra_april_2 = {arr:-1}
228    Ifra_sept_2 = {arr:-1}
229  ENDELSE
230;
231   Isnow_1 = read_ncdf(getenv('VAR1_Isnow'), allrecords = allrec, filename = std_file1_I )
232   caldat, time, mm
233   febr = where(mm EQ 2, cnt)
234   Isnow_febr_1 = {arr:1./float(cnt) * total(reform(Isnow_1.arr[*, *, temporary(febr)],nxt,nyt,cnt), 3), unit:Isnow_1.unit}
235   march = where(mm EQ 3, cnt)
236   Isnow_march_1 = {arr:1./float(cnt) * total(reform(Isnow_1.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Isnow_1.unit}
237   sept = where(mm EQ 9, cnt)
238   Isnow_sept_1 = {arr:1./float(cnt) * total(reform(Isnow_1.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Isnow_1.unit}
239   ;undefine, Isnow_1
240 ;
241   IF std_file2_I NE std_file1_I THEN BEGIN
242     Isnow_2 = read_ncdf(getenv('VAR2_Isnow'), allrecords = allrec, filename = std_file2_I )
243     caldat, time, mm
244     febr = where(mm EQ 2, cnt)
245     Isnow_febr_2 = {arr:1./float(cnt) * total(reform(Isnow_2.arr[*, *, temporary(febr)],nxt,nyt,cnt), 3), unit:Isnow_2.unit}
246     march = where(mm EQ 3, cnt)
247     Isnow_march_2 = {arr:1./float(cnt) * total(reform(Isnow_2.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Isnow_2.unit}
248     sept = where(mm EQ 9, cnt)
249     Isnow_sept_2 = {arr:1./float(cnt) * total(reform(Isnow_2.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Isnow_2.unit}
250   ;  undefine, Isnow_2
251   ENDIF ELSE BEGIN
252     Isnow_febr_2 = {arr:-1}
253     Isnow_march_2 = {arr:-1}
254     Isnow_sept_2 = {arr:-1}
255   ENDELSE
256;
257   Isal_1 = read_ncdf(getenv('VAR1_Isal'), allrecords = allrec, filename = std_file1_I )
258;SF  ready for mask : remove 0.15% for observations
259;SF
260;SF   Ifra_1 = read_ncdf(getenv('VAR1_Ifrac'), allrecords = allrec, filename = std_file1_I )
261;SF   msk = Ifra_1.arr gt 0.15 ; remove 0.15% for observations
262   caldat, time, mm
263   march = where(mm EQ 3, cnt)
264   febr = where(mm EQ 2, cnt)
265   Isal_march_1 = {arr:1./float(cnt) * total(reform(Isal_1.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Isal_1.unit}
266;SF
267;SF   Isal_1.arr = Isal_1.arr * msk
268;SF
269;SF   Isal_march_1 = {arr:1./float(cnt) * total(reform(Isal_1.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Isal_1.unit}
270   Isal_febr_1 = {arr:1./float(cnt) * total(reform(Isal_1.arr[*, *, temporary(febr)],nxt,nyt,cnt), 3), unit:Isal_1.unit}
271   sept = where(mm EQ 9, cnt)
272   Isal_sept_1 = {arr:1./float(cnt) * total(reform(Isal_1.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Isal_1.unit}
273   ;undefine, Isal_1
274 ;
275   Isal_2 = read_ncdf(getenv('VAR2_Isal'), allrecords = allrec, filename = std_file2_I )
276   IF std_file2_I NE std_file1_I THEN BEGIN
277     caldat, time, mm
278     march = where(mm EQ 3, cnt)
279     febr = where(mm EQ 2, cnt)
280     Isal_march_2 = {arr:1./float(cnt) * total(reform(Isal_2.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Isal_2.unit}
281     Isal_febr_2 = {arr:1./float(cnt) * total(reform(Isal_2.arr[*, *, temporary(febr)],nxt,nyt,cnt), 3), unit:Isal_2.unit}
282     sept = where(mm EQ 9, cnt)
283     Isal_sept_2 = {arr:1./float(cnt) * total(reform(Isal_2.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Isal_2.unit}
284   ;  undefine, Isal_2
285   ENDIF ELSE BEGIN
286     Isal_febr_2 = {arr:-1}
287     Isal_march_2 = {arr:-1}
288     Isal_sept_2 = {arr:-1}
289   ENDELSE
290;
291  IvelU_1 = read_ncdf(getenv('VAR1_IvelU'), allrecords = allrec, filename = std_file1_I )
292  IvelV_1 = read_ncdf(getenv('VAR1_IvelV'), allrecords = allrec, filename = std_file1_I )
293  Ivelo_1 = read_ncdf(getenv('VAR1_Ivelo'), allrecords = allrec, filename = std_file1_I )
294  caldat, time, mm
295  febr = where(mm EQ 2, cnt)
296  IvelU_febr_1 = {arr:1./float(cnt) * total(reform(IvelU_1.arr[*, *, febr],nxt,nyt,cnt), 3), unit:IvelU_1.unit, g: 'T'}
297  IvelV_febr_1 = {arr:1./float(cnt) * total(reform(IvelV_1.arr[*, *, febr],nxt,nyt,cnt), 3), unit:IvelV_1.unit, g: 'T'}
298  Ivelo_febr_1 = {arr:1./float(cnt) * total(reform(Ivelo_1.arr[*, *, febr],nxt,nyt,cnt), 3), unit:Ivelo_1.unit, g: 'T'}
299  march = where(mm EQ 3, cnt)
300  IvelU_march_1 = {arr:1./float(cnt) * total(reform(IvelU_1.arr[*, *, march],nxt,nyt,cnt), 3), unit:IvelU_1.unit, g: 'T'}
301  IvelV_march_1 = {arr:1./float(cnt) * total(reform(IvelV_1.arr[*, *, march],nxt,nyt,cnt), 3), unit:IvelV_1.unit, g: 'T'}
302  Ivelo_march_1 = {arr:1./float(cnt) * total(reform(Ivelo_1.arr[*, *, march],nxt,nyt,cnt), 3), unit:Ivelo_1.unit, g: 'T'}
303  sept = where(mm EQ 9, cnt)
304  IvelU_sept_1 = {arr:1./float(cnt) * total(reform(IvelU_1.arr[*, *, sept],nxt,nyt,cnt), 3), unit:IvelU_1.unit, g: 'T'}
305  IvelV_sept_1 = {arr:1./float(cnt) * total(reform(IvelV_1.arr[*, *, sept],nxt,nyt,cnt), 3), unit:IvelV_1.unit, g: 'T'}
306  Ivelo_sept_1 = {arr:1./float(cnt) * total(reform(Ivelo_1.arr[*, *, sept],nxt,nyt,cnt), 3), unit:Ivelo_1.unit, g: 'T'}
307;
308  IF std_file2_I NE std_file1_I THEN BEGIN
309    IvelU_2 = read_ncdf(getenv('VAR2_IvelU'), allrecords = allrec, filename = std_file2_I )
310    IvelV_2 = read_ncdf(getenv('VAR2_IvelV'), allrecords = allrec, filename = std_file2_I )
311    Ivelo_2 = read_ncdf(getenv('VAR2_Ivelo'), allrecords = allrec, filename = std_file2_I )
312    caldat, time, mm
313    febr = where(mm EQ 2, cnt)
314    IvelU_febr_2 = {arr:1./float(cnt) * total(reform(IvelU_2.arr[*, *, febr],nxt,nyt,cnt), 3), unit:IvelU_2.unit}
315    IvelV_febr_2 = {arr:1./float(cnt) * total(reform(IvelV_2.arr[*, *, febr],nxt,nyt,cnt), 3), unit:IvelV_2.unit}
316    Ivelo_febr_2 = {arr:1./float(cnt) * total(reform(Ivelo_2.arr[*, *, febr],nxt,nyt,cnt), 3), unit:Ivelo_2.unit}
317    march = where(mm EQ 3, cnt)
318    IvelU_march_2 = {arr:1./float(cnt) * total(reform(IvelU_2.arr[*, *, march],nxt,nyt,cnt), 3), unit:IvelU_2.unit}
319    IvelV_march_2 = {arr:1./float(cnt) * total(reform(IvelV_2.arr[*, *, march],nxt,nyt,cnt), 3), unit:IvelV_2.unit}
320    Ivelo_march_2 = {arr:1./float(cnt) * total(reform(Ivelo_2.arr[*, *, march],nxt,nyt,cnt), 3), unit:Ivelo_2.unit}
321    sept = where(mm EQ 9, cnt)
322    IvelU_sept_2 = {arr:1./float(cnt) * total(reform(IvelU_2.arr[*, *, sept],nxt,nyt,cnt), 3), unit:IvelU_2.unit}
323    IvelV_sept_2 = {arr:1./float(cnt) * total(reform(IvelV_2.arr[*, *, sept],nxt,nyt,cnt), 3), unit:IvelV_2.unit}   
324    Ivelo_sept_2 = {arr:1./float(cnt) * total(reform(Ivelo_2.arr[*, *, sept],nxt,nyt,cnt), 3), unit:Ivelo_2.unit}
325    undefine, Ifra_2
326  ENDIF ELSE BEGIN
327    IvelU_febr_2 = {arr:-1}
328    IvelV_febr_2 = {arr:-1}
329    Ivelo_febr_2 = {arr:-1}
330    IvelU_march_2 = {arr:-1}
331    IvelV_march_2 = {arr:-1}
332    Ivelo_march_2 = {arr:-1}
333    IvelU_sept_2 = {arr:-1}
334    IvelV_sept_2 = {arr:-1}
335    Ivelo_sept_2 = {arr:-1}
336  ENDELSE
337;
338  Iage_1 = read_ncdf(getenv('VAR1_Iage'), allrecords = allrec, filename = std_file1_I )
339  caldat, time, mm
340  febr = where(mm EQ 2, cnt)
341  Iage_febr_1 = {arr:1./float(cnt) * total(reform(Iage_1.arr[*, *, temporary(febr)],nxt,nyt,cnt), 3), unit:Iage_1.unit}
342  march = where(mm EQ 3, cnt)
343  Iage_march_1 = {arr:1./float(cnt) * total(reform(Iage_1.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Iage_1.unit}
344  sept = where(mm EQ 9, cnt)
345  Iage_sept_1 = {arr:1./float(cnt) * total(reform(Iage_1.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Iage_1.unit}
346  undefine, Iage_1
347;
348  IF std_file2_I NE std_file1_I THEN BEGIN
349     Iage_2 = read_ncdf(getenv('VAR2_Iage'), allrecords = allrec, filename = std_file2_I )
350     caldat, time, mm
351     febr = where(mm EQ 2, cnt)
352     Iage_febr_2 = {arr:1./float(cnt) * total(reform(Iage_2.arr[*, *, temporary(febr)],nxt,nyt,cnt), 3), unit:Iage_2.unit}
353     march = where(mm EQ 3, cnt)
354     Iage_march_2 = {arr:1./float(cnt) * total(reform(Iage_2.arr[*, *, temporary(march)],nxt,nyt,cnt), 3), unit:Iage_2.unit}
355     sept = where(mm EQ 9, cnt)
356     Iage_sept_2 = {arr:1./float(cnt) * total(reform(Iage_2.arr[*, *, temporary(sept)],nxt,nyt,cnt), 3), unit:Iage_2.unit}
357     undefine, Iage_2
358  ENDIF ELSE BEGIN
359     Iage_febr_2 = {arr:-1}
360     Iage_march_2 = {arr:-1}
361     Iage_sept_2 = {arr:-1}
362  ENDELSE
363;
364  jpt = 1
365;
366; shorter file names for legends...
367;
368  std_file1_T = file_basename(std_file1_T,'.nc')
369  std_file1_T = (strsplit(std_file1_T,'_grid_T',/extract,/regex))[0]
370  std_file2_T = file_basename(std_file2_T,'.nc')
371  std_file2_T = (strsplit(std_file2_T,'_grid_T',/extract,/regex))[0]
372  std_file1_U = file_basename(std_file1_U,'.nc')
373  std_file1_U = (strsplit(std_file1_U,'_grid_U',/extract,/regex))[0]
374  std_file2_U = file_basename(std_file2_U,'.nc')
375  std_file2_U = (strsplit(std_file2_U,'_grid_U',/extract,/regex))[0]
376  std_file1_V = file_basename(std_file1_V,'.nc')
377  std_file1_V = (strsplit(std_file1_V,'_grid_V',/extract,/regex))[0]
378  std_file2_V = file_basename(std_file2_V,'.nc')
379  std_file2_V = (strsplit(std_file2_V,'_grid_V',/extract,/regex))[0]
380  std_file1_I = file_basename(std_file1_I,'.nc')
381  std_file1_I = (strsplit(std_file1_I,'_icemod',/extract,/regex))[0]
382  std_file2_I = file_basename(std_file2_I,'.nc')
383  std_file2_I = (strsplit(std_file2_I,'_icemod',/extract,/regex))[0]
384
385;#########################################################################
386;######################  STANDARD PLOTS   ################################
387;#########################################################################
388
389  IF keyword_set(doplot) EQ 0 THEN doplot = 0
390
391; fixed color tabled
392  lct, 64
393  cnt = 0
394  htmltxt = ''
395;
396  cnt = cnt+1   &   blabla = 'Erp salinity damping term'
397  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_erp, ERP1, ERP2, _extra = ex
398;
399  cnt = cnt+1   &   blabla = 'Evaporation - Precipitation - Runoff term'
400  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_emp, EMP1, EMP2, _extra = ex
401;
402  cnt = cnt+1   &   blabla = 'Net heat flux'
403  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_qnet, Q1, Q2, QNET, _extra = ex
404;
405  cnt = cnt+1   &   blabla = 'Meridionnal Heat Transport'
406  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_mht, Q1.arr+geo, Q2.arr+geo, masknp, std_file_msksub, _extra = ex
407;
408  cnt = cnt+1   &   blabla = 'Global Barotropic stream Function'
409  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_bsf, U1, U2, _extra = ex
410;
411  cnt = cnt+1   &   blabla = 'mean Temperature diff with New Reynolds'
412  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_sst, T1, T2, TRey, _extra = ex
413;
414  cnt = cnt+1   &   blabla = 'mean Salinity diff with Levitus'
415  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_sss, S1, S2, SLev, _extra = ex
416;
417  cnt = cnt+1   &   blabla = 'Arctic mean Salinity diff with Levitus'
418  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_ArcSal, S1, SLev, _extra = ex       
419;
420  cnt = cnt+1   &   blabla = 'Arctic mean Salinity diff with Levitus and exp2'
421  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_ArcSal, S1, S2, SLev, _extra = ex
422;
423  cnt = cnt+1   &   blabla = 'Arctic mean Salinity diff with Levitus at z=100 meters'
424  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_ArcSal, S1, SLev, /z100, _extra = ex       
425;
426  cnt = cnt+1   &   blabla = 'Arctic mean Salinity diff with Levitus and exp2 at z=100 meters'
427  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_ArcSal, S1, S2, SLev, /z100, _extra = ex
428;
429  cnt = cnt+1   &   blabla = 'mean Temperature diff with Levitus at z=100 meters'
430  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_T100m, T1, T2, Tlev, _extra = ex
431; 
432  cnt = cnt+1   &   blabla = 'mean Salinity diff with Levitus at z=100 meters'
433  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_S100m, S1, S2, SLev, _extra = ex
434; 
435  cnt = cnt+1   &   blabla = 'Mixed layer depth'
436  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_mld, MLD1, MLD, _extra = ex
437;
438  cnt = cnt+1   &   blabla = 'Mixed layer depth differences'
439  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_mld, MLD1, MLD2, MLD, _extra = ex
440;
441  cnt = cnt+1   &   blabla = 'Zonal mean Mixed layer depth'
442  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_ZonMld, MLD1, MLD2, MLD, _extra = ex
443; 
444  cnt = cnt+1   &   blabla = 'Vertical Global mean T & S'
445  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_GlobMeanTS, T1, T2, TLev, S1, S2, SLev, _extra = ex
446;
447  cnt = cnt+1   &   blabla = 'Zonal mean Temperature diff with Levitus: Global'
448  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_zonal_T, T1, T2, TLev, _extra = ex
449;
450  cnt = cnt+1   &   blabla = 'Zonal mean Temperature diff with Levitus: Atlantic'
451  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_zonal_T, T1, T2, TLev, SUBBASIN = 'Atl', _extra = ex
452; 
453  cnt = cnt+1   &   blabla = 'Zonal mean Temperature diff with Levitus: Indian'
454  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_zonal_T, T1, T2, TLev, SUBBASIN = 'Ind', _extra = ex
455; 
456  cnt = cnt+1   &   blabla = 'Zonal mean Temperature diff with Levitus: Pacific'
457  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_zonal_T, T1, T2, TLev, SUBBASIN = 'Pac', _extra = ex
458;
459  cnt = cnt+1   &   blabla = 'Zonal mean Salinity diff with Levitus: Global'
460  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_zonal_S, S1, S2, SLev, _extra = ex
461;
462  cnt = cnt+1   &   blabla = 'Zonal mean Salinity diff with Levitus: Atlantic'
463  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_zonal_S, S1, S2, SLev, SUBBASIN = 'Atl', _extra = ex
464;
465  cnt = cnt+1   &   blabla = 'Zonal mean Salinity diff with Levitus: Indian'
466  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_zonal_S, S1, S2, SLev, SUBBASIN = 'Ind', _extra = ex
467;
468  cnt = cnt+1   &   blabla = 'Zonal mean Salinity diff with Levitus: Pacific'
469  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_zonal_S, S1, S2, SLev, SUBBASIN = 'Pac', _extra = ex
470; 
471  cnt = cnt+1   &   blabla = 'Meridional stream Function: Global (no Med)'
472  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_msf, V1, V2, SUBBASIN = 'GloNoMed', _extra = ex
473; 
474  cnt = cnt+1   &   blabla = 'Meridional stream Function: Atlantic'
475  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_msf, V1, V2, SUBBASIN = 'Atl', _extra = ex
476 ;
477  cnt = cnt+1   &   blabla = 'Meridional stream Function: Indian'
478  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_msf, V1, V2, SUBBASIN = 'Ind', _extra = ex
479 ;
480  cnt = cnt+1   &   blabla = 'Meridional stream Function: Indo-Pacific'
481  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_msf, V1, V2, SUBBASIN = 'IndoPac', _extra = ex
482;
483  cnt = cnt+1   &   blabla = 'Equatorial Temperature'
484  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_EqT, T1, T2, Tlev, _extra = ex
485;
486  cnt = cnt+1   &   blabla = 'Equatorial Salinity'
487  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_EqS, S1, S2, SLev, _extra = ex
488;
489  cnt = cnt+1   &   blabla = 'Equatorial zonal velocity'
490  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_EqU, U1, U2, _extra = ex
491;
492  cnt = cnt+1   &   blabla = 'Mediterranean salt tongue at depth=700'
493  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_Med_Sspread, S1, S2, SLev, 700, _extra = ex
494;
495  cnt = cnt+1   &   blabla = 'Mediterranean salt tongue at depth=1000'
496  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_Med_Sspread, S1, S2, SLev, 1000, _extra = ex
497;
498  cnt = cnt+1   &   blabla = 'Mediterranean water at lat=40¡ N'
499  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_Med_Sdepth, S1, S2, SLev, 40, _extra = ex
500;
501  cnt = cnt+1   &   blabla = 'Mediterranean water at lat=38¡ N'
502  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_Med_Sdepth, S1, S2, SLev, 38, _extra = ex
503;
504;;
505; all plot are done for ice fraction > 0.15%
506;;
507  cnt = cnt+1   &   blabla = 'Arctic Ice Fraction: MARCH'
508  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceFrac, Ifra_march_1, Ifra_march_2, /ARC, /MARCH, _extra = ex
509; 
510  cnt = cnt+1   &   blabla = 'Arctic Ice Fraction: SEPT'
511  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceFrac, Ifra_sept_1, Ifra_sept_2, /ARC, /SEPT, _extra = ex
512; 
513  cnt = cnt+1   &   blabla = 'Antarctic Ice Fraction: MARCH'
514  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceFrac, Ifra_febr_1, Ifra_febr_2, /ANT, /FEBR, _extra = ex
515;
516  cnt = cnt+1   &   blabla = 'Antarctic Ice Fraction: SEPT'
517  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceFrac, Ifra_sept_1, Ifra_sept_2, /ANT, /SEPT, _extra = ex
518;
519  cnt = cnt+1   &   blabla = 'Arctic Ice Thickness: JAN'
520  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceThick, Ithi_jan_1, Ithi_jan_2, Ifra_jan_1, Ifra_jan_2, /ARC, /JAN, _extra = ex
521;
522  cnt = cnt+1   &   blabla = 'Arctic Ice Thickness: APRIL'
523  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceThick, Ithi_april_1, Ithi_april_2, Ifra_april_1, Ifra_april_2, /ARC, /APRIL, _extra = ex
524; 
525  cnt = cnt+1   &   blabla = 'Antarctic Ice Thickness: APRIL'
526  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceThick, Ithi_april_1, Ithi_april_2, Ifra_april_1, Ifra_april_2, /ANT, /APRIL, _extra = ex
527;
528  cnt = cnt+1   &   blabla = 'Antarctic Ice Thickness: SEPT'
529  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceThick, Ithi_sept_1, Ithi_sept_2, Ifra_sept_1, Ifra_sept_2, /ANT, /SEPT, _extra = ex
530;
531  cnt = cnt+1   &   blabla = 'Arctic SNOW Thickness: MARCH'
532  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_SnowThick, Isnow_march_1, Isnow_march_2, Ifra_march_1, Ifra_march_2, /ARC, /MARCH, _extra = ex
533;
534  cnt = cnt+1   &   blabla = 'Arctic SNOW Thickness: SEPT'
535  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_SnowThick, Isnow_sept_1, Isnow_sept_2, Ifra_sept_1, Ifra_sept_2, /ARC, /SEPT, _extra = ex
536; 
537  cnt = cnt+1   &   blabla = 'Antarctic SNOW Thickness: MARCH'
538  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_SnowThick, Isnow_febr_1, Isnow_febr_2, Ifra_febr_1, Ifra_febr_2, /ANT, /FEBR, _extra = ex
539;
540  cnt = cnt+1   &   blabla = 'Antarctic SNOW Thickness: SEPT'
541  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_SnowThick, Isnow_sept_1, Isnow_sept_2, Ifra_sept_1, Ifra_sept_2, /ANT, /SEPT, _extra = ex
542;
543  cnt = cnt+1   &   blabla = 'Arctic Ice Salinity: MARCH'
544  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceSal, Isal_march_1, Isal_march_2, Ifra_march_1, Ifra_march_2, /ARC, /MARCH, _extra = ex
545;
546  cnt = cnt+1   &   blabla = 'Arctic Ice Salinity: SEPT'
547  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceSal, Isal_sept_1, Isal_sept_2, Ifra_sept_1, Ifra_sept_2, /ARC, /SEPT, _extra = ex
548; 
549  cnt = cnt+1   &   blabla = 'Antarctic Ice Salinity: FEBRUARY'
550  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceSal, Isal_febr_1, Isal_febr_2, Ifra_febr_1, Ifra_febr_2, /ANT, /FEBR, _extra = ex
551;
552  cnt = cnt+1   &   blabla = 'Antarctic Ice Salinity: SEPT'
553  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceSal, Isal_sept_1, Isal_sept_2, Ifra_sept_1, Ifra_sept_2, /ANT, /SEPT, _extra = ex
554;
555  cnt = cnt+1   &   blabla = 'Arctic Ice Velocity: FEBRUARY'
556  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceVel, IvelU_febr_1, IvelU_febr_2, IvelV_febr_1, IvelV_febr_2, Ivelo_febr_1, Ivelo_febr_2, /ARC, /FEBR, _extra = ex
557;
558  cnt = cnt+1   &   blabla = 'Arctic Ice Velocity: SEPT'
559  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceVel, IvelU_sept_1, IvelU_sept_2, IvelV_sept_1, IvelV_sept_2, Ivelo_sept_1, Ivelo_sept_2, /ARC,  /SEPT, _extra = ex
560;
561  cnt = cnt+1   &   blabla = 'Antartic Ice Velocity: MARCH'
562  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceVel, IvelU_march_1, IvelU_march_2, IvelV_march_1, IvelV_march_2, Ivelo_march_1, Ivelo_march_2, /ANT,  /MARCH, _extra = ex
563;
564  cnt = cnt+1   &   blabla = 'Antartic Ice Velocity: SEPT'
565  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceVel, IvelU_sept_1, IvelU_sept_2, IvelV_sept_1, IvelV_sept_2, Ivelo_sept_1, Ivelo_sept_2, /ANT,  /SEPT, _extra = ex
566;
567  cnt = cnt+1   &   blabla = 'Arctic Ice Age: MARCH'
568  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceAge, Iage_march_1, Iage_march_2, Ifra_march_1, Ifra_march_2, /ARC, /MARCH, _extra = ex
569; 
570  cnt = cnt+1   &   blabla = 'Arctic Ice Age: SEPT'
571  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceAge, Iage_sept_1, Iage_sept_2,  Ifra_sept_1, Ifra_sept_2, /ARC, /SEPT, _extra = ex
572; 
573  cnt = cnt+1   &   blabla = 'Antarctic Ice Age: FEBR'
574  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceAge, Iage_febr_1, Iage_febr_2, Ifra_febr_1, Ifra_febr_2, /ANT, /FEBR, _extra = ex
575;
576  cnt = cnt+1   &   blabla = 'Antarctic Ice Age: SEPT'
577  IF doplot EQ cnt OR doplot EQ 0 THEN std_plot_IceAge, Iage_sept_1, Iage_sept_2,  Ifra_sept_1, Ifra_sept_2, /ANT, /SEPT, _extra = ex
578;
579 
580  IF n_elements(htmltxt) GT 1 THEN putfile, psdir+'std_plot_html_body.txt', htmltxt[1:*]
581 
582  return
583END
Note: See TracBrowser for help on using the repository browser.