tower='wcreek'


yyyy='2002'
load,'/eddy/s1/cheas/flux/wcreek/wcreek'+yyyy+'_met.txt',met
load,'/eddy/s1/cheas/flux/wcreek/wcreek'+yyyy+'_flux.txt',flux

badval=-999

mon=met(1,*)
day=met(2,*)
hour=met(3,*)
jday=met(4,*)
fjday=met(5,*)
ustar=flux(11,*)
heat=flux(10,*)

wdir=met(20,*) ;wind diorection at 97ft
temp=met(30,*) ;temperature at 97ft

mo=temp*0.0+badval
r=where(heat NE badval AND ustar ge 0 AND temp NE badval,count)
IF count GE 1 THEN mo(r)=(-(temp(r)+273.15)*(ustar(r)^3)/(0.4*9.8*heat(r)/1004.0/1.29))


n=10000 ; number of points in x-direction
dx=0.5
x=findgen(n)*dx  ; m distance from the sensor
f=fltarr(48,n)+badval   ; footprint array (half hourly)
fcum=f
z0=2.3          ; roughtness  2.3m
zm=(29.6-16.4*0)        ; z-d m ; Bruce's estimate based on LAI 

dayf=fltarr(n)+0.0

nightf=fltarr(n)+0.0

wdirf0=fltarr(n)+0.0  ; 90-180 wind direction
wdirf1=wdirf0



!p.multi=[0,3,4]
!p.ticklen=1
!x.gridstyle=1
!y.gridstyle=1
!y.style=1
!y.minor=3
!y.ticks=5

printps1,'footprint.ps','portrait'


dayn=0.
nightn=0.
nwdirf0=0.
nwdirf1=0.

FOR i=222,223 DO begin  ; day loop

r=where(jday EQ i)
rmo=1.0D*mo(r)
dir=wdir(r)
  f=f*0+badval
  fcum=f         ; set intial values
  FOR j=0,len(r)-1,2+2 DO BEGIN  ; hour loop
    
    print,rmo(j)


    IF rmo(j) EQ -999 THEN GOTO,skip
    z_x,z0,1d*rmo(j),xplume,zplume
    fcum(j,*)=0d
    f(j,0)=0d
    FOR ip=0+1,n-1 DO BEGIN ; x point
       f(j,ip)=fy(rmo(j),z0,xplume,zplume,x(ip),zm)
       fcum(j,ip)=fcum(j,ip-1)+f(j,ip)  ;  * dx
    ENDFOR ; x point
     totalf=total(f(j,*))
      if totalf*dx gt 1 then f(j,*)=f(j,*)/(totalf*dx)

;     if rmo(j) gt 0 and totalf*dx gt 1 then   f(j,*)=f(j,*)/(totalf*dx)
;     if rmo(j) lt 0 and totalf*dx lt 1 then   f(j,*)=f(j,*)/(totalf*dx)

       fcum(j,*)=fcum(j,*);/(total(f(j,*))*dx)

    IF j Gt 6*2 AND j Lt 18*2 THEN BEGIN ;day
      dayf(*)=dayf(*)+f(j,*)
      dayn=dayn+1
    ENDIF else begin ;day
      nightf(*)=nightf(*)+f(j,*)
      nightn=nightn+1
   ENDelse ;night

    IF dir(j) GE 90 AND dir(j) LE 180 THEN BEGIN ; 90-180
    wdirf0(*)=wdirf0(*)+f(j,*)
    nwdirf0=nwdirf0+1
    ENDIF ELSE BEGIN
    wdirf1(*)=wdirf1(*)+f(j,*)
    nwdirf1=nwdirf1+1
    ENDELSE
     print,i,j,total(f(j,*))*dx


skip:
ymax=0.015
xmin=1e-3 & xmax=5
    plot,x/1000,f(j,*),title='day='+strtrim(string(i,'(I3)'),2)+' hh='+strtrim(string(1.0*j/2,'(F4.1)'),2)+'LST', xtitle='distance (km)', ytitle='Ft (1/m)',yrange=[0,ymax],/xlog,xrange=[xmin,xmax]

;    plot,x/1000,dx*fcum(j,*),title='day='+strtrim(string(i,'(I3)'),2)+' hh='+strtrim(string(1.0*j/2,'(F4.1)'),2)+'LST', xtitle='distance (km)', ytitle='cum Ft ',yrange=[0,1],yminor=2;,/xlog
    legend,['MO(m)='+string(rmo(j),'(I6)'),'Winddir='+string(dir(j),'(I4)')],psym=[3,3],/right,box=0,charsize=0.6,symsize=[0,0],margin=0.2



  ENDFOR  ; hour loop

ENDFOR ; day loop


!p.multi=[0,2,2]

  plot,x/1000,dayf(*)/dayn,title='DAY', ytitle='mean footprint(1/m)',xtitle='km',yrange=[0,ymax],ystyle=1,/xlog,xrange=[xmin,xmax]
    plot,x/1000,nightf(*)/nightn,title='night', ytitle='mean footprint (1/m)',xtitle='km',yrange=[0,ymax],/xlog,xrange=[xmin,xmax]
  plot,x/1000,wdirf0(*)/nwdirf0,title='90-180', ytitle='Mean footprint(1/m)',xtitle='km',yrange=[0,ymax],xrange=[xmin,xmax],/xlog
  plot,x/1000,wdirf1(*)/nwdirf1,title='other directions', ytitle='Mean footprint(1/m)',xtitle='km',yrange=[0,ymax],xrange=[xmin,xmax],/xlog
printps2


!y.style=0
END