source: TOOLS/CMIP6_FORCING/AER_STRAT/ECRAD/volc_ECRAD.sh @ 6267

dirpwd=`pwd`

#original CMIP6 files
#vv='v3'
#corrected CMIP6 files + extended in time
vv='v4_ECRAD'

#--choose the resolution
#--only LR and zoom works fully for now !
#lmdz='VVLR_l19'
#lmdz='VLR_l39'
#lmdz='VLR_l79'
lmdz='LR_l79'
#lmdz='MR_l79'
#--added for Frederique Cheruy
#lmdz='zoom_128x89_l79'

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_l19') then begin
dimlonlmdz=48
dimlatlmdz=37
dimz=19
output='/data/${USER}/CMIP6/VOLC/${lmdz}_${vv}/'
latfirst=90.
latinc=-180./float(dimlatlmdz-1)
latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
endif

if (lmdz eq 'VLR_l39') then begin
dimlonlmdz=96
dimlatlmdz=96
dimz=39
output='/data/${USER}/CMIP6/VOLC/${lmdz}_${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/${lmdz}_${vv}/'
latfirst=90.
latinc=-180./float(dimlatlmdz-1)
latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
endif

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

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

if (lmdz eq 'zoom_128x89_l79') 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

; Modif ThL: for ECRAD, we have 14 SW bands and 16 LW bands.
NSW=14
NLW=16
;
;--script only works on ciclad
dir='/bdd/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_ECRAD_radiation_v4_1850-2018.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)
;
;---exact altitudes of LMDZ -- L79
;--only exists for VLR, LR, MR and zoom_128x89 at the moment
;
filename='/bdd/input4MIPs/ZALT/zalt_zonmean_${lmdz}_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
;
;
;--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_ear_lmdz=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
ome_ear_lmdz=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
ggg_ear_lmdz=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
;
tau_ear_lmdz_ave=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
ome_ear_lmdz_ave=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
ggg_ear_lmdz_ave=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
;
tau_ear_lmdz_tr=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
ome_ear_lmdz_tr=fltarr(dimlatlmdz,dimz,NLW,month_in_year)
ggg_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
ome_ear_lmdz_ave(*,*,*,0)=1.e-15
ggg_ear_lmdz_ave(*,*,*,0)=0.0
;
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
ome_ear_lmdz(*,*,*,*)=1.e-15
ggg_ear_lmdz(*,*,*,*)=0.0
;
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)+ext_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*dz*frac
  ggg_sun_lmdz(j,k,nl,mth)=ggg_sun_lmdz(j,k,nl,mth)+ext_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*g_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)+ext_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*dz*frac
  ggg_sun_lmdz_ave(j,k,nl,0)=ggg_sun_lmdz_ave(j,k,nl,0)+ext_sun(l,kori,jj,nl)*omega_sun(l,kori,jj,nl)*g_sun(l,kori,jj,nl)*dz*frac
endfor
;
; For LW, script differs from previous versions in that
; (1.-omega_earth) factor becomes (omega_earth)
;
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)*dz*frac
  ome_ear_lmdz(j,k,nl,mth)=ome_ear_lmdz(j,k,nl,mth)+ext_earth(l,kori,jj,nl)*omega_earth(l,kori,jj,nl)*dz*frac
  ggg_ear_lmdz(j,k,nl,mth)=ggg_ear_lmdz(j,k,nl,mth)+ext_earth(l,kori,jj,nl)*omega_earth(l,kori,jj,nl)*g_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)*dz*frac
  ome_ear_lmdz_ave(j,k,nl,0)=ome_ear_lmdz_ave(j,k,nl,0)+ext_earth(l,kori,jj,nl)*omega_earth(l,kori,jj,nl)*dz*frac
  ggg_ear_lmdz_ave(j,k,nl,0)=ggg_ear_lmdz_ave(j,k,nl,0)+ext_earth(l,kori,jj,nl)*omega_earth(l,kori,jj,nl)*g_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)
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 ome_ear lmdz=', min(ome_ear_lmdz), max(ome_ear_lmdz)
print ,'min max ggg_ear lmdz=', min(ggg_ear_lmdz), max(ggg_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))
print ,'min max ome_ear lmdz l=1=', min(ome_ear_lmdz(*,*,*,0)), max(ome_ear_lmdz(*,*,*,0))
print ,'min max ggg_ear lmdz l=1=', min(ggg_ear_lmdz(*,*,*,0)), max(ggg_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','µm','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
;
; NOW FOR LW
;
opticstruct={lat:fltarr(dimlatlmdz),lev:fltarr(dimz),             $ 
             wav:fltarr(NLW),time:fltarr(month_in_year),          $
             tau_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
             ome_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
             ggg_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
opticstruct.ome_ear=ome_sun_lmdz
opticstruct.ggg_ear=ggg_sun_lmdz

;
attributes = {units:strarr(7),long_name:strarr(7)}
attributes.units = ['degrees_north','level','µm','month','-','-','-']
attributes.long_name = ['latitude','level','wavelength','time','tau_ear','ome_ear','g_ear']
;
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+'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)
;
ggg_ear_lmdz_ave(*,*,*,0)=ggg_ear_lmdz_ave(*,*,*,0)/ome_ear_lmdz_ave(*,*,*,0)
ome_ear_lmdz_ave(*,*,*,0)=ome_ear_lmdz_ave(*,*,*,0)/tau_ear_lmdz_ave(*,*,*,0)
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)
  ggg_ear_lmdz_ave(*,*,*,mth)=ggg_ear_lmdz_ave(*,*,*,0)
  ome_ear_lmdz_ave(*,*,*,mth)=ome_ear_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','µm','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),   $
             ome_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
             ggg_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
opticstruct.ome_ear=tau_ome_lmdz_ave
opticstruct.ggg_ear=tau_ggg_lmdz_ave
;
attributes = {units:strarr(7),long_name:strarr(7)}
attributes.units = ['degrees_north','level','µm','month','-','-','-']
attributes.long_name = ['latitude','level','wavelength','time','tau_ear','ome_ear','g_ear']
;
dimensions = {isdim:intarr(7), links:intarr(4,7)}
       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],[0,1,2,3],[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
ome_ear_lmdz_tr=w1*tau_ear_lmdz*ome_ear_lmdz+w2*tau_ear_lmdz_ave*ome_ear_lmdz_ave
ggg_ear_lmdz_tr=w1*tau_ear_lmdz*ome_ear_lmdz*ggg_ear_lmdz+ $ 
                w2*tau_ear_lmdz_ave*ome_ear_lmdz_ave*ggg_ear_lmdz_ave
ggg_ear_lmdz_tr=ggg_ear_lmdz_tr/ome_ear_lmdz_tr
ome_ear_lmdz_tr=ome_ear_lmdz_tr/tau_ear_lmdz_tr
;
;--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','µm','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),   $
             ome_ear:fltarr(dimlatlmdz,dimz,NLW,month_in_year),   $
             ggg_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
opticstruct.ome_ear=ome_ear_lmdz_tr
opticstruct.ggg_ear=ggg_ear_lmdz_tr
;
attributes = {units:strarr(7),long_name:strarr(7)}
attributes.units = ['degrees_north','level','µm','month','-','-','-']
attributes.long_name = ['latitude','level','wavelength','time','tau_ear','ome_ear','g_ear']
;
dimensions = {isdim:intarr(7), links:intarr(4,7)}
       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],[0,1,2,3],[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

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

qsub volc_${lmdz}.job