pro concatoxf

;------------------------------------
file1='input_totoptdep1.nc'
file2='input_totoptdep2.nc'
;------------------------------------

cdfid1 = ncdf_open(file1)
cdfid2 = ncdf_open(file2)

; coordinates
varid=ncdf_varid(cdfid1,'lat')
        ncdf_varget, cdfid1, varid, lati
lat=n_elements(lati)

varid=ncdf_varid(cdfid1,'lon')
        ncdf_varget, cdfid1, varid, longi
lon=n_elements(longi)

varid=ncdf_varid(cdfid1,'time')
        ncdf_varget, cdfid1, varid, time1
varid=ncdf_varid(cdfid2,'time')
        ncdf_varget, cdfid2, varid, time2
t1=n_elements(time1)
t2=n_elements(time2)

; concat time coordinates
ntime=t1+t2
timear=fltarr(ntime)
timear(0:t1-1)=time1(0:t1-1)
timear(t1:ntime-1)=time2(0:t2-1)+time1(t1-1)


; field
varid=ncdf_varid(cdfid1,'totoptdep')
        ncdf_varget, cdfid1, varid, totoptdep1
varid=ncdf_varid(cdfid2,'totoptdep')
        ncdf_varget, cdfid2, varid, totoptdep2


;loadct, 3
;!p.multi=[0,1,3]
;ref=totoptdep1(*,*,tata)
;smo=smooth(ref,toto)
;err=100.*abs(ref-smo)/ref
;err=(ref-smo)
;print, max(err), ref(where(err eq max(err)))
;contour, err, nlevels=40, /cell_fill
;contour, ref, nlevels=40, /cell_fill
;contour, smo, nlevels=40, /cell_fill
;stop



; concat field
totoptdep=fltarr(lon,lat,ntime)
totoptdep(*,*,0:t1-1)=totoptdep1(*,*,*)
totoptdep(*,*,t1:ntime-1)=totoptdep2(*,*,*)
fieldname1='totoptdep'
outfieldstring1='Total dust  optical depth (tauref/700Pa)'
unit1='Opacity/Pa'



;;;*********************************************************************

;
; SMOOTH wrt SPACE
;
totoptdep=smooth(totoptdep,2)


;
; INTERPOLATE wrt TIME
;

which_utc_is_TES=3
;; for Hellas 8AM is 02PM LT, i.e. TES measurements
;; IDL convention, starts 0 !!!!


nday=floor(float(ntime-1)/12.)-1

ref=fltarr(nday)
subsarr=fltarr(nday)

for i=0, n_elements(totoptdep(*,0,0))-1 do begin
	for j=0, n_elements(totoptdep(0,*,0))-1 do begin
		for each=0,nday-1 do begin
  		        subs=each*12 + which_utc_is_TES
			ref(each)=totoptdep(i,j,subs)
			subsarr(each)=subs
		endfor
		;;yeah=interpol(ref,subsarr,findgen(ntime), /spline)
		yeah=interpol(ref,subsarr,findgen(ntime) )
		totoptdep(i,j,*)=yeah
	endfor
endfor

;;;*********************************************************************





;--------------------------------
print,'writing netcdf file'

filencdf='input_totoptdep_tot.nc'
cdfid=ncdf_create(filencdf,/clobber)

londimid=ncdf_dimdef(cdfid,'lon',lon)
latdimid=ncdf_dimdef(cdfid,'lat',lat)
tdimid=ncdf_dimdef(cdfid,'time',ntime)
lonid=ncdf_vardef(cdfid,'lon',[londimid],/float)
ncdf_attput,cdfid,lonid,'long_name','longitude'
ncdf_attput,cdfid,lonid,'units','degrees_east'
latid=ncdf_vardef(cdfid,'lat',[latdimid],/float)
ncdf_attput,cdfid,latid,'long_name','latitude'
ncdf_attput,cdfid,latid,'units','degrees_north'
timeid=ncdf_vardef(cdfid,'time',[tdimid],/float)
ncdf_attput,cdfid,timeid,'long_name','model time'
ncdf_attput,cdfid,timeid,'units','hours since 00:00:00'
varid1=ncdf_vardef(cdfid,fieldname1,[londimid,latdimid,tdimid],/float)
ncdf_attput,cdfid,varid1,'Physics_diagnostic',outfieldstring1
ncdf_attput,cdfid,varid1,'units',unit1
ncdf_control,cdfid,/endef

ncdf_varput,cdfid,lonid,longi
ncdf_varput,cdfid,latid,lati
ncdf_varput,cdfid,timeid,timear
ncdf_varput,cdfid,varid1,totoptdep

ncdf_close,cdfid

end