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

dirpwd=`pwd`

vv='v4'

#--choose the resolution
#--only LR and zoom works fully for now !
#lmdz='VVLR'
#lmdz='VLR_L79'
#lmdz='VLR'
lmdz='LR'
#lmdz='MR'
#--added for Frederique Cheruy
#lmdz='zoom_128x89'

dirout='/data/'${USER}'/CMIP6/VOLC/'${lmdz}'_'${vv}'/'
echo $dirout
if [ ! -d ${dirout} ] ; then mkdir -p ${dirout} ; fi

cat > process_volc_${lmdz}.pro << EOF
pro regrid
;--version sans longitude
;--latitude - level - wavelength - time
;--sep 2017: add lat and month z variations of model layers
;--Olivier Boucher
;
lmdz='${lmdz}'
;
if (lmdz eq 'VVLR') then begin
dimlonlmdz=48
dimlatlmdz=37
dimz=19
output='/data/${USER}/CMIP6/VOLC/VVLR_${vv}/'
latfirst=90.
latinc=-180./float(dimlatlmdz-1)
latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
endif

if (lmdz eq 'VLR') then begin
dimlonlmdz=96
dimlatlmdz=96
dimz=39
output='/data/${USER}/CMIP6/VOLC/VLR_${vv}/'
latfirst=90.
latinc=-180./float(dimlatlmdz-1)
latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
endif

if (lmdz eq 'VLR_L79') then begin
dimlonlmdz=96
dimlatlmdz=96
dimz=79
output='/data/${USER}/CMIP6/VOLC/VLR_L79_${vv}/'
latfirst=90.
latinc=-180./float(dimlatlmdz-1)
latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
endif

if (lmdz eq 'LR') then begin
dimlonlmdz=144
dimlatlmdz=143
dimz=79
output='/data/${USER}/CMIP6/VOLC/LR_${vv}/'
latfirst=90.
latinc=-180./float(dimlatlmdz-1)
latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
endif

if (lmdz eq 'MR') then begin
dimlonlmdz=256
dimlatlmdz=257
dimz=79
output='/data/${USER}/CMIP6/VOLC/MR_${vv}/'
latfirst=90.
latinc=-180./float(dimlatlmdz-1)
latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
endif

if (lmdz eq 'zoom_128x89') then begin
dimlonlmdz=128
dimlatlmdz=89
dimz=79
output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
latitudelmdz=[90., 86.85495, 83.8430634, 80.9792252, 78.2779922, 75.7526321, 73.4137878, 71.2680206, 69.3164978, $
              67.5542145,   65.9702301, 64.5487747, 63.2709846, 62.1167107, 61.0660896, 60.1005592, 59.2034454, $ 
              58.360157,    57.5581551, 56.7868576, 56.0374718, 55.302887,  54.5775719, 53.8574829, 53.1399384, $
              52.4233971,   51.7071419, 50.9909439, 50.2747459, 49.5585518, 48.8423538, 48.1261597, 47.4099655, $
              46.6937675,   45.9775734, 45.2613792, 44.5451813, 43.8289604, 43.1125679, 42.395504,  41.6765518, $ 
              40.9534645,   40.2226982, 39.4792557, 38.7165909, 37.9265213, 37.099102,  36.2225075, 35.2828789, $ 
              34.2642632,   33.1487236, 31.9167614, 30.5481758, 29.0234509, 27.32551,   25.4415455, 23.3643932, $ 
              21.0930786,   18.6323872, 15.9917059, 13.1835403, 10.2221146, 7.12228632, 3.89880991, 0.565921485,$
              -2.86286402, -6.37483072, -9.95808029, -13.6015415, -17.2949867, -21.029068, -24.7953758, -28.5864811, $ 
             -32.3959961, -36.2186012, -40.0500298, -43.8870316, -47.7272758, -51.5692215, -55.4119453, -59.2549667, $
             -63.0980797, -66.9412079, -70.7843399, -74.6274719, -78.4706039, -82.313736, -86.156868, -90.]
endif

