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