source: trunk/hftp/iasi/aerosol/prog/opac.f @ 1

Last change on this file since 1 was 1, checked in by cbipsl, 18 years ago

Geisa inital import
File size: 80.4 KB
Line 
1      program opac
2
3ccccc ------------------------------------------------------------------c
4c     extracts and calculates optical properties of aerosols and clouds c
5c     from the data stored in directory ../optdat/.                     c
6c                                                                       c
7c     Version 3.0                                                       c
8c     ------------                                                      c
9c     01.10.97 new version based on Version 2.0 of 17.04.97             c
10c                                                                       c
11c     Version 3.1                                                       c
12c     ------------                                                      c
13c     21.10.97 new data format in ../optdat/                            c
14c     30.10.97 corrected: units [1/m] -> [1/km]                         c
15c     06.11.97 opt.depth for clouds corrected.                          c
16c     14.11.97 tested with Linux, SunOS, Dec Ultrix, HP Unix, DOS       c
17c                                                                       c
18c     14.11.97                                                  M. Hess c
19ccccc ------------------------------------------------------------------c
20
21      character*2  chum
22      character*4  comnam
23      character*6  inp
24      character*7  res
25      character*8  file,optnam,optun,opanam
26      character*10 opt
27      character*18 infile
28      character*20 rname
29      character*30 typnam,tname,tnew
30
31      logical wa,wb
32
33      real mixrat,numden
34     
35      integer nlt(10)
36
37      common /prog/   nprog
38      common /file/   lfile,file,rname,infile,inp,res,opt
39      common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8),
40     *                mixrat(10),dens(10,8)
41      common /mipaco/ sigmac(20),rminc(20,8),rmaxc(20,8),rmodc(20,8),
42     *                densc(20,8)
43      common /mipaty/ nco(20),mco(20,10),dnumb(20,10),nlay(20),
44     *                hmint(20,4),hmaxt(20,4),part(20,4),tname(20)
45      common /input/  mtyp,nuco(10),dnum(10),mcom,
46     *                hmin(5),hmax(5),hpar(5),
47     *                nwavea,indwa(100),wavea(100),
48     *                nwaveb,indwb(100),waveb(100),
49     *                nhc,indhum(8),nopt,indopt(20),tnew
50      common /atyp/   natyp,mcomp,ncomp(10),numden,
51     *                typnam,comnam(20)
52      common /layer/  mlay,nltyp(10),parlay(10,5),boundl(10),
53     *                boundu(10)
54      common /opar/   mopar,jnopar(20),nop,opanam(20),optnam(20),
55     *                optun(20)
56      common /wavel/  mlamb,alamb(61),niw,wlamb(61),
57     *                il035,il05,il055,il08,alambp(61),niwp
58      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
59      common /norm/   norm,mixnor
60
61      data norm  /1/
62      data nprog /0/
63      data nlt   /1,2,3,4,5,0,0,0,0,0/
64      data nhum  /0,50,70,80,90,95,98,99/,mhum/8/
65      data chum  /'00','50','70','80','90','95','98','99'/
66
67      data optnam /'ext.coef','sca.coef','abs.coef','sisc.alb',
68     *             'asym.par','op.depth',
69     *             '        ','turb.fac','li.ratio','pha.func',
70     *             'ext.rat ','abs.rat ',
71     *             '        ','ext.norm','        ','        ',
72     *             'visibil ','ref.real',
73     *             'ref.imag','        '/
74
75      data optun  /' [1/km] ',' [1/km] ',' [1/km] ','        ',
76     *             '        ','        ',
77     *             '        ','        ','  [sr]  ',' [1/km] ',
78     *             '[m^2/g] ','[m^2/g] ',
79     *             '        ','        ','        ','        ',
80     *             '  [km]  ','        ',
81     *             '        ','        '/
82
83      data comnam /'inso','waso','soot','ssam','sscm','minm','miam',
84     *             'micm','mitr','suso','stco','stma','cucc','cucp',
85     *             'cuma','fogr','cir1','cir2','cir3','    '/
86
87
88      print*,'********************************************************'
89      print*,'*  Optical Properties of Aerosols and Clouds OPAC 3.1  *'
90      print*,'*  --------------------------------------------------  *'
91      print*,'*                                                      *'
92      print*,'*  The contents of OPAC has been described in:         *'
93      print*,'*  M. Hess, P. Koepke and I. Schult (1997):            *'
94      print*,'*       "Optical Properties of Aerosols and Clouds:    *'
95      print*,'*        The software package OPAC"                    *'
96      print*,'*  submitted to Bull. Am. Meteor. Soc.                 *'     
97      print*,'*                                                      *'
98      print*,'*  last update: 14.11.97                     M. Hess   *'   
99      print*,'********************************************************'
100      print*,' '
101      print*,' '
102           
103ccccc ------------------------------------------------------------------c
104c     read data from configuration file opac.cfg                        c
105ccccc ------------------------------------------------------------------c
106
107      call readcfg
108
109ccccc -----------------------------------------------------------------c
110c     read name of input file                                          c
111ccccc -----------------------------------------------------------------c
112
113 4321 print*,' name of input file (without extension .inp),',
114     *            ' (max. 8 characters)'
115      print*,' default: opac.inp'
116      read (*,'(a)') file
117      call lang2 (file,8,lfile)
118      if (lfile.gt.8) then
119         print*,' lfile= ',lfile
120         print*,' file name too long! Try again!'
121         goto 4321
122      end if
123      if (file.ne.' ') then
124         infile=inp//file(1:lfile)//'.inp'
125      else
126         infile=inp//'opac.inp'
127         file='opac'
128         lfile=4
129      end if
130      write (*,*) ' input file ',infile,' used'
131
132ccccc ------------------------------------------------------------------c
133c     read data from input file                                         c
134ccccc ------------------------------------------------------------------c
135
136      call readinp
137
138ccccc ------------------------------------------------------------------c
139c     prepare data selection + check input data                         c
140c                                                                       c
141c     1. mixture of type                                                c
142ccccc ------------------------------------------------------------------c
143
144c      print*,mtyp
145      if(mtyp.eq.0) then
146         typnam=tnew
147
148         numden=0.
149         mcomp=0
150         do ic=1,5
151            if (dnum(ic).ne.0.) then
152               mcomp=mcomp+1
153               numden=numden+dnum(ic)
154            end if
155         end do
156         do ic=1,mcomp
157            ncomp(ic)=nuco(ic)
158            mixrat(ic)=dnum(ic)/numden
159            sigma(ic)=sigmac(nuco(ic))
160         end do
161      else
162         typnam=tname(mtyp)
163         mcomp=nco(mtyp)
164
165         numden=0.
166         do ic=1,mcomp
167            numden=numden+dnumb(mtyp,ic)
168         end do
169         do ic=1,mcomp
170            ncomp(ic)=mco(mtyp,ic)
171            mixrat(ic)=dnumb(mtyp,ic)/numden
172            sigma(ic)=sigmac(mco(mtyp,ic))
173         end do
174      end if
175
176c      print*,typnam,numden
177
178      if (mcomp.gt.1) then
179         do ic=1,mcomp
180            if (ncomp(ic).ge.11) then
181               print*,' '
182               print*,'    !! ATTENTION !!'
183               print*,' '
184               print*,' in the present version, clouds may not',
185     *                ' be mixed with other particles'
186               print*,' '
187               print*, ' Please change input file!'
188               stop
189            end if
190         end do
191      end if
192
193ccccc ------------------------------------------------------------------c
194c     2. relative humidity                                              c
195ccccc ------------------------------------------------------------------c
196
197      nih=0
198      do ih=1,nhc
199         if(indhum(ih).ne.0) then
200            nih=nih+1
201            khum(nih)=ih
202            ahum(nih)=nhum(ih)
203         end if
204      end do
205
206      if (nih.ne.0) then
207         do ih=1,nih
208            do ic=1,mcomp
209               dens(ic,ih)=densc(ncomp(ic),khum(ih))
210               rmin(ic,ih)=rminc(ncomp(ic),khum(ih))
211               rmax(ic,ih)=rmaxc(ncomp(ic),khum(ih))
212               rmod(ic,ih)=rmodc(ncomp(ic),khum(ih))
213c               print*,khum(ih),rmin(ic,ih),rmax(ic,ih),rmod(ic,ih)
214            end do
215         end do
216      else
217         print*,' '
218         print*,'    !! ATTENTION !!'
219         print*,' '
220         print*,' no relative humidity selected!'
221         print*,' '
222         print*, ' Please change input file!'
223         stop
224      end if
225
226ccccc ------------------------------------------------------------------c
227c     3. height profile                                                 c
228ccccc ------------------------------------------------------------------c
229
230      ht=0.
231      do il=1,5
232         ht=ht+(hmax(il)-hmin(il))
233      end do
234
235      if (mtyp.eq.0.or.ht.ne.0.) then
236         mlay=0
237         do il=1,5
238            if (hmax(il).ne.hmin(il)) then
239               mlay=mlay+1
240               nltyp(mlay)=nlt(il)
241               boundl(mlay)=hmin(il)
242               boundu(mlay)=hmax(il)
243               parlay(mlay,1)=hpar(il)
244            end if   
245         end do
246      else
247         mlay=nlay(mtyp)
248         do il=1,mlay
249            if (mlay.eq.1) then           ! for clouds
250               nltyp(1)=5
251            else
252               nltyp(il)=nlt(il)
253            end if   
254            boundl(il)=hmint(mtyp,il)
255            boundu(il)=hmaxt(mtyp,il)
256            parlay(il,1)=part(mtyp,il)
257         end do
258      end if
259
260ccccc ------------------------------------------------------------------c
261c     4. optical parameters                                             c
262ccccc ------------------------------------------------------------------c
263
264      mopar=nopt
265
266      nop=0
267      do iop=1,mopar
268         jnopar(iop)=indopt(iop)
269         if (jnopar(iop).eq.1) then
270            nop=nop+1
271            opanam(nop)=optnam(iop)
272         end if
273      end do
274      if (jnopar(15).eq.1) nop=nop-1   ! not wavelength dependent
275      if (jnopar(16).eq.1) nop=nop-1   !           "
276      if (jnopar(17).eq.1) nop=nop-1   !           "
277      if (jnopar(18).eq.1) nop=nop+1   ! refr.index has 2 values
278
279      do ic=1,mcomp
280         if (jnopar(11).eq.1.and.ncomp(ic).ge.11.or.
281     *    jnopar(12).eq.1.and.ncomp(ic).ge.11) then
282            print*,' '
283            print*,'    !! ATTENTION !!'
284            print*,' '
285            print*,' in the present version, no mass related',
286     *             ' parameters may be calculated'
287            print*,' for clouds'
288            print*,' '
289            print*, ' Please change input file!'
290            stop
291         end if
292      end do
293
294      if (jnopar(15).eq.1.and.jnopar(1).eq.0.or.
295     *    jnopar(15).eq.1.and.jnopar(4).eq.0.or.
296     *    jnopar(15).eq.1.and.jnopar(5).eq.0) then
297         print*,' '
298         print*,'    !! ATTENTION !!'
299         print*,' '
300         print*,' extinction coefficient, single scattering albedo and'
301         print*,' asymmetry parameter must be selected too'
302         print*,' for calculation of solar and terrestrial integrated'
303         print*,' parameters'
304         print*,' '
305         print*, ' Please change input file!'
306         stop
307      end if
308
309      if (jnopar(16).eq.1.and.jnopar(1).eq.0) then
310         print*,' '
311         print*,'    !! ATTENTION !!'
312         print*,' '
313         print*,' extinction coefficient must be selected, too'
314         print*,' for calculation of Angstrom coefficients'
315         print*,' '
316         print*, ' Please change input file!'
317         stop
318      end if
319
320      if (jnopar(18).eq.1.and.mcomp.gt.1) then
321         print*,' '
322         print*,'    !! ATTENTION !!'
323         print*,' '
324         print*,' refractive indices can only be printed for '
325         print*,' single components'
326         print*,' '
327         print*, ' Please change input file!'
328         stop
329      end if
330
331      if (jnopar(7).eq.1.or.jnopar(13).eq.1) then
332         print*,' '
333         print*,'    !! ATTENTION !!'
334         print*,' '
335         print*,' you selected an invalid optical parameter'
336         print*,' '
337         print*, ' Please change input file!'
338         stop
339      end if
340
341ccccc ------------------------------------------------------------------c
342c     5. wavelengths                                                    c
343ccccc ------------------------------------------------------------------c
344
345      wa=.false.
346      wb=.false.
347      do iw=1,nwavea
348         if (indwa(iw).ne.0) wa=.true.
349      end do
350      do iw=1,nwaveb
351         if (indwb(iw).ne.0) wb=.true.
352      end do
353      if (wa.and.wb) then
354         print*,' '
355         print*,'    !! ATTENTION !!'
356         print*,' '
357         print*,' you selected wavelengths in both wavelength tables!'
358         print*,' '
359         print*, ' Please change input file!'
360         stop
361      else if (wa) then
362         niw=0
363         mlamb=nwavea
364         do il=1,mlamb
365            wlamb(il)=wavea(il)
366            if (indwa(il).ne.0) then
367               niw=niw+1
368               alamb(niw)=wlamb(il)
369            end if
370         end do
371      else if (wb) then
372         niw=0
373         mlamb=nwaveb
374         do il=1,mlamb
375            wlamb(il)=waveb(il)
376            if (indwb(il).ne.0) then
377               niw=niw+1
378               alamb(niw)=wlamb(il)
379            end if
380         end do
381      end if
382      il05=0
383      il055=0
384      il035=0
385      il08=0
386      do iw=1,niw
387         if (alamb(iw).eq.0.35) il035=iw
388         if (alamb(iw).eq.0.5)   il05=iw
389         if (alamb(iw).eq.0.55) il055=iw
390         if (alamb(iw).eq.0.8)   il08=iw
391      end do
392      niwp=niw
393
394      if (jnopar(15).eq.1) then    ! bei solar/terrestr. integrierten Groessen
395         niwp=niw                  ! werden alle Wellenl. eingelesen, aber nicht
396         do il=1,niwp              ! gedruckt
397            alambp(il)=alamb(il)
398         end do
399         niw=mlamb
400         do il=1,niw
401            alamb(il)=wlamb(il)
402            if (alamb(il).eq.0.35) il035=il
403            if (alamb(il).eq.0.5)  il05=il
404            if (alamb(il).eq.0.55)  il055=il
405            if (alamb(il).eq.0.8)  il08=il
406         end do
407      end if
408
409      if (jnopar(14).eq.1.and.il055.eq.0) then
410         print*,' '
411         print*,'    !! ATTENTION !!'
412         print*,' '
413         print*,' wavelength 0.55 micron must be selected'
414         print*,' for calculation of normalized extinction coefficients'
415         print*,' '
416         print*, ' Please change input file!'
417         stop
418      end if
419
420      if (jnopar(16).eq.1.and.il035.eq.0.or.
421     *    jnopar(16).eq.1.and.il05.eq.0.or.
422     *    jnopar(16).eq.1.and.il08.eq.0) then
423         print*,' '
424         print*,'       !! ATTENTION !!'
425         print*,' '
426         print*,' wavelengths 0.35, 0.5 and 0.8 micron must be selected'
427         print*,' for calculation of Angstrom coefficients'
428         print*,' '
429         print*, ' Please change input file!'
430         stop
431      end if
432
433ccccc ------------------------------------------------------------------c
434c     Verzweigung ins Unterprogramm RAWOPT zur Berechnung der           c
435c     optischen Eigenschaften                                           c
436ccccc ------------------------------------------------------------------c
437
438      call rawopt
439
440      print*,' READY.'
441
442ccccc ------------------------------------------------------------------c
443c     Formate                                                           c
444ccccc ------------------------------------------------------------------c
445
446  100 format(i1)
447  110 format(e17.10)
448  120 format(a1)
449  130 format(e17.10)
450  150 format(i1,3x,f6.3)
451  200 format(i2)
452  300 format(a50)
453
454      stop
455      end
456      subroutine comment (com,ifile)
457
458      character*1 com,dum
459
460 1011 read (ifile,'(a1)') dum
461      if (dum.eq.com) then
462         goto 1011
463      else
464         backspace(ifile)
465      end if
466
467      return
468      end
469      subroutine readcfg
470
471ccccc ------------------------------------------------------------------c
472c     read microphysical parameters of components and types from        c
473c     file opac.cfg.                                                    c
474c                                                                       c
475c     30.09.97                                                M. Hess   c
476ccccc ------------------------------------------------------------------c
477
478      character*1  dum
479      character*4  os
480      character*6  inp
481      character*7  res
482      character*8  file
483      character*10 opt
484      character*18 infile
485      character*20 rname
486      character*30 tname
487
488      common /file/   lfile,file,rname,infile,inp,res,opt
489      common /mipaco/ sigmac(20),rminc(20,8),rmaxc(20,8),rmodc(20,8),
490     *                densc(20,8)
491      common /mipaty/ nco(20),mco(20,10),dnumb(20,10),nlay(20),
492     *                hmint(20,4),hmaxt(20,4),part(20,4),tname(20)
493
494      open (8,file='opac.cfg')
495
496      call comment('#',8)
497
498      read (8,'(a4)') os
499
500      if (os.eq.'dos ') then
501c         inp='input\'
502c         res='result\'
503c         opt='..\optdat\'
504      else if (os.eq.'unix') then
505         inp='input/'
506         res='result/'
507         opt='../optdat/'
508      else
509         stop ' wrong operating system!'
510      end if
511
512      call comment('#',8)
513
514      read (8,*) ncom
515      read (8,*) ntyp
516
517      call comment('#',8)
518
519      do icom=1,ncom
520         read(8,'(a1)') dum
521         read(8,'(a1)') dum
522         read(8,'(a1)') dum
523         do ih=1,8
524            if (icom.le.10) then
525            read (8,*) nh,rminc(icom,ih),rmaxc(icom,ih),rmodc(icom,ih),
526     *                 densc(icom,ih),sigmac(icom)
527            else
528            read (8,*) nh,rminc(icom,ih),rmaxc(icom,ih),rmodc(icom,ih),
529     *                 densc(icom,ih)
530            end if
531            read (8,'(a1)') dum
532            if (ih.eq.1.and.dum.eq.'#') then
533               do ihd=2,8
534                  rminc(icom,ihd)=rminc(icom,1)
535                  rmaxc(icom,ihd)=rmaxc(icom,1)
536                  rmodc(icom,ihd)=rmodc(icom,1)
537                  densc(icom,ihd)=densc(icom,1)
538               end do
539               goto 1111
540            else
541               backspace(8)
542            end if
543         end do
544         call comment('#',8)
545 1111    continue
546      end do
547
548      do icom=1,ncom
549      do ih=1,8
550c            print*,rminc(icom,ih),rmaxc(icom,ih),rmodc(icom,ih),
551c     *             densc(icom,ih)
552      end do
553      end do
554
555      call comment('#',8)
556
557      do ityp=1,ntyp
558         call comment('#',8)
559         read (8,200) nt,tname(ityp)
560  200    format(i2,2x,a30)
561         read (8,*) nco(ityp)
562         do ico=1,nco(ityp)
563            read (8,100) mco(ityp,ico),dnumb(ityp,ico)
564  100       format(i2,8x,e10.3)
565         end do
566         read (8,*) nlay(ityp)
567         do ilay=1,nlay(ityp)
568            read (8,*) hmint(ityp,ilay),hmaxt(ityp,ilay),part(ityp,ilay)
569         end do
570      end do
571
572      close(8)
573
574      return
575      end
576      subroutine readinp
577
578ccccc ------------------------------------------------------------------c
579c     read from input file *.inp which data to extract and calculate    c
580c     from the database.                                                c
581c                                                                       c
582c     21.09.97                                                M. Hess   c
583ccccc ------------------------------------------------------------------c
584
585      character*6  inp
586      character*7  res
587      character*8  file
588      character*10 opt
589      character*18 infile
590      character*20 rname
591      character*30 tnew
592
593      common /file/   lfile,file,rname,infile,inp,res,opt
594      common /input/  mtyp,nuco(10),dnum(10),mcom,
595     *                hmin(5),hmax(5),hpar(5),
596     *                nwavea,indwa(100),wavea(100),
597     *                nwaveb,indwb(100),waveb(100),
598     *                nhc,indhum(8),nopt,indopt(20),tnew
599
600      open (7,file=infile)
601      rname(1:7)=res
602      rname(8:(8+(lfile-1)))=file(1:lfile)
603      rname((9+(lfile-1)):(13+(lfile-2)))='.out'
604
605      call comment('#',7)
606
607ccccc ------------------------------------------------------------------c
608c     read no. of type                                                  c
609ccccc ------------------------------------------------------------------c
610
611      read (7,*) mtyp
612      call comment('#',7)
613      read (7,'(a30)') tnew
614      call comment('#',7)
615
616ccccc ------------------------------------------------------------------c
617c     read mixture of type                                              c
618ccccc ------------------------------------------------------------------c
619
620      mcom=0
621      do ic=1,5
622         read (7,*) nuco(ic),dnum(ic)
623         if (dnum(ic).gt.0.) mcom=mcom+1
624      end do
625
626      call comment('#',7)
627
628ccccc ------------------------------------------------------------------c
629c     read height profile                                               c
630ccccc ------------------------------------------------------------------c
631
632      do il=1,5
633         read (7,*) hmin(il),hmax(il),hpar(il)
634      end do
635
636      call comment('#',7)
637
638ccccc ------------------------------------------------------------------c
639c     read wavelengths                                                  c
640ccccc ------------------------------------------------------------------c
641
642      read (7,*) nwavea
643      do iw=1,nwavea
644         read (7,100) indwa(iw),wavea(iw)
645  100    format(i1,32x,f10.3)
646      end do
647
648      call comment('#',7)
649
650      read (7,*) nwaveb
651      do iw=1,nwaveb
652         read (7,100) indwb(iw),waveb(iw)
653      end do
654
655      call comment('#',7)
656
657ccccc ------------------------------------------------------------------c
658c     read relative humidity classes                                    c
659ccccc ------------------------------------------------------------------c
660
661      read (7,*) nhc
662      do ih=1,nhc
663         read (7,*) indhum(ih)
664      end do
665
666      call comment('#',7)
667
668ccccc ------------------------------------------------------------------c
669c     read selected optical parameters                                  c
670ccccc ------------------------------------------------------------------c
671
672      read (7,*) nopt
673      do iopt=1,nopt
674         read(7,*) indopt(iopt)
675      end do
676
677      close(7)
678
679      return
680      end
681      subroutine rawopt
682
683ccccc -----------------------------------------------------------------c
684c     Berechnung der optischen Parameter eines Aerosoltyps fr         c
685c     mehrere Wellenl„ngen und Feuchten                                c
686c                                                                      c
687c     Folgende Unterprogramme werden verwendet:                        c
688c                                                                      c
689c     - RAWMIC                                                         c
690c     - HEAD2                                                          c
691c     - OPTRAW                                                         c
692c     - OPTPAR                                                         c
693c     - OUT2                                                           c
694c                                                                      c
695c     Diese Version gehoert zur eingeschraenkten Datenbank OPAC und    c
696c     beruht auf RAWOPT vom 04.11.93.                                  c
697c                                                                      c
698c     13.05.94 Einlesen der Background-Extinktion.                     c
699c     12.01.96 Massen werden immer berechnet.                          c
700c     01.10.97 changed reading 'extback.dat'                           c
701c                                                                      c
702c     Stand: 01.10.97                                          M. Hess c
703ccccc -----------------------------------------------------------------c
704
705      real        mixrat,n,numden
706      integer     prnr,acnr,rht
707      character*1  dum
708      character*2  chum
709      character*3  atn,pat,typnam*30
710      character*4  comnam
711
712      common /prog/   nprog
713      common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8),
714     *                mixrat(10),dens(10,8)
715      COMMON /FTASTR/ EXTFTA(61),EXTSTR(61),extmit(61)
716      common /atyp/   natyp,mcomp,ncomp(10),numden,
717     *                typnam,comnam(20)
718      common /wavel/  mlamb,alamb(61),niw,wlamb(61),
719     *                il035,il05,il055,il08,alambp(61),niwp
720      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
721      common /mipoi/  latx,lonx,nl,prnr,rht(2),n(2),
722     *                njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2)
723      common /out/    oparam(16,2),phaf(167,2),indop(16)
724
725ccccc -----------------------------------------------------------------c
726c     Aufruf von RAWMIC fr kombinierte opt./mikrophys. Groessen       c
727ccccc -----------------------------------------------------------------c
728
729      call rawmic
730
731ccccc -----------------------------------------------------------------c
732c     Definition der Gr”áen, die in den Unterprogrammen gebraucht      c
733c     werden                                                           c
734ccccc -----------------------------------------------------------------c
735
736      n(1)=numden
737      nl=1
738      njc(1)=mcomp
739      do ic=1,mcomp
740         acnr(ic,1)=ncomp(ic)
741         acmr(ic,1)=mixrat(ic)
742      end do
743
744ccccc -----------------------------------------------------------------c
745c     Schleife ber alle verlangten Feuchteklassen                     c                                  c
746ccccc -----------------------------------------------------------------c
747
748      if (nih.eq.0) then
749         njh=1
750      else
751         njh=nih
752      end if
753
754      do ih=1,njh
755
756ccccc -----------------------------------------------------------------c
757c     Einlesen der Extinktionskoeffizienten des tropos. backgr.        c
758c     und des strat. backgr. und des Min. transported.                 c
759ccccc -----------------------------------------------------------------c
760
761c      print*,'read extback.dat' 
762      open (9,file='extback.dat')
763      IL=1
764      READ(9,'(a1)') dum
765      READ(9,'(a1)') dum     
766      DO IWL=1,mlamb
767         READ(9,*) WAVE,EXTFT,EXTST,extmi
768         do ila=1,niw
769            IF (WAVE.EQ.alamb(ila)) THEN
770               EXTFTA(IL)=EXTFT
771               EXTSTR(IL)=EXTST
772               extmit(il)=extmi
773               IL=IL+1
774            END IF
775         end do
776      end do
777
778      close (9)
779
780ccccc -----------------------------------------------------------------c
781c     Einlesen der optischen Rohdaten im Verzeichnis OPTDAT            c
782ccccc -----------------------------------------------------------------c
783
784c      print*,'start optraw' 
785      call optraw (ih)
786
787ccccc -----------------------------------------------------------------c
788c      Beschriftung des Output-Files                                   c
789ccccc -----------------------------------------------------------------c
790
791c      print*,'start head2'
792      call head2 (ih)
793
794ccccc -----------------------------------------------------------------c
795c     Schleife ber alle verlangten Wellenl„ngen                       c
796ccccc -----------------------------------------------------------------c
797
798         do il=1,niw
799
800ccccc -----------------------------------------------------------------c
801c     Auswahl der benoetigten Daten aus den eingelesenen               c
802ccccc -----------------------------------------------------------------c
803
804c      print*,'start optdat'
805            call optdat(il)
806
807ccccc -----------------------------------------------------------------c
808c     Berechnung der optischen Parameter                               c
809ccccc -----------------------------------------------------------------c
810
811c       print*,'start optpar'
812            call optpar(il,ih)
813
814         end do                ! Wellenlaengen-Schleife
815
816ccccc -----------------------------------------------------------------c
817c     Berechnung der solar und terrestr. integrierten Groessen und     c
818c     der Angstrom-Koeffizienten und der Sichtweite                    c
819ccccc -----------------------------------------------------------------c
820       
821c      print*,'start intco'
822
823         call intco
824
825      end do                   ! Feuchte-Schleife
826
827      return
828      end
829
830ccccc *****************************************************************c
831      subroutine head2 (ih)
832c     *****************************************************************c
833c                                                                      c
834c     -----------------------------------------------------------------c
835c     version for OPAC 3.0                                             c
836c                                                                      c
837c     29.09.97 new header                                              c
838c                                                                      c
839c     29.09.97                                                 M. Hess c
840ccccc -----------------------------------------------------------------c
841
842      real         mixrat,numden
843
844      character*2  chum
845      character*4  comnam
846      character*6  inp
847      character*7  res
848      character*8  file,optun,optnam,opanam
849      character*10 opt
850      character*18 infile
851      character*20 rname
852      character*30 typnam
853
854      common /file/   lfile,file,rname,infile,inp,res,opt
855      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
856      common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8),
857     *                mixrat(10),dens(10,8)
858      common /masse/  smas(10,8),smag(8),vsum(10,8),vges(8),xabmax
859      common /atyp/   natyp,mcomp,ncomp(10),numden,
860     *                typnam,comnam(20)
861      common /opar/   mopar,jnopar(20),nop,opanam(20),optnam(20),
862     *                optun(20)
863      common /out/    oparam(16,2),phaf(167,2),indop(16)
864      common /angle/  jnangle(167),angle(167),nia,ntheta
865      common /wavel/  mlamb,alamb(61),niw,wlamb(61),
866     *                il035,il05,il055,il08,alambp(61),niwp
867
868      open (10,file=rname(1:(13+(lfile-2))))
869
870CCCCC -----------------------------------------------------------------C
871C     Kopf des output-files                                            C
872CCCCC -----------------------------------------------------------------C
873
874      if (ih.eq.1) then
875         write(10,100) typnam,mcomp,nih,niwp
876         write (10,108)
877         write(10,112)
878         write (10,108)
879  100 format('#',12x,
880     *       'Optical Properties of Aerosols and Clouds (OPAC) 3.1'/
881     *       '#',11x,
882     *       '------------------------------------------------------'/
883     *        '#',1x,a30,/
884     *        '#',i3,' components'/
885     *        '#',i3,' relative humidities'/
886     *        '#',i3,' wavelengths')
887  112 format('# RH  COMP NUMBER    VOLUME    MASS      DENSITY  ',
888     *       'NUM.MIX  VOL.MIX   MAS.MIX'/
889     *       '# [%]     [1/cm^3]  [um^3/m^3] [ug/m^3]  [g/cm^3] ',
890     *       'RATIO    RATIO     RATIO')
891      do ihu=1,nih
892         do ic=1,mcomp
893           if (ic.eq.1) then
894                 write(10,124) ahum(ihu),comnam(ncomp(ic)),
895     *                         numden*mixrat(ic),
896     *                         vges(ihu)*vsum(ic,ihu),
897     *                         smag(ihu)*smas(ic,ihu),
898     *                         dens(ic,ihu),mixrat(ic),vsum(ic,ihu),
899     *                         smas(ic,ihu)
900           else
901                 write(10,121) comnam(ncomp(ic)),
902     *                         numden*mixrat(ic),
903     *                         vges(ihu)*vsum(ic,ihu),
904     *                         smag(ihu)*smas(ic,ihu),
905     *                         dens(ic,ihu),mixrat(ic),vsum(ic,ihu),
906     *                         smas(ic,ihu)
907           end if
908         end do
909      end do
910      write (10,108)
911  121 format ('#',5x,a4,1p3e10.3,1x,0pf4.2,3x,1p3e10.3)
912  124 format ('#',f4.0,1x,a4,1p3e10.3,1x,0pf4.2,3x,1p3e10.3)
913
914      write (10,110) xabmax
915  110 format ('#',1x,'only particles up to ',f5.1,' micrometer are ',
916     *               'considered for mass')
917
918      if (jnopar(10).eq.1) then
919         write (10,102) ntheta
920         write (10,103) (angle(ia),ia=1,ntheta)
921  102    format('#',i4,' scattering angles for phase functions below:')
922  103    format('#',8f10.3)
923      end if
924      write (10,107)
925  107 format('#===================================================',
926     *       '=========================')
927  108 format('#---------------------------------------------------',
928     *       '-------------------------')
929      end if
930
931CCCCC -----------------------------------------------------------------C
932C     Beschriftung fr verschiedene relative Feuchten                  c
933CCCCC -----------------------------------------------------------------C
934
935      if (ih.ne.1) then
936         write (10,108)
937      end if
938      WRITE(10,4000) ahum(ih)
939 4000 FORMAT('#',F4.0,'% relative humidity ')
940
941      return
942
943      entry names
944
945      if (jnopar(10).eq.1) then
946         kop=nop-1
947      else
948         kop=nop
949      end if
950
951      if (kop.le.7) then
952         write(10,4001) (optnam(indop(in)),in=1,kop)
953         write(10,4003) (optun(indop(in)),in=1,kop)
954 4001    format('# wavel.   ',7(1x,a8,1x))
955 4003    format('#  [um]    ',7(1x,a8,1x))
956      else
957         write(10,4001) (optnam(indop(in)),in=1,7)
958         write(10,4002) (optun(indop(in)),in=1,7)
959         write(10,4002) (optnam(indop(in)),in=8,kop)
960         write(10,4002) (optun(indop(in)),in=8,kop)
961 4002    format('#          ',7(1x,a8,1x))
962      end if
963
964      return
965      end
966
967ccccc *****************************************************************c
968      subroutine out2(ilp,il,iop)
969c     *****************************************************************c
970c                                                                      c
971c     output                                                           c
972c                                                                      c
973c     29.09.97                                                 M. Hess c
974ccccc -----------------------------------------------------------------c
975
976      character*8 optun,optnam,opanam
977
978      common /wavel/  mlamb,alamb(61),niw,wlamb(61),
979     *                il035,il05,il055,il08,alambp(61),niwp
980      common /out/    oparam(16,2),phaf(167,2),indop(16)
981      common /opar/   mopar,jnopar(20),nop,opanam(20),optnam(20),
982     *                optun(20)
983      common /angle/  jnangle(167),angle(167),nia,ntheta
984
985      if ((jnopar(10).eq.1.and.ilp.gt.1).or.ilp.eq.1) then
986         call names
987      end if
988
989      if (jnopar(10).eq.1) then
990         iop=nop-1
991      else
992         iop=nop
993      end if
994
995      if (iop.le.7) then
996         WRITE(10,4010) alamb(il),(oparam(ip,1),ip=1,iop)
997 4010    FORMAT(1p8e10.3)
998      else
999         WRITE(10,4010) alamb(il),(oparam(ip,1),ip=1,7)     ! 1. Zeile
1000         WRITE(10,4020) (oparam(ip,1),ip=8,iop)             ! 2. Zeile
1001 4020    FORMAT(10x,1p7e10.3)
1002      end if
1003
1004      if (jnopar(10).eq.1) then
1005         write(10,4002)
1006 4002    format('   phase function [1/km]')
1007         write(10,4010) (phaf(it,1),it=1,ntheta)
1008      end if
1009
1010      return
1011      end
1012      subroutine rawmic
1013
1014ccccc ------------------------------------------------------------------c
1015c     Berechnung der Gr”áenverteilung an Stuetzpunkten innerhalb der    c
1016c     Intervallgrenzen.                                                 c
1017c     Ausgabe dieser Wertepaare und der mikrophysikalischen Parameter   c
1018c     der Verteilungen aller Komponenten des Aerosoltyps.               c
1019c     Die Summe ueber die Verteilungen aller Komponenten wird ebenfalls c
1020c     berechnet.                                                        c
1021c                                                                       c
1022c                                                                       c
1023c     19.11.92 Berechnung der Volumenverteilung und des Gesamtvolumens  c
1024c     03.11.93 Masse der Komponenten und Gesamtmasse                    c
1025c              interaktives Einlesen geaenderter Radiusgrenzen          c
1026c                                                                       c
1027c     ab jetzt Version fuer OPAC.                                       c
1028c                                                                       c
1029c     08.11.93 interaktives Einlesen in dieser Version ausgeschaltet    c
1030c     04.12.95 Berechnung der Gesamtmasse nur bis zum max. Radius 7.5   c
1031c     15.01.96 Radiusfeld mit SORT1 erzeugt.                            c
1032c     17.01.96 Abscheideradius nicht groesser als Rmax                  c
1033c     22.01.96 keine Massenberechnung fuer Wolken                       c
1034c                                                                       c
1035c     Stand: 22.01.96                                           M. Hess c
1036ccccc ------------------------------------------------------------------c
1037
1038      real xr(220),xr2(220,8),dnsum(220,8),dn(220,10,8),numden,mixrat
1039      real xr3(220),dvsum(220,8),dv(220,10,8),r(220),v(220)
1040      real xr4(220),xab(10)
1041
1042      integer   nx2(8),indx(220)
1043      integer   nxs(8),nxc(10,8)
1044
1045      character*1  ent
1046      character*2  chum
1047      character*4  comnam
1048      character*6  inp
1049      character*7  res
1050      character*8  file
1051      character*10 opt
1052      character*18 infile
1053      character*20 rname
1054      character*30 typnam
1055
1056      common /file/   lfile,file,rname,infile,inp,res,opt
1057      common /prog/   nprog
1058      common /numdis/ sigma(10),rmin(10,8),rmax(10,8),rmod(10,8),
1059     *                mixrat(10),dens(10,8)
1060      common /masse/  smas(10,8),smag(8),vsum(10,8),vges(8),xabmax
1061      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
1062      common /atyp/   natyp,mcomp,ncomp(10),numden,
1063     *                typnam,comnam(20)
1064
1065      data deltar /0.015/
1066      data xrmin  /0.01/,xrmax /10./
1067      data xabmax /7.5/
1068
1069      open (10,file=rname)
1070
1071      pi=4.*atan(1.)
1072
1073ccccc ------------------------------------------------------------------c
1074c     Aendern der Radiusgrenzen bei Aufruf ueber RAWOPT                 c
1075ccccc ------------------------------------------------------------------c
1076
1077      if (nprog.eq.1) then
1078         write (*,*) ' Folgende Radiusgrenzen sind voreingestellt:'
1079         do ih=1,nih
1080            do ic=1,mcomp
1081               write (*,111) ahum(ih),ncomp(ic),rmin(ic,ih),rmax(ic,ih)
1082  111          format(' f= ',f3.0,'%',' Komponente ',i2,' rmin=',f6.3,
1083     *                ' rmax= ',f6.3)
1084            end do
1085         end do
1086
1087         write (*,112)
1088  112    format (/' sollen sie geaendert werden? (j/n) ')
1089
1090         read (*,'(a1)') ent
1091
1092         if (ent.eq.'j') then
1093            do ih=1,nih
1094               do ic=1,mcomp
1095                  write (*,113) ahum(ih),ncomp(ic)
1096  113             format(/' f= ',f3.0,'%',' Komponente ',i2,' rmin= ')
1097                  read(*,*) rmin(ic,ih)
1098                  write (*,114) ahum(ih),ncomp(ic)
1099  114             format(' f= ',f3.0,'%',' Komponente ',i2,' rmax= ')
1100                  read(*,*) rmax(ic,ih)
1101               end do
1102            end do
1103         end if
1104      end if
1105
1106ccccc ------------------------------------------------------------------c
1107c     Berechnung des Radiusgitters                                      c
1108ccccc ------------------------------------------------------------------c
1109
1110      xr(1)=xrmin
1111      ix=2
1112      xranf=alog10(xrmin)
1113      do while (xr(ix-1).lt.xrmax)
1114      xrl=xranf+deltar*(ix-1.)
1115      xr(ix)=10.**xrl
1116      ix=ix+1
1117      end do
1118      nx=ix-1
1119
1120ccccc ------------------------------------------------------------------c
1121c     Anfang Feuchteschleife                                            c
1122ccccc ------------------------------------------------------------------c
1123
1124      do ih=1,nih
1125
1126      do ix=1,nx
1127         xr2(ix,ih)=xr(ix)
1128      end do
1129
1130ccccc ------------------------------------------------------------------c
1131c     Einfuegen der Radiusgrenzen der Komponenten und des               c
1132c     Abscheideradius in das Gitter                                     c
1133ccccc ------------------------------------------------------------------c
1134
1135      kx=nx
1136      do ic=1,mcomp
1137         xab(ic)=xabmax
1138         kx=kx+1
1139         xr2(kx,ih)=rmin(ic,ih)
1140         kx=kx+1
1141         xr2(kx,ih)=rmax(ic,ih)
1142         if (rmax(ic,ih).lt.xabmax) then
1143            xab(ic)=rmax(ic,ih)
1144         end if
1145         kx=kx+1
1146         xr2(kx,ih)=xab(ic)
1147      end do
1148
1149      do ix=1,kx
1150         xr3(ix)=xr2(ix,ih)
1151      end do
1152
1153      call sort1(xr3,kx,xr4,indx)
1154
1155      xr2(1,ih)=xr4(1)
1156      lx=0
1157      do ix=2,kx
1158         if (xr4(ix).ne.xr4(ix-1)) then
1159            xr2(ix-lx,ih)=xr4(ix)
1160         else
1161            lx=lx+1
1162         end if
1163      end do
1164
1165      nx2(ih)=kx-lx
1166
1167ccccc ------------------------------------------------------------------c
1168c     Schleife ueber die Komponenten                                    c
1169ccccc ------------------------------------------------------------------c
1170
1171      do ic=1,mcomp
1172
1173         if (ncomp(ic).gt.10) goto 1101      ! keine Massenberechnung fuer Wolken
1174
1175ccccc ------------------------------------------------------------------c
1176c     Berechnung der Gr”áenverteilung an den Stuetzpunkten              c
1177ccccc ------------------------------------------------------------------c
1178
1179         do ix=1,nx2(ih)
1180
1181            if (xr2(ix,ih).ge.rmin(ic,ih).and.
1182     *          xr2(ix,ih).le.rmax(ic,ih)) then
1183
1184               dn(ix,ic,ih)=rlogn(sigma(ic),rmod(ic,ih),
1185     *                      numden*mixrat(ic),d,e,xr2(ix,ih))
1186               dv(ix,ic,ih)=vlogn(sigma(ic),rmod(ic,ih),
1187     *                      numden*mixrat(ic),pi,e,xr2(ix,ih))
1188               dv(ix,ic,ih)=dv(ix,ic,ih)*10.**(6)
1189               nxc(ic,ih)=nxc(ic,ih)+1
1190            else
1191               dn(ix,ic,ih)=0.
1192               dv(ix,ic,ih)=0.
1193            end if
1194
1195ccccc ------------------------------------------------------------------c
1196c        Berechnung der Summe ueber alle Komponenten                    c
1197ccccc ------------------------------------------------------------------c
1198
1199            dnsum(ix,ih)=dnsum(ix,ih)+dn(ix,ic,ih)
1200            dvsum(ix,ih)=dvsum(ix,ih)+dv(ix,ic,ih)
1201
1202         end do
1203
1204ccccc ------------------------------------------------------------------c
1205c     Berechnung des Aerosolvolumens und der Masse der Komponente       c
1206c     Volumen in (mikrometer**3/meter**3)                               c
1207c     Masse in (mikrogramm/m**3)                                        c
1208c                                                                       c
1209c     Volumen und Masse werden nur bis zum Radius XABMAX berechnet.     c
1210ccccc ------------------------------------------------------------------c
1211
1212         do ix=1,nx2(ih)
1213            r(ix)=xr2(ix,ih)
1214            if (r(ix).le.xab(ic)) then
1215               v(ix)=dv(ix,ic,ih)/(r(ix)*alog(10.))
1216c              v(ix)=dv(ix,ic,ih)
1217            else
1218               v(ix)=0.
1219            end if
1220         end do
1221         call gerin(r,v,erg,1,nx2(ih),ier)
1222         vsum(ic,ih)=erg
1223         smas(ic,ih)=vsum(ic,ih)*dens(ic,ih)*10.**(-6)
1224         vges(ih)=vges(ih)+vsum(ic,ih)
1225         smag(ih)=smag(ih)+smas(ic,ih)
1226
1227 1101    continue
1228      end do                                ! Komponentenschleife
1229      end do                                ! Feuchteschleife
1230
1231ccccc ------------------------------------------------------------------c
1232c     Berechnung der Volumen- und Massenmischungsverh„ltnisse           c
1233ccccc ------------------------------------------------------------------c
1234
1235      do ih=1,nih
1236         do ic=1,mcomp
1237            if (ncomp(ic).gt.10) goto 1102   ! keine Massenberechnung fuer Wolken
1238            vsum(ic,ih)=vsum(ic,ih)/vges(ih)
1239            smas(ic,ih)=smas(ic,ih)/smag(ih)
1240 1102       continue           
1241         end do
1242      end do
1243
1244ccccc -----------------------------------------------------------------c
1245c     Kopf des output-files (nicht bei Aufruf ueber RAWOPT)            c
1246ccccc -----------------------------------------------------------------c
1247
1248      if (nprog.eq.1) then
1249
1250      write(10,100) typnam,natyp,numden,mcomp,nih
1251  100 format(' 1 ; output of subroutine rawmic',10x,
1252     *       'calculation date: ',i2,'.',i2,'.',i4,
1253     *        2x,i2,':',i2,':',i2/
1254     *       ' season: ',a7,/
1255     *        1x,a30,' (',i2,') ', ' with the following microphysical',
1256     *             ' parameters:'/
1257     *        15x,'         number density: ',1pe10.3,' per cm**3',/,
1258     *        9x,'         number of components: ',i2,/,
1259     *        9x,'number of relative humidities: ',i2,/,2x,'RH    ',
1260     *        'NC     SIGMA     RMOD      RMIN      RMAX      MIXRAT',
1261     *        '    VMIXRAT')
1262      do ihu=1,nih
1263         do ic=1,mcomp
1264           if (ic.eq.1) then
1265              write(10,104) ahum(ihu),ncomp(ic),sigma(ic),rmod(ic,ihu),
1266     *                  rmin(ic,ihu),rmax(ic,ihu),
1267     *                          mixrat(ic),vsum(ic,ihu)
1268           else
1269              write(10,101) ncomp(ic),sigma(ic),rmod(ic,ihu),
1270     *                  rmin(ic,ihu),rmax(ic,ihu),
1271     *                          mixrat(ic),vsum(ic,ihu)
1272           end if
1273         end do
1274      end do
1275  101 format (i10,2x,6e10.3,i5)
1276  104 format (f5.0,i5,2x,6e10.3,i5)
1277      write (10,107) deltar
1278  107 format(' particle radii are given in micrometers, dN/dlogr in',
1279     *       ' particles/cm**3'/
1280     *       ' dV/dlogr in æm**3/m**3',5x,
1281     *       ' the radius interval is dlogr= ',f7.5/
1282     *       '====================================================',
1283     *       '============================')
1284
1285ccccc ------------------------------------------------------------------c
1286c     Ausgabe der berechneten Wertepaare der Komponenten                c
1287ccccc ------------------------------------------------------------------c
1288
1289      do ih=1,nih
1290         write (10,106) ahum(ih),vges(ih)
1291 106     format(' relative humidity: ',f3.0,'%',
1292     *          '    total volume: ',1pe10.3,' [æm**3/m**3]')
1293
1294         do ic=1,mcomp
1295            write (10,102) ncomp(ic),nxc(ic,ih)
1296  102       format('  RADIUS ','   dN/dlogr',' dV/dlogr component ',i2,
1297     *             ', ',i4,' values')
1298
1299            do ix=1,nx2(ih)
1300               if (dn(ix,ic,ih).ne.0.) then
1301                  write (10,105) xr2(ix,ih),dn(ix,ic,ih),dv(ix,ic,ih)
1302  105             format(1p3e10.3)
1303               end if
1304            end do
1305         end do
1306
1307ccccc ------------------------------------------------------------------c
1308c     Ausgabe der berechneten Wertepaare fuer die gesamte Verteilung    c
1309ccccc ------------------------------------------------------------------c
1310
1311         do ix=1,nx2(ih)
1312            if (dnsum(ix,ih).ne.0.) then
1313               nxs(ih)=nxs(ih)+1
1314            end if
1315         end do
1316         write (10,103) natyp,nxs(ih)
1317  103    format('  RADIUS ','   dN/dlogr','       aerosol type ',i2,
1318     *          ', ',i4,' values')
1319         do ix=1,nx2(ih)
1320            if (dnsum(ix,ih).ne.0.) then
1321               write (10,105) xr2(ix,ih),dnsum(ix,ih),dvsum(ix,ih)
1322            end if
1323         end do
1324
1325ccccc ------------------------------------------------------------------c
1326c     Ende der Feuchteschleife                                          c
1327ccccc ------------------------------------------------------------------c
1328
1329      end do
1330
1331      end if
1332
1333      return
1334      end
1335      FUNCTION RLOGN(SIG,RO,VN,D,E,X)
1336
1337CCCCC -----------------------------------------------------------------C
1338C     LOG-NORMAL-VERTEILUNG  DN/DLOGR IN cm**-3                        C
1339c     DIE RADIEN X UND RO MUESSEN IN um ANGEGEBEN WERDEN,              C
1340c     VN in cm**-3                                                     c
1341CCCCC -----------------------------------------------------------------C
1342
1343      PI=4.*ATAN(1.)
1344      A=VN/(SQRT(2.*PI)*ALOG10(SIG))
1345      B=ALOG10(RO)
1346      C=-2.*(ALOG10(SIG))**2
1347      RLOGN=A*EXP((ALOG10(X)-B)**2/C)
1348      RETURN
1349      END
1350      FUNCTION VLOGN(SIG,RO,VN,PI,E,X)
1351
1352CCCCC -----------------------------------------------------------------C
1353C     LOGNORMALE VOLUMENVERTEILUNG DV/DLOGR IN um**3 cm**-3            C
1354c     DIE RADIEN X UND RO MUESSEN IN um ANGEGEBEN WERDEN,              C
1355c     VN in cm**-3                                                     c
1356CCCCC -----------------------------------------------------------------C
1357
1358      A=VN/(SQRT(2.*PI)*ALOG10(SIG))
1359      B=ALOG10(RO)
1360      C=-2.*(ALOG10(SIG))**2
1361      RLOGN=A*EXP((ALOG10(X)-B)**2/C)
1362
1363      VLOGN=4./3.*PI*X**3.*RLOGN
1364
1365      RETURN
1366      END
1367
1368CCCCC *****************************************************************C
1369       subroutine optraw (ihum)
1370C     *****************************************************************C
1371C                                                                      c
1372c     neue Einleseroutine fuer Aerosoltypen (RAWOPT), nicht mehr       c
1373c     gleich der Version OPTCOM fuer ATLOPT.                           c
1374c                                                                      c
1375c     04.05.94 Ausschluss der Quellung bei Russ (Komponente 3)         c
1376c     15.05.94 Einlesen der Cirrus-Daten                               c
1377c     04.12.95 Einlesen der Brechungsindizes                           c
1378c     07.12.95 neue Namen fuer Dateien mit Komponentendaten            c
1379c                                                                      c
1380c     23.09.97 small changes for OPAC 3.0                              c
1381c     21.10.97 read new data format in ../optdat/ for OPAC 3.1         c
1382c                                                                      c
1383c     Stand: 21.10.97                                         M. Hess  c
1384CCCCC -----------------------------------------------------------------C
1385
1386      logical ende
1387      integer prnr,acnr,njc,rht
1388
1389      real n,numden
1390      real ex(61,6),sc(61,6),ab(61,6),si(61,6),as(61,6),ba(61,6)
1391      real ph(61,6,167),br(61,6),bi(61,6)
1392
1393      character*1  dum
1394      character*2  chum
1395      character*3  atn,pat
1396      character*4  comnam
1397      character*6  inp
1398      character*7  res
1399      character*8  file
1400      character*10 opt,dum2
1401      character*16 tap
1402      character*18 infile
1403      character*20 rname
1404      character*30 typnam
1405
1406      common /file/   lfile,file,rname,infile,inp,res,opt
1407      common /atyp/   natyp,mcomp,ncomp(10),numden,
1408     *                typnam,comnam(20)
1409      common /wavel/  mlamb,alamb(61),niw,wlamb(61),
1410     *                il035,il05,il055,il08,alambp(61),niwp
1411      common /hum/    khum(8),ahum(8),nih,nhum(8),mhum,chum(8)
1412      common /mipoi/  latx,lonx,nl,prnr,rht(2),n(2),
1413     *                njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2)
1414      common /oppoi/  ext(1,6),sca(1,6),abs(1,6),sis(1,6),asy(1,6),
1415     *                bac(1,6),pha(1,6,167),ext05,bre(1,6),bim(1,6)
1416      common /angle/  jnangle(167),angle(167),nia,ntheta
1417     
1418      save
1419
1420ccccc -----------------------------------------------------------------c
1421c     Schleife ber alle vorkommenden Komponenten                      c
1422ccccc -----------------------------------------------------------------c
1423
1424      do il=1,nl
1425         if (nih.eq.0) then
1426            do ihu=1,mhum
1427               if (nh(il).eq.ihu) then
1428                  khum(ihum)=ihu
1429               end if
1430            end do
1431         end if
1432
1433         ext05=0.
1434         do ic=1,njc(il)
1435            jc=acnr(ic,il)
1436
1437ccccc -----------------------------------------------------------------c
1438c     Ausschluá der Quellung bei insoluble, Russ und den               c
1439c     mineralischen Komponenten und bei den Wolken                     c                          c
1440ccccc -----------------------------------------------------------------c
1441
1442            if ( jc.eq.1.or.jc.eq.3.or.(jc.ge.6.and.jc.le.9).or.
1443     *           jc.gt.10 ) then
1444               iht=1
1445            else
1446               iht=khum(ihum)
1447            end if
1448
1449ccccc -----------------------------------------------------------------c
1450c     Bestimmung des Filenamens der gesuchten Komponente aus           c
1451c     Komponentennummer und Feuchteklasse                              c
1452ccccc -----------------------------------------------------------------c
1453
1454            tap(1:10)=opt
1455            tap(11:14)=comnam(jc)
1456            tap(15:16)=chum(iht)
1457c            print*,' tap= ',tap,jc,khum(ihum)
1458            ntap=70
1459
1460            open (ntap,file=tap,iostat=ios)
1461c            if (ios.ne.0) then
1462c               print*,' error while opening file ',tap
1463c               print*,'ios=',ios
1464c               stop
1465c            end if
1466
1467            do iline=1,100
1468               read (ntap,220) dum2
1469               if (dum2.eq.'# optical ') then
1470                  goto 2002
1471               end if
1472            end do
1473 2002       continue
1474            do iline=1,5
1475               read (ntap,200) dum
1476            end do                 
1477               
1478               do ilam=1,mlamb
1479                  read (ntap,500) rlamb,extco,scaco,absco,sisca,asymf,
1480     *                            exn,refr,refi
1481  500             format(2x,7e10.3,2e11.3)
1482                  ex(ilam,ic)=extco
1483                  sc(ilam,ic)=scaco       
1484                  ab(ilam,ic)=absco
1485                  if (rlamb.eq.0.55) then
1486                     ext05=ext05+ex(ilam,ic)*acmr(ic,1)*numden
1487                  end if
1488                  si(ilam,ic)=sisca
1489                  as(ilam,ic)=asymf
1490                  br(ilam,ic)=refr
1491                  bi(ilam,ic)=refi
1492                  if (rlamb.ne.wlamb(ilam)) then
1493                  print*,' '
1494                  print*,'        !! ATTENTION !!'
1495                  print*,' '
1496                  print*,' you have selected the wrong wavelength table'
1497                  print*,' Please change selection!'
1498                  print*,' '
1499                  print*,' Press <ENTER> to return to Input Menu'
1500                  read (*,'(a1)') dum
1501                  stop
1502                  end if
1503               end do
1504               read (ntap,'(7(/))')
1505               it=1
1506               ende=.false.
1507               do while (.not.ende)
1508                  read (ntap,510,end=511)
1509     *                  angle(it),(ph(ilam,ic,it),ilam=1,mlamb)
1510  510             format(e11.3,1x,70e10.3)
1511                  it=it+1
1512               end do
1513  511          ntheta=it-1
1514               do ilam=1,mlamb
1515                  do it=1,ntheta
1516                     ph(ilam,ic,it)=ph(ilam,ic,it)
1517                  end do
1518                  ba(ilam,ic)=ph(ilam,ic,ntheta)
1519               end do
1520               close (ntap)
1521
1522         end do
1523      end do
1524
1525ccccc -----------------------------------------------------------------c
1526c     Formate                                                          c
1527ccccc -----------------------------------------------------------------c
1528
1529  100 format(8e10.3)
1530  200 format(a1)
1531  220 format(a10)
1532  300 format(15x,f6.3,/,8x,e10.4,7x,e10.4,7x,e10.4,5x,f7.4,7x,f6.4,/)
1533  301 format(53x,f6.3,8x,e10.3,/,8x,f6.3//)
1534  302 format(8x,f6.3,8x,e10.3,9x,f6.3//)
1535  304 format(28x,f7.3//)
1536  303 format(7x,e10.3,7x,e10.3,7x,e10.3,7x,f7.4/)
1537  305 format(6x,a4)
1538  400 format(12(/))
1539  404 format(21(/))
1540 1010 format(70X,e10.3)
1541
1542      return
1543
1544      entry optdat (ilamb)
1545
1546ccccc -----------------------------------------------------------------c
1547c     Auswahl der gewuenschten optischen Daten aus den eingelesenen    c
1548ccccc -----------------------------------------------------------------c
1549
1550      do ic=1,njc(1)
1551         do il=1,mlamb
1552            if (alamb(ilamb).eq.wlamb(il)) then
1553               ext(1,ic)=ex(il,ic)
1554               sca(1,ic)=sc(il,ic)
1555               abs(1,ic)=ab(il,ic)
1556               sis(1,ic)=si(il,ic)
1557               asy(1,ic)=as(il,ic)
1558               bac(1,ic)=ba(il,ic)
1559               bre(1,ic)=br(il,ic)
1560               bim(1,ic)=bi(il,ic)
1561               do it=1,ntheta
1562                  pha(1,ic,it)=ph(il,ic,it)
1563               end do
1564            end if
1565         end do
1566      end do
1567
1568      return
1569      end
1570
1571CCCCC *****************************************************************C
1572      SUBROUTINE OPTPAR (ilamb,ihum)
1573C     *****************************************************************C
1574C                                                                      c
1575C     -----------------------------------------------------------------C
1576C     Berechnung und Ausdruck der gewnschten optischen Parameter      C
1577C                                                                      c
1578C     04.11.93 Parameter SCARA, ABSRA, OMERA eingefuehrt.              c
1579C     09.11.93 Version fuer OPAC (Aufruf OUT4 entfernt)                c
1580C     13.05.94 optische Dicke fuer alle Wellenlaengen                  c
1581C     14.05.94 Indizierung der opt. Param. fuer Ausgabe korrigiert.    c
1582C     28.07.94 normierter Ext.koeff. eingebaut                         c
1583C     04.12.95 Brechungsindizes eingebaut                              c
1584C     02.01.96 Berechnung der opt. Dicke ueber Schichttypen.           c
1585C     10.01.96 neue Berechnung der opt. Dicke (ohne Heff)              c
1586C     22.01.96 EXTRA statt SCARA                                       c
1587C                                                                      c
1588C     26.09.97 normalized extinction corrected.                        c
1589C                                                                      c
1590C     Stand: 26.09.97                                          M. Hess c
1591CCCCC -----------------------------------------------------------------C
1592
1593      integer prnr,acnr,rht
1594
1595      real n,numden
1596      REAL EXTN(2),ABSN(2),SCAN(2),PF18N(2),supf(167),phafu(167,2)
1597      REAL EXTA(2),ABSA(2),SCAA(2),SSA(2),ASF(2),PF18A(2)
1598      real scar(2),absr(2),omer(2)
1599
1600      character*3  atn,pat
1601      character*8  optnam,optun,opanam
1602      character*4  comnam
1603      character*30 typnam
1604
1605      common /atyp/   natyp,mcomp,ncomp(10),numden,
1606     *                typnam,comnam(20)
1607      common /layer/  mlay,nltyp(10),parlay(10,5),boundl(10),
1608     *                boundu(10)
1609      common /norm/   norm,mixnor
1610      common /opar/   mopar,jnopar(20),nop,opanam(20),optnam(20),
1611     *                optun(20)
1612      common /wavel/  mlamb,alamb(61),niw,wlamb(61),
1613     *                il035,il05,il055,il08,alambp(61),niwp
1614      common /mipoi/  latx,lonx,nl,prnr,rht(2),n(2),
1615     *                njc(2),acnr(5,2),acmr(5,2),nh(2),atn(2),pat(2)
1616      common /oppoi/  ext(1,6),sca(1,6),abs(1,6),sis(1,6),asy(1,6),
1617     *                bac(1,6),pha(1,6,167),ext05,bre(1,6),bim(1,6)
1618      COMMON /FTASTR/ EXTFTA(61),EXTSTR(61),extmit(61)
1619      common /out/    oparam(16,2),phaf(167,2),indop(16)
1620      common /sotin/  exi(61),ssi(61),asi(61)
1621      common /prog/   nprog
1622      common /angle/  jnangle(167),angle(167),nia,ntheta
1623      common /masse/  smas(10,8),smag(8),vsum(10,8),vges(8),xabmax
1624
1625CCCCC ------------------------------------------------------------------C
1626C     MISCHEN DES AEROSOL-TYPS                                          C
1627C     SUMM(E,A,S) : SUMME EXTINCTION, ABSORPTION, SCATTERING            C
1628C     SUPF18      : SUMME DES RUECKSTREUKOEFFIZIENTEN                   C
1629C     SUMASF      : ZWISCHENSUMME DES ASYMMETRIEFAKTORS (ASF)           C
1630C     SUMASF      : ZWISCHENSUMME DER SINGLE SCAT. ALB. (SSA)           C
1631CCCCC ------------------------------------------------------------------C
1632
1633      DO 10 L=1,NL
1634
1635      SUMME = 0.
1636      SUMMA = 0.
1637      SUMMS = 0.
1638      SUMSSA = 0.
1639      SUMASF = 0.
1640      SUPF18 = 0.
1641c     if (jnopar(10).eq.1) then
1642      do it=1,ntheta
1643         supf(it)=0.
1644      end do
1645c     end if
1646
1647      DO 20 JC=1,NJC(L)
1648
1649      SUMME = SUMME + ACMR(JC,L)*EXT(l,jc)
1650      SUMMA = SUMMA + ACMR(JC,L)*ABS(l,jc)
1651      SUMMS = SUMMS + ACMR(JC,L)*SCA(l,jc)
1652      SUMSSA = SUMSSA + ACMR(JC,L)*sis(l,jc)
1653     +                  *EXT(l,jc)
1654      SUMASF = SUMASF + ACMR(JC,L)*asy(l,jc)
1655     +                  *SCA(l,jc)
1656      SUPF18 = SUPF18 + ACMR(JC,L)*bac(l,jc)
1657c     if (jnopar(10).eq.1) then
1658         do it=1,ntheta
1659            supf(it)=supf(it)+acmr(jc,l)*pha(l,jc,it)
1660         end do
1661c     end if
1662   20 CONTINUE
1663
1664CCCCC -----------------------------------------------------------------C
1665C     Normierte  optische Parameter                                    c
1666CCCCC -----------------------------------------------------------------C
1667
1668      EXTN(L) = SUMME
1669      ABSN(L) = SUMMA
1670      SCAN(L) = SUMMS
1671      PF18N(L) = SUPF18
1672      if (jnopar(10).eq.1) then
1673      do it=1,ntheta
1674         phafu(it,l)=supf(it)
1675      end do
1676      end if
1677
1678      SSA(L) = SUMSSA/SUMME
1679      ASF(L) = SUMASF/SUMMS
1680
1681CCCCC -----------------------------------------------------------------C
1682C     ABSOLUTE OPTISCHE PARAMETER                                      C
1683CCCCC -----------------------------------------------------------------C
1684
1685      EXTA(L)= EXTN(L) * N(L)
1686      ABSA(L)= ABSN(L) * N(L)
1687      SCAA(L)= SCAN(L) * N(L)
1688      PF18A(L) = PF18N(L)* N(L)
1689      if (jnopar(10).eq.1.and.norm.eq.1) then
1690         do it=1,ntheta
1691            phafu(it,l)=phafu(it,l)*n(l)
1692         end do
1693      end if
1694      if (norm.eq.1) then
1695         EXTN(L)= EXTA(L)
1696         ABSN(L)= ABSA(L)
1697         SCAN(L)= SCAA(L)
1698         PF18N(L) = PF18A(L)
1699      end if
1700
1701c     if (jnopar(15).eq.1.or.jnopar(16).eq.1) then
1702         exi(ilamb)=extn(1)
1703         ssi(ilamb)=ssa(1)
1704         asi(ilamb)=asf(1)
1705c     end if
1706
1707c     if (jnopar(10).eq.1) then
1708         itp=1
1709         do it=1,ntheta
1710c           if (jnangle(it).eq.1) then
1711               phaf(itp,l)=phafu(it,l)
1712               itp=itp+1
1713c           end if
1714         end do
1715c     end if
1716
1717      if (jnopar(11).eq.1) then
1718         scar(l)=exta(l)/smag(ihum)*1000.  ! Einheit m**2/g
1719      end if
1720
1721      if (jnopar(12).eq.1) then
1722         absr(l)=absa(l)/smag(ihum)*1000.
1723      end if
1724
1725c     if (jnopar(13).eq.1) then
1726         kc=0
1727         do jc=1,njc(l)
1728            if (ncomp(jc).eq.3) kc=jc
1729         end do
1730         if (kc.ne.0) then
1731            omer(l)=smas(kc,ihum)/ssa(l)
1732         else
1733            omer(l)=99.
1734         end if
1735c     end if
1736
1737CCCCC -----------------------------------------------------------------C
1738C     AUSGABE DER DATEN                                                c
1739CCCCC -----------------------------------------------------------------C
1740
1741      iop=0
1742      kop=0
1743      if (jnopar(1).eq.1) then
1744         iop=iop+1
1745         indop(iop)=1
1746         oparam(iop,l)=extn(l)
1747      end if
1748      if (jnopar(2).eq.1) then
1749         iop=iop+1
1750         indop(iop)=2
1751         oparam(iop,l)=scan(l)
1752      end if
1753      if (jnopar(3).eq.1) then
1754         iop=iop+1
1755         indop(iop)=3
1756         oparam(iop,l)=absn(l)
1757      end if
1758      if (jnopar(4).eq.1) then
1759         iop=iop+1
1760         indop(iop)=4
1761         oparam(iop,l)=ssa(l)
1762      end if
1763      if (jnopar(5).eq.1) then
1764         iop=iop+1
1765         indop(iop)=5
1766         oparam(iop,l)=asf(l)
1767      end if
1768      if (jnopar(9).eq.1) then
1769         iop=iop+1
1770         indop(iop)=9
1771         oparam(iop,l)=exta(l)/pf18a(l)
1772      end if
1773      if (jnopar(11).eq.1) then
1774         iop=iop+1
1775         indop(iop)=11
1776         oparam(iop,l)=scar(l)
1777      end if
1778      if (jnopar(12).eq.1) then
1779         iop=iop+1
1780         indop(iop)=12
1781         oparam(iop,l)=absr(l)
1782      end if
1783      if (jnopar(13).eq.1) then
1784         iop=iop+1
1785         indop(iop)=13
1786         oparam(iop,l)=omer(l)
1787      end if
1788      if (jnopar(14).eq.1) then
1789         iop=iop+1
1790         indop(iop)=14
1791         oparam(iop,l)=extn(l)/ext05
1792         if (norm.eq.0) oparam(iop,l)=oparam(iop,l)*n(l)
1793      end if
1794      if (jnopar(18).eq.1) then         ! Brechungsindizes
1795         iop=iop+1
1796         indop(iop)=18
1797         oparam(iop,l)=bre(1,1)
1798         iop=iop+1
1799         indop(iop)=19
1800         oparam(iop,l)=bim(1,1)
1801      end if
1802   10 CONTINUE
1803
1804CCCCC -----------------------------------------------------------------C
1805C     OPTISCHE DICKE (nur bei absoluten Werten)                        c
1806CCCCC -----------------------------------------------------------------C
1807
1808      if (jnopar(6).eq.1.or.jnopar(7).eq.1.or.jnopar(8).eq.1) then
1809      if (norm.eq.1) then             ! nur fuer Absolutwerte
1810
1811CCCCC -----------------------------------------------------------------C
1812c     Bestimmung von HM, HFTA, HSTR, EXTFTA, EXTSTR aus den            c
1813c     eingelesenen Werten in /layer/ fuer RAWOPT                       c
1814CCCCC -----------------------------------------------------------------C
1815
1816         odepth=0.
1817         do il=1,mlay
1818            if (nltyp(il).eq.1) then                    ! mixing layer
1819               heff = parlay(il,1) *
1820     *                  (exp(-boundl(il)/parlay(il,1))-
1821     *                   exp(-boundu(il)/parlay(il,1)) )
1822               odepth = odepth + extn(1) * heff
1823c          print*,'OD(1)= ',odepth,' Heff= ',heff,' ext= ',extn(1)
1824            else if (nltyp(il).eq.2) then               ! mineral transported
1825               heff=(boundu(il)-boundl(il))
1826               odepth = odepth + extmit(ilamb) * parlay(il,1) * heff
1827c          print*,'OD(2)= ',odepth,' Heff= ',heff,' ext= ',extmit(ilamb)
1828            else if (nltyp(il).eq.3) then               ! free troposphere
1829               heff = parlay(il,1) *
1830     *                  (exp(-boundl(il)/parlay(il,1))-
1831     *                   exp(-boundu(il)/parlay(il,1)) )
1832               odepth = odepth + extfta(ilamb) * heff
1833c          print*,'OD(3)= ',odepth,' Heff= ',heff,' ext= ',extfta(ilamb)
1834            else if (nltyp(il).eq.4) then               ! stratosphere
1835               heff= (boundu(il)-boundl(il))
1836               odepth = odepth + extstr(ilamb) * heff
1837c          print*,'OD(4)= ',odepth,' heff= ',heff,' ext= ',extstr(ilamb)
1838            else if (nltyp(il).eq.5) then               ! cloud
1839               heff= (boundu(il)-boundl(il))
1840               odepth = odepth + extn(1) * heff
1841c          print*,'OD(5)= ',odepth,' heff= ',heff,' ext=',extn(1)
1842            end if
1843         end do
1844
1845         odeptha=odepth/alog(10.)
1846
1847         turbr=0.008569*alamb(ilamb)**(-4)*(1.+0.0113*alamb(ilamb)**
1848     *            (-2)+0.00013*alamb(ilamb)**(-4))
1849
1850         turbf=(odepth+turbr)/turbr
1851
1852         if (jnopar(6).eq.1) then
1853            iop=iop+1
1854            indop(iop)=6
1855            oparam(iop,1)=odepth
1856         end if
1857
1858         if (jnopar(7).eq.1) then
1859            iop=iop+1
1860            indop(iop)=7
1861            oparam(iop,1)=odeptha
1862         end if
1863
1864         if (jnopar(8).eq.1) then
1865            iop=iop+1
1866            indop(iop)=8
1867            oparam(iop,1)=turbf
1868         end if
1869
1870      else
1871         print*,' '
1872         print*, '          !! ATTENTION !!'
1873         print*,' '
1874         print*, ' optical depths may not be calculated together'
1875         print*, ' with normalized coefficients!'
1876         print*, ' Please change setting to absolute! '
1877         print*,' '
1878         print*, ' Press <ENTER> to return to Input Menu'
1879         read (*,'(a1)') dum
1880         stop
1881      end if
1882      end if
1883
1884      if (jnopar(15).eq.1) then
1885         do ilp=1,niwp
1886            if (alambp(ilp).eq.alamb(ilamb)) call out2(ilp,ilamb,iop)
1887         end do
1888      else
1889         call out2(ilamb,ilamb,iop)
1890      end if
1891
1892      RETURN
1893      END
1894      subroutine intco
1895
1896ccccc ------------------------------------------------------------------c
1897c     Integration von berechneten optischen Groessen mit Wichtung       c
1898c     mit der Solarkonstanten im solaren Spektralbereich und Wichtung   c
1899c     mit Planck(300 K) im terrestrischen.                              c
1900c                                                                       c
1901c     uebernommen von SOLINT.FOR vom 19.04.94                           c
1902c                                                                       c
1903c     ergaenzt um Berechnung der Angstrom-Koeffizienten bei 0.35/0.5    c
1904c     und 0.5/0.8 Mikrometer.                                           c
1905c                                                                       c
1906c     25.01.95                                                          c
1907c                                                                       c
1908c     21.01.96 Berechnung der Sichtweite                                c
1909c                                                                       c
1910c     Stand: 21.01.96                                           M. Hess c
1911ccccc ------------------------------------------------------------------c
1912
1913      real welsk(87),exis(87),ssis(87),asis(87)
1914      real sgrenz(100)
1915      real sol(100),terr(100),wter(100),fluss(100)
1916      real oint1(100),oint2(100),gint1(100),gint2(100),aint1(100)
1917      real atint1(100),atint2(100),aint2(100),eint1(100)
1918
1919      character*8 optnam,optun,opanam
1920
1921      common /opar/   mopar,jnopar(20),nop,opanam(20),optnam(20),
1922     *                optun(20)
1923      common /wavel/  mlamb,alamb(61),niw,wlamb(61),
1924     *                il035,il05,il055,il08,alambp(61),niwp
1925      common /oppoi/  ext(1,6),sca(1,6),abs(1,6),sis(1,6),asy(1,6),
1926     *                bac(1,6),pha(1,6,167),ext05,bre(1,6),bim(1,6)
1927      common /out/    oparam(16,2),phaf(167,2),indop(16)
1928      common /sotin/  exi(61),ssi(61),asi(61)
1929
1930      if(jnopar(15).eq.0.and.jnopar(16).eq.0.and.jnopar(17).eq.0) return
1931
1932ccccc ------------------------------------------------------------------c
1933c     solar und terrestr. integrierte Groessen                          c
1934ccccc ------------------------------------------------------------------c
1935
1936      open (5,file='wel.dat')
1937      read (5,*) nwelsk
1938      do iw=1,nwelsk
1939         read(5,*) welsk(iw)
1940      end do
1941      close(5)
1942
1943      if (jnopar(15).eq.1) then
1944         call trapez(alamb,exi,1,mlamb,welsk,exis,11,nwelsk-5,ier)
1945         call trapez(alamb,ssi,1,mlamb,welsk,ssis,11,nwelsk-5,ier)
1946         call trapez(alamb,asi,1,mlamb,welsk,asis,11,nwelsk-5,ier)
1947
1948         do il=1,nwelsk
1949            if (welsk(il).eq.0.3) i03=il
1950            if (welsk(il).eq.3.28) i3=il
1951            if (welsk(il).eq.7.992) i8=il
1952            if (welsk(il).eq.15.331) i15=il
1953         end do
1954
1955ccccc ------------------------------------------------------------------c
1956c        Einlesen der Solarkonstanten                                   c
1957ccccc ------------------------------------------------------------------c
1958
1959         open (5,file='solar.dat')
1960
1961         read (5,'(a1)') dum
1962         read (5,*) (sgrenz(ig),ig=1,38)
1963         read (5,'(a1)') dum
1964         read (5,*) (fluss(iv),iv=1,37)
1965
1966         do il=1,37
1967            sol(il)=fluss(il)/(sgrenz(il+1)-sgrenz(il))/4.
1968         end do
1969
1970         close (5)
1971
1972ccccc ------------------------------------------------------------------c
1973c        Einlesen der terrestrischen Strahlung fuer T=300 K             c
1974ccccc ------------------------------------------------------------------c
1975
1976         open (5,file='terr.dat')
1977
1978         do il=1,20
1979            read (5,*) wter(il),terr(il)
1980         end do
1981
1982         close (5)
1983
1984ccccc ------------------------------------------------------------------c
1985c        Integration                                                    c
1986ccccc ------------------------------------------------------------------c
1987
1988         do il=1,37
1989            oint1(il)=ssis(il)*exis(il)*sol(il)
1990            oint2(il)=exis(il)*sol(il)
1991            gint1(il)=asis(il)*exis(il)*ssis(il)*sol(il)
1992            gint2(il)=exis(il)*ssis(il)*sol(il)
1993            aint1(il)=exis(il)*(1.-ssis(il))*sol(il)
1994            aint2(il)=sol(il)
1995            eint1(il)=exis(il)*sol(il)
1996         end do
1997
1998         do il=1,20
1999            atint1(il)=exis(il+i8-1)*(1.-ssis(il+i8-1))*terr(il)
2000            atint2(il)=terr(il)
2001         end do
2002
2003         call gerin (welsk,oint1,o1erg,11,37,ier)
2004         call gerin (welsk,oint2,o2erg,11,37,ier)
2005         call gerin (welsk,gint1,g1erg,11,37,ier)
2006         call gerin (welsk,gint2,g2erg,11,37,ier)
2007         call gerin (welsk,aint1,a1erg,11,37,ier)
2008         call gerin (welsk,aint2,a2erg,11,37,ier)
2009         call gerin (welsk,eint1,eerg,11,37,ier)
2010
2011         call gerin (wter,atint1,at1erg,1,20,ier)
2012         call gerin (wter,atint2,at2erg,1,20,ier)
2013
2014         oint=o1erg/o2erg
2015         gint=g1erg/g2erg
2016         aint=a1erg/a2erg
2017         eint=eerg/a2erg
2018
2019         atint=at1erg/at2erg
2020      end if
2021
2022ccccc ------------------------------------------------------------------c
2023c     Angstrom-Koeffizienten                                            c
2024ccccc ------------------------------------------------------------------c
2025
2026      if (jnopar(16).eq.1) then
2027         al1=(alog10(exi(il05))-alog10(exi(il035)))/
2028     *       (alog10(alamb(il035))-alog10(alamb(il05)))
2029         al2=(alog10(exi(il08))-alog10(exi(il05)))/
2030     *       (alog10(alamb(il05))-alog10(alamb(il08)))
2031         be1=exi(il035)*alamb(il035)**(al1)
2032         be2=exi(il05)*alamb(il05)**(al2)
2033      end if
2034
2035ccccc ------------------------------------------------------------------c
2036c     Sichtweite                                                        c
2037ccccc ------------------------------------------------------------------c
2038
2039      if (jnopar(17).eq.1) then
2040         visib=3.0/(ext05+0.01159)
2041      end if
2042
2043ccccc ------------------------------------------------------------------c
2044c     Ausgabe in Output-File                                            c
2045ccccc ------------------------------------------------------------------c
2046
2047      write (10,107)
2048  107 format(' ---------------------------------------------------',
2049     *       '-------------------------')
2050      if (jnopar(15).eq.1) then
2051         write (10,120) eint,oint,gint,aint,atint
2052  120    format(' values weighted with solar spectrum',
2053     *          ' (0.3-3.3 micron):'/
2054     *          e10.3,' : extinction coefficient'/
2055     *          e10.3,' : single scattering albedo'/
2056     *          e10.3,' : asymmetry parameter'/
2057     *          e10.3,' : absorption coefficient'/
2058     *          ' values weighted with terrestrial spectrum at 300K',
2059     *          ' (8-15 micron):'/
2060     *          e10.3,' : absorption coefficient')
2061      end if
2062      if (jnopar(16).eq.1) then
2063         write (10,130) al1,be1,al2,be2
2064  130    format(' Angstrom coefficients calculated at ',
2065     *          '350 nm and 500 nm:'/
2066     *          e10.3,' : alpha'/
2067     *          e10.3,' : beta'/
2068     *          ' Angstrom coefficients calculated at ',
2069     *          '500 nm and 800 nm:'/
2070     *          e10.3,' : alpha'/
2071     *          e10.3,' : beta')
2072      end if
2073      if (jnopar(17).eq.1) then
2074         write (10,110) visib
2075  110    format(' Visibility: ',f6.2,' [km]')
2076      end if
2077
2078      return
2079      end
2080      SUBROUTINE LANG2(STEXT,mtext,LTEXT)
2081
2082CCCCC -----------------------------------------------------------------C
2083C     LAENGE DES TEXTES WIRD BERECHNET (ENDE SOBALD 2 BLANKS           C
2084C     HINTEREINANDER AUFTAUCHEN)                                       C
2085c                                                                      c
2086c     19.06.94 maximale Laenge mtext wird als Input mit uebergeben     c
2087c     03.03.95 Fehler korrigiert, falls ltext=mtext-1                  c
2088c     16.04.95 ganz neue Version                                       c
2089c                                                                      c
2090c     16.04.95                                               M. Hess   c
2091CCCCC -----------------------------------------------------------------C
2092
2093      CHARACTER*(*) STEXT
2094
2095      if (stext(mtext:mtext).ne.' ') then
2096         ltext=mtext
2097      else if (stext(mtext-1:mtext-1).eq.' ') then
2098         DO J=mtext,2,-1
2099            IF(STEXT(J:J).EQ.' '.AND.STEXT(J-1:J-1).EQ.' ') THEN
2100               LTEXT=J-2
2101            END IF
2102         end do
2103      else
2104         LTEXT=mtext-1
2105      end if
2106
2107      RETURN
2108      END
2109ccccc *****************************************************************c
2110      subroutine sort1(werte,nx,wsort,index)
2111c     *****************************************************************c
2112c                                                                      c
2113c     -----------------------------------------------------------------c
2114c     Sortieren eines 1-dimensionalen Feldes unter Verwendung von      c
2115c     MINMAX1.                                                         c
2116c                                                                      c
2117c     wsort : nach der Groesse sortiertes Feld                         c
2118c     index : zugehoeriges Feld der urspruenglichen Indizes            c
2119c                                                                      c
2120c     nx darf maximal 1000 sein!                                       c
2121c                                                                      c
2122c     Stand: 17.09.93                                          M. Hess c
2123ccccc -----------------------------------------------------------------c
2124
2125      integer index(nx)
2126      real werte(nx),wsort(nx),whilf(1000)
2127      logical ende
2128
2129      do ix=1,nx
2130         whilf(ix)=werte(ix)
2131      end do
2132
2133      do ix=1,nx-1
2134         call minmax1(whilf,nx,wmin,wmax)
2135         if (ix.eq.1) then                 ! Bestimmung des groessten Wertes
2136            wsort(nx)=wmax
2137            ende=.false.
2138            i=0
2139            do while (.not.ende)
2140               i=i+1
2141               if (whilf(i).eq.wmax) then
2142                  index(nx)=i
2143                  ende=.true.
2144               end if
2145            end do
2146         end if
2147         wsort(ix)=wmin                    ! Sortieren des restlichen Feldes
2148         ende=.false.
2149         i=0
2150         do while (.not.ende)
2151            i=i+1
2152            if (whilf(i).eq.wmin) then
2153               whilf(i)=wmax
2154               index(ix)=i
2155               ende=.true.
2156            end if
2157         end do
2158      end do
2159
2160      return
2161      end
2162ccccc *****************************************************************c
2163      subroutine minmax1(werte,nx,wmin,wmax)
2164c     *****************************************************************c
2165c                                                                      c
2166c     -----------------------------------------------------------------c
2167c     Bestimmung des Maximums und des Minimums aus einem               c
2168c     1-dimensionalen REAL-Feld                                        c
2169c                                                                      c
2170c     Stand: 14.02.92                                          M. Hess c
2171ccccc -----------------------------------------------------------------c
2172
2173      real werte(nx)
2174
2175      wmin=werte(1)
2176      wmax=werte(1)
2177
2178      do ix=1,nx
2179         if (werte(ix).ge.wmax) wmax=werte(ix)
2180         if (werte(ix).le.wmin) wmin=werte(ix)
2181      end do
2182
2183      return
2184      end
2185      SUBROUTINE TRAPEZ(X,Y,IANFA,IENDA,U,V,IANFN,IENDN,IERROR)
2186C
2187C     STAND VOM 15. JULI 1975
2188C     LINEARE INTERPOLATION EINES ORDINATENFELDES V(I), FUER EIN
2189C     ABSZISSENFELD U(I) FUER I=IANFN,IENDN AUS EINEM AUSGANGSFELD
2190C     (REFERENZFELD) Y(I),X(I), FUER I=IANFA,IENDA
2191C      U-FELD UND X-FELD  MUESSEN  GEORDNET SEIN
2192C
2193C     IERROR=0 IM NORMALFALL
2194C     IERROR=2 WENN DAS AUSGANGSFELD KLEINER IST ALS DAS ZU INTERPOLIERENDE
2195C              FELD
2196C
2197      DIMENSION X(IENDA),Y(IENDA),U(IENDN),V(IENDN)
2198      IERROR=0
2199C
2200      IDREHX=0
2201      IF (X(IANFA).LT.X(IENDA)) GOTO 2
2202C     DAS X-FELD WIRD ANSTEIGEND GEORDNET, DAS Y-FELD ENTSPRECHEND
2203      IORDX=(IENDA-IANFA+1)/2
2204      DO 14 I=IANFA,IORDX
2205      XUMORD=X(I)
2206      X(I)=X(IENDA+1-I)
2207      X(IENDA+1-I)=XUMORD
2208      YUMORD=Y(I)
2209      Y(I)=Y(IENDA+1-I)
2210   14 Y(IENDA+1-I)=YUMORD
2211      IDREHX=1
2212C
2213    2 IDREHU=0
2214      IF (U(IANFN).LT.U(IENDN)) GOTO 12
2215C     DAS U-FELD WIRD ANSTEIGEND GEORDNET
2216      IORDU=(IENDN-IANFN+1)/2
2217      DO 16 I=IANFN,IORDU
2218      UUMORD=U(I)
2219      U(I)=U(IENDN+1-I)
2220   16 U(IENDN+1-I)=UUMORD
2221      IDREHU=1
2222C
2223   12 IEXPO=0
2224      IANFNN=IANFN
2225    3 IF (U(IANFNN).GE.X(IANFA)) GOTO 13
2226C   EXTRAPOLATION BEI X(IANFA) VERHINDERN.
2227      IEXPO=1
2228      IANFNN=IANFNN+1
2229      GOTO 3
2230   13 IENDNN=IENDN
2231    4 IF (U(IENDNN).LE.X(IENDA)) GOTO 5
2232C   EXTRAPOLATION BEI X(IENDA) VERHINDERN.
2233      IEXPO=1
2234      IENDNN=IENDNN-1
2235      GOTO 4
2236    5 IF(IEXPO.EQ.1) WRITE(6,98) X(IANFA),X(IENDA),U(IANFN),U(IENDN),
2237     1 U(IANFNN),U(IENDNN)
2238   98 FORMAT(1H0,5X,'DAS REFERENZFELD GEHT VON X(IANFA) =',1PE10.3,1X,
2239     1'BIS X(IENDA) =',E10.3/1H ,'INTERPOLIERT WERDEN SOLLTE VON U(IANFN
2240     2) =',E10.3,' BIS U(IENDN) =',E10.3/1H ,5X,'WERTE AUSSERHALB VON',
2241     3E10.3,' UND',E10.3,' SIND DAHER NICHT BERECHNET.'/1H )
2242C
2243      IF (IEXPO.EQ.1) IERROR=2
2244      IANF=IANFA
2245      DO 20 J=IANFNN,IENDNN
2246      DO 10 I=IANF,IENDA
2247      IF (X(I)-U(J)) 10,7,8
2248    7 V(J)=Y(I)
2249      GOTO 22
2250    8 V(J)=Y(I-1)+(Y(I)-Y(I-1))/(X(I)-X(I-1))*(U(J)-X(I-1))
2251      GOTO 22
2252   10 CONTINUE
2253      GOTO 20
2254   22 IANF=I
2255   20 CONTINUE
2256      IF (IDREHX.EQ.0) GOTO 21
2257C     DAS X- UND Y-FELD WERDEN IN DIE AUSGANGS POSITION ZURUECKGESETZT
2258      DO 18 I=IANFA,IORDX
2259      XUMORD=X(I)
2260      X(I)=X(IENDA+1-I)
2261      X(IENDA+1-I)=XUMORD
2262      YUMORD=Y(I)
2263      Y(I)=Y(IENDA+1-I)
2264   18 Y(IENDA+1-I)=YUMORD
2265C
2266   21 IF (IDREHU.EQ.0) RETURN
2267      DO 19 I=IANFN,IORDU
2268      UUMORD=U(I)
2269      U(I)=U(IENDN+1-I)
2270      U(IENDN+1-I)=UUMORD
2271      VUMORD=V(I)
2272      V(I)=V(IENDN+1-I)
2273   19 V(IENDN+1-I)=VUMORD
2274C
2275      RETURN
2276      END
2277      SUBROUTINE GERIN (H,F,ERG,IA,IE,IERROR)
2278C     STAND VOM 20. FEB. 1974
2279C     TRAPEZINTEGRATION VON NICHT AEQUIDISTANTEN TABELL. FUNKTIONEN
2280C     ANZAHL DER ABSZISSENPUNKTE DARF GERADE UND UNGERADE SEIN
2281C     H=ABSZISSE,F=ORDINATE,ERG=INTEGRAL
2282C     INTEGRIERT WIRD VON H(IA) BIS H(IE)
2283C     FALLS NUR EIN STUETZWERT VORHANDEN (IA=IE) , WIRD DEM INTEGRAL
2284C     DIE ORDINATE DES STUETZWERTES ZUGEORDNET.
2285C     IERROR =0 IM NORMALFALL
2286C     IERROR=1, WENN IA GROESSER IE IST
2287C     IERROR=2 WENN IA=IE IST
2288C     IERROR=3 WENN ABSZISSENDIFFERENZEN ZU KLEIN WERDEN
2289C
2290      DIMENSION H(IE),F(IE)                                                   
2291C
2292      IERROR=0
2293      ERG=0.0
2294      IF (IE-IA.GE.0) GOTO 6
2295      WRITE(6,5) IA,IE
2296    5 FORMAT('0 IN DER SUBROUTINE GERIN  IST IA=',I3,'GROESSER IE=',I3)
2297      IERROR=1
2298      RETURN
2299C
2300    6 IF (IE-IA.NE.0) GOTO 3
2301      WRITE (6,15) IA
2302   15 FORMAT('0 IN DER SUBROUTINE GERIN IST IA=IE=',I3,' UND ERG=F(IA)')
2303      ERG = F(IA)
2304      IERROR=2
2305      RETURN
2306C
2307    3 IAP=IA+1
2308      DO 1 I=IAP,IE
2309      DIFH=H(I)-H(I-1)
2310      IF (ABS(DIFH/H(I)).GT.1.E-4) GOTO 4                                     
2311      WRITE (6,25) I
2312   25 FORMAT ('0 IN DER SUBROUTINE GERIN WIRD DIFH KLEINER ALS 1.E-4 BEI
2313     1 I=',I4)
2314C     DIE WARNUNG ENTSPRICHT EINEM FEHLER VON MAXIMAL 1 PROZENT
2315      IERROR=3
2316C
2317    4 ERG=ERG+(F(I)+F(I-1))*DIFH
2318    1 CONTINUE
2319      ERG=ERG/2.                                                               
2320C
2321      RETURN
2322      END
Note: See TracBrowser for help on using the repository browser.