program reformat_wf c****************************************************************************** c c Reformat wildfire from the txt file. c c Written by: Hongliang Zhang(Feb.,2011) c c****************************************************************************** implicit none c ioapi headers INCLUDE 'PARMS3.EXT' ! I/O parameters definitions INCLUDE 'FDESC3.EXT' ! file header data structure INCLUDE 'IODECL3.EXT' ! I/O definitions and declarations c name of the program CHARACTER*16 pname DATA pname/'extract_data'/ INTEGER LOGUNIT CHARACTER*16 Sfenv, sfenv2,sname(100),sfenv3,sfenv4 CHARACTER*256 Sfile, sfout,sftatout,sfile1,sfile2,msg,sfpbl, + sfht LOGICAL, external :: SETENVVAR integer n_out_col,n_out_row,ncf_day real,allocatable :: rdata1(:,:,:,:,:) real,allocatable :: rdata2(:,:,:,:) real, allocatable :: rdpbl(:,:) real, allocatable :: rdataht(:,:,:,:) real, allocatable :: rdht(:,:,:) real, allocatable :: rdatapbl(:,:,:) character*40 cdata(100) integer nx1,nx2,ny1,ny2,ipos(100),ncol,ncol1,ihr,ihr_mm integer t, itime, idate, ierr,iday,dx,dy,step,nlay_met integer i,j,k,n,iv,i1,j1,k1, itime1, idate1,i2,i3,m,iloc real xorg, yorg real vtot,area real ptopmax,pbotmax,ptophour,pbothour integer idx,idy,npositive,out_nlays3d character*16 svname(100),s16 character*1 s1 character*2 ss CHARACTER*2048 stemp real lati,longi,rtemp1,rtemp2,dayprof(24),rtemp3 real ht_temp, r_ht(0:20),frac_ht(20,24),rec_pbl real dayprof1(24),dayprof2(24),dayprof3(24),dayprof4(24) real behour(24),besize,behour0(24),dayprof0(24) logical ifasia c data dayprof / 0.0043, 0.0043, 0.0043, 0.03, 0.06, 0.1, c + 0.14, 0.17, 0.14, 0.12, 0.09, 0.06, c + 0.03, 0.0043, 0.0043,0.0043, 0.0043, 0.0043, c + 0.0043,0.0043, 0.0043,0.0043, 0.0043, 0.0043/ c for china c data behour for local data behour0 / + 0.03, 0.03, 0.03, 0.03,0.03,0.03, + 0.06, 0.10,0.2, 0.4, 0.7,0.8, + 0.9, 0.95, 0.99, 0.8, 0.7, 0.4, + 0.2, 0.1, 0.06, 0.03, 0.03, 0.03 / c + 0.06, 0.03, 0.03, 0.03, 0.03, 0.03 / c data behour / c + 0.06, 0.10, 0.2, 0.4, c + 0.7,0.8, 0.9, 0.95, 0.99, 0.8, c + 0.7, 0.4, 0.06, 0.03, 0.03, 0.03, c + 0.03, 0.03, 0.03, 0.03, 0.03, 0.03,0.03,0.03 / c from Table 4.1 of FINAL REPORT-1996 Fire EMissions Inventory for the c WRAP Region-Methodology and Emission Estimates data dayprof0 /0.0053, 0.0053, 0.0053, 0.0053, 0.0053, 0.0053, + 0.0053, 0.0053, 0.0053, 0.02, 0.04, 0.07, + 0.1, 0.13, 0.16, 0.2047, 0.12, 0.07, + 0.0053, 0.0053, 0.0053, 0.0053, 0.0053,0.0053 / c + 0.1, 0.13, 0.16, 0.17, 0.12, 0.07, c + 0.04, 0.0053, 0.0053, 0.0053, 0.0053,0.0053 / c domain and file descriptions INTEGER NCOLS INTEGER NROWS INTEGER NLAYS INTEGER NTHIK INTEGER NVARS INTEGER FTYPE INTEGER FSIZE REAL*8 P_ALP ! first, second, third map REAL*8 P_BET ! projection descriptive REAL*8 P_GAM ! parameters. REAL*8 XCENT ! lon for coord-system X=0 REAL*8 YCENT ! lat for coord-system Y=0 REAL*8 XORIG ! X-coordinate origin of grid (map units) REAL*8 YORIG ! Y-coordinate origin of grid REAL*8 XCELL ! X-coordinate cell dimension REAL*8 YCELL ! Y-coordinate cell dimension INTEGER VGTYP ! vertical coordinate type (VGSIGP3, ...) REAL VGTOP ! model-top, for sigma coord types. REAL VGLVS( MXLAYS3 + 1 ) ! vertical coord values. CHARACTER*16 GDNAM ! grid name (length NAMLEN3=16) INTEGER SDATE, STIME, TSTEP, DURATN, NSTEPS INTEGER JDATE, JTIME, JDATE2 real sp1, sp2, phi0,lon0 c variables to setup the output file c output variables integer ii,dday,is double precision x,y character*16 aero6(16) data aero6 /'PNCOM', 'PNA', 'PMG', 'PAL', 'PSI', 'PCL', + 'PK', 'PCA', 'PTI', 'PMN', 'PFE', 'PNH4', + 'PMOTHR','PH2O','PNO3', 'PSO4' / real faero6(16) data faero6 /0.2026, 0.004093, 0.0003733,0.0202,0.03069,0.04421, + 0.03481,0.0129,0.001461,0.00009967,0.009911,0.01009, + 0.207, 0.00,0.001797, 0.0437/ c sp1=30. c sp2=60. c phi0=45. c lon0=-97. c rtemp1=0. c read in read(*,*) sp1,sp2,phi0,lon0 print*,sp1,sp2,phi0,lon0 read(*,*) xorg, yorg read(*,*) idx,idy read(*,*)dday read(*,*)ncf_day read(*,'(A)')sfile1 read(*,'(A)')sfpbl read(*,'(A)')sfht read(*,'(A)')sfile read(*,'(A)')sfout c**************need modify if for other cases other than North American and China ifasia=.false. if(lon0.gt. 80.and. lon0.lt.130) ifasia=.true. print*,'ifasia',ifasia c c start IOAPI LOGUNIT = INIT3() c set input filename to enviromental variable SFENV='SFILE' SFENV2='SFOUT' sfenv3='sfpbl' sfenv4='sfht' c read in PBL from netcdf file IF ( .NOT.SETENVVAR(Sfenv3,sfpbl) ) & CALL M3EXIT(pname,0,0,'FAIL TO SET INPUT FILE',-1) IF ( .NOT.OPEN3(SFENV3,FSREAD3,pname) ) & CALL M3EXIT(pname,0,0,'FAIL TO OPEN INPUT FILE',-1) IF ( .NOT.DESC3(SFENV3) ) & CALL M3EXIT(pname,0,0,'FAIL TO GET FILE DESC',-1) c allocate memory for the average variable allocate(rdatapbl(ncols3d,nrows3d,mxrec3d), + stat=ierr) if (ierr.ne.0) & CALL M3EXIT(pname,0,0,'FAIL TO ALLOCATE MEMORY FOR rdata',-1) allocate(rdpbl(ncols3d,nrows3d), + stat=ierr) if (ierr.ne.0) & CALL M3EXIT(pname,0,0,'FAIL TO ALLOCATE MEMORY FOR rdata',-1) rdatapbl=0. jdate=sdate3d jtime=stime3d tstep=tstep3d c .. loop over all time steps do step=1, mxrec3d print*,'time step',jdate,jtime,tstep IF ( .NOT.READ3(SFENV3,'PBL',1,jdate,jtime,rdpbl)) + CALL M3EXIT (pname,0,0,'FAIL TO READ DATA',-1) do i=1, ncols3d do j=1, nrows3d rdatapbl(i,j,step)=rdpbl(i,j) enddo enddo c .. advance time call nextime( jdate, jtime, tstep) print*,'after nextime',jdate,jtime enddo if (.not.close3(sfenv3)) + call m3exit(pname,idate,itime,'Cannot write data', -1) c closed met file after reading all PBL data to rdatapbl IF ( .NOT.SETENVVAR(Sfenv4,sfht) ) & CALL M3EXIT(pname,0,0,'FAIL TO SET INPUT FILE',-1) IF ( .NOT.OPEN3(SFENV4,FSREAD3,pname) ) & CALL M3EXIT(pname,0,0,'FAIL TO OPEN INPUT FILE',-1) IF ( .NOT.DESC3(SFENV4) ) & CALL M3EXIT(pname,0,0,'FAIL TO GET FILE DESC',-1) c allocate memory for the average variable allocate(rdht(ncols3d,nrows3d,mxrec3d), + stat=ierr) if (ierr.ne.0) & CALL M3EXIT(pname,0,0,'FAIL TO ALLOCATE MEMORY FOR rdata',-1) allocate(rdataht(ncols3d,nrows3d,nlays3d,mxrec3d), + stat=ierr) if (ierr.ne.0) & CALL M3EXIT(pname,0,0,'FAIL TO ALLOCATE MEMORY FOR rdata',-1) rdataht=0. jdate=sdate3d jtime=stime3d tstep=tstep3d out_nlays3d=min(12,nlays3d) c .. loop over all time steps do step=1, mxrec3d print*,'time step',jdate,jtime,tstep IF ( .NOT.READ3(SFENV4,'ZF',ALLAYS3,jdate,jtime,rdht)) + CALL M3EXIT (pname,0,0,'FAIL TO READ DATA',-1) do i=1, ncols3d do j=1, nrows3d do k=1,nlays3d rdataht(i,j,k,step)=rdht(i,j,k) enddo enddo enddo c .. advance time call nextime( jdate, jtime, tstep) print*,'after nextime',jdate,jtime enddo if (.not.close3(sfenv4)) + call m3exit(pname,idate,itime,'Cannot write data', -1) c closed met file after reading all layer height data to rdataht IF ( .NOT.SETENVVAR(Sfenv,Sfile) ) & CALL M3EXIT(pname,0,0,'FAIL TO SET INPUT FILE',-1) IF ( .NOT.OPEN3(SFENV,FSREAD3,pname) ) & CALL M3EXIT(pname,0,0,'FAIL TO OPEN INPUT FILE',-1) IF ( .NOT.DESC3(SFENV) ) & CALL M3EXIT(pname,0,0,'FAIL TO GET FILE DESC',-1) c c nlays3d=nlays3d+2 ! increase the layers we want to calculate c VGLVS3d(nlays3d)=0.85 c VGLVS3d(nlays3d+1)=0.775 sdate3d=ncf_day nlays3d=out_nlays3d c adding species if not exists iloc=0 call slocate(vname3d, nvars3d, 'CCO_OH', iloc) if (iloc.eq.0) then nvars3d=nvars3d+1 units3d(nvars3d)='moles/s' vname3d(nvars3d) ='CCO_OH' vtype3d(nvars3d)=M3REAL VDEsc3d(nvars3d) = 'Model species CCO_OH' endif call slocate(vname3d, nvars3d, 'HCOOH', iloc) if (iloc.eq.0) then nvars3d=nvars3d+1 units3d(nvars3d)='moles/s' vname3d(nvars3d) ='HCOOH' vtype3d(nvars3d)=M3REAL VDEsc3d(nvars3d) = 'Model species HCOOH' endif call slocate(vname3d, nvars3d, 'RNO3', iloc) if (iloc.eq.0) then nvars3d=nvars3d+1 units3d(nvars3d)='moles/s' vname3d(nvars3d) ='RNO3' vtype3d(nvars3d)=M3REAL VDEsc3d(nvars3d) = 'Model species RNO3' endif call slocate(vname3d, nvars3d, 'PH2O', iloc) if (iloc.eq.0) then nvars3d=nvars3d+1 units3d(nvars3d)='g/s' vname3d(nvars3d) ='PH2O' vtype3d(nvars3d)=M3REAL VDEsc3d(nvars3d) = 'Model species PH2O' endif c allocate memory for the average variable allocate(rdata1(ncols3d,nrows3d,nlays3d,nvars3d,mxrec3d), + stat=ierr) if (ierr.ne.0) & CALL M3EXIT(pname,0,0,'FAIL TO ALLOCATE MEMORY FOR rdata',-1) rdata1=0. allocate(rdata2(ncols3d,nrows3d,nlays3d,nvars3d), + stat=ierr) if (ierr.ne.0) & CALL M3EXIT(pname,0,0,'FAIL TO ALLOCATE MEMORY FOR rdata',-1) rdata1=0. print*,'memory allocated for rdata1,rdata2', + ncols3d,nrows3d,nlays3d,nvars3d,mxrec3d IF ( .NOT.SETENVVAR(sfenv2,sfout) ) & CALL M3EXIT(pname,0,0,'FAIL TO SET INPUT FILE',-1) c open output file IF ( .NOT.OPEN3(SFENV2,FSUNKN3,pname) ) & CALL M3EXIT(pname,0,0,'FAIL TO OPEN INPUT FILE',-1) c get file description IF ( .NOT.DESC3(SFENV2) ) & CALL M3EXIT(pname,0,0,'FAIL TO GET FILE DESC',-1) c close template file if (.not.close3(sfenv)) + call m3exit(pname,idate,itime,'Cannot write data', -1) print*,'output file opened' open(unit=1,file=sfile1,status='old') call GETNCOL(1,ncol1,ipos,ierr) READ (1,'(A)') stemp print*,'first line:',trim(stemp) do i=1,ncol1 i1=ipos(i)+1 i2=ipos(i+1)-1 if (i2.ge.i1) then read (stemp(i1:i2),'(a)') sname(i) else print*,'error',i,i1,i2 stop endif enddo c adjust variables name c sname(34)='MACR' c sname(44)='POC' c sname(45)='PEC' sname(35)='MACR' sname(45)='POC' sname(46)='PEC' call GETNCOL(1,ncol,ipos,ierr) if(ncol.ne.ncol1.and.ncol.ne.46)then print*,'not consistent',ncol,ncol1 stop endif 20 READ (1,'(A)',end=90,iostat=ierr) stemp do i=1,ncol i1=ipos(i)+1 i2=ipos(i+1)-1 if (i2.ge.i1) then read (stemp(i1:i2),'(a)') cdata(i) else print*,'error',i,i1,i2 stop endif enddo c write cdata to integer and real,if the day we need? read(cdata(1),'(a)')iday if(iday.eq.dday)then c judge the domain read(cdata(4),'(F30.10)')lati read(cdata(5),'(F30.10)')longi call ll2lam(sp1,sp2,phi0,lon0,lati,longi,x,y) c dx=int(x+2898000)/36000+1 c dy=int(y+2787771)/36000+1 dx=int(x-xorg)/idx+1 dy=int(y-yorg)/idy+1 if(dx.gt.0.and.dx.lt.ncols3d.and. + dy.gt.0.and.dy.lt.nrows3d)then ! in the domain print*,'find recond for',iday,dx,dy c decide the hourly, no hour information needed c read(cdata(2),'(I)')ihr_mm c ihr=int(ihr_mm/100)+1 c decide the time zone if(ifasia)then do j=1,24 behour(j)=behour0(mod(j+7,24)+1) dayprof(j)=dayprof0(mod(j+7,24)+1) enddo else if(longi.lt.-127.5)then ! West US, +7 do j=1,24 behour(j)=behour0(mod(24+j-6,24)+1) dayprof(j)=dayprof0(mod(24+j-6,24)+1) enddo else if(longi.ge.-127.5.and.longi.lt.-112.5)then ! Mon US, +6 do j=1,24 behour(j)=behour0(mod(24+j-5,24)+1) dayprof(j)=dayprof0(mod(24+j-5,24)+1) enddo else if(longi.ge.-112.5.and.longi.lt.-97.5)then ! Mid, +5 do j=1,24 behour(j)=behour0(mod(24+j-4,24)+1) dayprof(j)=dayprof0(mod(24+j-4,24)+1) enddo else ! east US, +4 do j=1,24 behour(j)=behour0(mod(24+j-3,24)+1) dayprof(j)=dayprof0(mod(24+j-3,24)+1) enddo endif endif c********************************** c decide the fraction for each layers c read in area read(cdata(6),'(F30.10)')area area=area/4046. ! converted m2 to arce if(area.lt.10.) then pbotmax=10. ptopmax=160. besize=0.4 else if(area.lt.100.)then pbotmax=900. ptopmax=2400. besize=0.6 else if(area.lt.1000.)then pbotmax=2200. ptopmax=6400. besize=0.75 else if(area.lt.5000.)then pbotmax=3000. ptopmax=7200. besize=0.85 else pbotmax=3000. ptopmax=8000. besize=0.90 endif do j=1,24 c rec_pbl=rdatapbl(dx,dy,j) ptophour=(behour(j)**2)*(besize**2)*ptopmax pbothour=(behour(j)**2)*(besize**2)*pbotmax print*,j,ptopmax,pbotmax,ptophour,pbothour r_ht(0)=0. do i=1,nlays3d r_ht(i)=rdataht(dx,dy,i,j) enddo do i=1,nlays3d if(r_ht(i).lt.pbothour)then ! lower than fire frac_ht(i,j)=(1.-besize)*(r_ht(i)-r_ht(i-1))/pbothour c print*,'lower than fire' else if(r_ht(i).gt.ptophour.and.r_ht(i-1).lt.pbothour)then ! within one layer frac_ht(i,j)=besize+(1.-besize)*(pbothour-r_ht(i-1))/pbothour c print*,'include the fire' else if(r_ht(i).lt.ptophour.and.r_ht(i-1).gt.pbothour)then ! within file frac_ht(i,j)=besize*(r_ht(i)-r_ht(i-1))/(ptophour-pbothour) else if (r_ht(i-1).ge.ptophour) then ! above fire frac_ht(i,j)=0. else if (r_ht(i).lt.ptophour.and.r_ht(i).gt.pbothour .and. + r_ht(i-1).lt.pbothour) then !! include lower part frac_ht(i,j)=besize*(r_ht(i)-pbothour)/(ptophour-pbothour) + + (1.-besize)*(pbothour-r_ht(i-1))/pbothour else if (r_ht(i-1).lt.ptophour.and.r_ht(i).gt.ptophour.and. + r_ht(i-1).gt.pbothour) then !! frac_ht(i,j)=besize*(ptophour-r_ht(i-1))/(ptophour-pbothour) else print*,'wrong in vertical distribution',i,j frac_ht(i,j)=0. c stop endif if(i.eq.nlays3d) then if(r_ht(i).lt.ptophour) then if(r_ht(i).lt.pbothour)then frac_ht(i,j)=1. else frac_ht(i,j)=frac_ht(i,j)+besize* + (ptophour-r_ht(i))/(ptophour-pbothour) endif endif endif c print*,i,j,r_ht(i),r_ht(i-1) c print*,pbothour,ptophour,besize,frac_ht(i,j) enddo vtot=0. do i=1,nlays3d c print*,'frac_ht',i,frac_ht(i,j) vtot=vtot+frac_ht(i,j) enddo if( vtot.lt.0.95 .and. vtot.gt.1.01) then print*,'wrong in vertical distribution2',j,vtot stop endif enddo c print*,frac_ht c decide species m/day to m/s c do ii=8,43 do ii=9,44 call slocate(vname3d, nvars3d, sname(ii), iloc) if (iloc.gt.0) then read(cdata(ii),'(F30.10)')rtemp1 rtemp1=rtemp1/3600 do i=1,nlays3d do j=1,24 rdata1(dx,dy,i,iloc,j)=rdata1(dx,dy,i,iloc,j)+ + rtemp1*frac_ht(i,j)*dayprof(j) c print*,'frac_ht',i,j,frac_ht(i,j),dayprof(j), c + rdata1(dx,dy,i,iloc,j) enddo if(i.eq.nlays3d) then c stop endif enddo else print*,'gas species not found',ii,sname(ii) endif enddo c for oc ec kg/d to g/s do ii=45,46 call slocate(vname3d, nvars3d, sname(ii), iloc) if (iloc.gt.0) then read(cdata(ii),'(F30.10)')rtemp1 endif rtemp1=rtemp1*1000/3600 do i=1,nlays3d do j=1,24 rdata1(dx,dy,i,iloc,j)=rdata1(dx,dy,i,iloc,j)+ + rtemp1*frac_ht(i,j)*dayprof(j) enddo enddo enddo c*************************************************************** c PNCOM, PNA, PMG, PAL,PSI, PCL, PK, PCA, PTI, PMN, PFE, PNH4, PMOTHR, c PH2O,PNO3,PSO4 do is=1,16 print*,aero6(is) call slocate(vname3d,nvars3d,trim(aero6(is)),iloc) if(iloc.gt.0)then read(cdata(47),'(F30.10)')rtemp1 do i=1,nlays3d do j=1,24 rdata1(dx,dy,i,iloc,j)=rdata1(dx,dy,i,iloc,j)+ + rtemp1*faero6(is)*frac_ht(i,j)*dayprof(j)/3.6 ! fso4 enddo enddo else print*,'species not found',is,aero6(is) do i=1,nvars3d print*,i,vname3d(i) enddo stop endif enddo call slocate(vname3d, nvars3d,'PMFINE', iloc) if(iloc.gt.0)then read(cdata(47),'(F30.10)')rtemp1 read(cdata(45),'(F30.10)')rtemp2 read(cdata(46),'(F30.10)')rtemp3 print*,rtemp1,rtemp2,rtemp3 print*,rtemp1*(1.-faero6(15))-rtemp2-rtemp3 do i=1,nlays3d do j=1,24 rdata1(dx,dy,i,iloc,j)=rdata1(dx,dy,i,iloc,j)+ + max(rtemp1*(1.-faero6(15)-faero6(16)-faero6(12)) + -rtemp2-rtemp3,0.) + *frac_ht(i,j)*dayprof(j)/3.6 enddo enddo else print*,'no PMFINE is found in vanme3d', vname3d stop endif c************************************* call slocate(vname3d, nvars3d,'PMC', iloc) if(iloc.gt.0)then read(cdata(47),'(F30.10)')rtemp1 read(cdata(48),'(F30.10)')rtemp2 do i=1,nlays3d do j=1,24 rdata1(dx,dy,i,iloc,j)=rdata1(dx,dy,i,iloc,j)+ + (rtemp2-rtemp1)*frac_ht(i,j)*dayprof(j)/3.6 enddo enddo endif endif endif go to 20 90 continue close(1) c get file description IF ( .NOT.DESC3(SFENV2) ) & CALL M3EXIT(pname,0,0,'FAIL TO GET FILE DESC',-1) jdate=sdate3d jtime=stime3d tstep=10000 c .. loop over all time steps do step=1, 25 print*,'looping, step is',tstep, 'of nsteps=',25 print*,'time step',jdate,jtime,tstep npositive=0 do i=1, ncols3d do j=1,nrows3d do k=1,nvars3d do ii=1,nlays3d if(isnan(rdata1(i,j,ii,k,step)))rdata1(i,j,ii,k,step)=0. rdata2(i,j,ii,k)=rdata1(i,j,ii,k,step) if(ii.gt.1)then if(rdata1(i,j,ii,k,step).gt.0.)npositive=npositive+1 endif if(step.eq.25)then rdata2(i,j,ii,k)=rdata1(i,j,ii,k,1) endif enddo enddo enddo enddo if(npositive.gt.0.)print*,'npositive great than 0',npositive if (.not.write3(sfenv2,'ALL',jdate,jtime,rdata2) ) & call m3exit(pname,jdate,jtime,'Cannot write data sf2', -1) c .. advance time call nextime( jdate, jtime, tstep) print*,'after nextime',jdate,jtime enddo if (.not.close3(sfenv2)) + call m3exit(pname,idate,itime,'Cannot write data', -1) end