program opac ccccc ------------------------------------------------------------------c c extracts and calculates optical properties of aerosols and clouds c c from the data stored in directory ../optdat/. c c c c Version 3.0 c c ------------ c c 01.10.97 new version based on Version 2.0 of 17.04.97 c c c c Version 3.1 c c ------------ c c 21.10.97 new data format in ../optdat/ c c 30.10.97 corrected: units [1/m] -> [1/km] c c 06.11.97 opt.depth for clouds corrected. c c 14.11.97 tested with Linux, SunOS, Dec Ultrix, HP Unix, DOS c c c c 14.11.97 M. Hess c ccccc ------------------------------------------------------------------c character*2 chum character*4 comnam character*6 inp character*7 res character*8 file,optnam,optun,opanam character*10 opt character*18 infile character*20 rname character*30 typnam,tname,tnew logical wa,wb real mixrat,numden integer nlt(10) common /prog/ nprog common /file/ lfile,file,rname,infile,inp,res,opt common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8), * mixrat(10),dens(10,8) common /mipaco/ sigmac(20),rminc(20,8),rmaxc(20,8),rmodc(20,8), * densc(20,8) common /mipaty/ nco(20),mco(20,10),dnumb(20,10),nlay(20), * hmint(20,4),hmaxt(20,4),part(20,4),tname(20) common /input/ mtyp,nuco(10),dnum(10),mcom, * hmin(5),hmax(5),hpar(5), * nwavea,indwa(100),wavea(100), * nwaveb,indwb(100),waveb(100), * nhc,indhum(8),nopt,indopt(20),tnew common /atyp/ natyp,mcomp,ncomp(10),numden, * typnam,comnam(20) common /layer/ mlay,nltyp(10),parlay(10,5),boundl(10), * boundu(10) common /opar/ mopar,jnopar(20),nop,opanam(20),optnam(20), * optun(20) common /wavel/ mlamb,alamb(61),niw,wlamb(61), * il035,il05,il055,il08,alambp(61),niwp common /hum/ khum(8),ahum(8),nih,nhum(8),mhum,chum(8) common /norm/ norm,mixnor data norm /1/ data nprog /0/ data nlt /1,2,3,4,5,0,0,0,0,0/ data nhum /0,50,70,80,90,95,98,99/,mhum/8/ data chum /'00','50','70','80','90','95','98','99'/ data optnam /'ext.coef','sca.coef','abs.coef','sisc.alb', * 'asym.par','op.depth', * ' ','turb.fac','li.ratio','pha.func', * 'ext.rat ','abs.rat ', * ' ','ext.norm',' ',' ', * 'visibil ','ref.real', * 'ref.imag',' '/ data optun /' [1/km] ',' [1/km] ',' [1/km] ',' ', * ' ',' ', * ' ',' ',' [sr] ',' [1/km] ', * '[m^2/g] ','[m^2/g] ', * ' ',' ',' ',' ', * ' [km] ',' ', * ' ',' '/ data comnam /'inso','waso','soot','ssam','sscm','minm','miam', * 'micm','mitr','suso','stco','stma','cucc','cucp', * 'cuma','fogr','cir1','cir2','cir3',' '/ print*,'********************************************************' print*,'* Optical Properties of Aerosols and Clouds OPAC 3.1 *' print*,'* -------------------------------------------------- *' print*,'* *' print*,'* The contents of OPAC has been described in: *' print*,'* M. Hess, P. Koepke and I. Schult (1997): *' print*,'* "Optical Properties of Aerosols and Clouds: *' print*,'* The software package OPAC" *' print*,'* submitted to Bull. Am. Meteor. Soc. *' print*,'* *' print*,'* last update: 14.11.97 M. Hess *' print*,'********************************************************' print*,' ' print*,' ' ccccc ------------------------------------------------------------------c c read data from configuration file opac.cfg c ccccc ------------------------------------------------------------------c call readcfg ccccc -----------------------------------------------------------------c c read name of input file c ccccc -----------------------------------------------------------------c 4321 print*,' name of input file (without extension .inp),', * ' (max. 8 characters)' print*,' default: opac.inp' read (*,'(a)') file call lang2 (file,8,lfile) if (lfile.gt.8) then print*,' lfile= ',lfile print*,' file name too long! Try again!' goto 4321 end if if (file.ne.' ') then infile=inp//file(1:lfile)//'.inp' else infile=inp//'opac.inp' file='opac' lfile=4 end if write (*,*) ' input file ',infile,' used' ccccc ------------------------------------------------------------------c c read data from input file c ccccc ------------------------------------------------------------------c call readinp ccccc ------------------------------------------------------------------c c prepare data selection + check input data c c c c 1. mixture of type c ccccc ------------------------------------------------------------------c c print*,mtyp if(mtyp.eq.0) then typnam=tnew numden=0. mcomp=0 do ic=1,5 if (dnum(ic).ne.0.) then mcomp=mcomp+1 numden=numden+dnum(ic) end if end do do ic=1,mcomp ncomp(ic)=nuco(ic) mixrat(ic)=dnum(ic)/numden sigma(ic)=sigmac(nuco(ic)) end do else typnam=tname(mtyp) mcomp=nco(mtyp) numden=0. do ic=1,mcomp numden=numden+dnumb(mtyp,ic) end do do ic=1,mcomp ncomp(ic)=mco(mtyp,ic) mixrat(ic)=dnumb(mtyp,ic)/numden sigma(ic)=sigmac(mco(mtyp,ic)) end do end if c print*,typnam,numden if (mcomp.gt.1) then do ic=1,mcomp if (ncomp(ic).ge.11) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' in the present version, clouds may not', * ' be mixed with other particles' print*,' ' print*, ' Please change input file!' stop end if end do end if ccccc ------------------------------------------------------------------c c 2. relative humidity c ccccc ------------------------------------------------------------------c nih=0 do ih=1,nhc if(indhum(ih).ne.0) then nih=nih+1 khum(nih)=ih ahum(nih)=nhum(ih) end if end do if (nih.ne.0) then do ih=1,nih do ic=1,mcomp dens(ic,ih)=densc(ncomp(ic),khum(ih)) rmin(ic,ih)=rminc(ncomp(ic),khum(ih)) rmax(ic,ih)=rmaxc(ncomp(ic),khum(ih)) rmod(ic,ih)=rmodc(ncomp(ic),khum(ih)) c print*,khum(ih),rmin(ic,ih),rmax(ic,ih),rmod(ic,ih) end do end do else print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' no relative humidity selected!' print*,' ' print*, ' Please change input file!' stop end if ccccc ------------------------------------------------------------------c c 3. height profile c ccccc ------------------------------------------------------------------c ht=0. do il=1,5 ht=ht+(hmax(il)-hmin(il)) end do if (mtyp.eq.0.or.ht.ne.0.) then mlay=0 do il=1,5 if (hmax(il).ne.hmin(il)) then mlay=mlay+1 nltyp(mlay)=nlt(il) boundl(mlay)=hmin(il) boundu(mlay)=hmax(il) parlay(mlay,1)=hpar(il) end if end do else mlay=nlay(mtyp) do il=1,mlay if (mlay.eq.1) then ! for clouds nltyp(1)=5 else nltyp(il)=nlt(il) end if boundl(il)=hmint(mtyp,il) boundu(il)=hmaxt(mtyp,il) parlay(il,1)=part(mtyp,il) end do end if ccccc ------------------------------------------------------------------c c 4. optical parameters c ccccc ------------------------------------------------------------------c mopar=nopt nop=0 do iop=1,mopar jnopar(iop)=indopt(iop) if (jnopar(iop).eq.1) then nop=nop+1 opanam(nop)=optnam(iop) end if end do if (jnopar(15).eq.1) nop=nop-1 ! not wavelength dependent if (jnopar(16).eq.1) nop=nop-1 ! " if (jnopar(17).eq.1) nop=nop-1 ! " if (jnopar(18).eq.1) nop=nop+1 ! refr.index has 2 values do ic=1,mcomp if (jnopar(11).eq.1.and.ncomp(ic).ge.11.or. * jnopar(12).eq.1.and.ncomp(ic).ge.11) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' in the present version, no mass related', * ' parameters may be calculated' print*,' for clouds' print*,' ' print*, ' Please change input file!' stop end if end do if (jnopar(15).eq.1.and.jnopar(1).eq.0.or. * jnopar(15).eq.1.and.jnopar(4).eq.0.or. * jnopar(15).eq.1.and.jnopar(5).eq.0) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' extinction coefficient, single scattering albedo and' print*,' asymmetry parameter must be selected too' print*,' for calculation of solar and terrestrial integrated' print*,' parameters' print*,' ' print*, ' Please change input file!' stop end if if (jnopar(16).eq.1.and.jnopar(1).eq.0) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' extinction coefficient must be selected, too' print*,' for calculation of Angstrom coefficients' print*,' ' print*, ' Please change input file!' stop end if if (jnopar(18).eq.1.and.mcomp.gt.1) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' refractive indices can only be printed for ' print*,' single components' print*,' ' print*, ' Please change input file!' stop end if if (jnopar(7).eq.1.or.jnopar(13).eq.1) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' you selected an invalid optical parameter' print*,' ' print*, ' Please change input file!' stop end if ccccc ------------------------------------------------------------------c c 5. wavelengths c ccccc ------------------------------------------------------------------c wa=.false. wb=.false. do iw=1,nwavea if (indwa(iw).ne.0) wa=.true. end do do iw=1,nwaveb if (indwb(iw).ne.0) wb=.true. end do if (wa.and.wb) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' you selected wavelengths in both wavelength tables!' print*,' ' print*, ' Please change input file!' stop else if (wa) then niw=0 mlamb=nwavea do il=1,mlamb wlamb(il)=wavea(il) if (indwa(il).ne.0) then niw=niw+1 alamb(niw)=wlamb(il) end if end do else if (wb) then niw=0 mlamb=nwaveb do il=1,mlamb wlamb(il)=waveb(il) if (indwb(il).ne.0) then niw=niw+1 alamb(niw)=wlamb(il) end if end do end if il05=0 il055=0 il035=0 il08=0 do iw=1,niw if (alamb(iw).eq.0.35) il035=iw if (alamb(iw).eq.0.5) il05=iw if (alamb(iw).eq.0.55) il055=iw if (alamb(iw).eq.0.8) il08=iw end do niwp=niw if (jnopar(15).eq.1) then ! bei solar/terrestr. integrierten Groessen niwp=niw ! werden alle Wellenl. eingelesen, aber nicht do il=1,niwp ! gedruckt alambp(il)=alamb(il) end do niw=mlamb do il=1,niw alamb(il)=wlamb(il) if (alamb(il).eq.0.35) il035=il if (alamb(il).eq.0.5) il05=il if (alamb(il).eq.0.55) il055=il if (alamb(il).eq.0.8) il08=il end do end if if (jnopar(14).eq.1.and.il055.eq.0) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' wavelength 0.55 micron must be selected' print*,' for calculation of normalized extinction coefficients' print*,' ' print*, ' Please change input file!' stop end if if (jnopar(16).eq.1.and.il035.eq.0.or. * jnopar(16).eq.1.and.il05.eq.0.or. * jnopar(16).eq.1.and.il08.eq.0) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' wavelengths 0.35, 0.5 and 0.8 micron must be selected' print*,' for calculation of Angstrom coefficients' print*,' ' print*, ' Please change input file!' stop end if ccccc ------------------------------------------------------------------c c Verzweigung ins Unterprogramm RAWOPT zur Berechnung der c c optischen Eigenschaften c ccccc ------------------------------------------------------------------c call rawopt print*,' READY.' ccccc ------------------------------------------------------------------c c Formate c ccccc ------------------------------------------------------------------c 100 format(i1) 110 format(e17.10) 120 format(a1) 130 format(e17.10) 150 format(i1,3x,f6.3) 200 format(i2) 300 format(a50) stop end subroutine comment (com,ifile) character*1 com,dum 1011 read (ifile,'(a1)') dum if (dum.eq.com) then goto 1011 else backspace(ifile) end if return end subroutine readcfg ccccc ------------------------------------------------------------------c c read microphysical parameters of components and types from c c file opac.cfg. c c c c 30.09.97 M. Hess c ccccc ------------------------------------------------------------------c character*1 dum character*4 os character*6 inp character*7 res character*8 file character*10 opt character*18 infile character*20 rname character*30 tname common /file/ lfile,file,rname,infile,inp,res,opt common /mipaco/ sigmac(20),rminc(20,8),rmaxc(20,8),rmodc(20,8), * densc(20,8) common /mipaty/ nco(20),mco(20,10),dnumb(20,10),nlay(20), * hmint(20,4),hmaxt(20,4),part(20,4),tname(20) open (8,file='opac.cfg') call comment('#',8) read (8,'(a4)') os if (os.eq.'dos ') then c inp='input\' c res='result\' c opt='..\optdat\' else if (os.eq.'unix') then inp='input/' res='result/' opt='../optdat/' else stop ' wrong operating system!' end if call comment('#',8) read (8,*) ncom read (8,*) ntyp call comment('#',8) do icom=1,ncom read(8,'(a1)') dum read(8,'(a1)') dum read(8,'(a1)') dum do ih=1,8 if (icom.le.10) then read (8,*) nh,rminc(icom,ih),rmaxc(icom,ih),rmodc(icom,ih), * densc(icom,ih),sigmac(icom) else read (8,*) nh,rminc(icom,ih),rmaxc(icom,ih),rmodc(icom,ih), * densc(icom,ih) end if read (8,'(a1)') dum if (ih.eq.1.and.dum.eq.'#') then do ihd=2,8 rminc(icom,ihd)=rminc(icom,1) rmaxc(icom,ihd)=rmaxc(icom,1) rmodc(icom,ihd)=rmodc(icom,1) densc(icom,ihd)=densc(icom,1) end do goto 1111 else backspace(8) end if end do call comment('#',8) 1111 continue end do do icom=1,ncom do ih=1,8 c print*,rminc(icom,ih),rmaxc(icom,ih),rmodc(icom,ih), c * densc(icom,ih) end do end do call comment('#',8) do ityp=1,ntyp call comment('#',8) read (8,200) nt,tname(ityp) 200 format(i2,2x,a30) read (8,*) nco(ityp) do ico=1,nco(ityp) read (8,100) mco(ityp,ico),dnumb(ityp,ico) 100 format(i2,8x,e10.3) end do read (8,*) nlay(ityp) do ilay=1,nlay(ityp) read (8,*) hmint(ityp,ilay),hmaxt(ityp,ilay),part(ityp,ilay) end do end do close(8) return end subroutine readinp ccccc ------------------------------------------------------------------c c read from input file *.inp which data to extract and calculate c c from the database. c c c c 21.09.97 M. Hess c ccccc ------------------------------------------------------------------c character*6 inp character*7 res character*8 file character*10 opt character*18 infile character*20 rname character*30 tnew common /file/ lfile,file,rname,infile,inp,res,opt common /input/ mtyp,nuco(10),dnum(10),mcom, * hmin(5),hmax(5),hpar(5), * nwavea,indwa(100),wavea(100), * nwaveb,indwb(100),waveb(100), * nhc,indhum(8),nopt,indopt(20),tnew open (7,file=infile) rname(1:7)=res rname(8:(8+(lfile-1)))=file(1:lfile) rname((9+(lfile-1)):(13+(lfile-2)))='.out' call comment('#',7) ccccc ------------------------------------------------------------------c c read no. of type c ccccc ------------------------------------------------------------------c read (7,*) mtyp call comment('#',7) read (7,'(a30)') tnew call comment('#',7) ccccc ------------------------------------------------------------------c c read mixture of type c ccccc ------------------------------------------------------------------c mcom=0 do ic=1,5 read (7,*) nuco(ic),dnum(ic) if (dnum(ic).gt.0.) mcom=mcom+1 end do call comment('#',7) ccccc ------------------------------------------------------------------c c read height profile c ccccc ------------------------------------------------------------------c do il=1,5 read (7,*) hmin(il),hmax(il),hpar(il) end do call comment('#',7) ccccc ------------------------------------------------------------------c c read wavelengths c ccccc ------------------------------------------------------------------c read (7,*) nwavea do iw=1,nwavea read (7,100) indwa(iw),wavea(iw) 100 format(i1,32x,f10.3) end do call comment('#',7) read (7,*) nwaveb do iw=1,nwaveb read (7,100) indwb(iw),waveb(iw) end do call comment('#',7) ccccc ------------------------------------------------------------------c c read relative humidity classes c ccccc ------------------------------------------------------------------c read (7,*) nhc do ih=1,nhc read (7,*) indhum(ih) end do call comment('#',7) ccccc ------------------------------------------------------------------c c read selected optical parameters c ccccc ------------------------------------------------------------------c read (7,*) nopt do iopt=1,nopt read(7,*) indopt(iopt) end do close(7) return end subroutine rawopt ccccc -----------------------------------------------------------------c c Berechnung der optischen Parameter eines Aerosoltyps f�r c c mehrere Wellenl�ngen und Feuchten c c c c Folgende Unterprogramme werden verwendet: c c c c - RAWMIC c c - HEAD2 c c - OPTRAW c c - OPTPAR c c - OUT2 c c c c Diese Version gehoert zur eingeschraenkten Datenbank OPAC und c c beruht auf RAWOPT vom 04.11.93. c c c c 13.05.94 Einlesen der Background-Extinktion. c c 12.01.96 Massen werden immer berechnet. c c 01.10.97 changed reading 'extback.dat' c c c c Stand: 01.10.97 M. Hess c ccccc -----------------------------------------------------------------c real mixrat,n,numden integer prnr,acnr,rht character*1 dum character*2 chum character*3 atn,pat,typnam*30 character*4 comnam common /prog/ nprog common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8), * mixrat(10),dens(10,8) COMMON /FTASTR/ EXTFTA(61),EXTSTR(61),extmit(61) common /atyp/ natyp,mcomp,ncomp(10),numden, * typnam,comnam(20) common /wavel/ mlamb,alamb(61),niw,wlamb(61), * il035,il05,il055,il08,alambp(61),niwp common /hum/ khum(8),ahum(8),nih,nhum(8),mhum,chum(8) common /mipoi/ latx,lonx,nl,prnr,rht(2),n(2), * njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2) common /out/ oparam(16,2),phaf(167,2),indop(16) ccccc -----------------------------------------------------------------c c Aufruf von RAWMIC f�r kombinierte opt./mikrophys. Groessen c ccccc -----------------------------------------------------------------c call rawmic ccccc -----------------------------------------------------------------c c Definition der Gr��en, die in den Unterprogrammen gebraucht c c werden c ccccc -----------------------------------------------------------------c n(1)=numden nl=1 njc(1)=mcomp do ic=1,mcomp acnr(ic,1)=ncomp(ic) acmr(ic,1)=mixrat(ic) end do ccccc -----------------------------------------------------------------c c Schleife �ber alle verlangten Feuchteklassen c c ccccc -----------------------------------------------------------------c if (nih.eq.0) then njh=1 else njh=nih end if do ih=1,njh ccccc -----------------------------------------------------------------c c Einlesen der Extinktionskoeffizienten des tropos. backgr. c c und des strat. backgr. und des Min. transported. c ccccc -----------------------------------------------------------------c c print*,'read extback.dat' open (9,file='extback.dat') IL=1 READ(9,'(a1)') dum READ(9,'(a1)') dum DO IWL=1,mlamb READ(9,*) WAVE,EXTFT,EXTST,extmi do ila=1,niw IF (WAVE.EQ.alamb(ila)) THEN EXTFTA(IL)=EXTFT EXTSTR(IL)=EXTST extmit(il)=extmi IL=IL+1 END IF end do end do close (9) ccccc -----------------------------------------------------------------c c Einlesen der optischen Rohdaten im Verzeichnis OPTDAT c ccccc -----------------------------------------------------------------c c print*,'start optraw' call optraw (ih) ccccc -----------------------------------------------------------------c c Beschriftung des Output-Files c ccccc -----------------------------------------------------------------c c print*,'start head2' call head2 (ih) ccccc -----------------------------------------------------------------c c Schleife �ber alle verlangten Wellenl�ngen c ccccc -----------------------------------------------------------------c do il=1,niw ccccc -----------------------------------------------------------------c c Auswahl der benoetigten Daten aus den eingelesenen c ccccc -----------------------------------------------------------------c c print*,'start optdat' call optdat(il) ccccc -----------------------------------------------------------------c c Berechnung der optischen Parameter c ccccc -----------------------------------------------------------------c c print*,'start optpar' call optpar(il,ih) end do ! Wellenlaengen-Schleife ccccc -----------------------------------------------------------------c c Berechnung der solar und terrestr. integrierten Groessen und c c der Angstrom-Koeffizienten und der Sichtweite c ccccc -----------------------------------------------------------------c c print*,'start intco' call intco end do ! Feuchte-Schleife return end ccccc *****************************************************************c subroutine head2 (ih) c *****************************************************************c c c c -----------------------------------------------------------------c c version for OPAC 3.0 c c c c 29.09.97 new header c c c c 29.09.97 M. Hess c ccccc -----------------------------------------------------------------c real mixrat,numden character*2 chum character*4 comnam character*6 inp character*7 res character*8 file,optun,optnam,opanam character*10 opt character*18 infile character*20 rname character*30 typnam common /file/ lfile,file,rname,infile,inp,res,opt common /hum/ khum(8),ahum(8),nih,nhum(8),mhum,chum(8) common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8), * mixrat(10),dens(10,8) common /masse/ smas(10,8),smag(8),vsum(10,8),vges(8),xabmax common /atyp/ natyp,mcomp,ncomp(10),numden, * typnam,comnam(20) common /opar/ mopar,jnopar(20),nop,opanam(20),optnam(20), * optun(20) common /out/ oparam(16,2),phaf(167,2),indop(16) common /angle/ jnangle(167),angle(167),nia,ntheta common /wavel/ mlamb,alamb(61),niw,wlamb(61), * il035,il05,il055,il08,alambp(61),niwp open (10,file=rname(1:(13+(lfile-2)))) CCCCC -----------------------------------------------------------------C C Kopf des output-files C CCCCC -----------------------------------------------------------------C if (ih.eq.1) then write(10,100) typnam,mcomp,nih,niwp write (10,108) write(10,112) write (10,108) 100 format('#',12x, * 'Optical Properties of Aerosols and Clouds (OPAC) 3.1'/ * '#',11x, * '------------------------------------------------------'/ * '#',1x,a30,/ * '#',i3,' components'/ * '#',i3,' relative humidities'/ * '#',i3,' wavelengths') 112 format('# RH COMP NUMBER VOLUME MASS DENSITY ', * 'NUM.MIX VOL.MIX MAS.MIX'/ * '# [%] [1/cm^3] [um^3/m^3] [ug/m^3] [g/cm^3] ', * 'RATIO RATIO RATIO') do ihu=1,nih do ic=1,mcomp if (ic.eq.1) then write(10,124) ahum(ihu),comnam(ncomp(ic)), * numden*mixrat(ic), * vges(ihu)*vsum(ic,ihu), * smag(ihu)*smas(ic,ihu), * dens(ic,ihu),mixrat(ic),vsum(ic,ihu), * smas(ic,ihu) else write(10,121) comnam(ncomp(ic)), * numden*mixrat(ic), * vges(ihu)*vsum(ic,ihu), * smag(ihu)*smas(ic,ihu), * dens(ic,ihu),mixrat(ic),vsum(ic,ihu), * smas(ic,ihu) end if end do end do write (10,108) 121 format ('#',5x,a4,1p3e10.3,1x,0pf4.2,3x,1p3e10.3) 124 format ('#',f4.0,1x,a4,1p3e10.3,1x,0pf4.2,3x,1p3e10.3) write (10,110) xabmax 110 format ('#',1x,'only particles up to ',f5.1,' micrometer are ', * 'considered for mass') if (jnopar(10).eq.1) then write (10,102) ntheta write (10,103) (angle(ia),ia=1,ntheta) 102 format('#',i4,' scattering angles for phase functions below:') 103 format('#',8f10.3) end if write (10,107) 107 format('#===================================================', * '=========================') 108 format('#---------------------------------------------------', * '-------------------------') end if CCCCC -----------------------------------------------------------------C C Beschriftung f�r verschiedene relative Feuchten c CCCCC -----------------------------------------------------------------C if (ih.ne.1) then write (10,108) end if WRITE(10,4000) ahum(ih) 4000 FORMAT('#',F4.0,'% relative humidity ') return entry names if (jnopar(10).eq.1) then kop=nop-1 else kop=nop end if if (kop.le.7) then write(10,4001) (optnam(indop(in)),in=1,kop) write(10,4003) (optun(indop(in)),in=1,kop) 4001 format('# wavel. ',7(1x,a8,1x)) 4003 format('# [um] ',7(1x,a8,1x)) else write(10,4001) (optnam(indop(in)),in=1,7) write(10,4002) (optun(indop(in)),in=1,7) write(10,4002) (optnam(indop(in)),in=8,kop) write(10,4002) (optun(indop(in)),in=8,kop) 4002 format('# ',7(1x,a8,1x)) end if return end ccccc *****************************************************************c subroutine out2(ilp,il,iop) c *****************************************************************c c c c output c c c c 29.09.97 M. Hess c ccccc -----------------------------------------------------------------c character*8 optun,optnam,opanam common /wavel/ mlamb,alamb(61),niw,wlamb(61), * il035,il05,il055,il08,alambp(61),niwp common /out/ oparam(16,2),phaf(167,2),indop(16) common /opar/ mopar,jnopar(20),nop,opanam(20),optnam(20), * optun(20) common /angle/ jnangle(167),angle(167),nia,ntheta if ((jnopar(10).eq.1.and.ilp.gt.1).or.ilp.eq.1) then call names end if if (jnopar(10).eq.1) then iop=nop-1 else iop=nop end if if (iop.le.7) then WRITE(10,4010) alamb(il),(oparam(ip,1),ip=1,iop) 4010 FORMAT(1p8e10.3) else WRITE(10,4010) alamb(il),(oparam(ip,1),ip=1,7) ! 1. Zeile WRITE(10,4020) (oparam(ip,1),ip=8,iop) ! 2. Zeile 4020 FORMAT(10x,1p7e10.3) end if if (jnopar(10).eq.1) then write(10,4002) 4002 format(' phase function [1/km]') write(10,4010) (phaf(it,1),it=1,ntheta) end if return end subroutine rawmic ccccc ------------------------------------------------------------------c c Berechnung der Gr��enverteilung an Stuetzpunkten innerhalb der c c Intervallgrenzen. c c Ausgabe dieser Wertepaare und der mikrophysikalischen Parameter c c der Verteilungen aller Komponenten des Aerosoltyps. c c Die Summe ueber die Verteilungen aller Komponenten wird ebenfalls c c berechnet. c c c c c c 19.11.92 Berechnung der Volumenverteilung und des Gesamtvolumens c c 03.11.93 Masse der Komponenten und Gesamtmasse c c interaktives Einlesen geaenderter Radiusgrenzen c c c c ab jetzt Version fuer OPAC. c c c c 08.11.93 interaktives Einlesen in dieser Version ausgeschaltet c c 04.12.95 Berechnung der Gesamtmasse nur bis zum max. Radius 7.5 c c 15.01.96 Radiusfeld mit SORT1 erzeugt. c c 17.01.96 Abscheideradius nicht groesser als Rmax c c 22.01.96 keine Massenberechnung fuer Wolken c c c c Stand: 22.01.96 M. Hess c ccccc ------------------------------------------------------------------c real xr(220),xr2(220,8),dnsum(220,8),dn(220,10,8),numden,mixrat real xr3(220),dvsum(220,8),dv(220,10,8),r(220),v(220) real xr4(220),xab(10) integer nx2(8),indx(220) integer nxs(8),nxc(10,8) character*1 ent character*2 chum character*4 comnam character*6 inp character*7 res character*8 file character*10 opt character*18 infile character*20 rname character*30 typnam common /file/ lfile,file,rname,infile,inp,res,opt common /prog/ nprog common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8), * mixrat(10),dens(10,8) common /masse/ smas(10,8),smag(8),vsum(10,8),vges(8),xabmax common /hum/ khum(8),ahum(8),nih,nhum(8),mhum,chum(8) common /atyp/ natyp,mcomp,ncomp(10),numden, * typnam,comnam(20) data deltar /0.015/ data xrmin /0.01/,xrmax /10./ data xabmax /7.5/ open (10,file=rname) pi=4.*atan(1.) ccccc ------------------------------------------------------------------c c Aendern der Radiusgrenzen bei Aufruf ueber RAWOPT c ccccc ------------------------------------------------------------------c if (nprog.eq.1) then write (*,*) ' Folgende Radiusgrenzen sind voreingestellt:' do ih=1,nih do ic=1,mcomp write (*,111) ahum(ih),ncomp(ic),rmin(ic,ih),rmax(ic,ih) 111 format(' f= ',f3.0,'%',' Komponente ',i2,' rmin=',f6.3, * ' rmax= ',f6.3) end do end do write (*,112) 112 format (/' sollen sie geaendert werden? (j/n) ') read (*,'(a1)') ent if (ent.eq.'j') then do ih=1,nih do ic=1,mcomp write (*,113) ahum(ih),ncomp(ic) 113 format(/' f= ',f3.0,'%',' Komponente ',i2,' rmin= ') read(*,*) rmin(ic,ih) write (*,114) ahum(ih),ncomp(ic) 114 format(' f= ',f3.0,'%',' Komponente ',i2,' rmax= ') read(*,*) rmax(ic,ih) end do end do end if end if ccccc ------------------------------------------------------------------c c Berechnung des Radiusgitters c ccccc ------------------------------------------------------------------c xr(1)=xrmin ix=2 xranf=alog10(xrmin) do while (xr(ix-1).lt.xrmax) xrl=xranf+deltar*(ix-1.) xr(ix)=10.**xrl ix=ix+1 end do nx=ix-1 ccccc ------------------------------------------------------------------c c Anfang Feuchteschleife c ccccc ------------------------------------------------------------------c do ih=1,nih do ix=1,nx xr2(ix,ih)=xr(ix) end do ccccc ------------------------------------------------------------------c c Einfuegen der Radiusgrenzen der Komponenten und des c c Abscheideradius in das Gitter c ccccc ------------------------------------------------------------------c kx=nx do ic=1,mcomp xab(ic)=xabmax kx=kx+1 xr2(kx,ih)=rmin(ic,ih) kx=kx+1 xr2(kx,ih)=rmax(ic,ih) if (rmax(ic,ih).lt.xabmax) then xab(ic)=rmax(ic,ih) end if kx=kx+1 xr2(kx,ih)=xab(ic) end do do ix=1,kx xr3(ix)=xr2(ix,ih) end do call sort1(xr3,kx,xr4,indx) xr2(1,ih)=xr4(1) lx=0 do ix=2,kx if (xr4(ix).ne.xr4(ix-1)) then xr2(ix-lx,ih)=xr4(ix) else lx=lx+1 end if end do nx2(ih)=kx-lx ccccc ------------------------------------------------------------------c c Schleife ueber die Komponenten c ccccc ------------------------------------------------------------------c do ic=1,mcomp if (ncomp(ic).gt.10) goto 1101 ! keine Massenberechnung fuer Wolken ccccc ------------------------------------------------------------------c c Berechnung der Gr��enverteilung an den Stuetzpunkten c ccccc ------------------------------------------------------------------c do ix=1,nx2(ih) if (xr2(ix,ih).ge.rmin(ic,ih).and. * xr2(ix,ih).le.rmax(ic,ih)) then dn(ix,ic,ih)=rlogn(sigma(ic),rmod(ic,ih), * numden*mixrat(ic),d,e,xr2(ix,ih)) dv(ix,ic,ih)=vlogn(sigma(ic),rmod(ic,ih), * numden*mixrat(ic),pi,e,xr2(ix,ih)) dv(ix,ic,ih)=dv(ix,ic,ih)*10.**(6) nxc(ic,ih)=nxc(ic,ih)+1 else dn(ix,ic,ih)=0. dv(ix,ic,ih)=0. end if ccccc ------------------------------------------------------------------c c Berechnung der Summe ueber alle Komponenten c ccccc ------------------------------------------------------------------c dnsum(ix,ih)=dnsum(ix,ih)+dn(ix,ic,ih) dvsum(ix,ih)=dvsum(ix,ih)+dv(ix,ic,ih) end do ccccc ------------------------------------------------------------------c c Berechnung des Aerosolvolumens und der Masse der Komponente c c Volumen in (mikrometer**3/meter**3) c c Masse in (mikrogramm/m**3) c c c c Volumen und Masse werden nur bis zum Radius XABMAX berechnet. c ccccc ------------------------------------------------------------------c do ix=1,nx2(ih) r(ix)=xr2(ix,ih) if (r(ix).le.xab(ic)) then v(ix)=dv(ix,ic,ih)/(r(ix)*alog(10.)) c v(ix)=dv(ix,ic,ih) else v(ix)=0. end if end do call gerin(r,v,erg,1,nx2(ih),ier) vsum(ic,ih)=erg smas(ic,ih)=vsum(ic,ih)*dens(ic,ih)*10.**(-6) vges(ih)=vges(ih)+vsum(ic,ih) smag(ih)=smag(ih)+smas(ic,ih) 1101 continue end do ! Komponentenschleife end do ! Feuchteschleife ccccc ------------------------------------------------------------------c c Berechnung der Volumen- und Massenmischungsverh�ltnisse c ccccc ------------------------------------------------------------------c do ih=1,nih do ic=1,mcomp if (ncomp(ic).gt.10) goto 1102 ! keine Massenberechnung fuer Wolken vsum(ic,ih)=vsum(ic,ih)/vges(ih) smas(ic,ih)=smas(ic,ih)/smag(ih) 1102 continue end do end do ccccc -----------------------------------------------------------------c c Kopf des output-files (nicht bei Aufruf ueber RAWOPT) c ccccc -----------------------------------------------------------------c if (nprog.eq.1) then write(10,100) typnam,natyp,numden,mcomp,nih 100 format(' 1 ; output of subroutine rawmic',10x, * 'calculation date: ',i2,'.',i2,'.',i4, * 2x,i2,':',i2,':',i2/ * ' season: ',a7,/ * 1x,a30,' (',i2,') ', ' with the following microphysical', * ' parameters:'/ * 15x,' number density: ',1pe10.3,' per cm**3',/, * 9x,' number of components: ',i2,/, * 9x,'number of relative humidities: ',i2,/,2x,'RH ', * 'NC SIGMA RMOD RMIN RMAX MIXRAT', * ' VMIXRAT') do ihu=1,nih do ic=1,mcomp if (ic.eq.1) then write(10,104) ahum(ihu),ncomp(ic),sigma(ic),rmod(ic,ihu), * rmin(ic,ihu),rmax(ic,ihu), * mixrat(ic),vsum(ic,ihu) else write(10,101) ncomp(ic),sigma(ic),rmod(ic,ihu), * rmin(ic,ihu),rmax(ic,ihu), * mixrat(ic),vsum(ic,ihu) end if end do end do 101 format (i10,2x,6e10.3,i5) 104 format (f5.0,i5,2x,6e10.3,i5) write (10,107) deltar 107 format(' particle radii are given in micrometers, dN/dlogr in', * ' particles/cm**3'/ * ' dV/dlogr in �m**3/m**3',5x, * ' the radius interval is dlogr= ',f7.5/ * '====================================================', * '============================') ccccc ------------------------------------------------------------------c c Ausgabe der berechneten Wertepaare der Komponenten c ccccc ------------------------------------------------------------------c do ih=1,nih write (10,106) ahum(ih),vges(ih) 106 format(' relative humidity: ',f3.0,'%', * ' total volume: ',1pe10.3,' [�m**3/m**3]') do ic=1,mcomp write (10,102) ncomp(ic),nxc(ic,ih) 102 format(' RADIUS ',' dN/dlogr',' dV/dlogr component ',i2, * ', ',i4,' values') do ix=1,nx2(ih) if (dn(ix,ic,ih).ne.0.) then write (10,105) xr2(ix,ih),dn(ix,ic,ih),dv(ix,ic,ih) 105 format(1p3e10.3) end if end do end do ccccc ------------------------------------------------------------------c c Ausgabe der berechneten Wertepaare fuer die gesamte Verteilung c ccccc ------------------------------------------------------------------c do ix=1,nx2(ih) if (dnsum(ix,ih).ne.0.) then nxs(ih)=nxs(ih)+1 end if end do write (10,103) natyp,nxs(ih) 103 format(' RADIUS ',' dN/dlogr',' aerosol type ',i2, * ', ',i4,' values') do ix=1,nx2(ih) if (dnsum(ix,ih).ne.0.) then write (10,105) xr2(ix,ih),dnsum(ix,ih),dvsum(ix,ih) end if end do ccccc ------------------------------------------------------------------c c Ende der Feuchteschleife c ccccc ------------------------------------------------------------------c end do end if return end FUNCTION RLOGN(SIG,RO,VN,D,E,X) CCCCC -----------------------------------------------------------------C C LOG-NORMAL-VERTEILUNG DN/DLOGR IN cm**-3 C c DIE RADIEN X UND RO MUESSEN IN um ANGEGEBEN WERDEN, C c VN in cm**-3 c CCCCC -----------------------------------------------------------------C PI=4.*ATAN(1.) A=VN/(SQRT(2.*PI)*ALOG10(SIG)) B=ALOG10(RO) C=-2.*(ALOG10(SIG))**2 RLOGN=A*EXP((ALOG10(X)-B)**2/C) RETURN END FUNCTION VLOGN(SIG,RO,VN,PI,E,X) CCCCC -----------------------------------------------------------------C C LOGNORMALE VOLUMENVERTEILUNG DV/DLOGR IN um**3 cm**-3 C c DIE RADIEN X UND RO MUESSEN IN um ANGEGEBEN WERDEN, C c VN in cm**-3 c CCCCC -----------------------------------------------------------------C A=VN/(SQRT(2.*PI)*ALOG10(SIG)) B=ALOG10(RO) C=-2.*(ALOG10(SIG))**2 RLOGN=A*EXP((ALOG10(X)-B)**2/C) VLOGN=4./3.*PI*X**3.*RLOGN RETURN END CCCCC *****************************************************************C subroutine optraw (ihum) C *****************************************************************C C c c neue Einleseroutine fuer Aerosoltypen (RAWOPT), nicht mehr c c gleich der Version OPTCOM fuer ATLOPT. c c c c 04.05.94 Ausschluss der Quellung bei Russ (Komponente 3) c c 15.05.94 Einlesen der Cirrus-Daten c c 04.12.95 Einlesen der Brechungsindizes c c 07.12.95 neue Namen fuer Dateien mit Komponentendaten c c c c 23.09.97 small changes for OPAC 3.0 c c 21.10.97 read new data format in ../optdat/ for OPAC 3.1 c c c c Stand: 21.10.97 M. Hess c CCCCC -----------------------------------------------------------------C logical ende integer prnr,acnr,njc,rht real n,numden real ex(61,6),sc(61,6),ab(61,6),si(61,6),as(61,6),ba(61,6) real ph(61,6,167),br(61,6),bi(61,6) character*1 dum character*2 chum character*3 atn,pat character*4 comnam character*6 inp character*7 res character*8 file character*10 opt,dum2 character*16 tap character*18 infile character*20 rname character*30 typnam common /file/ lfile,file,rname,infile,inp,res,opt common /atyp/ natyp,mcomp,ncomp(10),numden, * typnam,comnam(20) common /wavel/ mlamb,alamb(61),niw,wlamb(61), * il035,il05,il055,il08,alambp(61),niwp common /hum/ khum(8),ahum(8),nih,nhum(8),mhum,chum(8) common /mipoi/ latx,lonx,nl,prnr,rht(2),n(2), * njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2) common /oppoi/ ext(1,6),sca(1,6),abs(1,6),sis(1,6),asy(1,6), * bac(1,6),pha(1,6,167),ext05,bre(1,6),bim(1,6) common /angle/ jnangle(167),angle(167),nia,ntheta save ccccc -----------------------------------------------------------------c c Schleife �ber alle vorkommenden Komponenten c ccccc -----------------------------------------------------------------c do il=1,nl if (nih.eq.0) then do ihu=1,mhum if (nh(il).eq.ihu) then khum(ihum)=ihu end if end do end if ext05=0. do ic=1,njc(il) jc=acnr(ic,il) ccccc -----------------------------------------------------------------c c Ausschlu� der Quellung bei insoluble, Russ und den c c mineralischen Komponenten und bei den Wolken c c ccccc -----------------------------------------------------------------c if ( jc.eq.1.or.jc.eq.3.or.(jc.ge.6.and.jc.le.9).or. * jc.gt.10 ) then iht=1 else iht=khum(ihum) end if ccccc -----------------------------------------------------------------c c Bestimmung des Filenamens der gesuchten Komponente aus c c Komponentennummer und Feuchteklasse c ccccc -----------------------------------------------------------------c tap(1:10)=opt tap(11:14)=comnam(jc) tap(15:16)=chum(iht) c print*,' tap= ',tap,jc,khum(ihum) ntap=70 open (ntap,file=tap,iostat=ios) c if (ios.ne.0) then c print*,' error while opening file ',tap c print*,'ios=',ios c stop c end if do iline=1,100 read (ntap,220) dum2 if (dum2.eq.'# optical ') then goto 2002 end if end do 2002 continue do iline=1,5 read (ntap,200) dum end do do ilam=1,mlamb read (ntap,500) rlamb,extco,scaco,absco,sisca,asymf, * exn,refr,refi 500 format(2x,7e10.3,2e11.3) ex(ilam,ic)=extco sc(ilam,ic)=scaco ab(ilam,ic)=absco if (rlamb.eq.0.55) then ext05=ext05+ex(ilam,ic)*acmr(ic,1)*numden end if si(ilam,ic)=sisca as(ilam,ic)=asymf br(ilam,ic)=refr bi(ilam,ic)=refi if (rlamb.ne.wlamb(ilam)) then print*,' ' print*,' !! ATTENTION !!' print*,' ' print*,' you have selected the wrong wavelength table' print*,' Please change selection!' print*,' ' print*,' Press <ENTER> to return to Input Menu' read (*,'(a1)') dum stop end if end do read (ntap,'(7(/))') it=1 ende=.false. do while (.not.ende) read (ntap,510,end=511) * angle(it),(ph(ilam,ic,it),ilam=1,mlamb) 510 format(e11.3,1x,70e10.3) it=it+1 end do 511 ntheta=it-1 do ilam=1,mlamb do it=1,ntheta ph(ilam,ic,it)=ph(ilam,ic,it) end do ba(ilam,ic)=ph(ilam,ic,ntheta) end do close (ntap) end do end do ccccc -----------------------------------------------------------------c c Formate c ccccc -----------------------------------------------------------------c 100 format(8e10.3) 200 format(a1) 220 format(a10) 300 format(15x,f6.3,/,8x,e10.4,7x,e10.4,7x,e10.4,5x,f7.4,7x,f6.4,/) 301 format(53x,f6.3,8x,e10.3,/,8x,f6.3//) 302 format(8x,f6.3,8x,e10.3,9x,f6.3//) 304 format(28x,f7.3//) 303 format(7x,e10.3,7x,e10.3,7x,e10.3,7x,f7.4/) 305 format(6x,a4) 400 format(12(/)) 404 format(21(/)) 1010 format(70X,e10.3) return entry optdat (ilamb) ccccc -----------------------------------------------------------------c c Auswahl der gewuenschten optischen Daten aus den eingelesenen c ccccc -----------------------------------------------------------------c do ic=1,njc(1) do il=1,mlamb if (alamb(ilamb).eq.wlamb(il)) then ext(1,ic)=ex(il,ic) sca(1,ic)=sc(il,ic) abs(1,ic)=ab(il,ic) sis(1,ic)=si(il,ic) asy(1,ic)=as(il,ic) bac(1,ic)=ba(il,ic) bre(1,ic)=br(il,ic) bim(1,ic)=bi(il,ic) do it=1,ntheta pha(1,ic,it)=ph(il,ic,it) end do end if end do end do return end CCCCC *****************************************************************C SUBROUTINE OPTPAR (ilamb,ihum) C *****************************************************************C C c C -----------------------------------------------------------------C C Berechnung und Ausdruck der gew�nschten optischen Parameter C C c C 04.11.93 Parameter SCARA, ABSRA, OMERA eingefuehrt. c C 09.11.93 Version fuer OPAC (Aufruf OUT4 entfernt) c C 13.05.94 optische Dicke fuer alle Wellenlaengen c C 14.05.94 Indizierung der opt. Param. fuer Ausgabe korrigiert. c C 28.07.94 normierter Ext.koeff. eingebaut c C 04.12.95 Brechungsindizes eingebaut c C 02.01.96 Berechnung der opt. Dicke ueber Schichttypen. c C 10.01.96 neue Berechnung der opt. Dicke (ohne Heff) c C 22.01.96 EXTRA statt SCARA c C c C 26.09.97 normalized extinction corrected. c C c C Stand: 26.09.97 M. Hess c CCCCC -----------------------------------------------------------------C integer prnr,acnr,rht real n,numden REAL EXTN(2),ABSN(2),SCAN(2),PF18N(2),supf(167),phafu(167,2) REAL EXTA(2),ABSA(2),SCAA(2),SSA(2),ASF(2),PF18A(2) real scar(2),absr(2),omer(2) character*3 atn,pat character*8 optnam,optun,opanam character*4 comnam character*30 typnam common /atyp/ natyp,mcomp,ncomp(10),numden, * typnam,comnam(20) common /layer/ mlay,nltyp(10),parlay(10,5),boundl(10), * boundu(10) common /norm/ norm,mixnor common /opar/ mopar,jnopar(20),nop,opanam(20),optnam(20), * optun(20) common /wavel/ mlamb,alamb(61),niw,wlamb(61), * il035,il05,il055,il08,alambp(61),niwp common /mipoi/ latx,lonx,nl,prnr,rht(2),n(2), * njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2) common /oppoi/ ext(1,6),sca(1,6),abs(1,6),sis(1,6),asy(1,6), * bac(1,6),pha(1,6,167),ext05,bre(1,6),bim(1,6) COMMON /FTASTR/ EXTFTA(61),EXTSTR(61),extmit(61) common /out/ oparam(16,2),phaf(167,2),indop(16) common /sotin/ exi(61),ssi(61),asi(61) common /prog/ nprog common /angle/ jnangle(167),angle(167),nia,ntheta common /masse/ smas(10,8),smag(8),vsum(10,8),vges(8),xabmax CCCCC ------------------------------------------------------------------C C MISCHEN DES AEROSOL-TYPS C C SUMM(E,A,S) : SUMME EXTINCTION, ABSORPTION, SCATTERING C C SUPF18 : SUMME DES RUECKSTREUKOEFFIZIENTEN C C SUMASF : ZWISCHENSUMME DES ASYMMETRIEFAKTORS (ASF) C C SUMASF : ZWISCHENSUMME DER SINGLE SCAT. ALB. (SSA) C CCCCC ------------------------------------------------------------------C DO 10 L=1,NL SUMME = 0. SUMMA = 0. SUMMS = 0. SUMSSA = 0. SUMASF = 0. SUPF18 = 0. c if (jnopar(10).eq.1) then do it=1,ntheta supf(it)=0. end do c end if DO 20 JC=1,NJC(L) SUMME = SUMME + ACMR(JC,L)*EXT(l,jc) SUMMA = SUMMA + ACMR(JC,L)*ABS(l,jc) SUMMS = SUMMS + ACMR(JC,L)*SCA(l,jc) SUMSSA = SUMSSA + ACMR(JC,L)*sis(l,jc) + *EXT(l,jc) SUMASF = SUMASF + ACMR(JC,L)*asy(l,jc) + *SCA(l,jc) SUPF18 = SUPF18 + ACMR(JC,L)*bac(l,jc) c if (jnopar(10).eq.1) then do it=1,ntheta supf(it)=supf(it)+acmr(jc,l)*pha(l,jc,it) end do c end if 20 CONTINUE CCCCC -----------------------------------------------------------------C C Normierte optische Parameter c CCCCC -----------------------------------------------------------------C EXTN(L) = SUMME ABSN(L) = SUMMA SCAN(L) = SUMMS PF18N(L) = SUPF18 if (jnopar(10).eq.1) then do it=1,ntheta phafu(it,l)=supf(it) end do end if SSA(L) = SUMSSA/SUMME ASF(L) = SUMASF/SUMMS CCCCC -----------------------------------------------------------------C C ABSOLUTE OPTISCHE PARAMETER C CCCCC -----------------------------------------------------------------C EXTA(L)= EXTN(L) * N(L) ABSA(L)= ABSN(L) * N(L) SCAA(L)= SCAN(L) * N(L) PF18A(L) = PF18N(L)* N(L) if (jnopar(10).eq.1.and.norm.eq.1) then do it=1,ntheta phafu(it,l)=phafu(it,l)*n(l) end do end if if (norm.eq.1) then EXTN(L)= EXTA(L) ABSN(L)= ABSA(L) SCAN(L)= SCAA(L) PF18N(L) = PF18A(L) end if c if (jnopar(15).eq.1.or.jnopar(16).eq.1) then exi(ilamb)=extn(1) ssi(ilamb)=ssa(1) asi(ilamb)=asf(1) c end if c if (jnopar(10).eq.1) then itp=1 do it=1,ntheta c if (jnangle(it).eq.1) then phaf(itp,l)=phafu(it,l) itp=itp+1 c end if end do c end if if (jnopar(11).eq.1) then scar(l)=exta(l)/smag(ihum)*1000. ! Einheit m**2/g end if if (jnopar(12).eq.1) then absr(l)=absa(l)/smag(ihum)*1000. end if c if (jnopar(13).eq.1) then kc=0 do jc=1,njc(l) if (ncomp(jc).eq.3) kc=jc end do if (kc.ne.0) then omer(l)=smas(kc,ihum)/ssa(l) else omer(l)=99. end if c end if CCCCC -----------------------------------------------------------------C C AUSGABE DER DATEN c CCCCC -----------------------------------------------------------------C iop=0 kop=0 if (jnopar(1).eq.1) then iop=iop+1 indop(iop)=1 oparam(iop,l)=extn(l) end if if (jnopar(2).eq.1) then iop=iop+1 indop(iop)=2 oparam(iop,l)=scan(l) end if if (jnopar(3).eq.1) then iop=iop+1 indop(iop)=3 oparam(iop,l)=absn(l) end if if (jnopar(4).eq.1) then iop=iop+1 indop(iop)=4 oparam(iop,l)=ssa(l) end if if (jnopar(5).eq.1) then iop=iop+1 indop(iop)=5 oparam(iop,l)=asf(l) end if if (jnopar(9).eq.1) then iop=iop+1 indop(iop)=9 oparam(iop,l)=exta(l)/pf18a(l) end if if (jnopar(11).eq.1) then iop=iop+1 indop(iop)=11 oparam(iop,l)=scar(l) end if if (jnopar(12).eq.1) then iop=iop+1 indop(iop)=12 oparam(iop,l)=absr(l) end if if (jnopar(13).eq.1) then iop=iop+1 indop(iop)=13 oparam(iop,l)=omer(l) end if if (jnopar(14).eq.1) then iop=iop+1 indop(iop)=14 oparam(iop,l)=extn(l)/ext05 if (norm.eq.0) oparam(iop,l)=oparam(iop,l)*n(l) end if if (jnopar(18).eq.1) then ! Brechungsindizes iop=iop+1 indop(iop)=18 oparam(iop,l)=bre(1,1) iop=iop+1 indop(iop)=19 oparam(iop,l)=bim(1,1) end if 10 CONTINUE CCCCC -----------------------------------------------------------------C C OPTISCHE DICKE (nur bei absoluten Werten) c CCCCC -----------------------------------------------------------------C if (jnopar(6).eq.1.or.jnopar(7).eq.1.or.jnopar(8).eq.1) then if (norm.eq.1) then ! nur fuer Absolutwerte CCCCC -----------------------------------------------------------------C c Bestimmung von HM, HFTA, HSTR, EXTFTA, EXTSTR aus den c c eingelesenen Werten in /layer/ fuer RAWOPT c CCCCC -----------------------------------------------------------------C odepth=0. do il=1,mlay if (nltyp(il).eq.1) then ! mixing layer heff = parlay(il,1) * * (exp(-boundl(il)/parlay(il,1))- * exp(-boundu(il)/parlay(il,1)) ) odepth = odepth + extn(1) * heff c print*,'OD(1)= ',odepth,' Heff= ',heff,' ext= ',extn(1) else if (nltyp(il).eq.2) then ! mineral transported heff=(boundu(il)-boundl(il)) odepth = odepth + extmit(ilamb) * parlay(il,1) * heff c print*,'OD(2)= ',odepth,' Heff= ',heff,' ext= ',extmit(ilamb) else if (nltyp(il).eq.3) then ! free troposphere heff = parlay(il,1) * * (exp(-boundl(il)/parlay(il,1))- * exp(-boundu(il)/parlay(il,1)) ) odepth = odepth + extfta(ilamb) * heff c print*,'OD(3)= ',odepth,' Heff= ',heff,' ext= ',extfta(ilamb) else if (nltyp(il).eq.4) then ! stratosphere heff= (boundu(il)-boundl(il)) odepth = odepth + extstr(ilamb) * heff c print*,'OD(4)= ',odepth,' heff= ',heff,' ext= ',extstr(ilamb) else if (nltyp(il).eq.5) then ! cloud heff= (boundu(il)-boundl(il)) odepth = odepth + extn(1) * heff c print*,'OD(5)= ',odepth,' heff= ',heff,' ext=',extn(1) end if end do odeptha=odepth/alog(10.) turbr=0.008569*alamb(ilamb)**(-4)*(1.+0.0113*alamb(ilamb)** * (-2)+0.00013*alamb(ilamb)**(-4)) turbf=(odepth+turbr)/turbr if (jnopar(6).eq.1) then iop=iop+1 indop(iop)=6 oparam(iop,1)=odepth end if if (jnopar(7).eq.1) then iop=iop+1 indop(iop)=7 oparam(iop,1)=odeptha end if if (jnopar(8).eq.1) then iop=iop+1 indop(iop)=8 oparam(iop,1)=turbf end if else print*,' ' print*, ' !! ATTENTION !!' print*,' ' print*, ' optical depths may not be calculated together' print*, ' with normalized coefficients!' print*, ' Please change setting to absolute! ' print*,' ' print*, ' Press <ENTER> to return to Input Menu' read (*,'(a1)') dum stop end if end if if (jnopar(15).eq.1) then do ilp=1,niwp if (alambp(ilp).eq.alamb(ilamb)) call out2(ilp,ilamb,iop) end do else call out2(ilamb,ilamb,iop) end if RETURN END subroutine intco ccccc ------------------------------------------------------------------c c Integration von berechneten optischen Groessen mit Wichtung c c mit der Solarkonstanten im solaren Spektralbereich und Wichtung c c mit Planck(300 K) im terrestrischen. c c c c uebernommen von SOLINT.FOR vom 19.04.94 c c c c ergaenzt um Berechnung der Angstrom-Koeffizienten bei 0.35/0.5 c c und 0.5/0.8 Mikrometer. c c c c 25.01.95 c c c c 21.01.96 Berechnung der Sichtweite c c c c Stand: 21.01.96 M. Hess c ccccc ------------------------------------------------------------------c real welsk(87),exis(87),ssis(87),asis(87) real sgrenz(100) real sol(100),terr(100),wter(100),fluss(100) real oint1(100),oint2(100),gint1(100),gint2(100),aint1(100) real atint1(100),atint2(100),aint2(100),eint1(100) character*8 optnam,optun,opanam common /opar/ mopar,jnopar(20),nop,opanam(20),optnam(20), * optun(20) common /wavel/ mlamb,alamb(61),niw,wlamb(61), * il035,il05,il055,il08,alambp(61),niwp common /oppoi/ ext(1,6),sca(1,6),abs(1,6),sis(1,6),asy(1,6), * bac(1,6),pha(1,6,167),ext05,bre(1,6),bim(1,6) common /out/ oparam(16,2),phaf(167,2),indop(16) common /sotin/ exi(61),ssi(61),asi(61) if(jnopar(15).eq.0.and.jnopar(16).eq.0.and.jnopar(17).eq.0) return ccccc ------------------------------------------------------------------c c solar und terrestr. integrierte Groessen c ccccc ------------------------------------------------------------------c open (5,file='wel.dat') read (5,*) nwelsk do iw=1,nwelsk read(5,*) welsk(iw) end do close(5) if (jnopar(15).eq.1) then call trapez(alamb,exi,1,mlamb,welsk,exis,11,nwelsk-5,ier) call trapez(alamb,ssi,1,mlamb,welsk,ssis,11,nwelsk-5,ier) call trapez(alamb,asi,1,mlamb,welsk,asis,11,nwelsk-5,ier) do il=1,nwelsk if (welsk(il).eq.0.3) i03=il if (welsk(il).eq.3.28) i3=il if (welsk(il).eq.7.992) i8=il if (welsk(il).eq.15.331) i15=il end do ccccc ------------------------------------------------------------------c c Einlesen der Solarkonstanten c ccccc ------------------------------------------------------------------c open (5,file='solar.dat') read (5,'(a1)') dum read (5,*) (sgrenz(ig),ig=1,38) read (5,'(a1)') dum read (5,*) (fluss(iv),iv=1,37) do il=1,37 sol(il)=fluss(il)/(sgrenz(il+1)-sgrenz(il))/4. end do close (5) ccccc ------------------------------------------------------------------c c Einlesen der terrestrischen Strahlung fuer T=300 K c ccccc ------------------------------------------------------------------c open (5,file='terr.dat') do il=1,20 read (5,*) wter(il),terr(il) end do close (5) ccccc ------------------------------------------------------------------c c Integration c ccccc ------------------------------------------------------------------c do il=1,37 oint1(il)=ssis(il)*exis(il)*sol(il) oint2(il)=exis(il)*sol(il) gint1(il)=asis(il)*exis(il)*ssis(il)*sol(il) gint2(il)=exis(il)*ssis(il)*sol(il) aint1(il)=exis(il)*(1.-ssis(il))*sol(il) aint2(il)=sol(il) eint1(il)=exis(il)*sol(il) end do do il=1,20 atint1(il)=exis(il+i8-1)*(1.-ssis(il+i8-1))*terr(il) atint2(il)=terr(il) end do call gerin (welsk,oint1,o1erg,11,37,ier) call gerin (welsk,oint2,o2erg,11,37,ier) call gerin (welsk,gint1,g1erg,11,37,ier) call gerin (welsk,gint2,g2erg,11,37,ier) call gerin (welsk,aint1,a1erg,11,37,ier) call gerin (welsk,aint2,a2erg,11,37,ier) call gerin (welsk,eint1,eerg,11,37,ier) call gerin (wter,atint1,at1erg,1,20,ier) call gerin (wter,atint2,at2erg,1,20,ier) oint=o1erg/o2erg gint=g1erg/g2erg aint=a1erg/a2erg eint=eerg/a2erg atint=at1erg/at2erg end if ccccc ------------------------------------------------------------------c c Angstrom-Koeffizienten c ccccc ------------------------------------------------------------------c if (jnopar(16).eq.1) then al1=(alog10(exi(il05))-alog10(exi(il035)))/ * (alog10(alamb(il035))-alog10(alamb(il05))) al2=(alog10(exi(il08))-alog10(exi(il05)))/ * (alog10(alamb(il05))-alog10(alamb(il08))) be1=exi(il035)*alamb(il035)**(al1) be2=exi(il05)*alamb(il05)**(al2) end if ccccc ------------------------------------------------------------------c c Sichtweite c ccccc ------------------------------------------------------------------c if (jnopar(17).eq.1) then visib=3.0/(ext05+0.01159) end if ccccc ------------------------------------------------------------------c c Ausgabe in Output-File c ccccc ------------------------------------------------------------------c write (10,107) 107 format(' ---------------------------------------------------', * '-------------------------') if (jnopar(15).eq.1) then write (10,120) eint,oint,gint,aint,atint 120 format(' values weighted with solar spectrum', * ' (0.3-3.3 micron):'/ * e10.3,' : extinction coefficient'/ * e10.3,' : single scattering albedo'/ * e10.3,' : asymmetry parameter'/ * e10.3,' : absorption coefficient'/ * ' values weighted with terrestrial spectrum at 300K', * ' (8-15 micron):'/ * e10.3,' : absorption coefficient') end if if (jnopar(16).eq.1) then write (10,130) al1,be1,al2,be2 130 format(' Angstrom coefficients calculated at ', * '350 nm and 500 nm:'/ * e10.3,' : alpha'/ * e10.3,' : beta'/ * ' Angstrom coefficients calculated at ', * '500 nm and 800 nm:'/ * e10.3,' : alpha'/ * e10.3,' : beta') end if if (jnopar(17).eq.1) then write (10,110) visib 110 format(' Visibility: ',f6.2,' [km]') end if return end SUBROUTINE LANG2(STEXT,mtext,LTEXT) CCCCC -----------------------------------------------------------------C C LAENGE DES TEXTES WIRD BERECHNET (ENDE SOBALD 2 BLANKS C C HINTEREINANDER AUFTAUCHEN) C c c c 19.06.94 maximale Laenge mtext wird als Input mit uebergeben c c 03.03.95 Fehler korrigiert, falls ltext=mtext-1 c c 16.04.95 ganz neue Version c c c c 16.04.95 M. Hess c CCCCC -----------------------------------------------------------------C CHARACTER*(*) STEXT if (stext(mtext:mtext).ne.' ') then ltext=mtext else if (stext(mtext-1:mtext-1).eq.' ') then DO J=mtext,2,-1 IF(STEXT(J:J).EQ.' '.AND.STEXT(J-1:J-1).EQ.' ') THEN LTEXT=J-2 END IF end do else LTEXT=mtext-1 end if RETURN END ccccc *****************************************************************c subroutine sort1(werte,nx,wsort,index) c *****************************************************************c c c c -----------------------------------------------------------------c c Sortieren eines 1-dimensionalen Feldes unter Verwendung von c c MINMAX1. c c c c wsort : nach der Groesse sortiertes Feld c c index : zugehoeriges Feld der urspruenglichen Indizes c c c c nx darf maximal 1000 sein! c c c c Stand: 17.09.93 M. Hess c ccccc -----------------------------------------------------------------c integer index(nx) real werte(nx),wsort(nx),whilf(1000) logical ende do ix=1,nx whilf(ix)=werte(ix) end do do ix=1,nx-1 call minmax1(whilf,nx,wmin,wmax) if (ix.eq.1) then ! Bestimmung des groessten Wertes wsort(nx)=wmax ende=.false. i=0 do while (.not.ende) i=i+1 if (whilf(i).eq.wmax) then index(nx)=i ende=.true. end if end do end if wsort(ix)=wmin ! Sortieren des restlichen Feldes ende=.false. i=0 do while (.not.ende) i=i+1 if (whilf(i).eq.wmin) then whilf(i)=wmax index(ix)=i ende=.true. end if end do end do return end ccccc *****************************************************************c subroutine minmax1(werte,nx,wmin,wmax) c *****************************************************************c c c c -----------------------------------------------------------------c c Bestimmung des Maximums und des Minimums aus einem c c 1-dimensionalen REAL-Feld c c c c Stand: 14.02.92 M. Hess c ccccc -----------------------------------------------------------------c real werte(nx) wmin=werte(1) wmax=werte(1) do ix=1,nx if (werte(ix).ge.wmax) wmax=werte(ix) if (werte(ix).le.wmin) wmin=werte(ix) end do return end SUBROUTINE TRAPEZ(X,Y,IANFA,IENDA,U,V,IANFN,IENDN,IERROR) C C STAND VOM 15. JULI 1975 C LINEARE INTERPOLATION EINES ORDINATENFELDES V(I), FUER EIN C ABSZISSENFELD U(I) FUER I=IANFN,IENDN AUS EINEM AUSGANGSFELD C (REFERENZFELD) Y(I),X(I), FUER I=IANFA,IENDA C U-FELD UND X-FELD MUESSEN GEORDNET SEIN C C IERROR=0 IM NORMALFALL C IERROR=2 WENN DAS AUSGANGSFELD KLEINER IST ALS DAS ZU INTERPOLIERENDE C FELD C DIMENSION X(IENDA),Y(IENDA),U(IENDN),V(IENDN) IERROR=0 C IDREHX=0 IF (X(IANFA).LT.X(IENDA)) GOTO 2 C DAS X-FELD WIRD ANSTEIGEND GEORDNET, DAS Y-FELD ENTSPRECHEND IORDX=(IENDA-IANFA+1)/2 DO 14 I=IANFA,IORDX XUMORD=X(I) X(I)=X(IENDA+1-I) X(IENDA+1-I)=XUMORD YUMORD=Y(I) Y(I)=Y(IENDA+1-I) 14 Y(IENDA+1-I)=YUMORD IDREHX=1 C 2 IDREHU=0 IF (U(IANFN).LT.U(IENDN)) GOTO 12 C DAS U-FELD WIRD ANSTEIGEND GEORDNET IORDU=(IENDN-IANFN+1)/2 DO 16 I=IANFN,IORDU UUMORD=U(I) U(I)=U(IENDN+1-I) 16 U(IENDN+1-I)=UUMORD IDREHU=1 C 12 IEXPO=0 IANFNN=IANFN 3 IF (U(IANFNN).GE.X(IANFA)) GOTO 13 C EXTRAPOLATION BEI X(IANFA) VERHINDERN. IEXPO=1 IANFNN=IANFNN+1 GOTO 3 13 IENDNN=IENDN 4 IF (U(IENDNN).LE.X(IENDA)) GOTO 5 C EXTRAPOLATION BEI X(IENDA) VERHINDERN. IEXPO=1 IENDNN=IENDNN-1 GOTO 4 5 IF(IEXPO.EQ.1) WRITE(6,98) X(IANFA),X(IENDA),U(IANFN),U(IENDN), 1 U(IANFNN),U(IENDNN) 98 FORMAT(1H0,5X,'DAS REFERENZFELD GEHT VON X(IANFA) =',1PE10.3,1X, 1'BIS X(IENDA) =',E10.3/1H ,'INTERPOLIERT WERDEN SOLLTE VON U(IANFN 2) =',E10.3,' BIS U(IENDN) =',E10.3/1H ,5X,'WERTE AUSSERHALB VON', 3E10.3,' UND',E10.3,' SIND DAHER NICHT BERECHNET.'/1H ) C IF (IEXPO.EQ.1) IERROR=2 IANF=IANFA DO 20 J=IANFNN,IENDNN DO 10 I=IANF,IENDA IF (X(I)-U(J)) 10,7,8 7 V(J)=Y(I) GOTO 22 8 V(J)=Y(I-1)+(Y(I)-Y(I-1))/(X(I)-X(I-1))*(U(J)-X(I-1)) GOTO 22 10 CONTINUE GOTO 20 22 IANF=I 20 CONTINUE IF (IDREHX.EQ.0) GOTO 21 C DAS X- UND Y-FELD WERDEN IN DIE AUSGANGS POSITION ZURUECKGESETZT DO 18 I=IANFA,IORDX XUMORD=X(I) X(I)=X(IENDA+1-I) X(IENDA+1-I)=XUMORD YUMORD=Y(I) Y(I)=Y(IENDA+1-I) 18 Y(IENDA+1-I)=YUMORD C 21 IF (IDREHU.EQ.0) RETURN DO 19 I=IANFN,IORDU UUMORD=U(I) U(I)=U(IENDN+1-I) U(IENDN+1-I)=UUMORD VUMORD=V(I) V(I)=V(IENDN+1-I) 19 V(IENDN+1-I)=VUMORD C RETURN END SUBROUTINE GERIN (H,F,ERG,IA,IE,IERROR) C STAND VOM 20. FEB. 1974 C TRAPEZINTEGRATION VON NICHT AEQUIDISTANTEN TABELL. FUNKTIONEN C ANZAHL DER ABSZISSENPUNKTE DARF GERADE UND UNGERADE SEIN C H=ABSZISSE,F=ORDINATE,ERG=INTEGRAL C INTEGRIERT WIRD VON H(IA) BIS H(IE) C FALLS NUR EIN STUETZWERT VORHANDEN (IA=IE) , WIRD DEM INTEGRAL C DIE ORDINATE DES STUETZWERTES ZUGEORDNET. C IERROR =0 IM NORMALFALL C IERROR=1, WENN IA GROESSER IE IST C IERROR=2 WENN IA=IE IST C IERROR=3 WENN ABSZISSENDIFFERENZEN ZU KLEIN WERDEN C DIMENSION H(IE),F(IE) C IERROR=0 ERG=0.0 IF (IE-IA.GE.0) GOTO 6 WRITE(6,5) IA,IE 5 FORMAT('0 IN DER SUBROUTINE GERIN IST IA=',I3,'GROESSER IE=',I3) IERROR=1 RETURN C 6 IF (IE-IA.NE.0) GOTO 3 WRITE (6,15) IA 15 FORMAT('0 IN DER SUBROUTINE GERIN IST IA=IE=',I3,' UND ERG=F(IA)') ERG = F(IA) IERROR=2 RETURN C 3 IAP=IA+1 DO 1 I=IAP,IE DIFH=H(I)-H(I-1) IF (ABS(DIFH/H(I)).GT.1.E-4) GOTO 4 WRITE (6,25) I 25 FORMAT ('0 IN DER SUBROUTINE GERIN WIRD DIFH KLEINER ALS 1.E-4 BEI 1 I=',I4) C DIE WARNUNG ENTSPRICHT EINEM FEHLER VON MAXIMAL 1 PROZENT IERROR=3 C 4 ERG=ERG+(F(I)+F(I-1))*DIFH 1 CONTINUE ERG=ERG/2. C RETURN END