NSW=6
NLW=16
;
;--script only works on ciclad
dir='/prodigfs/project/input4MIPs/VOLC/${vv}/'
;;;not available for v4
;;;filename=dir+'CMIP_1850_2014_extinction_550nm_${vv}.nc'
;;;NETCDFREAD,filename,'altitude',altitude,dimaltitude
;;;NETCDFREAD,filename,'ext550',ext550,dimext
;
filename=dir+'CMIP_IPSL-CM6_radiation_${vv}.nc'
;
NETCDFREAD,filename,'altitude',altitude,dimaltitude
NETCDFREAD,filename,'latitude',latitude,dimlatitude
NETCDFREAD,filename,'wl1_sun',wl1_sun,solar_bands
NETCDFREAD,filename,'wl2_sun',wl2_sun,solar_bands
NETCDFREAD,filename,'wl1_earth',wl1_earth,terrestrial_bands
NETCDFREAD,filename,'wl2_earth',wl2_earth,terrestrial_bands
NETCDFREAD,filename,'month',month,dimmonth
NETCDFREAD,filename,'ext_sun',ext_sun,dimext
NETCDFREAD,filename,'omega_sun',omega_sun,dimext
NETCDFREAD,filename,'g_sun',g_sun,dimext
NETCDFREAD,filename,'ext_earth',ext_earth,dimext
NETCDFREAD,filename,'omega_earth',omega_earth,dimext
NETCDFREAD,filename,'g_earth',g_earth,dimext
;
print ,'wl1_sun=', wl1_sun
print ,'wl2_sun=', wl2_sun
print ,'wl1_ear=', wl1_earth
print ,'wl2_ear=', wl2_earth
;
print ,'min max ext_sun=', min(ext_sun),   max(ext_sun)
print ,'min max ome_sun=', min(omega_sun), max(omega_sun)
;for nl=0, NSW-1 do begin
;print ,'min max ome_sun nl=', nl, min(omega_sun(*,*,*,nl)), max(omega_sun(*,*,*,nl))
;endfor
print ,'min max ggg_sun=', min(g_sun),     max(g_sun)
print ,'min max ext_ear=', min(ext_earth), max(ext_earth)
;
dimlat=dimlatitude(0)
dimalt=dimaltitude(0)
dimtime=dimmonth(0)
;
;--vertical resolution in Beiping's code
dz=0.5
;
month_in_year=12
;;month_in_year=1
;
;--compute optical depth of input data
;--dimension ext_sun(l,k,j,nl)
;ext_sun0=ext_sun(0,*,*,*)
;print , 'size ext_sun0=', size(ext_sun0)
;tau_sun=TOTAL(ext_sun0,2)*dz
;print ,'min max tau_sun l=0 =', min(tau_sun),   max(tau_sun)
;
zz=fltarr(dimz)
;
;---approximate altitudes of LMDZ - L19 and L39
;---NOT RECOMMENDED
;
if (dimz eq 19) then begin
zz=[ 0.07125416, 0.2402302, 0.4865242, 0.8627044, 1.426484, 2.244528,    $
    3.394085, 4.952366, 6.950302, 9.26943, 11.59338, 13.71214, 15.83651, $
    18.27338, 21.16397, 24.63156, 28.9039, 34.57184, 44.51565] ;
endif
;
if (dimz eq 39) then begin
zz=[0.0338988, 0.1106602, 0.2139233, 0.3609663, 0.5685416, 0.8534464,     $
    1.233232, 1.726868, 2.355076, 3.139813, 4.101937, 5.25541, 6.596076,  $
    8.085103, 9.636482, 11.13441, 12.50023, 13.75765, 15.00424, 16.32207, $
    17.74403, 19.27899, 20.93219, 22.70941, 24.61736, 26.66389, 28.85831, $
    31.21178, 33.73778, 36.45278, 39.37709, 42.53607, 45.96209, 49.69786, $
    53.80379, 58.37727, 63.61497, 70.07092, 80.52708]
endif
;
;---exact altitudes of LMDZ -- L79
;--only exists for LR, MR and zoom_128x89 at the moment
if (dimz eq 79) then begin
  filename='./zalt_zonmean_${lmdz}_l79_rev.nc'
  NETCDFREAD,filename,'GEOP',zz,dimzz
  NETCDFREAD,filename,'LAT',zzlat,dimzzlat
  NETCDFREAD,filename,'TIME_COUNTER',zztime,dimzztime
  dimzzilat=dimzzlat[0]
  dimzzitime=dimzztime[0]
  ;--becareful zz comes with four dimensions
  ;--lon lat k time 
  print, 'GEOP size=', size(zz)
  ;--becareful zzlat comes with South Pole first
  print, 'LAT from zalt field=', zzlat
  if (dimzzilat ne dimlatlmdz) then begin
    print , 'PB dimension latitude'
  endif
