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

dirpwd=`pwd`

vv='v3'

#--choose the resolution
#--only LR works for now !
#lmdz='VVLR'
#lmdz='VLR_L79'
#lmdz='VLR'
lmdz='LR'
#lmdz='MR'

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

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

if (lmdz eq 'VLR') then begin
dimlonlmdz=96
dimlatlmdz=96
dimz=39
output='/data/${USER}/CMIP6/VOLC/VLR_${vv}/'
endif

if (lmdz eq 'VLR_L79') then begin
dimlonlmdz=96
dimlatlmdz=96
dimz=79
output='/data/${USER}/CMIP6/VOLC/VLR_L79_${vv}/'
endif

if (lmdz eq 'LR') then begin
dimlonlmdz=144
dimlatlmdz=143
dimz=79
output='/data/${USER}/CMIP6/VOLC/LR_${vv}/'
endif

if (lmdz eq 'MR') then begin
dimlonlmdz=256
dimlatlmdz=257
dimz=79
output='/data/${USER}/CMIP6/VOLC/MR_${vv}/'
endif

latfirst=90.
latinc=-180./float(dimlatlmdz-1)
lonfirst=-180.
loninc=360./float(dimlonlmdz-1)
latitudelmdz=latfirst+latinc*indgen(dimlatlmdz)
longitudelmdz=lonfirst+loninc*indgen(dimlonlmdz)
NSW=6
NLW=16
;
;--script only works on ciclad
dir='/prodigfs/project/input4MIPs/VOLC/${vv}/'
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 at the moment
if (dimz eq 79) then begin
  filename='./zalt_zonmean_LR_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
;
;
;--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
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)
;
;--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_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
;
for year=0, dimtime/month_in_year-1 do begin
;for year=141, 142 do begin ;--Pinatubo
;
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))
;
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
;
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
;
end
EOF

cat > volc.job << EOF2
#PBS -N process_volc
#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
regrid
eof
EOF2

#qsub volc.job