source: TOOLS/CMIP6_FORCING/AER_STRAT/volc.sh @ 6048

Last change on this file since 6048 was 6048, checked in by oboucher, 3 years ago

new path for input data

  • Property svn:executable set to *
File size: 19.3 KB
Line 
1dirpwd=`pwd`
2
3#original CMIP6 files
4vv='v3'
5#corrected CMIP6 files + extended in time
6vv='v4'
7
8#--choose the resolution
9#--only LR and zoom works fully for now !
10#lmdz='VVLR_l19'
11lmdz='VLR_l39'
12#lmdz='VLR_l79'
13#lmdz='LR_l79'
14#lmdz='MR_l79'
15#--added for Frederique Cheruy
16#lmdz='zoom_128x89_l79'
17
18dirout='/data/'${USER}'/CMIP6/VOLC/'${lmdz}'_'${vv}'/'
19echo $dirout
20if [ ! -d ${dirout} ] ; then mkdir -p ${dirout} ; fi
21
22cat > process_volc_${lmdz}.pro << EOF
23pro regrid
24;--version sans longitude
25;--latitude - level - wavelength - time
26;--sep 2017: add lat and month z variations of model layers
27;--Olivier Boucher
28;
29lmdz='${lmdz}'
30;
31if (lmdz eq 'VVLR_l19') then begin
32dimlonlmdz=48
33dimlatlmdz=37
34dimz=19
35output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
36latfirst=90.
37latinc=-180./float(dimlatlmdz-1)
38latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
39endif
40
41if (lmdz eq 'VLR_l39') then begin
42dimlonlmdz=96
43dimlatlmdz=96
44dimz=39
45output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
46latfirst=90.
47latinc=-180./float(dimlatlmdz-1)
48latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
49endif
50
51if (lmdz eq 'VLR_l79') then begin
52dimlonlmdz=96
53dimlatlmdz=96
54dimz=79
55output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
56latfirst=90.
57latinc=-180./float(dimlatlmdz-1)
58latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
59endif
60
61if (lmdz eq 'LR_l79') then begin
62dimlonlmdz=144
63dimlatlmdz=143
64dimz=79
65output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
66latfirst=90.
67latinc=-180./float(dimlatlmdz-1)
68latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
69endif
70
71if (lmdz eq 'MR_l79') then begin
72dimlonlmdz=256
73dimlatlmdz=257
74dimz=79
75output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
76latfirst=90.
77latinc=-180./float(dimlatlmdz-1)
78latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
79endif
80
81if (lmdz eq 'zoom_128x89_l79') then begin
82dimlonlmdz=128
83dimlatlmdz=89
84dimz=79
85output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
86latitudelmdz=[90., 86.85495, 83.8430634, 80.9792252, 78.2779922, 75.7526321, 73.4137878, 71.2680206, 69.3164978, $
87              67.5542145,   65.9702301, 64.5487747, 63.2709846, 62.1167107, 61.0660896, 60.1005592, 59.2034454, $
88              58.360157,    57.5581551, 56.7868576, 56.0374718, 55.302887,  54.5775719, 53.8574829, 53.1399384, $
89              52.4233971,   51.7071419, 50.9909439, 50.2747459, 49.5585518, 48.8423538, 48.1261597, 47.4099655, $
90              46.6937675,   45.9775734, 45.2613792, 44.5451813, 43.8289604, 43.1125679, 42.395504,  41.6765518, $
91              40.9534645,   40.2226982, 39.4792557, 38.7165909, 37.9265213, 37.099102,  36.2225075, 35.2828789, $
92              34.2642632,   33.1487236, 31.9167614, 30.5481758, 29.0234509, 27.32551,   25.4415455, 23.3643932, $
93              21.0930786,   18.6323872, 15.9917059, 13.1835403, 10.2221146, 7.12228632, 3.89880991, 0.565921485,$
94              -2.86286402, -6.37483072, -9.95808029, -13.6015415, -17.2949867, -21.029068, -24.7953758, -28.5864811, $
95             -32.3959961, -36.2186012, -40.0500298, -43.8870316, -47.7272758, -51.5692215, -55.4119453, -59.2549667, $
96             -63.0980797, -66.9412079, -70.7843399, -74.6274719, -78.4706039, -82.313736, -86.156868, -90.]
97endif
98
99NSW=6
100NLW=16
101;
102;--script only works on ciclad
103dir='/bdd/input4MIPs/VOLC/${vv}/'
104;;;not available for v4
105;;;filename=dir+'CMIP_1850_2014_extinction_550nm_${vv}.nc'
106;;;NETCDFREAD,filename,'altitude',altitude,dimaltitude
107;;;NETCDFREAD,filename,'ext550',ext550,dimext
108;
109filename=dir+'CMIP_IPSL-CM6_radiation_${vv}.nc'
110;
111NETCDFREAD,filename,'altitude',altitude,dimaltitude
112NETCDFREAD,filename,'latitude',latitude,dimlatitude
113NETCDFREAD,filename,'wl1_sun',wl1_sun,solar_bands
114NETCDFREAD,filename,'wl2_sun',wl2_sun,solar_bands
115NETCDFREAD,filename,'wl1_earth',wl1_earth,terrestrial_bands
116NETCDFREAD,filename,'wl2_earth',wl2_earth,terrestrial_bands
117NETCDFREAD,filename,'month',month,dimmonth
118NETCDFREAD,filename,'ext_sun',ext_sun,dimext
119NETCDFREAD,filename,'omega_sun',omega_sun,dimext
120NETCDFREAD,filename,'g_sun',g_sun,dimext
121NETCDFREAD,filename,'ext_earth',ext_earth,dimext
122NETCDFREAD,filename,'omega_earth',omega_earth,dimext
123NETCDFREAD,filename,'g_earth',g_earth,dimext
124;
125print ,'wl1_sun=', wl1_sun
126print ,'wl2_sun=', wl2_sun
127print ,'wl1_ear=', wl1_earth
128print ,'wl2_ear=', wl2_earth
129;
130print ,'min max ext_sun=', min(ext_sun),   max(ext_sun)
131print ,'min max ome_sun=', min(omega_sun), max(omega_sun)
132;for nl=0, NSW-1 do begin
133;print ,'min max ome_sun nl=', nl, min(omega_sun(*,*,*,nl)), max(omega_sun(*,*,*,nl))
134;endfor
135print ,'min max ggg_sun=', min(g_sun),     max(g_sun)
136print ,'min max ext_ear=', min(ext_earth), max(ext_earth)
137;
138dimlat=dimlatitude(0)
139dimalt=dimaltitude(0)
140dimtime=dimmonth(0)
141;
142;--vertical resolution in Beiping's code
143dz=0.5
144;
145month_in_year=12
146;;month_in_year=1
147;
148;--compute optical depth of input data
149;--dimension ext_sun(l,k,j,nl)
150;ext_sun0=ext_sun(0,*,*,*)
151;print , 'size ext_sun0=', size(ext_sun0)
152;tau_sun=TOTAL(ext_sun0,2)*dz
153;print ,'min max tau_sun l=0 =', min(tau_sun),   max(tau_sun)
154;
155zz=fltarr(dimz)
156;
157;---exact altitudes of LMDZ -- L79
158;--only exists for VLR, LR, MR and zoom_128x89 at the moment
159;
160filename='/bdd/input4MIPs/ZALT/zalt_zonmean_${lmdz}_rev.nc'
161NETCDFREAD,filename,'GEOP',zz,dimzz
162NETCDFREAD,filename,'LAT',zzlat,dimzzlat
163NETCDFREAD,filename,'TIME_COUNTER',zztime,dimzztime
164dimzzilat=dimzzlat[0]
165dimzzitime=dimzztime[0]
166;--becareful zz comes with four dimensions
167;--lon lat k time
168print, 'GEOP size=', size(zz)
169;--becareful zzlat comes with South Pole first
170print, 'LAT from zalt field=', zzlat
171if (dimzzilat ne dimlatlmdz) then begin
172  print , 'PB dimension latitude'
173endif
174;
175;
176;--reconstructing the vertical coordinate at interfaces (in unit km)
177;--and reverse lat axis for zzi  i==>ii
178;--and reverse zz axis for zzi   dimz-1-k==>k
179;--and forget about lon axis : index 0 in zz
180zzi=fltarr(dimzzilat,dimzzitime,dimz+1)
181for i=0, dimzzilat-1 do begin
182 ii=dimzzilat-1-i
183 for t=0, dimzzitime-1 do begin
184  zzi(ii,t,0)=0.0
185  for k=1, dimz-1 do begin
186   zzi(ii,t,k)=(zz(0,i,dimz-1-(k-1),t)+zz(0,i,dimz-1-k,t))/2.0
187  endfor
188  zzi(ii,t,dimz)=100.00
189 endfor
190endfor
191print, 'zzi=',zzi
192;
193lev=indgen(dimz)+1
194;
195dimzori=dimalt
196zziori=fltarr(dimzori+1)
197zziori(0)=altitude(0)-0.25
198for k=1, dimzori do begin
199 zziori(k)=altitude(k-1)+0.25
200endfor
201;
202;--550 nm properties
203;;;not available for v4
204;;;tau_550_lmdz=fltarr(dimlatlmdz,dimz,month_in_year)
205;
206;--SW properties
207tau_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
208ome_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
209ggg_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
210;
211tau_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
212ome_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
213ggg_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
214;
215tau_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
216ome_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
217ggg_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
218;
219;--LW properties tau_abs
220tau_ear_lmdz=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
221;
222tau_ear_lmdz_ave=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
223;
224tau_ear_lmdz_tr=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
225;
226tau_sun_lmdz_ave(*,*,*,0)=1.e-15
227ome_sun_lmdz_ave(*,*,*,0)=1.e-15
228ggg_sun_lmdz_ave(*,*,*,0)=0.0
229tau_ear_lmdz_ave(*,*,*,0)=1.e-15
230;
231yearend=dimtime/month_in_year-1
232;
233for year=0, yearend do begin
234;
235;;;not available for v4
236;;;tau_550_lmdz(*,*,*)=1.e-15
237tau_sun_lmdz(*,*,*,*)=1.e-15
238ome_sun_lmdz(*,*,*,*)=1.e-15
239ggg_sun_lmdz(*,*,*,*)=0.0
240tau_ear_lmdz(*,*,*,*)=1.e-15
241;
242chyr=strcompress(1850+year,/rem)
243;
244for mth=0,month_in_year-1 do begin
245;
246;--timestep
247l=mth+month_in_year*year
248print,'year mth l=',chyr, mth, l
249;
250;regridding
251for j=0, dimlatlmdz-1  do begin
252;
253;--finding latitude in beiping luo data
254jj=0
255for jluo=0,dimlat-2 do begin
256  if (latitudelmdz(j) gt (latitude(jluo)+latitude(jluo+1))/2. ) then begin
257    jj=jluo+1
258  endif
259endfor
260;
261for k=0, dimz-1 do begin
262for kori=0, dimzori-1 do begin
263;
264;fraction de la maille kori qui se trouve dans la maille k
265frac= max([0.0,min([zzi(j,mth,k+1),zziori(kori+1)])-max([zzi(j,mth,k),zziori(kori)])])/(zziori(kori+1)-zziori(kori))
266;
267;;;not available for v4
268;;;tau_550_lmdz(j,k,mth)=tau_550_lmdz(j,k,mth)+ext550(l,kori,jj)*dz*frac
269;
270for nl=0, NSW-1 do begin
271  tau_sun_lmdz(j,k,nl,mth)=tau_sun_lmdz(j,k,nl,mth)+ext_sun(l,kori,jj,nl)*dz*frac
272  ome_sun_lmdz(j,k,nl,mth)=ome_sun_lmdz(j,k,nl,mth)+omega_sun(l,kori,jj,nl)*ext_sun(l,kori,jj,nl)*dz*frac
273  ggg_sun_lmdz(j,k,nl,mth)=ggg_sun_lmdz(j,k,nl,mth)+g_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*ext_sun(l,kori,jj,nl)*dz*frac
274;
275  tau_sun_lmdz_ave(j,k,nl,0)=tau_sun_lmdz_ave(j,k,nl,0)+ext_sun(l,kori,jj,nl)*dz*frac
276  ome_sun_lmdz_ave(j,k,nl,0)=ome_sun_lmdz_ave(j,k,nl,0)+omega_sun(l,kori,jj,nl)*ext_sun(l,kori,jj,nl)*dz*frac
277  ggg_sun_lmdz_ave(j,k,nl,0)=ggg_sun_lmdz_ave(j,k,nl,0)+g_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*ext_sun(l,kori,jj,nl)*dz*frac
278endfor
279;
280for nl=0, NLW-1 do begin
281  tau_ear_lmdz(j,k,nl,mth)=tau_ear_lmdz(j,k,nl,mth)+ext_earth(l,kori,jj,nl)*(1.-omega_earth(l,kori,jj,nl))*dz*frac
282  tau_ear_lmdz_ave(j,k,nl,0)=tau_ear_lmdz_ave(j,k,nl,0)+ext_earth(l,kori,jj,nl)*(1.-omega_earth(l,kori,jj,nl))*dz*frac
283endfor
284;
285endfor
286;--end lat loop
287;
288endfor
289endfor
290;--end k loops
291;
292endfor
293;--end month loop
294;
295;renormalizing intensive SW properties
296ggg_sun_lmdz(*,*,*,*)=ggg_sun_lmdz(*,*,*,*)/ome_sun_lmdz(*,*,*,*)
297ome_sun_lmdz(*,*,*,*)=ome_sun_lmdz(*,*,*,*)/tau_sun_lmdz(*,*,*,*)
298;
299;saving netcdf file
300;
301print ,'min max ext_sun lmdz=', min(tau_sun_lmdz), max(tau_sun_lmdz)
302print ,'min max ome_sun lmdz=', min(ome_sun_lmdz), max(ome_sun_lmdz)
303;for nl=0, NSW-1 do begin
304;print ,'min max ome_sun lmdz nl =', nl, min(ome_sun_lmdz(*,*,nl,*)), max(ome_sun_lmdz(*,*,nl,*))
305;endfor
306print ,'min max ggg_sun lmdz=', min(ggg_sun_lmdz), max(ggg_sun_lmdz)
307print ,'min max ext_ear lmdz=', min(tau_ear_lmdz), max(tau_ear_lmdz)
308;
309print ,'min max ext_sun lmdz l=1=', min(tau_sun_lmdz(*,*,*,0)), max(tau_sun_lmdz(*,*,*,0))
310print ,'min max ome_sun lmdz l=1=', min(ome_sun_lmdz(*,*,*,0)), max(ome_sun_lmdz(*,*,*,0))
311print ,'min max ggg_sun lmdz l=1=', min(ggg_sun_lmdz(*,*,*,0)), max(ggg_sun_lmdz(*,*,*,0))
312print ,'min max ext_ear lmdz l=1=', min(tau_ear_lmdz(*,*,*,0)), max(tau_ear_lmdz(*,*,*,0))
313;
314;--compute optical depth of output data
315;tau_sun=TOTAL(tau_sun_lmdz,2)
316;print ,'min max tau_sun_lmdz vert=', min(tau_sun),   max(tau_sun)
317;
318opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $
319             wav:fltarr(NSW),time:fltarr(month_in_year),          $
320             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
321             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
322             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
323;
324opticstruct.lat=latitudelmdz
325opticstruct.lev=lev
326opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
327opticstruct.time=float(indgen(month_in_year)+1)
328opticstruct.tau_sun=tau_sun_lmdz
329opticstruct.ome_sun=ome_sun_lmdz
330opticstruct.ggg_sun=ggg_sun_lmdz
331;
332attributes = {units:strarr(7),long_name:strarr(7)}
333attributes.units =     ['degrees_north','level','meters','month','-','-','-']
334attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
335;
336dimensions = {isdim:intarr(7), links:intarr(4,7)}
337       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
338       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
339                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
340                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
341;
342netcdfwrite,output+'tauswstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
343            attributes=attributes, dimensions=dimensions
344;
345opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $
346             wav:fltarr(NLW),time:fltarr(month_in_year),          $
347             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }
348;
349opticstruct.lat=latitudelmdz
350opticstruct.lev=lev
351opticstruct.wav=(wl1_earth+wl2_earth)/2.
352opticstruct.time=float(indgen(month_in_year)+1)
353opticstruct.tau_ear=tau_ear_lmdz
354;
355attributes = {units:strarr(5),long_name:strarr(5)}
356attributes.units = ['degrees_north','level','cm-1','month','-']
357attributes.long_name = ['latitude','level','wavenumber','time','tau_ear']
358;
359dimensions = {isdim:intarr(5), links:intarr(4,5)}
360       dimensions.isdim =  [1,1,1,1,0]  ; (1=dimension, 0=variable)
361       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
362                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
363                            [0,1,2,3]  ]
364;
365netcdfwrite,output+'taulwstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
366            attributes=attributes, dimensions=dimensions
367;
368;;;not available in v4
369;;;opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $
370;;;             time:fltarr(month_in_year),                          $
371;;;             tau550:fltarr(dimlatlmdz,dimz,month_in_year) }
372;
373;;;opticstruct.lat=latitudelmdz
374;;;opticstruct.lev=lev
375;;;opticstruct.time=float(indgen(month_in_year)+1)
376;;;opticstruct.tau550=tau_550_lmdz
377;
378;;;attributes = {units:strarr(4),long_name:strarr(4)}
379;;;attributes.units = ['degrees_north','level','month','-']
380;;;attributes.long_name = ['latitude','level','time','tau550']
381;
382;;;dimensions = {isdim:intarr(4), links:intarr(3,4)}
383;;;       dimensions.isdim =  [1,1,1,0]  ; (1=dimension, 0=variable)
384;;;       dimensions.links = [ [-1,-1,-1],[-1,-1,-1],   $
385;;;                            [-1,-1,-1],[0,1,2]  ]
386;
387;;;netcdfwrite,output+'tau550strat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
388;;;            attributes=attributes, dimensions=dimensions
389;
390endfor
391;--end loop on years
392;
393;now deal with average conditions
394ggg_sun_lmdz_ave(*,*,*,0)=ggg_sun_lmdz_ave(*,*,*,0)/ome_sun_lmdz_ave(*,*,*,0)
395ome_sun_lmdz_ave(*,*,*,0)=ome_sun_lmdz_ave(*,*,*,0)/tau_sun_lmdz_ave(*,*,*,0)
396tau_sun_lmdz_ave(*,*,*,0)=tau_sun_lmdz_ave(*,*,*,0)/float(dimtime)
397tau_ear_lmdz_ave(*,*,*,0)=tau_ear_lmdz_ave(*,*,*,0)/float(dimtime)
398;
399for mth=1, month_in_year-1 do begin
400  ggg_sun_lmdz_ave(*,*,*,mth)=ggg_sun_lmdz_ave(*,*,*,0)
401  ome_sun_lmdz_ave(*,*,*,mth)=ome_sun_lmdz_ave(*,*,*,0)
402  tau_sun_lmdz_ave(*,*,*,mth)=tau_sun_lmdz_ave(*,*,*,0)
403  tau_ear_lmdz_ave(*,*,*,mth)=tau_ear_lmdz_ave(*,*,*,0)
404endfor
405;
406opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $
407             wav:fltarr(NSW),time:fltarr(month_in_year),          $
408             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
409             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
410             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
411;
412opticstruct.lat=latitudelmdz
413opticstruct.lev=lev
414opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
415opticstruct.time=float(indgen(month_in_year))
416opticstruct.tau_sun=tau_sun_lmdz_ave
417opticstruct.ome_sun=ome_sun_lmdz_ave
418opticstruct.ggg_sun=ggg_sun_lmdz_ave
419;
420attributes = {units:strarr(7),long_name:strarr(7)}
421attributes.units =     ['degrees_north','level','meters','month','-','-','-']
422attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
423;
424dimensions = {isdim:intarr(7), links:intarr(4,7)}
425       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
426       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
427                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
428                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
429;
430netcdfwrite,output+'tauswstrat.2D.ave.nc',opticstruct,clobber=1, $
431            attributes=attributes, dimensions=dimensions
432;
433opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $
434             wav:fltarr(NLW),time:fltarr(month_in_year),          $
435             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }
436;
437opticstruct.lat=latitudelmdz
438opticstruct.lev=lev
439opticstruct.wav=(wl1_earth+wl2_earth)/2.
440opticstruct.time=float(indgen(month_in_year))
441opticstruct.tau_ear=tau_ear_lmdz_ave
442;
443attributes = {units:strarr(5),long_name:strarr(5)}
444attributes.units = ['degrees_north','level','cm-1','month','-']
445attributes.long_name = ['latitude','level','wavenumber','time','tau_ear']
446;
447dimensions = {isdim:intarr(5), links:intarr(4,5)}
448       dimensions.isdim =  [1,1,1,1,0]  ; (1=dimension, 0=variable)
449       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
450                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
451                            [0,1,2,3]  ]
452;
453netcdfwrite,output+'taulwstrat.2D.ave.nc',opticstruct,clobber=1, $
454            attributes=attributes, dimensions=dimensions
455;
456;--finally prepare the 10-year transition period
457;--that is 2015-2023 for v3 or 2017-2025 for v4
458;--we mix the last year of the period with the climatology
459;
460for year=1850+yearend+1,1850+yearend+9 do begin
461;
462chyr=strcompress(year,/rem)
463print,'year =',chyr
464;--compute weights
465w1=float(1850+yearend+10-year)/10.
466w2=1.0-w1
467;--SW properties
468tau_sun_lmdz_tr=w1*tau_sun_lmdz+w2*tau_sun_lmdz_ave
469ome_sun_lmdz_tr=w1*tau_sun_lmdz*ome_sun_lmdz+w2*tau_sun_lmdz_ave*ome_sun_lmdz_ave
470ggg_sun_lmdz_tr=w1*tau_sun_lmdz*ome_sun_lmdz*ggg_sun_lmdz+ $
471                w2*tau_sun_lmdz_ave*ome_sun_lmdz_ave*ggg_sun_lmdz_ave
472ggg_sun_lmdz_tr=ggg_sun_lmdz_tr/ome_sun_lmdz_tr
473ome_sun_lmdz_tr=ome_sun_lmdz_tr/tau_sun_lmdz_tr
474;
475;--LW properties
476tau_ear_lmdz_tr=w1*tau_ear_lmdz+w2*tau_ear_lmdz_ave
477;
478;--prepare SW output
479opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $
480             wav:fltarr(NSW),time:fltarr(month_in_year),          $
481             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
482             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
483             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
484;
485opticstruct.lat=latitudelmdz
486opticstruct.lev=lev
487opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
488opticstruct.time=float(indgen(month_in_year)+1)
489opticstruct.tau_sun=tau_sun_lmdz_tr
490opticstruct.ome_sun=ome_sun_lmdz_tr
491opticstruct.ggg_sun=ggg_sun_lmdz_tr
492;
493attributes = {units:strarr(7),long_name:strarr(7)}
494attributes.units =     ['degrees_north','level','meters','month','-','-','-']
495attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
496;
497dimensions = {isdim:intarr(7), links:intarr(4,7)}
498       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
499       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
500                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
501                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
502;
503netcdfwrite,output+'tauswstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
504            attributes=attributes, dimensions=dimensions
505;
506;--prepare LW output
507opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $
508             wav:fltarr(NLW),time:fltarr(month_in_year),          $
509             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }
510;
511opticstruct.lat=latitudelmdz
512opticstruct.lev=lev
513opticstruct.wav=(wl1_earth+wl2_earth)/2.
514opticstruct.time=float(indgen(month_in_year)+1)
515opticstruct.tau_ear=tau_ear_lmdz_tr
516;
517attributes = {units:strarr(5),long_name:strarr(5)}
518attributes.units = ['degrees_north','level','cm-1','month','-']
519attributes.long_name = ['latitude','level','wavenumber','time','tau_ear']
520;
521dimensions = {isdim:intarr(5), links:intarr(4,5)}
522       dimensions.isdim =  [1,1,1,1,0]  ; (1=dimension, 0=variable)
523       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
524                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
525                            [0,1,2,3]  ]
526;
527netcdfwrite,output+'taulwstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
528            attributes=attributes, dimensions=dimensions
529;
530endfor
531;--end loop on years
532;
533end
534EOF
535
536cat > volc_${lmdz}.job << EOF2
537#PBS -N process_volc_${lmdz}
538#PBS -S /bin/bash
539#PBS -q week # there exist: short, day, days3, week...
540#PBS -k eo  # to write the output of stdin
541### Max memory
542#PBS -l vmem=10gb   # virtual memory
543#PBS -l  mem=10gb
544
545cd $dirpwd
546
547idl << eof
548.r netcdf
549.r process_volc_${lmdz}
550regrid
551eof
552EOF2
553
554qsub volc_${lmdz}.job
Note: See TracBrowser for help on using the repository browser.