endif
;
if (dimz eq 39) then begin
  filename='./zalt_zonmean_${lmdz}_l39_rev.nc'
  NETCDFREAD,filename,'GEOP',zz,dimzz
  NETCDFREAD,filename,'LAT',zzlat,dimzzlat
  NETCDFREAD,filename,'TIME_COUNTER',zztime,dimzztime
  dimzzilat=dimzzlat[0]
  dimzzitime=dimzztime[0]
  ;--becareful zz comes with four dimensions
  ;--lon lat k time
  print, 'GEOP size=', size(zz)
  ;--becareful zzlat comes with South Pole first
  print, 'LAT from zalt field=', zzlat
  if (dimzzilat ne dimlatlmdz) then begin
    print , 'PB dimension latitude'
  endif
endif
;
;
;--reconstructing the vertical coordinate at interfaces (in unit km)
;--and reverse lat axis for zzi  i==>ii
;--and reverse zz axis for zzi   dimz-1-k==>k
;--and forget about lon axis : index 0 in zz
zzi=fltarr(dimzzilat,dimzzitime,dimz+1)
for i=0, dimzzilat-1 do begin
 ii=dimzzilat-1-i
 for t=0, dimzzitime-1 do begin
  zzi(ii,t,0)=0.0
  for k=1, dimz-1 do begin
   zzi(ii,t,k)=(zz(0,i,dimz-1-(k-1),t)+zz(0,i,dimz-1-k,t))/2.0
  endfor
  zzi(ii,t,dimz)=100.00
 endfor
endfor
print, 'zzi=',zzi
;
lev=indgen(dimz)+1
;
dimzori=dimalt
zziori=fltarr(dimzori+1)
zziori(0)=altitude(0)-0.25
for k=1, dimzori do begin
 zziori(k)=altitude(k-1)+0.25
endfor
;
;--550 nm properties
;;;not available for v4
;;;tau_550_lmdz=fltarr(dimlatlmdz,dimz,month_in_year)
;
;--SW properties
tau_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
ome_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
ggg_sun_lmdz=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
;
tau_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
ome_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
ggg_sun_lmdz_ave=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
;
tau_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
ome_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
ggg_sun_lmdz_tr=fltarr(dimlatlmdz,dimz,NSW,month_in_year)
;
;--LW properties tau_abs
tau_ear_lmdz=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
;
tau_ear_lmdz_ave=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
;
tau_ear_lmdz_tr=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
;
tau_sun_lmdz_ave(*,*,*,0)=1.e-15
ome_sun_lmdz_ave(*,*,*,0)=1.e-15
ggg_sun_lmdz_ave(*,*,*,0)=0.0
tau_ear_lmdz_ave(*,*,*,0)=1.e-15
;
yearend=dimtime/month_in_year-1
;
for year=0, yearend do begin
;
;;;not available for v4
;;;tau_550_lmdz(*,*,*)=1.e-15
tau_sun_lmdz(*,*,*,*)=1.e-15
ome_sun_lmdz(*,*,*,*)=1.e-15
ggg_sun_lmdz(*,*,*,*)=0.0
tau_ear_lmdz(*,*,*,*)=1.e-15
;
chyr=strcompress(1850+year,/rem)
;
for mth=0,month_in_year-1 do begin
;
;--timestep
l=mth+month_in_year*year
print,'year mth l=',chyr, mth, l
;
;regridding
for j=0, dimlatlmdz-1  do begin
;
;--finding latitude in beiping luo data
jj=0
for jluo=0,dimlat-2 do begin
  if (latitudelmdz(j) gt (latitude(jluo)+latitude(jluo+1))/2. ) then begin
    jj=jluo+1
  endif
