; this code is used to calculate the crosswind integrated footprint at 
; an upwind distance.
; changed on 7/01/2004 using latest parameters
; changed 0n 8/31/2004 using fit_all.pro results


FUNCTION cblfootprintz,xp,nz,rmo,z0

; zi=1000.0
 wstar=1.0; input from the dummny
; u=1.0D
; q=1.0D
 zs=0.0D  ; source height
 nzs=zs;/zi  nondimensional source height
 sigmw=0.5*wstar
 s=0.5D
; r=0.8D
  r=1.0 ;2/42004 

 g1=(1.0+r*r)/(1+3*r*r)  ; 2/4/2004 add a 'r'
 g2=1+r*r

 w1=g1*s/2.0+sqrt(g1*g1*s*s+4.0/g2)/2.0
 w1=w1*sigmw

 w2=g1*s/2.0-sqrt(g1*g1*s*s+4.0/g2)/2.0
 w2=w2*sigmw

 a1=w2/(w2-w1)
 a2=-w1/(w2-w1)

 sigm1=r*w1
 sigm2=-r*w2


a=fltarr(2)+0D
sigm=fltarr(2)+0D
w=fltarr(2)+0D

a(0)=a1 & a(1)=a2

sigm(0)=sigm1 & sigm(1)=sigm2

w(0)=w1 & w(1)=w2

;width=coef1*(z/zi)^(coef2)
;;;;;;  z0=1.0D-5 ; assume independent of z0
;  rmo=-0.05D
;;  coef1_pdf=1.675D
;;  coef2_pdf=0.968D


;   coef1_ran=(0.4575-4.989*rmo) 
;   coef2_ran=(0.7333+6.377*rmo) 
   coef1_ran=-17.569*rmo^2-5.8085*rmo+0.4506
   coef2_ran=34.45*rmo^2+8.5801*rmo+0.7602

; position =coef3*(z/zi)^(coef4)
;  coef3_ran=0.572-(3.231)*rmo 
;  coef4_ran=0.70D 
  coef3_ran=-10.701*rmo^2-4.1947*rmo+0.5565
  coef4_ran=-35.002*rmo^2-1.7016*rmo+0.6828

  coef3_pdf=1.2305D

  coef4_pdf=0.993D

;
; maximal value=coef5*(z/zi)^(coef6)
;   coef5_ran=1.2575+10.528*rmo 
;  coef6_ran=-0.6877-2.81*rmo 
   coef5_ran=314.35*rmo^3+147.31*rmo^2+18.212*rmo+1.3414
   coef6_ran=-3.0408*rmo-0.6937

    pw=0.0D
    factor=1.0

    width_pdf=-5.4377*nz^5+17.058*nz^4-17.822*nz^3+5.7192*nz^2+1.131*nz  ; S=0.5 
    width_ran=coef1_ran*nz^coef2_ran
; new parameters for width
     p0=0.4114 & p1=-7.0704 & p2=12.5299  & p3=959.6516 & p4=0.6666
     p5=5.0849 & p6=97.2557 & p7=-62.7665 & p8=-30.4509 & p9=0.0
     p10=0.0
     coef_aa=p0+p1*rmo+p2*z0+p3*rmo*z0+p8*rmo*rmo+p9*z0*z0
     coef_bb=p4+p5*rmo+p6*z0+p7*rmo*z0
     width_ran=1d*coef_aa*(nz^coef_bb)+p10*(nz^3)
;end of new for width

    width=width_pdf/width_ran

    position_pdf=coef3_pdf*(nz^coef4_pdf) ;s=0.5 
    position_ran=1d*coef3_ran*nz^coef4_ran 


;   position_pdf=0.4762*nz^4-1.9972*nz^3+1.2424*nz^2+1.3428*nz
; new parameters for position
     p0=0.6145 & p1=-1.9699 & p2=-9.7974  & p3=494.1766 & p4=0.7830
     p5=2.4452 & p6=13.4432 & p7=-915.30331 & p8=7.8899577  & p9=0.0
     p10=0.0
     coef_aa=p0+p1*rmo+p2*z0+p3*rmo*z0+p8*rmo*rmo+p9*z0*z0
     coef_bb=p4+p5*rmo+p6*z0+p7*rmo*z0
     position_ran=1d*coef_aa*(nz^coef_bb)+p10*(nz^3)
;end of new for position

   position=position_pdf/position_ran


   maximal_ran=coef5_ran*nz^coef6_ran 
   maximal_pdf=0.3766*nz^(-1.0028) ;s=0.5 
; new parameters for Maxi
     p0=1.3272 & p1=11.244429 & p2=11.534338  & p3=-500.2998 & p4=-0.694832
     p5=-4.7810256 & p6=-45.592302 & p7=679.84818 & p8=64.237277 & p9=0.0
     p10=-0.95424732
     coef_aa=p0+p1*rmo+p2*z0+p3*rmo*z0+p8*rmo*rmo+p9*z0*z0
     coef_bb=p4+p5*rmo+p6*z0+p7*rmo*z0
     maximal_ran=1d*coef_aa*(nz^coef_bb)+p10*(nz^3)
;end of new for Max



   maximal=maximal_ran/maximal_pdf 

   factor=0.5*width+0.5*maximal
;  factor(iz)=2.5
  adjustx=xp*factor+(position_pdf-factor*position_ran)
    IF adjustx LT 0 THEN BEGIN & ftz=0.0 & return, ftz & ENDIF
    FOR j=0,1 DO BEGIN
     FOR k=-10D,10 DO BEGIN

      ww2=(nz-nzs+2.0D*k)*wstar/adjustx
      ww1=(-nz-nzs+2.0D*k)*wstar/adjustx

       e1=-(ww1-w(j))^2/(2D*sigm(j)*sigm(j))
       e2=-(ww2-w(j))^2/(2D*sigm(j)*sigm(j))
       exp1=0.0D
       exp2=0.0D
       IF abs(e1) LE 150 THEN $
       exp1=exp(e1)
       IF abs(e2) LE 150 THEN $
         exp2=exp(e2)
        mid=a(j)/(sqrt(2.0D*3.14)*sigm(j))*(-ww1*exp1+ww2*exp2)*wstar/adjustx/wstar
      pw=pw+mid

     ENDFOR
    ENDFOR
contd:
    ftz=factor*pw ; dimensionless, wstar/ndx ~ u*zi/x

return,ftz




end