; this code is used to compute the footprint for e.g., a year at a
; tower using surface  footprint formula (HW1994)

; use another way to reduce the computation time.


tower='wlef'

print,'input yr, yyyy'
yr=''
read,yr

print,'input use a guessd z0 and D (y/n)'
guess=''
read,guess

if guess eq 'y' then begin 
guessnumber='0'
print,'input guess number= (0,1,2,3...)'
read,guessnumber
endif

badval=-999.
type='10'    ; type is an indicator of coefficient of horizontal diffusion , sigmy=type/10*(expression)
 



z0=0.9d ; roughness length
displace=0.6*15.0 ; 2/3 canopy height, at WLEF, canopy is roughly 15m
displace=15.0 ; m canopy =20m
width=1.0  ; 5/13/05

IF guess EQ 'y' THEN BEGIN
 IF guessnumber EQ '0' THEN BEGIN 
 z0=0.4d
 displace=10;m   ; use an arbitary z0 and displace depth to consider grassy area, since the respiration rates for aspen is too high if use z0=0.9 or 0.5
 subname='footprintcases/guess' 
ENDIF

  IF guessnumber EQ '1' THEN BEGIN
  z0=0.3d
  displace=8.0 ;m
 subname='footprintcases/guess'+guessnumber 
  ENDIF

  IF guessnumber EQ '2' THEN BEGIN
  z0=0.2d
  displace=6.0
 subname='footprintcases/guess'+guessnumber 
  ENDIF


  IF guessnumber EQ '3' THEN BEGIN
  z0=0.1d
  displace=4.0
 subname='footprintcases/guess'+guessnumber 
  ENDIF

; experiemnts about different horizontal diffusion coeffient using z0
; and D at guess2

    IF guessnumber EQ '4' THEN BEGIN
  z0=0.2d
  displace=6.0
 subname='footprintcases/guess'+guessnumber 
  width=0.5
  ENDIF

    IF guessnumber EQ '5' THEN BEGIN
  z0=0.2d
  displace=6.0
 subname='footprintcases/guess'+guessnumber 
  width=2.0
  ENDIF




ENDIF ELSE BEGIN

  z0=0.5d
  displace=15.0d ; control case  ; same as ft_long_surface_bruce.pro 
  subname='footprintcases/control'  
ENDELSE


;hor_coef=1d*fix(type)/10.0

hor_coef=1d*fix(type)/10.0*width

nx=401
ny=401
grid=30.0D 
nworkx=nx*1.5/2 ; size of the work array
nworky=ny*1.5

work=fltarr(nworkx,nworky)+0D ;
xw=findgen(nworkx)*grid
yw=findgen(nworky)*grid
; translate in x, y direction
;       ^ y 
;       ! 
;-------!---------->x wind direction
;       !source at (0,0)   
;       ! 

xo=xw(nworkx-1)
yo=yw((nworky+1)/2-1) ; origin coordinates  ; add -1 on 7/4/2004

xw=xw-xo
yw=yw-yo



ft=fltarr(nx,ny)+0.0D
veg=fltarr(nx,ny)
fttemp=ft
x=findgen(nx)*grid
y=findgen(ny)*grid
xc=x((nx+1)/2-1)
yc=y((ny+1)/2-1) ; tower coordinates ; add -1 on 7/4/2004, make sure tower is in the center of the map. 

;transform coordiant into one centered at the tower
x=x-xc
y=y-yc

; open output files
;outfilename='/eddy/s3/wang/thesis/ft/decompose_flux/results/vegratio_'+tower+yr+'_30.txt'+type
IF guess EQ 'y' THEN outfilename='/eddy/s3/wang/thesis/ft/decompose_flux/results/'+subname+'/vegratio_'+tower+yr+'_30.txt'+type else outfilename='/eddy/s3/wang/thesis/ft/decompose_flux/results/'+subname+'/vegratio_'+tower+yr+'_30.txt'+type


openw,14,outfilename


; load, veg data
temparr=fltarr(nx)

IF tower EQ 'wlef' THEN vegmapdataname='/eddy/s3/wang/thesis/ft/landmap/dnr_wi/croppedwlef-20k-wlc401.dat'

openr,12,vegmapdataname

FOR j=0,ny-1 DO BEGIN
 readf,12,temparr
 veg(*,ny-1-j)=temparr   ; y=0 (southmost), y=ny-1 (northmost) in veg
ENDFOR
close,12

;input data to compute wstar, 
; heat flux, m-o length, zi, wind speed, wind direction,


; load met file

get_data,'/eddy/s3/wang/thesis/ablbudget/data/WLEF'+yr+'_met.txt',met,column=44,rowmax=10000L ;LST

