PROGRAM readucddata !**************************************************************************** ! ! The purpose of this program is to extract a map of pollutants concentration ! from the UCD air quality model outputs ! ! Written by: Jianlin Hu (July 2013) ! University of Caifornia, Davis ! jjlhu@ucdavis.edu ! !**************************************************************************** implicit none c parameters for standard bins integer nsbins parameter (nsbins=33) real rsbins(nsbins) c parameters for domain dimensions integer ngmx,ngmy parameter (ngmx=110,ngmy=110) real map(ngmx,ngmy) ! data array for target species c data arrays integer nact,nalim,nblim,ntype,ntime parameter (nact=616,nalim=51,nblim=15,ntype=1,ntime=1) real gconc(ngmx,ngmy,nact) ! gas conc real aconc(ntype,ntime,nblim,nalim,ngmx,ngmy) ! pm conc real radp(ntype,ntime,nblim,ngmx,ngmy) ! radius real rnump(ntype,ntime,nblim,ngmx,ngmy) ! number conc c real conc0(nblim),conc1(nsbins),ddp(nsbins),dp0(nblim) c character*256 infile,outfile,sline character*16 polnam,gas_spec(nact),pm_spec(nalim) character*10 sunit,dummy,upcase c real radius1,radius2 real rsum,rsum2,rsumi c integer igoa,ierr,ipoln,idx,iwrite,locate integer ix,iy,it,itt,ib,inalim1,inalim2,ip c c -- standard bins data rsbins /1.29867E-08,2.36865E-08,4.22182E-08,7.48169E-08, + 1.33506E-07,2.1702E-07,2.88735E-07,3.41979E-07, + 3.99025E-07,4.68274E-07,5.52078E-07,6.39735E-07, + 7.13674E-07,7.65597E-07,8.09096E-07,8.56052E-07, + 9.05733E-07,9.69392E-07,1.06169E-06,1.17625E-06, + 1.30017E-06,1.4388E-06,1.60325E-06,1.79474E-06, + 2.01141E-06,2.25165E-06,2.51768E-06,2.89401E-06, + 3.5359E-06,4.45657E-06,5.61048E-06,7.07131E-06, + 8.91251E-06/ c read in information c-- input file READ (5,'(a)') infile c-- target PM size range READ (5,*) radius1 , radius2 c-- target species READ (5,'(a120)') sline IF ( sline(1:1).NE.'*' ) THEN READ (sline,'(a16,1x,a4,1x,a4)') polnam, sunit, dummy dummy = UPCASE(dummy) IF ( dummy(1:3).EQ.'GAS' ) THEN igoa = 1 ELSE igoa = 0 ENDIF ENDIF c-- output file READ (5,'(a)') outfile c init array DO ix = 1 , NGMX DO iy = 1 , NGMY map(ix,iy) = -999. ENDDO ENDDO c read in data CALL PREAD(infile,ngmx,ngmy,nact,nalim,nblim,ntype,ntime, + gconc,aconc,radp,rnump,gas_spec,pm_spec) IF ( igoa.EQ.1 ) THEN idx = locate(polnam,gas_spec,nact) c PRINT * , 'LOCATING GAS' , polnam(i) , ' AT' , idx ipoln = idx ELSEIF ( INDEX(polnam,'MASS').EQ.0 ) THEN idx = LOCATE(polnam,pm_spec,nalim) c PRINT * , 'LOCATING PM' , polnam(i) , ' AT ' , idx ipoln = idx ELSEIF ( INDEX(polnam,'MASS').NE.0 ) THEN c c print *, 'this step is for PM MASS' c print *, polnam(i) ipoln = -1 ELSE PRINT * , 'Unindentified species.' , polnam STOP ENDIF c loop over all grid cells DO ix = 1 , ngmx DO iy = 1 , ngmy c gas phase species IF(igoa.EQ.1 .AND. ipoln.GT.0) THEN map(ix,iy) = gconc(ix,iy,ipoln) ELSEIF(igoa.EQ.0 .AND. ipoln.NE.0) THEN c particle phase species rsum = 0 DO it = 1 , NTYPE DO itt = 1 , NTIME IF ( ipoln.EQ.-1 ) THEN inalim1 = 1 inalim2 = NALIM - 1 ! dry particle only ELSE inalim1 = ipoln inalim2 = ipoln ENDIF c DO ip = inalim1 , inalim2 c exclude the dummy mass variable for total PM if (index(pm_spec(ip),'TRACER').ne.0 + .and.ipoln.eq.-1) goto 1212 c rsumi = 0. DO ib = 1 , NBLIM conc0(ib)= aconc(it,itt,ib,ip,ix,iy) dp0(ib)= radp(it,itt,ib,ix,iy)*2. rsumi = rsumi + aconc(it,itt,ib,ip,ix,iy) ENDDO IF ( rsumi.gt.0. ) THEN CALL DMDDP(conc0,dp0,NBLIM, & conc1,rsbins,nsbins,ddp,1) DO ib = 1 , nsbins conc1(ib) = conc1(ib)*ddp(ib) ENDDO CALL TOTALMASS(conc1,rsbins, & nsbins,radius1,radius2,1,rsum2) c=================================================== c eliminate ultrafine nitrate, sulfate, and ammonium if(radius2.le.2.5e-7) then ! for particles smaller than 0.25um if(index(pm_spec(ip),'N(V)').ne.0 .or. $ index(pm_spec(ip),'S(VI)').ne.0 .or. $ index(pm_spec(ip),'N(-III)').ne.0 ) then rsum2=0. endif endif c=================================================== rsum = rsum + rsum2 ENDIF 1212 ENDDO ! ip ENDDO ! itt ENDDO ! it map(ix,iy) = rsum ENDIF ! igoa ENDDO ! ix ENDDO ! iy c=================================================== c--WRITE OUTPUT iwrite=12 OPEN (UNIT=iwrite,FILE=outfile,STATUS='unknown',iostat=ierr) if (ierr.ne.0) then print *,'Fail to open file' stop endif DO ix=1,ngmx DO iy=1,ngmy IF(map(ix,iy).gt.0.) then write(iwrite,'(i4,1a,i4,1a,e12.3)') + ix,',',iy,',',map(ix,iy) ENDIF ENDDO ENDDO CLOSE (iwrite) c=================================================== c-- END subroutine pread(filesv,ngmx,ngmy,nact,nalim,nblim,ntype,ntime, + gconc,aconc,radp,rnump,gas_spec,pm_spec) c**************************************************************************** c c Read the output created by UCD air quality model c c Written by: Jianlin Hu (July 2013) c University of Caifornia, Davis c jjlhu@ucdavis.edu c c inputs: c filesv: input file name c outputs: c gconc: gas concentrations c aconc: pm concentrations c gas_spec: gas species list c pm_spec: pm species list c c**************************************************************************** implicit none integer ngmx,ngmy,nact,nalim,nblim,ntype,ntime c-- real gconc(ngmx,ngmy,nact) ! gas conc real aconc(ntype,ntime,nblim,nalim,ngmx,ngmy) ! pm conc real radp(ntype,ntime,nblim,ngmx,ngmy) ! radius real rnump(ntype,ntime,nblim,ngmx,ngmy) ! number conc character filesv*256 ! input file name character*16 gas_spec(nact), pm_spec(nalim) ! gas and pm species lists c--other variables integer ngmz,maxtrj,ntrajsite parameter(ngmz=16,maxtrj=24,ntrajsite=10) real dzrho(ngmz) real timmin,tavg,height real tcut(ntype,ntime) integer istart(maxtrj,ntrajsite) character trajcd(10)*6 c integer i,j,jj,ix,iy,icolm,incm,ic,inz,i1,i2,i3 ! counter variables integer iread,iyr,imn,idy,nvert,nact1,nalim1,ntsite,maxtrj1 real rdum1,rdum2,rdum3,rdum(nalim) ! dummy variables c--initialize the conc array gconc=-999. aconc=-999. c--open input file-- 910 iread = 11 c open(unit=iread,file=filesv, + status='old',form='unformatted') c--read header information of the file-- read(iread) iyr, imn, idy, timmin, tavg read(iread) height, nvert,(dzrho(j),j=1,nvert) read(iread) nact1,(gas_spec(j),j=1,nact1) read(iread) nalim1,(pm_spec(j),j=1,nalim1) read(iread) ntsite, (trajcd(j),j=1,ntsite) read(iread) maxtrj1,((istart(j,jj),j=1,maxtrj1),jj=1,ntsite) c print *,'read header done.' c--read data c gas 980 read(iread,end=981)timmin,ix,iy,icolm, + (gconc(ix,iy,ic),ic=1,nact) if (ix.eq.0.and.iy.eq.0) then close(iread) return endif c particles 982 read(iread,end=981)inz,i1,i2,i3 c check if it is the end of the cell record if (inz.eq.0.and.i1.eq.0.and.i2.eq.0.and.i3.eq.0) then backspace(iread) read(iread,end=981)inz,i1,i2,i3, + rdum1,rdum2,rdum3, + (rdum(j),j=1,nalim) goto 980 endif backspace(iread) read(iread,end=981)inz,i1,i2,i3,tcut(i1,i2), + radp(i1,i2,i3,ix,iy), + rnump(i1,i2,i3,ix,iy), + (aconc(i1,i2,i3,j,ix,iy),j=1,nalim) goto 982 981 continue close(iread) return end FUNCTION LOCATE(Cdum,Slistp,Nalim) c**************************************************************************** c c written by: Mike Kleeman (Sept 1995) c c The purpose of this function is to find the variable cdum in the list c slistp of size nalim. c c*************************************************************************** c c inputs: c c cdum chracter variable being search for c slistp character list being search in c nalim size of the character list being searched c IMPLICIT REAL(A-H,O-Z) C*** Start of declarations inserted by SPAG INTEGER i , LOCATE , Nalim C*** End of declarations inserted by SPAG CHARACTER*16 Cdum , Slistp(Nalim) c--loop through the list and try to find cdum-- LOCATE = 0 DO i = 1 , Nalim IF ( Cdum.EQ.Slistp(i) ) THEN LOCATE = i RETURN ENDIF ENDDO c--return to the calling program-- END SUBROUTINE DMDDP(Conc,Dp0,Nb0,Concs,Dps,Nbs,Ddp,Ilog) c**************************************************************************** c c The purpose of this subroutine is to convert binned mass distribution c to a dm/ddp (or dm/dlogdp)distriubtion and regrid it into standard size bins c c Written by: Qi Ying (Feburary 2003) c University of California, Davis c Davis, CA 95616 c c Input: c conc concentration in the size bin (dm) c dp0 mean diameter of the bin c nb0 number of total bins c dps standard mean diameter of the regridded bins c nbs number of regridded bins c ilog using logdp (ilog=1) c c Output: c conc normalized concentration in the orignal bins (dm/ddp) c concs normalized concentration in the regridded bins (dm/dlogdp) c ddp bin width (ddp or dlogdp) c c***************************************************************************** IMPLICIT NONE INTEGER MAXB PARAMETER (MAXB=100) c input parameters INTEGER Nb0 , Nbs , Ilog REAL Conc(Nb0) , Dp0(Nb0) , Concs(MAXB) , Dps(MAXB) , Ddp(MAXB) c local variables INTEGER i , j REAL dp(MAXB) REAL dp1(3,MAXB) , dp2(3,MAXB) , conc1(MAXB) INTEGER ibin REAL br , bl , br1 , bl1 INTEGER iset , ip1 , ip2 , nb REAL rtemp c clear concentrations DO i = 1 , Nbs Concs(i) = 0. ENDDO c make a local copy of dp0 DO i = 1 , Nb0 dp(i) = Dp0(i) ENDDO nb = Nb0 c find fist and last non zero position iset = 0 DO i = 1 , nb IF ( dp(i).NE.0 .AND. iset.EQ.0 ) THEN ip1 = i iset = 1 ENDIF ENDDO IF ( iset.EQ.0 ) STOP 'dmddp: all zero diameter in size bins' iset = 0 DO i = nb , 1 , -1 IF ( dp(i).NE.0 .AND. iset.EQ.0 ) THEN ip2 = i iset = 1 ENDIF ENDDO c testing code c print *,'original nb=',nb c do i=1,nb c print *,i,dp(i),log10(dp(i)) c enddo c move contrations upward DO i = ip1 , ip2 dp(i-ip1+1) = dp(i) Conc(i-ip1+1) = Conc(i) ENDDO c correct nb nb = ip2 - ip1 + 1 c sorting DO i = 1 , nb DO j = i + 1 , nb IF ( dp(i).GT.dp(j) ) THEN rtemp = dp(i) dp(i) = dp(j) dp(j) = rtemp rtemp = Conc(i) Conc(i) = Conc(j) Conc(j) = rtemp ENDIF ENDDO ENDDO c checking diameter sequence DO i = 1 , nb - 1 IF ( dp(i).GT.dp(i+1) ) THEN PRINT * , 'bin sequence error.' STOP ENDIF ENDDO c convert input diameter to log scale IF ( Ilog.EQ.1 ) THEN DO i = 1 , nb dp1(2,i) = LOG10(dp(i)) ENDDO DO i = 1 , Nbs dp2(2,i) = LOG10(Dps(i)) ENDDO ELSE DO i = 1 , nb dp1(2,i) = dp(i) ENDDO DO i = 1 , Nbs dp2(2,i) = Dps(i) ENDDO ENDIF c debug c print *,'final, nb=',nb c do i=1,nb c print *,i,dp1(2,i) c enddo c compute upper and lower limit of each bin c -- mid points DO i = 2 , nb - 1 dp1(1,i) = (dp1(2,i)+dp1(2,i-1))/2. dp1(3,i) = (dp1(2,i)+dp1(2,i+1))/2. ENDDO DO i = 2 , Nbs - 1 dp2(1,i) = (dp2(2,i)+dp2(2,i-1))/2. dp2(3,i) = (dp2(2,i)+dp2(2,i+1))/2. ENDDO c -- first point dp1(3,1) = dp1(1,2) dp1(1,1) = dp1(2,1) - (dp1(3,1)-dp1(2,1)) dp2(3,1) = dp2(1,2) dp2(1,1) = dp2(2,1) - (dp2(3,1)-dp2(2,1)) c -- last point dp1(1,nb) = dp1(3,nb-1) dp1(3,nb) = dp1(2,nb) + (dp1(2,nb)-dp1(1,nb)) dp2(1,Nbs) = dp2(3,Nbs-1) dp2(3,Nbs) = dp2(2,Nbs) + (dp2(2,Nbs)-dp2(1,Nbs)) c save standard bin width DO i = 1 , Nbs Ddp(i) = dp2(3,i) - dp2(1,i) ENDDO c check bin c do i=1,nb c print *,i,dp1(1,i),dp1(2,i),dp1(3,i) c enddo c do j=1,nbs c print *,j,dp2(1,j),dp2(2,j),dp2(3,j) c enddo c normalize concentration in the orginal size bin DO i = 1 , nb IF ( dp1(3,i).EQ.dp1(1,i) ) dp1(3,i) = dp1(3,i)*1.1 Conc(i) = Conc(i)/(dp1(3,i)-dp1(1,i)) ENDDO c regrid into normalized bins DO i = 1 , Nbs c print *,'testing standard bin ',i bl = dp2(1,i) br = dp2(3,i) c -- find out the original bins that is in the regrided bin ibin = 0 Concs(i) = 0. DO j = 1 , nb c print *,' checking raw bins:',j bl1 = dp1(1,j) br1 = dp1(3,j) IF ( bl.GE.bl1 .AND. br.LE.br1 ) THEN ! easy case Concs(i) = Conc(j) ibin = ibin + 1 c print *,' complete conver the standard bin',concs(i) GOTO 100 ELSEIF ( bl1.GE.bl .AND. br1.LE.br ) THEN ! all in the bin Concs(i) = Concs(i) + Conc(j)*(br1-bl1)/(br-bl) ibin = ibin + 1 c print *,' all in the standard bin' ELSEIF ( bl1.GE.bl .AND. bl1.LT.br ) THEN ! left half in the bin Concs(i) = Concs(i) + Conc(j)*(br-bl1)/(br-bl) c print *,' left half in the standard bin',concs(i) ibin = ibin + 1 ELSEIF ( br1.GT.bl .AND. br1.LE.br ) THEN ! right half in the bin ibin = ibin + 1 c print *,' right half in the standard bin',concs(i) Concs(i) = Concs(i) + Conc(j)*(br1-bl)/(br-bl) ENDIF ENDDO c print *,ibin,' bins in standard size bin ',i,'concs=',concs(i) 100 ENDDO END subroutine totalmass(conc,dp0,nbs,rd10,rd20,ilog,rsum) c************************************************************************* c c The purpose of this subroutine is to calculate the total mass c within diameter range rd1 and rd2. c c Written by: Qi Ying (Feburary 2003) c University of California, Davis c Davis CA 95616 c c Input: conc concentration c dp0 mid-point diameter of size bins c nbs number of standard sizebins c rd10 lower range of diameter c rd20 upper range of diameter c iog logarithmic size distribution (=1) c c Output: c rsum total mass in the size range. c c************************************************************************* implicit real(a-h,o-z) integer maxb parameter (maxb=100) integer nbs real conc(nbs),dp0(nbs) real rd1,rd2 integer ilog real rsum real dp1(maxb),dp(3,maxb) integer nb c make a local copy of dp0 do i=1,nbs dp1(i)=dp0(i) enddo nb=nbs c find fist and last non zero position iset=0 do i=1,nb if (dp1(i).ne.0.and.iset.eq.0) then ip1=i iset=1 endif enddo if (iset.eq.0) stop 'tmass: all zero diameter in size bins' iset=0 do i=nb,1,-1 if (dp1(i).ne.0.and.iset.eq.0) then ip2=i iset=1 endif enddo c move contrations upward do i=ip1,ip2 dp1(i-ip1+1)=dp1(i) conc(i-ip1+1)=conc(i) enddo c correct nb nb=ip2-ip1+1 c sorting do i=1,nb do j=i+1,nb if (dp1(i).gt.dp1(j)) then rtemp=dp1(i) dp1(i)=dp1(j) dp1(j)=rtemp rtemp=conc(i) conc(i)=conc(j) conc(j)=rtemp endif enddo enddo c checking diameter sequence do i=1,nb-1 if (dp1(i).gt.dp1(i+1)) then print *,'bin sequence error.' stop endif enddo c convert input diameter to log scale if (ilog.eq.1) then do i=1,nb dp(2,i)=log10(dp1(i)) enddo rd1=log10(rd10) rd2=log10(rd20) else do i=1,nb dp(2,i)=dp1(i) enddo rd1=rd10 rd2=rd20 endif c compute upper and lower limit of each bin c -- mid points do i=2,nb-1 dp(1,i)=(dp(2,i)+dp(2,i-1))/2. dp(3,i)=(dp(2,i)+dp(2,i+1))/2. enddo c -- first point dp(3,1)=dp(1,2) dp(1,1)=dp(2,1)-(dp(3,1)-dp(2,1)) c -- last point dp(1,nb)=dp(3,nb-1) dp(3,nb)=dp(2,nb)+(dp(2,nb)-dp(1,nb)) c check nb c print *,'size bins:',nb c do i=1,nb c print *,i,dp(1,i),dp(2,i),dp(3,i),conc(i) c enddo c stop c summing up rsum=0. do i=1,nb d1=dp(1,i) d2=dp(3,i) if (d1.ge.rd1.and.d1.lt.rd2 .and. + d2.le.rd2.and.d2.gt.rd1) then rsum=rsum+conc(i) elseif(d1.lt.rd1.and.d2.gt.d1) then rsum=rsum+conc(i)*(d2-rd1)/(d2-d1) elseif(d1.lt.rd2.and.d2.gt.rd2) then rsum=rsum+conc(i)*(rd2-d1)/(d2-d1) elseif(d1.lt.rd1.and.d2.gt.rd2) then rsum=rsum+conc(i)*(rd2-rd1)/(d2-d1) else c print *,'not in the bin',d1,d2,rd1,rd2 c print *,'mass=',conc(i),rsum c print *,'error',d2,d1,rd1,rd2,rsum c istop=1 endif c print *,i,dp(2,i),conc(i),rsum enddo c if (istop.eq.1) then c print *,'rsum=',rsum c endif return end FUNCTION UPCASE(RUNNAM) CHARACTER*10 UPCASE, RUNNAM, NEWNAM INTEGER IDEL IDEL = ICHAR('a') - ICHAR('A') DO 100 I=1, 10 IF (RUNNAM(I:I) .LT. 'a' .or. RUNNAM(I:I) .GT. 'z') THEN NEWNAM(I:I) = RUNNAM(I:I) ELSE NEWNAM(I:I) = CHAR(ICHAR(RUNNAM(I:I)) - IDEL) ENDIF 100 CONTINUE UPCASE = NEWNAM RETURN END