; this code is used to compute footprint using HW 1992 and 1994
;
;  4/28/2004 ,  considering large roughness length z0, its effects should be considered in the calculation of wind profiles and concentration, e..g, use zm-z0 instead of zm,.....
;  Note that, If L is too small, (e.g., abs(L)= or < z0, the plume height and distance cannot be correctly calculated using the formulas. SO L should >> z0

FUNCTION fy,mo,z0,x,z,xp,zm

; input
; mo --- M-O length (m) ; should be smaller than z0, and if abs(mo)>1e5, then mo=1e5 or -1e5
; z0 ---- roughness length (m)
; x, z --- horizontal distance(m) and plume height (m) built by the
;          code of z_x.pro, should call this code before calculation
; xp   --- distance from the sensor (m)
; zm   ---- measurement height (m) should be smaller than 100m

;  output
;  fy ----- cross-wind integrated footprint (1/m)


b=1.5164 ; b=gamma(1/r)/gamma(2/r), r=1.5
p=1.55
a=0.73  ;  a=r*gamma(2/r)/gamma(1/r)^2 
c=0.63
beta=5.0
r=1.5
fy=0.0

; compute psim, phih

zmean=zbar(x,z,xp)

psimphih,p*zmean,mo,psim,phih
psimphih,z0,mo,psimz0,phihz0
dzdx=0.4*0.4/((alog(p*zmean/z0)-psim+0*psimz0)*phih)

; compute the wind speeds

IF mo GE 0 THEN BEGIN
 
 uzm=(alog(zm/z0)+beta*zm/mo-beta*z0/mo); coef, ustar/0.4 is ignored
; czbar=(c*zmean)>z0   ; take max(c*zmean,z0), 
                       ;if z0 is small enough, the two methods give the similar results
 czbar=z0+c*(zmean-z0)

 uzbar=(alog(czbar/z0)+beta*zmean/mo-beta*z0/mo); coef, ustar/0.4 is ignored, 
; note that the formula is different from wind profile in the surface layer

 if uzbar le 0.0  then return,fy
ENDIF ELSE BEGIN

 psimphih,zm,mo,psimzm,phih
 psimphih,z0,mo,psimz0,phihz0
 uzm=(alog(zm/z0)-psimzm+psimz0)
;
; czbar=(c*zmean)>z0   ; take max(c*zmean,z0)
  czbar=z0+c*(zmean-z0)
; see smae note for case of L >0

 psimphih,czbar,mo,psimbar,phih

 uzbar=(alog(czbar/z0)-psimbar+psimz0); coef, ustar/0.4 is ignored

if uzbar le 0 then return,fy
ENDELSE


expz=0.0
zm1=zm-z0  &  zmean1=zmean-z0
eee=(zm1/b/zmean1)^r
IF eee LE 20 THEN expz=exp(-eee)

fy=dzdx*(zm1/(zmean1*zmean1))*(uzm/uzbar)*a*expz
return,fy
end