mon=met(1,*)
day=met(2,*)
hr=met(3,*)
jday=met(4,*)
fjday=met(5,*)
nhr=len(hr)


u=met(15,*)  ; wind speed; m/s at 30m

wdir=met(18,*) ; wind direction at 30m

temp=met(21,*); temperature at 30m C

par1=met(36,*);par at wlef
par2=met(37,*);par at Willow spring clear-cut site
par3=met(38,*);par at Willow spring full forest cover
par=par1*0.0+badval

; wind directions at three levels

wdir30=met(18,*)
wdir122=met(17,*)
wdir396=met(16,*)
utemp=fltarr(3)+10.0; used to compute the mean wind direction
udtemp=fltarr(3)
wdirmean=fltarr(nhr)+badval

; compute mean wind dirction

FOR i=0,nhr-1 DO BEGIN
 udtemp(0)=wdir30(i)
 udtemp(1)=wdir122(i)
 udtemp(2)=wdir396(i)
 wdirmean,utemp, udtemp,umean,udmean
 wdirmean(i)=udmean

; print,i,wdir30(i),wdir122(i),wdir396(i),wdirmean(i)
ENDFOR



; load flux file

get_data,'/eddy/s3/wang/thesis/ablbudget/data/WLEF'+yr+'_flux.txt',flx,column=42,rowmax=10000L ;LST

heat=flx(37,*); preferred sensible heat, w/m2

ustar=flx(39,*); ustar at 30m

co2=flx(12,*); co2 flux at 30m

mo=ustar*0.0+badval

r=where(heat NE badval AND ustar NE badval 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)


;---------------------------------------------------------------------------------------------------
; compute the ft
number=0.0
totalco2=0.0

FOR k=0,len(hr)-1 DO BEGIN
 print,'jday=',jday(k),hr(k)
; set initial values
 a50=badval
 a90=badval
 s50=badval
 s90=badval
 ratio=fltarr(12)+badval



;judge the conditions to select which ft model is used.

;IF wdir30(k) NE badval  AND mo(k) NE badval  AND mo(k) LE 0 AND mon(k) GE 5 AND mon(k) LE 10 and hr(k) ge 9 and hr(k) le 17 THEN BEGIN 

;7/11/2004 including cases with unstable and weak stable stratisfication

IF wdir30(k) NE badval  AND mo(k) NE badval  AND (mo(k) LE 0 OR mo(k) GT 300) THEN BEGIN 


;compute the relationship between x and mean height of plume
; abs(mo) should >> z0

 ;mo(k)=mo(k)<(-z0*10.0) ;mo is negative here

  IF abs(mo(k)) LE z0*20.0 THEN mo(k)=z0*20.0*(abs(mo(k))/mo(k))
  IF mo(k) LT 0 AND (hr(k) LE 7 OR hr(k) GE 20) THEN mo(k)=5000.0  ; netural

z_x,z0,mo(k),xplume,zplume

 pztemp=fltarr(nworkx)
FOR i=0,nworkx-1 DO BEGIN
 

;calculate ft_z ; crosswind-integral ft 
  pz=fy(mo(k),z0,xplume,zplume,abs(xw(i)+1.0e-10),30.0-displace) ; 1/m unit
;   print,xw(i),pz
  pztemp(i)=pz
;  sy=hor_coef*(abs(1.3/mo(k))+0.1)*(abs(xw(i))+1e-10)/sqrt(1.0+0.0001*abs(xw(i))+1e-10)+1.0  ;min sigmY=1.0

