source: branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/ORCA2_LIM/IDL_scripts/std_ts_SNOW_Vol.pro @ 4320

Last change on this file since 4320 was 4320, checked in by flavoni, 7 years ago

update IDL_scripts for ORCA2_LIM and ORCA2_LIM3, add plot of max_mld output 5days, see ticket: #724

File size: 6.9 KB
Line 
1function read_arr2d, filename, varname, t1, t2
2;; function that read input file and return 2d array with monthly timecounter
3nyear = (t2-t1+1)/12
4arr2d = ncdf_lec(filename, VAR=varname)
5arr2d = arr2d[t1:t2]
6arr2d = reform(arr2d,12,nyear) ; put in 2D array
7arr2d = total(arr2d,2)/nyear ; total over 2th dimension (i.e.years)
8
9return, arr2d
10end
11
12;
13;
14pro std_ts_SNOW_Vol, masknp, s_iodir_data,  POSTSCRIPT = postscript, _extra = ex
15
16  compile_opt idl2, strictarrsubs
17 
18@common
19@std_common
20
21; get exp1 info
22  vICE1     = getenv('VAR1_SNOW')     &   prefix = getenv('V1SNOW_PREF')  &   suffix = getenv('V1SNOW_SUFF')
23
24;;;;SF ???? come faccio lo snow volume se ho solo snowthickness???
25;;;;SF ???? quale variabile in output nel file iodef? e moltiplicarla per snowthicness?
26
27
28
29  v1_Ithick = getenv('VAR1_Ithick')  &   prefix = getenv('V1It_PREF')   &   suffix = getenv('V1It_SUFF')
30; get exp2 info
31  vICE2     = getenv('VAR2_ICE')     &   prefix2 = getenv('V2ICE_PREF')   &   suffix2 = getenv('V2ICE_SUFF')
32  v2_Ithick = getenv('VAR2_Ithick')  &   prefix2 = getenv('V2It_PREF')    &   suffix2 = getenv('V2It_SUFF')
33
34; get ice climatology info
35;
36;; METTERE CLIMATOLOGIA VOLUME SE TROVATA
37;
38;  std_file_ice =  isafile(getenv('FILE_ICE'), title = 'ICE Extent Climatology', iodir = std_iodir_climato)
39;
40;  time_ice = ncdf_lec( std_file_ice, VAR='time' )
41;  time_ice = (time_ice - FLOOR(time_ice) ) * 12
42;  time_ice = round(time_ice) ; round to nearest integer
43;  t1 = where(time_ice eq 0)
44;  t1 = t1[0] ;  jannuary
45;  t2 = where(time_ice eq 11, count)
46;  t2 = t2[count-1] ; last day of december
47;  nyear = (t2-t1+1)/12
48;  vICE_Ithick_NH = read_arr2d(std_file_ice, getenv('VAR1_Ithick'), t1, t2 )
49;  vICE_Ithick_SH = read_arr2d(std_file_ice, getenv('VAR_ICE_EXT_SH'), t1, t2 )
50;
51; vICE_area_NH = read_arr2d(std_file_ice, getenv('VAR_ICE_area_NH'), t1, t2 )
52;  vICE_area_SH = read_arr2d(std_file_ice, getenv('VAR_ICE_area_SH'), t1, t2 )
53
54  cdti3 = string(cnt, format = '(i3.3)')
55  print, cdti3 + ') ' + blabla
56  filename = cdti3 + '_ts_ICE_'+prefix
57  if prefix NE prefix2 then filename = filename + '_'+prefix2
58  if KEYWORD_SET(postscript) then openps, filename+'.ps', portrait = 1
59  d1_d2 = '('+strtrim(date1, 1)+' - '+strtrim(date2, 1)+')'
60  iodir = std_iodir_data
61  ; ICE Area(=Surface) in NORTH Hemisphere
62  domdef, 0, jpi-1, 30, 90, /xindex
63  ICE_N = rseries_ncdf(vICE1, date1, date2, prefix, suffix, FIRSTONLY = 1 - allrec)
64  ICE_thick = rseries_ncdf(v1_Ithick, date1, date2, prefix, suffix, FIRSTONLY = 1 - allrec)
65  ; Volume = Area(=Surface) * Thickness
66  ICE_vol_N = (ICE_N.arr < 1.e10 ) * ( ICE_thick.arr < 1.e10)  ; limited mask value of 1.e20, because 1.e20 * 1.e20 = inf for idl
67  ICE_vol_N = grossemoyenne(ICE_vol_N, 'xy', /integration, mask2d = masknp)
68  if jpt mod 12 ne 0 then stop
69  nyr=jpt/12.
70  ICE_vol_N = reform(ICE_vol_N, 12, nyr)
71  ICE_vol_N = total(ICE_vol_N,2)/nyr
72 ; ICE_vol_N = {arr:ICE_vol_N * 1.e-3, unit : '10^3 Km^3'}
73  ICE_vol_N = {arr:ICE_vol_N * 1.e-9, unit : '10^9 Km^3'}
74  ;ICE Area(=Surface) in SOUTH Hemisphere
75  domdef, 0, jpi-1, -90, -30, /xindex
76  ICE_S = rseries_ncdf(vICE1, date1, date2, prefix, suffix, FIRSTONLY = 1 - allrec)
77  ICE_thick = rseries_ncdf(v1_Ithick, date1, date2, prefix, suffix, FIRSTONLY = 1 - allrec)
78  ; Volume = Area(=Surface) * Thickness
79  ICE_vol_S = (ICE_S.arr < 1.e10 ) * ( ICE_thick.arr < 1.e10)  ; limited mask value of 1.e20, because 1.e20 * 1.e20 = inf for idl
80  ICE_vol_S = grossemoyenne(ICE_vol_S, 'xy', /integration, mask2d = masknp)
81  if jpt mod 12 ne 0 then stop
82  nyr=jpt/12.
83  ICE_vol_S = reform(ICE_vol_S, 12, nyr)
84  ICE_vol_S = total(ICE_vol_S,2)/nyr
85;print, 'prima valore di ice vol s'
86;print, ICE_vol_S
87  ICE_vol_S = {arr:ICE_vol_S * 1.e-9, unit : '10^9 Km^3'}
88;print, 'DOPO valore di ice vol s'
89;print, ICE_vol_S
90  ;
91  title = 'Northern Hemisphere'+'!C'+prefix+' '+d1_d2+'!C'+'Global Annual Mean Ice Volume (Black SOLID simulation) + Observation (dashed)'
92  jpt=12
93  time=julday(1,15,1900)+30*lindgen(12)
94  pltt, ICE_vol_N, 't', 0., 15., 19000101, 19001231, /REMPLI, /PORTRAIT,MIN = 0., MAX = 30000. $
95        , small = [1, 2, 1], YTITLE = '10^9 Km^3 ', TITLE = title, _extra = ex
96;
97  title ='Southern Hemisphere' +'!C'+prefix+' '+d1_d2+' - '+'!C'+'Global Annual Mean Ice Volume (Black SOLID simulation) + Observation (dashed)'
98  pltt, ICE_vol_S, 't', 0., 15.,19000101 ,19001231 , /REMPLI, /NOERASE , MIN = 0., MAX = 11000. $
99        , small = [1, 2, 2], YTITLE = '10^9 Km^3 ', TITLE = title, _extra = ex
100;
101
102;;return
103  htmltxt = [ htmltxt, '<hr>'+blabla, '<br><img width="80%" src='+filename+'.png  />  ' ]
104  if KEYWORD_SET(postscript) then closeps
105
106  if prefix NE prefix2 then BEGIN
107
108    d1_d2_2 = '('+strtrim(date1_2, 1)+' - '+strtrim(date2_2, 1)+')'
109    tsave = time
110    domdef, 0, jpi-1, 30, 90, /xindex
111    ;ICE Extent in NORTH Hemisphere
112    ICE_N2 = rseries_ncdf(vICE2, date1, date2, prefix, suffix, FIRSTONLY = 1 - allrec)
113    ICE_thick2 = rseries_ncdf(v2_Ithick, date1, date2, prefix, suffix, FIRSTONLY = 1 - allrec)
114    ; Volume = Area(=Surface) * Thickness
115    ICE_vol_N2 = (ICE_N2.arr < 1.e10 ) * ( ICE_thick2.arr < 1.e10) ; limited mask value of 1.e20, because 1.e20 * 1.e20 = inf for idl
116    ICE_vol_N2 = grossemoyenne(ICE_vol_N2, 'xy', /integration, mask2d = masknp)
117    if jpt mod 12 ne 0 then stop
118    nyr=jpt/12.
119    ICE_vol_N2 = reform(ICE_vol_N2, 12, nyr)
120    ICE_vol_N2 = total(ICE_vol_N2,2)/nyr
121    ICE_vol_N2 = {arr:ICE_vol_N2 * 1.e-9, unit : '10^3 Km^3'}
122
123    ;ICE Area(=Surface) in SOUTH Hemisphere
124    domdef, 0, jpi-1, -90, -30, /xindex
125    ICE_S2 = rseries_ncdf(vICE2, date1, date2, prefix, suffix, FIRSTONLY = 1 - allrec)
126    ICE_thick2 = rseries_ncdf(v2_Ithick, date1, date2, prefix, suffix, FIRSTONLY = 1 - allrec)
127    ; Volume = Area(=Surface) * Thickness
128    ICE_vol_S2 = (ICE_S2.arr < 1.e10 ) * ( ICE_thick.arr < 1.e10) ; limited mask value of 1.e20, because 1.e20 * 1.e20 = inf for idl
129    ICE_vol_S2 = grossemoyenne(ICE_vol_S2, 'xy', /integration, mask2d = masknp)
130    if jpt mod 12 ne 0 then stop
131    nyr=jpt/12.
132    ICE_vol_S2 = reform(ICE_vol_S2, 12, nyr)
133    ICE_vol_S2 = total(ICE_vol_S2,2)/nyr
134    ICE_vol_S2 = {arr:ICE_vol_S2 * 1.e-9, unit : '10^3 Km^3'}
135
136    time = tsave   &   IF n_elements(time) NE jpt THEN stop
137
138    if KEYWORD_SET(postscript) then openps, filename+'.ps', portrait = 1
139
140    title = prefix+' '+d1_d2+' - '+prefix2+' '+d1_d2_2+'!C'+'Global Annual Mean Ice Volume (North. Hemisp.)'
141    pltt, ICE_vol_N.arr - ICE_vol_N2.arr, 't', -2., 2., date1, date2, /REMPLI, /PORTRAIT, window = 2 $
142          , COLOR = 250, small = [1, 2, 1], YTITLE = '10^9 Km^3 ', TITLE = title, _extra = ex
143   ;
144    title = prefix+' '+d1_d2+' - '+prefix2+' '+d1_d2_2+'!C'+'Global Annual Mean Ice Area (South. Hemisp.)'
145    pltt, ICE_vol_S.arr - ICE_vol_S2.arr, 't', -2., 2., date1, date2, /REMPLI, /NOERASE $
146          , COLOR = 250, small = [1, 2, 2], YTITLE = '10^9 Km^3 ', TITLE = title, _extra = ex   
147
148    htmltxt = [ htmltxt, '<hr>'+blabla, '<br><img width="80%" src='+filename+'.png  />  ' ]
149    if KEYWORD_SET(postscript) then closeps
150
151  endif
152
153  domdef
154
155  return
156end
Note: See TracBrowser for help on using the repository browser.