endfor
;
for k=0, dimz-1 do begin
for kori=0, dimzori-1 do begin
;
;fraction de la maille kori qui se trouve dans la maille k
frac= max([0.0,min([zzi(j,mth,k+1),zziori(kori+1)])-max([zzi(j,mth,k),zziori(kori)])])/(zziori(kori+1)-zziori(kori))
;
;;;not available for v4
;;;tau_550_lmdz(j,k,mth)=tau_550_lmdz(j,k,mth)+ext550(l,kori,jj)*dz*frac
;
for nl=0, NSW-1 do begin
  tau_sun_lmdz(j,k,nl,mth)=tau_sun_lmdz(j,k,nl,mth)+ext_sun(l,kori,jj,nl)*dz*frac
  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
  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
;
  tau_sun_lmdz_ave(j,k,nl,0)=tau_sun_lmdz_ave(j,k,nl,0)+ext_sun(l,kori,jj,nl)*dz*frac
  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
  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
endfor
;
for nl=0, NLW-1 do begin
  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
  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
endfor
;
endfor
;--end lat loop
;
endfor
endfor 
;--end k loops
;
endfor
;--end month loop
;
;renormalizing intensive SW properties
ggg_sun_lmdz(*,*,*,*)=ggg_sun_lmdz(*,*,*,*)/ome_sun_lmdz(*,*,*,*)
ome_sun_lmdz(*,*,*,*)=ome_sun_lmdz(*,*,*,*)/tau_sun_lmdz(*,*,*,*)
;
;saving netcdf file
;
print ,'min max ext_sun lmdz=', min(tau_sun_lmdz), max(tau_sun_lmdz)
print ,'min max ome_sun lmdz=', min(ome_sun_lmdz), max(ome_sun_lmdz)
;for nl=0, NSW-1 do begin
;print ,'min max ome_sun lmdz nl =', nl, min(ome_sun_lmdz(*,*,nl,*)), max(ome_sun_lmdz(*,*,nl,*))
;endfor
print ,'min max ggg_sun lmdz=', min(ggg_sun_lmdz), max(ggg_sun_lmdz)
print ,'min max ext_ear lmdz=', min(tau_ear_lmdz), max(tau_ear_lmdz)
;
print ,'min max ext_sun lmdz l=1=', min(tau_sun_lmdz(*,*,*,0)), max(tau_sun_lmdz(*,*,*,0))
print ,'min max ome_sun lmdz l=1=', min(ome_sun_lmdz(*,*,*,0)), max(ome_sun_lmdz(*,*,*,0))
print ,'min max ggg_sun lmdz l=1=', min(ggg_sun_lmdz(*,*,*,0)), max(ggg_sun_lmdz(*,*,*,0))
print ,'min max ext_ear lmdz l=1=', min(tau_ear_lmdz(*,*,*,0)), max(tau_ear_lmdz(*,*,*,0))
;
;--compute optical depth of output data
;tau_sun=TOTAL(tau_sun_lmdz,2)
;print ,'min max tau_sun_lmdz vert=', min(tau_sun),   max(tau_sun)
;
opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
             wav:fltarr(NSW),time:fltarr(month_in_year),          $ 
             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
;
opticstruct.lat=latitudelmdz
opticstruct.lev=lev
opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
opticstruct.time=float(indgen(month_in_year)+1)
opticstruct.tau_sun=tau_sun_lmdz
opticstruct.ome_sun=ome_sun_lmdz
opticstruct.ggg_sun=ggg_sun_lmdz
;
attributes = {units:strarr(7),long_name:strarr(7)}
attributes.units =     ['degrees_north','level','meters','month','-','-','-']
attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
;
dimensions = {isdim:intarr(7), links:intarr(4,7)}
       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
;
netcdfwrite,output+'tauswstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
            attributes=attributes, dimensions=dimensions
;
opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
             wav:fltarr(NLW),time:fltarr(month_in_year),          $
             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }
;
opticstruct.lat=latitudelmdz
opticstruct.lev=lev
opticstruct.wav=(wl1_earth+wl2_earth)/2.
opticstruct.time=float(indgen(month_in_year)+1)
opticstruct.tau_ear=tau_ear_lmdz
;
attributes = {units:strarr(5),long_name:strarr(5)}
attributes.units = ['degrees_north','level','cm-1','month','-']
attributes.long_name = ['latitude','level','wavenumber','time','tau_ear']
;
dimensions = {isdim:intarr(5), links:intarr(4,5)}
       dimensions.isdim =  [1,1,1,1,0]  ; (1=dimension, 0=variable)
       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [0,1,2,3]  ]
