; this code is used to compute footprint for a year using CBL footprit ; model tower='wlef' print,'input yr, yyyy' yr='' read,yr print,'input use guessed zo 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 print,'input height, 122, 396' level='' read,level badval=-999.0 type='10' z0=0.9d ; roughness length ;displace=0.6*15.0 ; 2/3 canopy height, at WLEF, canopy is roughly 15m displace=15.0 ; these two values for CBL modeling z0_forest=0.9d displace_forest=15.0 width=1.0 ; sigm=width*sigm ; default=1.0, 5/13/2005 IF guess EQ 'y' THEN BEGIN IF guessnumber EQ '0' THEN BEGIN z0_surf=0.4d displace_surf=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_surf=0.3d displace_surf=8d;m subname='footprintcases/guess'+guessnumber ENDIF IF guessnumber EQ '2' THEN BEGIN z0_surf=0.2d displace_surf=6d;m subname='footprintcases/guess'+guessnumber ENDIF IF guessnumber EQ '3' THEN BEGIN z0_surf=0.1d displace_surf=4d;m subname='footprintcases/guess'+guessnumber ENDIF ; experiments about the effects of horizontal diffusion coefficient ; use guess2 case IF guessnumber EQ '4' THEN BEGIN z0_surf=0.2d displace_surf=6d;m subname='footprintcases/guess'+guessnumber width=0.5 ENDIF IF guessnumber EQ '5' THEN BEGIN z0_surf=0.2d displace_surf=6d;m subname='footprintcases/guess'+guessnumber width=2.0 ENDIF ENDIF ELSE BEGIN z0_surf=0.5d displace_surf=15.0 ; control case subname='footprintcases/control' ENDELSE ;hor_coef=1d*fix(type)/10.0 hor_coef=1d*fix(type)/10.0*width ; 5/13/2005 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 to make sure tower is at the center of the map ;transform coordiant into one centered at the tower x=x-xc y=y-yc ; open output file outfilename='/eddy/s3/wang/thesis/ft/decompose_flux/results/vegratio_'+tower+yr+'_'+level+'.txt'+type IF guess EQ 'y' THEN outfilename='/eddy/s3/wang/thesis/ft/decompose_flux/results/'+subname+'/vegratio_'+tower+yr+'_'+level+'.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 zi, ;load,'/eddy/s3/wang/thesis/ablbudget/data/ruc/'+yr+'/'+yr+'wlefruc_pbl32.txt',pbl ; UTC load,'/eddy/s3/wang/thesis/ablbudget/data/ml/wlef_ml.dat'+yr,pbl ;LST zi1=pbl(12,*) zi=zi1*0.0+badval IF yr EQ '1998' OR yr EQ '1999' THEN BEGIN load,'/eddy/s3/wang/thesis/ablbudget/data/ml/'+yr+'ml_wlefobslst.txt',pbl2 ; radar obs ; lst zi=pbl2(5,*) ENDIF r=where(zi EQ -999) zi(r)=zi1(r) ; use directly obs plus emipirical computation ; load met file ;load,'/eddy/s3/wang/thesis/ablbudget/data/WLEF'+yr+'_met.txt',met ;LST 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) IF level EQ '396' THEN index=13 IF level EQ '122' THEN index=14 u=met(index,*) ; wind speed; m/s wdir=met(index+3,*) ; wind direction 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 ;load,'/eddy/s3/wang/thesis/ablbudget/data/WLEF'+yr+'_flux.txt',flx ;LST 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 IF level EQ '122' THEN co2=flx(13,*); co2 flux(NEE) at 122m IF level EQ '396' THEN co2=flx(14,*); co2 flux(NEE) at 396m mo=ustar*0.0+badval wstar=mo 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) ; 7/6/2004, change heat ne badval in to heat >0 r=where(heat GE 0.0 AND temp NE badval AND zi NE badval,count) IF count GE 1 THEN wstar(r)=(9.8*(heat(r)/1004.0/1.29)*zi(r)/(temp(r)+273.15))^(0.3333) ;compute wind speed at about 100m (or 0.1 zi) FOR i=0,len(hr)-1 DO BEGIN IF mo(i) NE badval AND mo(i) LE 0 THEN BEGIN ; abs(mo) should be << z0 mo(i)=mo(i)<(-20.0*z0) hh=(0.1*zi[i])>100 ya=(1-16*hh/mo(i))^(0.25) psim=2*alog((1+ya)/2.0)+alog((1+ya*ya)/2)+2*atan((1-ya)/(1+ya)) utheory=(ustar(i)/0.4)*(alog(hh/z0)-psim) ;print,hh,ustar(i),psim,wdir(i),u(i),utheory IF u(i) EQ badval OR u(i) EQ 0 THEN u(i)=utheory IF wdirmean(i) EQ -999 THEN u(i)=badval ENDIF ENDFOR ;-------------------------------------------------------------------- ; 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 a95=badval s50=badval s95=badval ratio=fltarr(12)+badval z0=z0_forest displace=displace_forest ; take vales for forest as default; if within surface layer, the values will change ;judge the conditions to select which ft model is used. IF wdirmean(k) NE badval AND u(k) NE badval AND wstar(k) NE badval AND mo(k) NE badval AND zi(k) GE 400.0 AND float(level)/zi(k) LT 0.5 AND abs(mo(k)/zi(k)) LE 0.125 AND mo(k) LE 0 THEN BEGIN ; only consider unstable condition rmo=mo(k)/zi(k) rz0=z0/zi(k) ; rz0 should be in the range of fiitting 3e-5 to 2e-3 rz0=rz0>3e-5 rz0=rz0<2e-3 zp=1.0D*(1.0D*fix(level)-displace)/zi(k) print,'sigmY width=',width print,jday(k),hr(k),zp,wdirmean(k),rmo,mo(k),zi(k),co2(k),'rz0=',rz0,wstar(k),u(k) ; use forest z0 and d to judge if the measurement can be within ; surface layer (approximation, no big deal) surface_flag=0 IF zp LE 0.1 THEN BEGIN ;use surface layer model for zp<=0.1 IF guess EQ 'y' THEN BEGIN ; guess z0=z0_surf displace=displace_surf ENDIF ;guess z_x,z0,mo(k),xplume,zplume surface_flag=1 print,'zp=',zp,' using surface model', 'zo=',z0 ENDIF ; use surface layer model pztemp=fltarr(nworkx) FOR i=0,nworkx-1 DO BEGIN ; dimensionless horiaontal distance nxw nxw=abs(xw(i))*wstar(k)/(u(k)*zi(k))+1.0e-10 ;calculate ft_z ; crosswind-integral ft pz=0.0D IF surface_flag EQ 1 THEN BEGIN ; surface model is used IF guess EQ 'y' THEN BEGIN ; guess z0=z0_surf displace=displace_surf ENDIF ;guess pz=fy(mo(k),z0,xplume,zplume,abs(xw(i)),[float(level)-displace]) ; 1/m unit ; pz=pz/(wstar(k)/(u(k)*zi(k))) ; dimensionless ;; sy=hor_coef*(abs(1.3/mo(k))+0.1)*abs(xw(i))/sqrt(1.0+0.0001*abs(xw(i)))+1.0 ;min sy=1 ; nsy=sy*wstar(k)/(u(k)*zi(k)) ; 8/31/04 to be consistent, use Briggs formula for urban conditions sy=hor_coef*(abs(4.0/mo(k))+0.16)*abs(xw(i))/sqrt(1.0+0.0004*abs(xw(i)))+1.0 ENDIF ELSE BEGIN ; CBL model IF nxw LE 4 THEN pz=(wstar(k)/(u(k)*zi(k)))*cblfootprintz(nxw,zp,rmo,rz0) ;1/m nsy=hor_coef*0.6*nxw*wstar(k)/u(k) ; sy=0.6*wstar*x/u p 88 of book Lecture on air pollution modeling; 11/18 experiments times 2 ; nsy=sy*wstar(k)/(u(k)*zi(k)) ; non-dimensional ; print,i,nxw,zp,rmo,pz sy=0.6*wstar(k)*(abs(xw(i)))/u(k)+1.0 ; 8/31/04 to be consistent, use Briggs formula for urban conditions ; sy=hor_coef*(abs(1.3/mo(k))+0.16)*(abs(xw(i)))/sqrt(1.0+0.0004*(abs(xw(i)))) ENDELSE ; CBL model ; sy=sy>1.0 pztemp(i)=pz FOR j=0,nworky-1 DO BEGIN ; crosswind distance, y -direction ; nyw=abs(yw(j))*wstar(k)/(u(k)*zi(k)) py=0.0D ; expnyw=1.0D*nyw*nyw/(2.0*nsy*nsy) expnyw=1.0d*yw(j)*yw(j)/(2*sy*sy) ; 7/6/2004 , change i ne 0 in to i ne nworkx-1 IF expnyw LE 150 AND i NE nworkx-1 THEN py=exp(-expnyw)/(sqrt(6.28)*sy) ; work(i,j)=py*(wstar(k)/(u(k)*zi(k)))*pz*(wstar(k)/(u(k)*zi(k))) ; ; 1/m2 unit work(i,j)=py*pz ;; printf,20,i,j,pz,py,work(i,j) ENDFOR ENDFOR ; IF surface_flag EQ 0 THEN plot,xw,pztemp,title=string(jday(k))+string(hr(k))+'CBL' ; IF surface_flag EQ 1 THEN plot,xw,pztemp,title=string(jday(k))+string(hr(k))+'SURFACE' ;contour,work,xw/1000.,yw/1000. IF surface_flag EQ 0 THEN print,'=========CBL model is used=========z0=,d=',z0,displace ELSE print,'=======Surface model is used=====z0=,d=',z0,displace print,'--',total(pztemp)*grid ; *wstar(k)/(u(k)*zi(k)) print,'--',total(work)*((grid)^2) IF total(work)*((grid)^2) GE 2 THEN BEGIN & print,'skip this hour' & GOTO,skip & ENDIF ; transform the results in the wind direction coordinate into real ; coordinate fttemp=fttemp*0.0 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,wdirmean(k),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-30 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,ix,iy)+1e-8 ; distance square d2=distance2(xwind,ywind,ix1,iy1)+1e-8 d3=distance2(xwind,ywind,ix2,iy2)+1e-8 d4=distance2(xwind,ywind,ix3,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 tot_temp=total(fttemp) IF tot_temp NE 0.0 THEN fttemp=fttemp/(tot_temp*grid*grid) ; a50=findcontour(fttemp,grid,0.5,s50) ; a90=findcontour(fttemp,grid,0.9,s90) ;; contour,fttemp,x/1000.,y/1000.,nlevels=50 ;,levels=[a90,a50] findvegarea_dnr,fttemp,veg,ratio ; print, 'ratio=',ratio ; print,'total fttemp=',total(fttemp) ; contour,fttemp ;; ft=ft+fttemp & number=number+1 If co2(k) ne badval then begin ; weighted ft ft=ft+fttemp; *abs(co2(k)) totalco2=totalco2+1.0; abs(co2(k)) endif ; weighted ft ENDIF skip: ; par(k)=meanbadval(met(36:38,k)) printf,14,met(0:5,k),a50,a95,s50,s95,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 footprint outgridname='/eddy/s3/wang/thesis/ft/decompose_flux/results/ft_long_cbl_grid'+yr+'_'+level+'.txt'+type IF guess EQ 'y' THEN outgridname='/eddy/s3/wang/thesis/ft/decompose_flux/results/'+subname+'/ft_long_cbl_grid'+yr+'_'+level+'.txt'+type openw,16,outgridname 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 '+outgridname end