#RES=48x36
#RES=96x95
RES=144x142
#enlarge Sato file into a lat-lon field
cdo enlarge,GRID/grid-giss-enlarged INPUT/tau_zm_NOAA-GISS.nc tmp0.nc

#remap enlarged lat-lon field onto LMDZ CMIP5 grid
cdo remapcon,GRID/grid-lmdz-lonlat_${RES} tmp0.nc tmp1.nc

#regridding on the vertical
rm -f regrid.pro
cat << eod >> regrid.pro
pro regrid
filename1='tmp1.nc'
NETCDFREAD,filename1,'tau',TAU
NETCDFREAD,filename1,'lat',lat,dimlat0
NETCDFREAD,filename1,'lon',lon,dimlon0
NETCDFREAD,filename1,'time',time,dimtime0
NETCDFREAD,filename1,'altitude',altitude,dimaltitude
dimlat=dimlat0(0)
dimlon=dimlon0(0)
dimtime=dimtime0(0)
;
;---approximate altitudes of LMDZ -L39
dimz=39
zz=fltarr(dimz)
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]
zzi=fltarr(dimz+1)
zzi(0)=0.0
for k=1, dimz-1 do begin
 zzi(k)=(zz(k-1)+zz(k))/2.0
endfor
zzi(dimz)=100.00
print, 'zzi=',zzi
lev=indgen(dimz)+1
;
dimzori=4
zziori=fltarr(dimzori+1)
zziori=[15.,20.,25.,30.,35.]
;
tau2=fltarr(dimlon,dimlat,dimz,dimtime)
tau2(*,*,*,*)=0.0
;
for k=0, dimz-1 do begin
for kori=0, dimzori-1 do begin
print, 'k=',k,' kori=',kori
;
;fraction de la maille kori qui se trouve dans la maille k
frac= max([0.0,min([zzi(k+1),zziori(kori+1)])-max([zzi(k),zziori(kori)])])/(zziori(kori+1)-zziori(kori))
;
;regridding
for i=0,dimlon-1  do begin
for j=0,dimlat-1  do begin
for l=0,dimtime-1 do begin
  tau2(i,j,k,l)=tau2(i,j,k,l)+tau(i,j,kori,l)*frac
endfor
endfor
endfor
;
endfor
endfor
;
;saving netcdf file
;
taustruct={lon:fltarr(dimlon),lat:fltarr(dimlat), $
     lev:fltarr(dimz),time:fltarr(dimtime), $
     taustrat:fltarr(dimlon,dimlat,dimz,dimtime) }
;
taustruct.lon=lon
taustruct.lat=lat
taustruct.lev=lev
taustruct.time=time
taustruct.taustrat=tau2
;
attributes = {units:strarr(5),long_name:strarr(5)}
attributes.units = ['degrees_east','degrees_north','level','month','taustrat']
attributes.long_name = ['longitude','latitude','level', 'time', 'TAUSTRAT']
;
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_${RES}/taustrat.nc',taustruct,clobber=1, $
            attributes=attributes, dimensions=dimensions
;
end
eod
#
#--calling IDL
#
/opt/idl-6.4/idl/bin/idl << eod
.r netcdf
.r regrid
regrid
exit
eod
#
#
rm -f tmp0.nc tmp1.nc regrid.pro
#
for yr in {1850..2010}
do
month=`echo $((${yr}-1850))*12+1 | bc`
for mth in {2..12}
do
month=${month}','`echo $((${yr}-1850))*12+${mth} | bc`
done
echo ${month}
cdo seltimestep,${month} OUTPUT_${RES}/taustrat.nc OUTPUT_${RES}/taustrat.${yr}.nc
done