;
netcdfwrite,output+'taulwstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
            attributes=attributes, dimensions=dimensions
;
;;;not available in v4
;;;opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
;;;             time:fltarr(month_in_year),                          $
;;;             tau550:fltarr(dimlatlmdz,dimz,month_in_year) }
;
;;;opticstruct.lat=latitudelmdz
;;;opticstruct.lev=lev
;;;opticstruct.time=float(indgen(month_in_year)+1)
;;;opticstruct.tau550=tau_550_lmdz
;
;;;attributes = {units:strarr(4),long_name:strarr(4)}
;;;attributes.units = ['degrees_north','level','month','-']
;;;attributes.long_name = ['latitude','level','time','tau550']
;
;;;dimensions = {isdim:intarr(4), links:intarr(3,4)}
;;;       dimensions.isdim =  [1,1,1,0]  ; (1=dimension, 0=variable)
;;;       dimensions.links = [ [-1,-1,-1],[-1,-1,-1],   $
;;;                            [-1,-1,-1],[0,1,2]  ]
;
;;;netcdfwrite,output+'tau550strat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
;;;            attributes=attributes, dimensions=dimensions
;
endfor
;--end loop on years
;
;now deal with average conditions
ggg_sun_lmdz_ave(*,*,*,0)=ggg_sun_lmdz_ave(*,*,*,0)/ome_sun_lmdz_ave(*,*,*,0)
ome_sun_lmdz_ave(*,*,*,0)=ome_sun_lmdz_ave(*,*,*,0)/tau_sun_lmdz_ave(*,*,*,0)
tau_sun_lmdz_ave(*,*,*,0)=tau_sun_lmdz_ave(*,*,*,0)/float(dimtime)
tau_ear_lmdz_ave(*,*,*,0)=tau_ear_lmdz_ave(*,*,*,0)/float(dimtime)
;
for mth=1, month_in_year-1 do begin
  ggg_sun_lmdz_ave(*,*,*,mth)=ggg_sun_lmdz_ave(*,*,*,0)
  ome_sun_lmdz_ave(*,*,*,mth)=ome_sun_lmdz_ave(*,*,*,0)
  tau_sun_lmdz_ave(*,*,*,mth)=tau_sun_lmdz_ave(*,*,*,0)
  tau_ear_lmdz_ave(*,*,*,mth)=tau_ear_lmdz_ave(*,*,*,0)
endfor
;
opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
             wav:fltarr(NSW),time:fltarr(month_in_year),          $ 
             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
;
opticstruct.lat=latitudelmdz
opticstruct.lev=lev
opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
opticstruct.time=float(indgen(month_in_year))
opticstruct.tau_sun=tau_sun_lmdz_ave
opticstruct.ome_sun=ome_sun_lmdz_ave
opticstruct.ggg_sun=ggg_sun_lmdz_ave
;
attributes = {units:strarr(7),long_name:strarr(7)}
attributes.units =     ['degrees_north','level','meters','month','-','-','-']
attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
;
dimensions = {isdim:intarr(7), links:intarr(4,7)}
       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
;
netcdfwrite,output+'tauswstrat.2D.ave.nc',opticstruct,clobber=1, $
            attributes=attributes, dimensions=dimensions
;
opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
             wav:fltarr(NLW),time:fltarr(month_in_year),          $
             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }
;
opticstruct.lat=latitudelmdz
opticstruct.lev=lev
opticstruct.wav=(wl1_earth+wl2_earth)/2.
opticstruct.time=float(indgen(month_in_year))
opticstruct.tau_ear=tau_ear_lmdz_ave
;
attributes = {units:strarr(5),long_name:strarr(5)}
attributes.units = ['degrees_north','level','cm-1','month','-']
attributes.long_name = ['latitude','level','wavenumber','time','tau_ear']
;
dimensions = {isdim:intarr(5), links:intarr(4,5)}
       dimensions.isdim =  [1,1,1,1,0]  ; (1=dimension, 0=variable)
       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [0,1,2,3]  ]
;
netcdfwrite,output+'taulwstrat.2D.ave.nc',opticstruct,clobber=1, $
            attributes=attributes, dimensions=dimensions