;  8/30/2004, use briggs(1973)'s fitting under urban condition 

  sy=hor_coef*(abs(4.0/mo(k))+0.16)*abs(xw(i))/sqrt(1.0+0.0004*abs(xw(i)))+1.0

  ;min sigmY=1.0 (note that 0.1 is changed to 0.16 to approximate the coefs in the Briggs' formula) 
  
   IF mo(k) GE 300 OR mo(k) LT -300 THEN sy=hor_coef*(4./300+0.16)*(abs(xw(i))>1.00)/sqrt(1+0.0004*(abs(xw(i))>1.00))  ; Briggs(1973) P-G stability type "D",sigmY, Akula Venkatram(1988),Dispersion in the stable boundary layer, in: Lectures on air pollution modeling, 229-265. edited by Akula Venkatram and John C Wyngarrd. 
; 8/30/04 , change abs(xw(i))>100 into abs(xw(i))>1.0

  sy=sy>1.0

;  print,i,nxw,zp,rmo,pz  
     FOR j=0,nworky-1 DO BEGIN
; crosswind distance, y -direction

  py=0.0D
  expnyw=1.0D*yw(j)*yw(j)/(2.0*sy*sy)
; 7/6/2004 , change i ne 0 into i ne nworkx-1, xw(nworkx-1)=0 

  IF expnyw LE 150 AND i NE 0 THEN py=exp(-expnyw)/(sqrt(6.28)*sy)
  work(i,j)=py*pz ; 1/m2 unit
;   print,i,j,pz,py

     ENDFOR

ENDFOR

  print,'total(work)=',total(work)*grid*grid
  print,'zo=',z0,'d=',displace,'sigmY width=',width
;contour,work,xw/1000.,yw/1000.


; transform the results in the wind direction coordinate into real
; coordinate
  fttemp=fttemp*0.0
; determine wind direction to be used for strong unstable consditions, 
; use mean wind direction, otherwise, use wdir30 , if wdir30 eq
; badval, then use mean value.

  winddirection=wdir30(k)
  IF winddirection EQ badval  THEN winddirection=wdirmean(k)

print,jday(k),hr(k),wdir30(k),winddirection,mo(k),co2(k)

  FOR i=0,nx-1 DO begin

  FOR j=0,ny-1 DO begin

; transform coordinate
; get x,y in real coordinates transformed into wind coordinates, find
; where it is, and get an estimate

xytowindcoordinates,winddirection,x(i),y(j),xwind,ywind

; for downwind direction, no need to calculate
 IF xwind GE 0.0 THEN GOTO, ccc
  
; find the index number 
  ix=min([fix((xwind-xw(0))/grid),nworkx])
  iy=min([fix((ywind-yw(0))/grid),nworky])
  IF abs(work(ix,iy)) LE 1.0e-15 THEN GOTO,ccc

; interpolate the value from near four points
  ix1=min([ix+1,nworkx])
  iy1=min([iy,nworky])
  ix2=min([ix+1,nworkx])
  iy2=min([iy+1,nworky])
  ix3=min([ix,nworkx])
  iy3=min([iy+1,nworky])

  d1=distance2(xwind,ywind,xw(ix),yw(iy))+1e-8 ; distance square
  d2=distance2(xwind,ywind,xw(ix1),yw(iy1))+1e-8
  d3=distance2(xwind,ywind,xw(ix2),yw(iy2))+1e-8
  d4=distance2(xwind,ywind,xw(ix3),yw(iy3))+1e-8

  d=1.0/d1+1.0/d2+1.0/d3+1.0/d4
  fttemp(i,j)=(work(ix,iy)/d1+work(ix1,iy1)/d2+work(ix2,iy2)/d3+work(ix3,iy3)/d4)/d


ccc:
  ENDFOR 

 ENDFOR

;normalize

 IF total(fttemp) NE 0 THEN fttemp=fttemp/(total(fttemp)*grid*grid)

 a50=findcontour(fttemp,grid,0.5,source_area=s50)
 a90=findcontour(fttemp,grid,0.9,source_area=s90)
; contour,fttemp,x/1000.,y/1000.,levels=[a95,a50+1.e-10]

 findvegarea_dnr,fttemp,veg,ratio
 print,total(ratio),'---'

 if co2(k) ne badval then begin ;weighted ft  
   ft=ft+fttemp;        *abs(co2(k))
   totalco2=totalco2+1.0    ;abs(co2(k))
; & number=number+1
 endif ;weighted ft

ENDIF     

 par(k)=meanbadval(met(36:38,k))
; print,'par=',par(k),par1(k),par2(k),par3(k)
 printf,14,met(0:5,k),a50,a90,s50,s90,par1(k),co2(k),ratio,format='(3I5,f6.1,I5,f10.4,6e15.4,12f11.5)'

ENDFOR
 
;ft=ft/number  ; mean footprint 
 ft=ft/totalco2 ; weighted mean footprint

;outgridfile='/eddy/s3/wang/thesis/ft/decompose_flux/results/ft_long_surface_grid'+yr+'.txt'+type
IF guess EQ 'y' THEN outgridfile='/eddy/s3/wang/thesis/ft/decompose_flux/results/'+subname+'/ft_long_surface_grid'+yr+'.txt'+type else outgridfile='/eddy/s3/wang/thesis/ft/decompose_flux/results/'+subname+'/ft_long_surface_grid'+yr+'.txt'+type

openw,16,outgridfile

printf,16,nx,'nx'
printf,16,ny,'ny'
printf,16,grid,'grid size (m)'
printf,16,ft
close,16

close,14

print,'saving '+outfilename
print,'saving '+outgridfile

end