;
;--finally prepare the 10-year transition period 
;--that is 2015-2023 for v3 or 2017-2025 for v4
;--we mix the last year of the period with the climatology
;
for year=1850+yearend+1,1850+yearend+9 do begin
;
chyr=strcompress(year,/rem)
print,'year =',chyr
;--compute weights
w1=float(1850+yearend+10-year)/10.
w2=1.0-w1
;--SW properties
tau_sun_lmdz_tr=w1*tau_sun_lmdz+w2*tau_sun_lmdz_ave
ome_sun_lmdz_tr=w1*tau_sun_lmdz*ome_sun_lmdz+w2*tau_sun_lmdz_ave*ome_sun_lmdz_ave
ggg_sun_lmdz_tr=w1*tau_sun_lmdz*ome_sun_lmdz*ggg_sun_lmdz+ $ 
                w2*tau_sun_lmdz_ave*ome_sun_lmdz_ave*ggg_sun_lmdz_ave
ggg_sun_lmdz_tr=ggg_sun_lmdz_tr/ome_sun_lmdz_tr
ome_sun_lmdz_tr=ome_sun_lmdz_tr/tau_sun_lmdz_tr
;
;--LW properties
tau_ear_lmdz_tr=w1*tau_ear_lmdz+w2*tau_ear_lmdz_ave
;
;--prepare SW output
opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
             wav:fltarr(NSW),time:fltarr(month_in_year),          $ 
             tau_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
             ome_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year),   $
             ggg_sun:fltarr(dimlatlmdz,dimz,NSW,month_in_year) }
;
opticstruct.lat=latitudelmdz
opticstruct.lev=lev
opticstruct.wav=(wl1_sun(0:NSW-1)+wl2_sun(0:NSW-1))/2.
opticstruct.time=float(indgen(month_in_year)+1)
opticstruct.tau_sun=tau_sun_lmdz_tr
opticstruct.ome_sun=ome_sun_lmdz_tr
opticstruct.ggg_sun=ggg_sun_lmdz_tr
;
attributes = {units:strarr(7),long_name:strarr(7)}
attributes.units =     ['degrees_north','level','meters','month','-','-','-']
attributes.long_name = ['latitude','level','wavelength','time','tau_sun','ome_sun','g_sun']
;
dimensions = {isdim:intarr(7), links:intarr(4,7)}
       dimensions.isdim =  [1,1,1,1,0,0,0]  ; (1=dimension, 0=variable)
       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [0,1,2,3],[0,1,2,3],[0,1,2,3]  ]
;
netcdfwrite,output+'tauswstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
            attributes=attributes, dimensions=dimensions
;
;--prepare LW output
opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
             wav:fltarr(NLW),time:fltarr(month_in_year),          $
             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year) }
;
opticstruct.lat=latitudelmdz
opticstruct.lev=lev
opticstruct.wav=(wl1_earth+wl2_earth)/2.
opticstruct.time=float(indgen(month_in_year)+1)
opticstruct.tau_ear=tau_ear_lmdz_tr
;
attributes = {units:strarr(5),long_name:strarr(5)}
attributes.units = ['degrees_north','level','cm-1','month','-']
attributes.long_name = ['latitude','level','wavenumber','time','tau_ear']
;
dimensions = {isdim:intarr(5), links:intarr(4,5)}
       dimensions.isdim =  [1,1,1,1,0]  ; (1=dimension, 0=variable)
       dimensions.links = [ [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [-1,-1,-1,-1],[-1,-1,-1,-1],   $
                            [0,1,2,3]  ]
;
netcdfwrite,output+'taulwstrat.2D.'+chyr+'.nc',opticstruct,clobber=1, $
            attributes=attributes, dimensions=dimensions
;
endfor
;--end loop on years
;
end
EOF

cat > volc_${lmdz}.job << EOF2
#PBS -N process_volc_${lmdz}
#PBS -S /bin/bash
#PBS -q week # there exist: short, day, days3, week...
#PBS -k eo  # to write the output of stdin
### Max memory
#PBS -l vmem=10gb   # virtual memory
#PBS -l  mem=10gb

cd $dirpwd

module load idl

idl << eof
.r netcdf 
.r process_volc_${lmdz}
regrid
eof
EOF2

qsub volc_${lmdz}.job