      program volpis

*
*               ^ ^   VOLPIS
*              <* *> -------~~~*
*                -      |  |
*              !   !   !  !
*
****************************************************************
* THIS CODE HAS BEEN SUBMITTED TO COMPUTERS AND GEOSCIENCES
* ON JUNE 29TH, 2006, TOGETHER WITH THE RESEARCH PAPER: 
* "A FREQUENCY DOMAIN INVERSION CODE TO RETRIEVE TIME-DEPENDENT
* PARAMETERS OF LONG PERIOD VOLCANIC EVENTS" 
* BY S.CESCA AND T.DAHM
****************************************************************
* S. Cesca, Univ. Hamburg, December 2005
* S. Cesca, version 4, UHH, March 2006
* S. Cesca, version 6, UHH, May 2006
* S. Cesca, version 8, UHH, June 2006
* S. Cesca, version 9, UHH, December 2006
*
* 1. Comments on VOLPIS original code
*
* VOLPIS, VOlcanic Long Period Inversion for the Source
* is an inversion program, which main structure is based on 
* mtinvers code (Dahm).
* It uses observed displacements and synthetic Green Function
* to determinate source mechanism of LP event as the sum
* of a moment tensor component (MT) and a single force (SF).
* No previous assumption are done on time history of this
* source components, except that this has the same time length 
* then used traces.
* Inversion can be done both in time domain, using the approach
* formulated by Ohminato et al. (1998) and Chouet et al. (2003),
* or in frequency domain. Note that this last approach differs
* from amplitude spectra moment tensor inversion technique, as
* in this code it is fit the full spectra (amplitude and phase).
*
* 2. Comments on later versions 
* 
* Version 4 (UHH, Feb 2006)
* Version 4 include determination of solution including 
* 2 possible constrains on the source history: common
* source history for all source components, and common MT +
* common SF.
* These are obtained using as input the unconstrained source
* inverted by VOLPIS (as done by the original version), and using
* singular value decomposition approach (Vasco,  1989, BSSA; 
* Numerical recipes in Fortran)
* Output include fits relative to constrained source.
*
* Version 6-7
* Version 6-7 include the possibility of forward modelling
* some new files must be included:
*   - volpis6.input, is similar to older version volpis.input
*     but has an extra last line.
*     the line as the following format
*        ex. 1, 0, 45, 60, 3, 6
*        1: run forward modeling (other values run the inversion)
*(VERS7) 0: indicates no stf convolution (other numbers represent
*                           length in seconds of Fz. Brustle-Muller)
*        45: strike angle of the crack (in degree)
*        60: dip angle of the crack (in degree)
*        3: mu       
*        6: lambda
*
* Version 8, May-June 2006
* Version 8 is summaryzing all changes introduced in the previous
* versions, including some new features:
*  - improvement in the management of SAC input files
*  - generalization of the data windows to include also North
*    and East direction. This is helpful to handle the problem
*    of the inversion for different source locations, as rotation
*    of the axis does not have to be considered in the preprocessing.
*  - added file volpis.erg to check correct response of the code
*  - generalization of forward modeling, including 3 possible cases:
*    generic source (0), pure crack (1), pure DC (2).
*  - improvement of the code script by adding subroutines
*    CALMIS (calculate misfit, L2 norm) and 
*    SFORW (solving forward problem, calculating synthetic 
*    seismograms for given MT and GF) 
* Version 8 is described in the paper Cesca and Dahm (2006) which 
* is under preparation and will be soon submitted to a scientific
* journal.
* Version 8 is accompanied by text files describing input/output 
* formats and how the code works.
*                          Simone Cesca, UHH, June 2006 
*
* Version 9
* This version of volpis code include an estimation of errors
* for the unconstrained solution.
* Errors traces are saved in files named:
*         emt*.dat, esf*.dat (error)
*         zmt*.dat, zsf*.dat (unconstrained solution + error)
*         ymt*.dat, ysf*.dat (unconstrained solution - error) 
*
* Simone Cesca, UHH, November/December 2006
*
*
*
* 3. I/O Files (updated for version 4)
*
* Input files
* volpis.input : inversion parameters 
* retard.dat   : traces parameters
* displ.c.n    : observed displacements, coordinate c station n
* MTx.c.n      : Green functions related to MT components
* SFx.c.n      : Green functions related to SF components
*
* Output files
* result.fit   : fits of the inversions
*                  line 1: fit of unconstrained source (case 1)
*                  line 2: fit of common history source (case 2)
*                  line 3: fit of common MT + common SF (case 3)
* dispi*       : filtered displacements (case 1) 
* synti*       : filtered synthetics (case 1)
* dispf*       : filtered displacements (case 2) 
* synth*       : filtered synthetics (case 2)
* dispc*       : filtered displacements (case 3) 
* syntc*       : filtered synthetics (case 3)
* imt*.dat     : mt(t) source solution (case 1)
* isf*.dat     : sf(t) source solution (case 1)
* mmt*.dat     : mt(t) source solution (case 2)
* fsf*.dat     : sf(t) source solution (case 2)
* cmt*.dat     : mt(t) source solution (case 3)
* csf*.dat     : sf(t) source solution (case 3)
* volpis.erg
*
*********************************************************************************
*
* Some comments on the version of the original code are below the text
*********************************************************************************

********************************************************************************
*     0. DECLARATIVE
********************************************************************************

      include 'volpis.inc'
      include 'codeco_common.f'

      ! --------- other local parameter ------------
      integer*4 iballa(wmaxIBALLA)
      character*1 cbuf(c_bufsize),simrisp
      integer*4 nchecksum,ierr
      integer i4,i5,i6
      integer igsein,igseout ! (1) output waveforms in gse2, else in 13bit ascii
      integer ngridmod
      real    xxx(wmax)

      !---------------------------------------------
      if (c_sigsize.ne.wmaxIBALLA) then
       write(*,*)"c_sigsize.ne.wmaxIBALLA"
       write(*,*)"change parameter in codeco_common.f or mtinvers.inc"
       stop
      end if
     
      igsein  = 0
      igseout = 0
      pi = 4.*atan(1.)
      rad = 180./pi
      ci = cmplx(0.,1.)
      iadd = 0    ! Daten werden um 0 Punkte verlaengert
      iadd = 20    ! Daten werden um 20 Punkte verlaengert



********************************************************************************
*     MAIN PROGRAM
********************************************************************************

********************************************************************************
*     1. READ PARAMETER FROM PROMPT OR USING VOLPIS.INPUT
********************************************************************************
      ! control output file
      open(unit=40,file='volpis.erg')

      ! plot VOLPIS drawing
      call introvolp      

      !-------------------------------------------
      ! name of control file
      !-------------------------------------------
      write(*,*)'Step 1: Reading volpis6.input input file'
      write(40,*)'Step 1: Reading volpis6.input input file'
c      write(*,*)'Name des Steuerungsfiles (retard.dat)'
      read(*,*)cinput

      !-------------------------------------------
      !             inversion mode               
      !isampf : root of 2, idicates number of frequency samples for FFT
      !tap :  Taper on time window in percentage (-0.25 < tap <0.25),
      !       negative values enlarge chosen window
      !iscal: 0 no scaling, 1 scaling proportional to r
      !igame = 0: 6 indep MT 
      !        1: Single Force
      !        2: SF + MT
      !ifreq = 0: time domain inversion 
      !        1: frequency domain inversion 
      !--------------------------------------------
c      write(*,*)'npts (isampf) Potenz 2,',
c     & ' Taper-length (-0.25 < tap < 0.25)',
c     & ' lopass fcu, fc',' iscal',' igame',' ifr'
      read(*,*)isampf,tap,fcu,fc,iscal,igame,ifr
      write(40,*)'isampf ',isampf
      write(40,*)'tap    ',tap
      write(40,*)'fcu    ',fcu
      write(40,*)'fc     ',fc
      write(40,*)'iscal  ',iscal
      write(40,*)'fgame  ',igame
      write(40,*)'ifr    ',ifr
      read(*,*)eta,mu,qp,qs
      
      wfikt = 0.0                ! must be zero for seismic sensors

      if(tap.lt.-0.25)tap = -0.25
      if(tap.gt.+0.25)tap = +0.25
      fa = log(real(isampf))/log(2.)
      write(40,*)'isampf=',isampf,' power',fa
      if(fa-int(fa).gt.0.000001)then
       write(*,*)' ERROR: isampf must be power of 2'
       write(40,*)' ERROR: isampf must be power of 2'
       stop
      end if
      if(real(isampf).le.1.6*isamp(1))then   ! not used anymore, f.k. 3.1.01
       write(*,*)'ERROR: isampf too small, or isamp too large'
       write(40,*)'ERROR: isampf too small, or isamp too large'
       stop
      end if
      lxx=isampf/2+1
      if(fc.le.0)then
       write(*,*)' ERROR: fc must be > 0 '
       write(40,*)' ERROR: fc must be > 0 '    
       stop
      end if

      !--------------------------------------------
      ! input output format
      ! igsein,igseout: 0=13bit ascii format,
      !                 1=gse format, 
      !                 2=sac-binary format
      !--------------------------------------------

      write(*,*)'igsein,igseout'
      read(*,*)igsein,igseout
      write(40,*)'I-mode ',igsein
      write(40,*)'O-mode ',igseout


      read(*,*)v6i,v6istf,v6fmod
      if (v6i.eq.1) then
       write(*,*)'Forward mode chosen!'     
       write(40,*)'Forward mode, parameters:'
       if (v6istf.eq.0) then
        write(40,*)'STF type : none, already included in GF'
       else
        write(40,*)'STF type :      Brustle-Muller'
        write(40,*)'STF length :    ',v6istf
       end if
       if (v6fmod.eq.0) then
        write(*,*)'Moment tensor : Generic Source'
        write(40,*)'Moment tensor : Generic Source'
	read(*,*)v8m(1),v8m(2),v8m(3),v8m(4),v8m(5),v8m(6)
	read(*,*)v6f1,v6f2,v6f3	
	write(40,*)'          M11 : ',v8m(1)
	write(40,*)'          M12 : ',v8m(2)
	write(40,*)'          M22 : ',v8m(3)
	write(40,*)'          M13 : ',v8m(4)
	write(40,*)'          M23 : ',v8m(5)
	write(40,*)'          M33 : ',v8m(6)
	write(40,*)'          F1  : ',v6f1
	write(40,*)'          F2  : ',v6f2
	write(40,*)'          F3  : ',v6f3
       else if (v6fmod.eq.1) then
        write(*,*)'Moment tensor : Pure crack'
        write(40,*)'Moment tensor : Pure crack'
        read(*,*)v6str,v6dip,v6mu,v6lam            
	write(40,*)'       Strike : ',v6str
	write(40,*)'       Dip    : ',v6dip
	write(40,*)'       Mu     : ',v6mu
	write(40,*)'       Lambda : ',v6lam
       else if (v6fmod.eq.2) then
        write(*,*)'Moment tensor : Pure DC'
        write(40,*)'Moment tensor : Pure DC'
        read(*,*)v6str,v6dip,v6rak            
	write(40,*)'       Strike : ',v6str
	write(40,*)'       Dip    : ',v6dip
	write(40,*)'       Rake   : ',v6rak
       else
        write(*,*)'ERROR: Source type, incorrect input'
        write(40,*)'ERROR: Source type, incorrect input'
	stop
       end if 
      else 
       write(*,*)'Inverse mode'
       if (v6istf.eq.0) then
        write(*,*)'STF type : none, already included in GF'
       else
        write(*,*)'STF type :      Brustle-Muller'
        write(*,*)'STF length :    ',v6istf
       end if
      end if
      

c     FORWARD MODELING
c     setting source term configuration
c      
      if (v6i.eq.1) then
c     Calculating moment tensor components

      if (v6fmod.eq.0) then
      ! case of generic MT

       xm1d(1,1)=v8m(1)
       xm1d(1,2)=v8m(2)
       xm1d(2,2)=v8m(3)
       xm1d(1,3)=v8m(4)
       xm1d(2,3)=v8m(5)
       xm1d(3,3)=v8m(6)        
       xf1d(1)=v6f1
       xf1d(2)=v6f2
       xf1d(3)=v6f3  

      else if (v6fmod.eq.1) then
      ! case of pure crack

      !original principal axis tensor
      !  lambda     0        0
      !    0     lambda      0
      !    0        0    lambda+2mu
      scalf=(3*v6lam)+(2*v6mu)
      xm0d(1,1) = v6lam/scalf
      xm0d(1,2) = 0.d0
      xm0d(2,1) = 0.d0 
      xm0d(1,3) = 0.d0
      xm0d(3,1) = 0.d0
      xm0d(2,2) = v6lam/scalf
      xm0d(2,3) = 0.d0
      xm0d(3,2) = 0.d0
      xm0d(3,3) = (v6lam+(2*v6mu))/scalf

      !Euler matrix
      
      do i=1,3
       do j=1,3
        rza(i,j)=0.d0
	rxb(i,j)=0.d0
	rzg(i,j)=0.d0
       end do
      end do

      v6str=v6str*PI/180
      v6dip=v6dip*PI/180
      rza(1,1)=dcos(v6str)
      rza(1,2)=-dsin(v6str)
      rza(2,1)=dsin(v6str)
      rza(2,2)=dcos(v6str)
      rza(3,3)=1.d0
      rxb(1,1)=1.d0
      rxb(2,2)=dcos(v6dip)
      rxb(2,3)=-dsin(v6dip)
      rxb(3,2)=dsin(v6dip)
      rxb(3,3)=dcos(v6dip)
   
      call matmul1(3,3,3,rza,3,rxb,3,rhilf,3)
      
      do k=1,3
       do m=1,3
        xm1d(k,m)=0.d0
	do i=1,3
	 do j=1,3
          xm1d(k,m)=xm1d(k,m)+xm0d(i,j)*rhilf(k,i)*rhilf(m,j) 
	 end do
	end do
       end do
      end do
      xf1d(1)=0.d0
      xf1d(2)=0.d0
      xf1d(3)=0.d0  

      else if (v6fmod.eq.2) then
      ! case of pure double couple
      ! originally coded in mtdekomp.f, T. Dahm

       v6str = pi*v6str/180.
       v6dip = pi*v6dip/180.
       v6rak = pi*v6rak/180.

       xm0d(1,1) = 1.d0
       xm0d(1,2) = 0.d0
       xm0d(2,1) = xm0d(1,2) 
       xm0d(1,3) = 0.d0
       xm0d(3,1) = xm0d(1,3)
       xm0d(2,2) = 0.d0
       xm0d(2,3) = 0.d0
       xm0d(3,2) = xm0d(2,3)
       xm0d(3,3) = -1.d0

c      Euler Matrix
       do i=1,3
        do j=1,3
         rza(i,j) = 0.d0
         rxb(i,j) = 0.d0
         rzg(i,j) = 0.d0
         ry45(i,j) = 0.d0
        end do
       end do

       rza(1,1) = dcos(v6str)
       rza(1,2) = -dsin(v6str)
       rza(2,1) = +dsin(v6str)
       rza(2,2) = dcos(v6str)
       rza(3,3) = 1.d0
       rxb(1,1) = 1.d0
       rxb(2,2) = dcos(v6dip)
       rxb(2,3) = -dsin(v6dip)
       rxb(3,2) = +dsin(v6dip)
       rxb(3,3) = dcos(v6dip)
       rzg(1,1) = dcos(v6rak)
       rzg(1,2) = +dsin(v6rak)
       rzg(2,1) = -dsin(v6rak)
       rzg(2,2) = dcos(v6rak)
       rzg(3,3) = 1.d0
       phi = pi*45.d0/180.
       ry45(1,1) = dcos(phi)
       ry45(1,3) = +dsin(phi)
       ry45(2,2) = 1.d0
       ry45(3,1) = -dsin(phi)
       ry45(3,3) = dcos(phi)

       call matmul1(3,3,3,rza,3,rxb,3,rhilf,3)
       call matmul1(3,3,3,rhilf,3,rzg,3,rbet,3)
       do k=1,3
        do m=1,3
         rhilf(k,m) = rbet(k,m)
        end do
       end do

       call matmul1(3,3,3,rhilf,3,ry45,3,rbet,3)

       do k=1,3
        do m=1,3
         xm1d(k,m) = 0.d0
         do i=1,3
          do j=1,3
           xm1d(k,m) = xm1d(k,m) + xm0d(i,j)*rbet(k,i)*rbet(m,j)
          end do
	 end do
	end do
       end do
      xf1d(1)=0.d0
      xf1d(2)=0.d0
      xf1d(3)=0.d0  
      
      else
       write(*,*)'ERROR: incorrect source type'
       write(40,*)'ERROR: incorrect source type'
       stop
      end if

      write(*,*)'M11 ',xm1d(1,1)
      write(*,*)'M12 ',xm1d(1,2),xm1d(2,1)
      write(*,*)'M13 ',xm1d(1,3),xm1d(3,1)
      write(*,*)'M22 ',xm1d(2,2)
      write(*,*)'M23 ',xm1d(2,3),xm1d(3,2)
      write(*,*)'M33 ',xm1d(3,3)
      write(*,*)'F1  ',xf1d(1)
      write(*,*)'M2  ',xf1d(2)
      write(*,*)'M3  ',xf1d(3)
      write(40,*)'M11 ',xm1d(1,1)
      write(40,*)'M12 ',xm1d(1,2),xm1d(2,1)
      write(40,*)'M13 ',xm1d(1,3),xm1d(3,1)
      write(40,*)'M22 ',xm1d(2,2)
      write(40,*)'M23 ',xm1d(2,3),xm1d(3,2)
      write(40,*)'M33 ',xm1d(3,3)
      write(40,*)'F1  ',xf1d(1)
      write(40,*)'M2  ',xf1d(2)
      write(40,*)'M3  ',xf1d(3)
      
      
      end if

********************************************************************************
*     2. READ CONTROL FILE RETARD.DAT 
********************************************************************************

      write(*,*)'Step 2: Reading retard.dat control file'
      write(40,*)'Step 2: Reading retard.dat control file'
      open(2,file=cinput)
      idummy = isamp(1)
      call comment(2)
      read(2,'(a25)')cdir
      call searchbl1(cdir,26,ibl)
      j=0
      nstatc(1)=0
      nstatc(2)=0
      nstatc(3)=0
      nstatc(4)=0
      nstatc(5)=0
      nstatc(6)=0
      nstatc(7)=0
      nstatc(8)=0
      
1     j=j+1
      call comment(2)
      read(2,123,end=2)
     &  comp(j),istat(j),cphas(j),vel(j),iretard(j),
     &  iretardg(j),toff(j),azi(j),isamp(j),wt(j)
      write(*,123)
     &  comp(j),istat(j),cphas(j),vel(j),iretard(j),
     &  iretardg(j),toff(j),azi(j),isamp(j),wt(j)
123   format(a1,1x,i3,1x,a8,1x,f7.3,1x,2(i5,1x),2(f8.3,1x),i3,1x,f6.3)
       if(isamp(j).le.0)isamp(j)=idummy

       ifreq(j)=ifr

       toff(j) = toff(j)/rad
       dazi = azi(j)/rad

       sa(j,1) = dsin(toff(j))
       sa(j,2) = 0.
       sa(j,3) = dcos(toff(j))

        vel(j) = vel(j)*1000.
        if(comp(j).eq.'r' .or. comp(j).eq.'R')then
         nstatc(1)=nstatc(1)+1
         ikomp(j)=1
         comp(j)='r'
        else if(comp(j).eq.'z' .or. comp(j).eq.'Z')then
         nstatc(2)=nstatc(2)+1
         ikomp(j)=2
         comp(j)='z'
        else if(comp(j).eq.'t' .or. comp(j).eq.'T')then
         nstatc(3)=nstatc(3)+1
         ikomp(j)=3
         comp(j)='t'
	else if (comp(j).eq.'n' .or. comp(j).eq.'N')then
         nstatc(7)=nstatc(7)+1
         ikomp(j)=7
         comp(j)='n'
        else if(comp(j).eq.'e' .or. comp(j).eq.'E')then
         nstatc(8)=nstatc(8)+1
         ikomp(j)=8
         comp(j)='e'
	else
         write(*,*)' component must be r, t, z, n or e '
         stop 
        end if

       ! changing window length depending on tap
       if(tap.lt.0)then
        isamp(j) = int( 1.0*isamp(j) / (1.+2.0*tap) )
        iretard(j) = iretard(j) + int(isamp(j)*tap) 
        iretardg(j) = iretardg(j) + int(isamp(j)*tap) 
c        write(*,*)iretard(j),iretardg(j)
       end if

       ! ------- ein paar tests -----------
       lx=isampf

       if(real(lx).le.1.0*isamp(j))then
        write(*,*)'in retard.dat: j=',j,' phas=',cphas(j)
        write(*,*)'isampf too small, or isamp too large'
        write(40,*)'in retard.dat: j=',j,' phas=',cphas(j)
        write(40,*)'isampf too small, or isamp too large'
        stop
       end if
       if(istat(j).gt.999)then
        write(*,*)'istat(j).gt.999, j=',j
        write(40,*)'istat(j).gt.999, j=',j
        stop
       end if

       go to 1

2     close(2)

      nstat = j-1
      if(nstat.gt.nmax)then
       write(*,*)'nstat.gt.nmax'
       write(40,*)'nstat.gt.nmax'
       stop
      end if
      if(tap.lt.0.0)tap = -tap      ! jetzt gilt 0 < tap < 0.25

      iwl=isamp(1)

********************************************************************************
*     3. READ GROUND DISPLACEMENT DATA, PROCESSING AND WRITE INTO DISPI.*
********************************************************************************

      write(*,*)'Step 3: Reading displ.* data files'    
      write(40,*)'Step 3: Reading displ.* data files'         
      
      if(igseout.eq.1)then
       cfil0(1)= 'dispf.r.gse'
       cfil0(2)= 'dispf.z.gse'
       cfil0(3)= 'dispf.t.gse'
       cfil0(7)= 'dispf.n.gse'
       cfil0(8)= 'dispf.e.gse'
      else
       cfil0(1)= 'dispf.r'
       cfil0(2)= 'dispf.z'
       cfil0(3)= 'dispf.t'
       cfil0(7)= 'dispf.z'
       cfil0(8)= 'dispf.t'     
      end if
      ifil(1) = 31
      ifil(2) = 32
      ifil(3) = 33
      ifil(7) = 37
      ifil(5) = 38
      do i=1,3
       if(nstatc(i).gt.0)open(ifil(i),file=cfil0(i))
      end do
      do i=7,8
       if(nstatc(i).gt.0)open(ifil(i),file=cfil0(i))
      end do
      ione = -1
      
      do i=1,nstat   !...loop over stations as given in retard.dat

        if(istat(i).lt.10)then
         write(char3,'(a1,i1)')'.',istat(i)
        else if(istat(i).lt.100)then
         write(char3,'(a1,i2)')'.',istat(i)
        else
         write(char3,'(a1,i3)')'.',istat(i)
        end if
        cfil(i) = cdir(1:ibl)//'displ.'//comp(i)//char3

        if(igsein.eq.0)then   ! 13 bit ascii format
         open(1,file=cfil(i))
         read(1,'(a1)')cdummy
         read(1,110)nlay,mdeck,iso,idummy
         do k=1,nlay
          read(1,'(5f10.4)')x1,x2,(nv(j),j=1,3)
         end do
         read(1,150) nstat1
         if(nstat1.ne.1)then
          write(*,*)'nstat1.ne.1 in ',cfil(i)
          write(40,*)'nstat1.ne.1 in ',cfil(i)
          stop
         end if
         read(1,160)x1
         read(1,170)dtr
        else if(igsein.eq.1) then
         write(*,*)'input in gse2 not yet implemented'
         write(40,*)'input in gse2 not yet implemented'
         stop
        else                  ! sac binary input


          !call rsac1(cfil(i),xxx,lx,abal,dtr,wmaxIBALLA,ierr)
          call rsac1(cfil(i),xxx,lx,abal,dtr,wmax,ierr)
          if (ierr.gt.0) then
           write(*,"('reading error , ierr= ',i4)") ierr
           write(40,"('reading error , ierr= ',i4)") ierr
           write(40,*)'lx  = ',lx
           write(40,*)'dtr = ',dtr
           write(40,*)'abal= ',abal
           if(ierr.eq.108) write(*,*)'file ',cfil(i),' not existing'
           ! 801 = 'file not evenly spaced'
           ! 803 = 'only idim points read'
           ! call getfhv('E',x1,ierr)
           ! write(*,"('ending value of indep. variable, e=',f15.3)")x1
           stop
          else if (ierr.lt.0) then
           write(40,*)'too many points, file ',cfil(i),' cutted',ierr
          end if
           write(40,*)'station i=',i,' read C'
          call getfhv('DIST',x1,ierr)
           re(i) = x1
          call getfhv('AZ',x1,ierr)
          az = x1
          call getfhv ('USER1',x1,ierr)
           komp = int(x1)
          call getkhv('KSTNM',cinput,ierr)

           write(40,*)'station i=',i,' read D'
           write(40,*)' lx    ',lx
           write(40,*)' abal  ',abal
           write(40,*)' dtr   ',dtr
           write(40,*)' ierr  ',ierr
           write(40,*)' dist  ',re(i)
           write(40,*)' az    ',az
           write(40,*)' kstnm ',cinput
           write(40,*)' nevid ',komp 
           write(40,*)'station i=',i,' read E'

        end if

        dt=dtr
        ione = ione+1
      
        if (fc.gt.0) then 
       !-------------------------------------------
       ! calculate cosine taper 
       ! letztlich wird isamp Bereich zum rechnen + Inversion verwendet.
       ! Der taper ist isamp*tap (0 < tap < 0.25) auf jeder Seite
       ! beginnt bei iretard.
       ! Wenn das Nutzsignal den isamp-Bereich voll ausfuellt (z.B. vorne
       ! zu knapp gelegt) dann kann es evtl. zu syst. Fehlern durch den 
       ! Taper kommen.
       !-------------------------------------------
         fscal = 1000000.

         ifenl = int(isamp(i)*tap+0.5)
         if (ifenl .eq. 0) then
          ifenl = 1
          write(40,*)'ifenl was 0, changed to 1!!!!'
         endif
         do k=1,isamp(i)
          if(k .le. +ifenl)then
           ifen(i,k) = fscal*0.5*(1.+cos(pi*(1.*k-ifenl)/ifenl))
          else if(k.ge. isamp(i)-ifenl)then
           ifen(i,k) = fscal*0.5*
     &           (1.+cos(pi*(1.*k-(isamp(i)-ifenl))/ifenl))
          else 
           ifen(i,k) = fscal*1.0
          end if
         end do

        end if
 
c        write(*,*)'Step 3.3'
 
       !-------------------------------------------
       !  diverse einmalige Befehle wenn ione == 0, ionel ==0
       !-------------------------------------------
        if(ione.eq.0) then
         ! ----------- header von dispi.* schreiben
         do j1=1,3
          if(nstatc(j1).gt.0)then

          if(igseout.eq.1)then
           write(ifil(j1),'(a50)')'Raumwellen'
           write(ifil(j1),'(a50)')
     & 'scale with: GSEcounts * calib *1.E-9 ym'
           write(ifil(j1),151) nstatc(j1)
           write(ifil(j1),161)(re(1),j=1,nstatc(j1))
	   
           write(ifil(j1),171)0.,0.,dt,0.,0.
          else
           write(ifil(j1),'(a50)')'Raumwellen'
           write(ifil(j1),111)1,0,0,j1
           write(ifil(j1),'(5f10.4)')x1,x2,(nv(j),j=1,3)
           write(ifil(j1),151) nstatc(j1)
           write(ifil(j1),161)(re(1),j=1,nstatc(j1))
           
           write(ifil(j1),171)0.,0.,dt,0.,0.	   
          end if

          end if
         end do
 
       !--------------------------------------------
       !  BP Filter in frequency                       
       !--------------------------------------------
         if(fc.gt.0)then 
          ifail=0
          write(40,*)'Traces Filter, f=',fc,'Hz'
          dtr = dt
          call filtkaus(dtr,fcu,fc,1,isampf,filt,ifail)
          if(ifail.ne.0)then
            write(40,*)'ifail filter .ne. 0, stopped'
            write(*,*)'ifail filter .ne. 0, stopped'
            stop
          end if
         end if

         ione = ione + 1
        end if

       !--------------------------------------------
       !     end diverse einmalige Befehle
       !--------------------------------------------

        if(igsein.eq.0)then
         read(1,5001)re(i),ABAL,komp,lx,BALMAX,az
         !if(iscal.gt.0)BALMAX = BALMAX*sqrt(1.e+4+re(i)**2)/1.e+3
         if(iscal.gt.0)BALMAX = BALMAX*re(i)/1000.
        else
         if(iscal.gt.0)then
          BALMAX = re(i)/1000.
         else
          BALMAX = 1.
         end if
        end if
        BALMAX = BALMAX*sign(1.0,wt(i))
        wt(i) = abs(wt(i))

        !..... Daten werden spaeter um iadd Punkte verlaengert
        abal = abal -iadd*dtr

        ibal = int(abal/dtr)

        if(abs(azi(i)-az).gt.0.1)then  ! test whether azi is correct
         write(40,*)'warning'
         write(40,*)'i=',i,'azi .ne. az, r=',re(i),azi(i),az
         az=azi(i)
	 write(40,*)'...corrected'
        end if
        itmp=isampf

        if((lx+iadd).lt.itmp)then
         write(40,*)'lx+iadd.lt.isamp, i=',i,' r=',re(i)
         write(40,*)'lx=',lx,' iadd=',iadd,' isamp=',itmp
        end if
        if((lx+iadd).gt.wmaxIBALLA)then
         write(*,*)'lx+iadd.gt.wmaxIBALLA , r=',re(i)
         write(40,*)'lx+iadd.gt.wmaxIBALLA , r=',re(i)
	 stop
        end if

        if(igsein.eq.0)then
         read(1,5002)(IBALLA(i1),i1=1,lx)
         !....... Daten um iadd punkte nach vorne verlaengern
         do i1=lx+iadd,iadd+1,-1
          iballa(i1)=iballa(i1-iadd)
         end do
         do i1=1,iadd
          iballa(i1)=iballa(iadd+1)*0.5*
     &              (1-cos(pi*(real(i1)-1.)/real(iadd)))
         end do
        else
         do i1=lx+iadd,iadd+1,-1
          xxx(i1)=xxx(i1-iadd)
         end do
         do i1=1,iadd
          xxx(i1)=xxx(iadd+1)*0.5*
     &              (1-cos(pi*(real(i1)-1.)/real(iadd)))
         end do
        end if

        lx = lx+iadd

       !--------------------------------------------
       !  Taper in time domain                                
       !--------------------------------------------
        i2 = iretard(i)-ibal
        ! kbeg = i2+int(isamp(i)/2) - int(itmp/2)-1
        kbeg = i2-1
        !write(*,*)abal,ibal,i2,kbeg
        if(kbeg.lt.0)then
         write(*,*)cfil(i),kbeg
         write(*,*)'ERROR!!! kbeg < 0 '
         write(40,*)cfil(i),kbeg
         write(40,*)'ERROR!!! kbeg < 0 '
         !kbeg = 0
         stop
        end if
        write(40,*)'kbeg',kbeg,isamp(i),lx,iadd
        do k=kbeg+1,kbeg+isamp(i)       ! Zeitfenster isamp
          if(igsein.eq.0)then
             x2=(iballa(k)/10000.-5.0001)*balmax/4.999
             ddd(i,k-kbeg)=iballa(k)
	     dddmax(i)=balmax
          else
             x2 = xxx(k)*balmax
          end if
          x2 = x2*ifen(i,k-kbeg) / fscal
c         if(kbeg+isamp(i).le.lx-iadd)then changed F.K. 26.11.02
          if(k.le.lx-iadd)then
            cx(k-kbeg) = cmplx(real(x2),0.0)
          else
            cx(k-kbeg) = cmplx(0.0,0.0)
          endif
        end do
c       do k=kbeg+isamp(i)+1,kbeg+itmp  ! Rest mit Nullen, alte Formulierung f.k.8.5.02
        do k=kbeg+isamp(i)+1,wmax+kbeg  ! Rest mit Nullen
          x2 = 0.
          cx(k-kbeg) = cmplx(real(x2),0.0)
        end do

c        write(*,*)'Step 3.10'
        if(cphas(i)(8:8).eq.'P')then
         q = qp
         !write(*,*)cphas(i)(8:8),', qp = ',q,' used'
        else if(cphas(i)(8:8).eq.'S')then
         q = qs
         !write(*,*)cphas(i)(8:8),', qs = ',q,' used'
        else
         q = qp
         !write(*,*)cphas(i)(8:8),', q  = ',q,' used, unknown phase'
        end if
        tstar = re(i) / (vel(i)*q)
c        write(*,*)'here below'
c        write(*,*)cphas(i)(8:8),re(i),vel(i),q,tstar

       !--------------------------------------------
       ! Filter in frequency                               
       !--------------------------------------------
        call fork(isampf,cx,-1.0)
        wnyq = 2.*pi / (2.*dt )
        !v6 add if
	if (v6i.ne.1) then
	 do k=1,isampf
          cx(k)=cx(k)*cmplx(filt(k),0.)
         end do
	end if
        do k=2,lxx
         w=2.*pi*real(k-1)/(real(isampf)*dt)
	 carg=cmplx(0.5*w*tstar,-w*tstar*log(w/wnyq)/pi)
	 cx(k)=cx(k)*cdexp(carg)
         cx(isampf-k+2)=conjg(cx(k))
        end do
        cx(lxx)=cmplx(real(cx(lxx)),0.)

        if(ifr.eq.1)then
	 do k=1,lxx
c          write(*,*)i,k,' fdat ',cx(k)
          fdat(i,k)=cx(k)
c          synth(k,i) = cabs(cx(k))
         end do
        end if
	
        call fork(isampf,cx,+1.0)
        do i1=1,isampf
         dat(i,i1)=dble(cx(i1))
         synth(i1,i) = real(cx(i1))
        end do

       !--------------------------------------------
       ! Write to dispi.*                          
       !--------------------------------------------

        abal = iretard(i)*dt

        lx = isamp(i)
        komp = ikomp(i)
         
        do i1=1,lx
          xxx(i1) = synth(i1,i)
        enddo

        BALMAX = AMAX(xxx,lx)
        BAL = abs(AMIN(xxx,lx))
        IF(BAL.GT.BALMAX) BALMAX = BAL
        IF(BALMAX.EQ.0.) then
         bal = 1.
        else
         bal = BALMAX
        end if

        if(igseout.eq.1)then  !-- 24 bit gse2 format

         call init_hdrvars_mtinvers
         hdr_debug    = 4   !  muss wohl groesser als 3 sein ?
         hdr_nsamp    = lx

         hdr_smprate = 1.0 / dt

         call add_time_shift(abal) !-add t-shift abal, calculate new date
         hdr_station  = cphas(i)
         hdr_chan     = '__'//comp(i)
         if(komp.eq.1)then          ! rad
           hdr_hang   = azi(i)
           hdr_vang   = 90.0
         else if(komp.eq.2)then     ! ver, but comp=s has also komp=2
           hdr_hang   = -1.
           hdr_vang   = 0.0
         else if(komp.eq.3)then     ! tra
           hdr_hang   = azi(i)-90.
           hdr_vang   = 90.0
         end if
         hdr_calfac  = BAL*59.6  ! bal * 1E9/2**24
         !write(*,*)cphas(i),' calfac= ',hdr_calfac

         do i1=1,lx
          synth(i1,i)=synth(i1,i)/hdr_calfac
          iballa(i1)=1.E+9 *synth(i1,i)
         end do
         
         do i1=1,10  ! check free space in character variables
          if(hdr_station(i1:i1).eq.' ')hdr_station(i1:i1)='_'
         end do

         iout = ifil(ikomp(i))
         call gseout(iout,cbuf,iballa,nchecksum)

        else     ! ---------- 13 bit ascii format

         do i1=1,lx
          synth(i1,i)=4.999*synth(i1,i)/BAL + 5.0001
          iballa(i1)=10000.*synth(i1,i)
         end do
         write(ifil(ikomp(i)),6001)
     &      re(i),ABAL,komp,lx,dddmax(i),azi(i)
         write(ifil(ikomp(i)),6002)(ddd(i,i1),i1=1,lx)

        end if

        if(igsein.eq.0)then
         read(1,5003)re(i),sn,se,sz
         close(1)
        end if


        if(comp(i).eq.'r') then     ! radial component

         sn = +cos(azi(i)/rad)*1.  -sin(azi(i)/rad)*0.
         se = +sin(azi(i)/rad)*1.  +cos(azi(i)/rad)*0.
         sz = 0.
         vkomp(i,1) = 1.
         vkomp(i,2) = 0.
         vkomp(i,3) = 0.
         vkomp(i,4) = 0.
         vkomp(i,5) = 0.	 

        else if(comp(i).eq.'t') then     ! transversal component

         sn = +cos(azi(i)/rad)*0.  -sin(azi(i)/rad)*1.
         se = +sin(azi(i)/rad)*0.  +cos(azi(i)/rad)*1.
         sz = 0.
         vkomp(i,1) = 0.
         vkomp(i,2) = 1.
         vkomp(i,3) = 0.
         vkomp(i,4) = 0.
         vkomp(i,5) = 0.	 

        else if(comp(i).eq.'z') then     ! vertical component

         sn = 0.
         se = 0.
         sz = 1.
         vkomp(i,1) = 0.
         vkomp(i,2) = 0.
         vkomp(i,3) = 1.
         vkomp(i,4) = 0.
         vkomp(i,5) = 0.	 

        else if(comp(i).eq.'n') then     ! north component

         sn = 1.
         se = 0.
         sz = 0.
         vkomp(i,1) = 0.
         vkomp(i,2) = 0.
         vkomp(i,3) = 0.
         vkomp(i,4) = 1.
         vkomp(i,5) = 0.	 

        else if(comp(i).eq.'e') then     ! east component

         sn = 0.
         se = 1.
         sz = 0.
         vkomp(i,1) = 0.
         vkomp(i,2) = 0.
         vkomp(i,3) = 0.
         vkomp(i,4) = 0.
         vkomp(i,5) = 1.	 	 

        end if
        if(igseout.ne.1)write(ifil(ikomp(i)),6003)re(i),sn,se,sz
      end do

      !--------------------------------------------
      ! Closing dispi.* 
      !--------------------------------------------
      do i=1,8
       if(nstatc(i).gt.0) close(ifil(i))
      end do

110   FORMAT(4I5)
140   FORMAT(20X,3F10.4)
150   FORMAT(I5)
160   FORMAT(7F10.3)
170   FORMAT(20x,F10.6,25X)
5001  FORMAT(2F15.5,I5/,I10,5X,E15.4,f15.5)
5002  FORMAT(16I5)
5003  FORMAT(4F15.5)


********************************************************************************
*     4. READ AND PROCESS GREEN FUNCTIONS
********************************************************************************

      if (v6i.ne.1) go to 2233
      
      !option use STF=Fz.Brustle-Muller
      do i=1,nstat
       winlen=dt*isamp(i)
       if (v6istf.ge.winlen) then
       write(40,*)'WARNING: stf too long or data too short at ',i
       end if
      end do
      if (v6istf.ne.0) then
       do ihfkt=1,isampf
        hfktlen=ihfkt*dt
        if (hfktlen.le.v6istf) then
         hfkt(ihfkt)=0.75*pi/v6istf*dsin(pi*ihfkt*dt/v6istf)**3
        else
         hfkt(ihfkt)=0.
        end if
        do igreen=1,9
         sfun(igreen,ihfkt)=cmplx(hfkt(ihfkt),0.e0)
         if (igreen.eq.1) then
  	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xm1d(1,1) 
         else if (igreen.eq.2) then
 	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xm1d(1,2)
         else if (igreen.eq.3) then
	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xm1d(2,2)
	 else if (igreen.eq.4) then
	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xm1d(1,3)
         else if (igreen.eq.5) then
	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xm1d(2,3)
         else if (igreen.eq.6) then
	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xm1d(3,3)
         else if (igreen.eq.7) then 
	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xf1d(1)
         else if (igreen.eq.8) then
	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xf1d(2)
         else if (igreen.eq.9) then
	  sfun(igreen,ihfkt)=sfun(igreen,ihfkt)*xf1d(3)	 
	 end if
        end do
       end do 
      else
       do igreen=1,9
        if (igreen.eq.1) then
  	  sfun(igreen,1)=cmplx(xm1d(1,1),0.) 
         else if (igreen.eq.2) then
  	  sfun(igreen,1)=cmplx(xm1d(1,2),0.) 
         else if (igreen.eq.3) then
  	  sfun(igreen,1)=cmplx(xm1d(2,2),0.) 
	 else if (igreen.eq.4) then
  	  sfun(igreen,1)=cmplx(xm1d(1,3),0.) 
         else if (igreen.eq.5) then
  	  sfun(igreen,1)=cmplx(xm1d(2,3),0.) 
         else if (igreen.eq.6) then
  	  sfun(igreen,1)=cmplx(xm1d(3,3),0.) 
         else if (igreen.eq.7) then 
  	  sfun(igreen,1)=cmplx(xf1d(1),0.) 
         else if (igreen.eq.8) then
  	  sfun(igreen,1)=cmplx(xf1d(2),0.) 
         else if (igreen.eq.9) then
  	  sfun(igreen,1)=cmplx(xf1d(3),0.)
	 end if
        end do
       
      end if    
      !end generation STF (at this point in f-domain)     

2233  continue

      write(*,*)'Step 4: Reading Greens Functions' 
      write(40,*)'Step 4: Reading Greens Functions' 

      cfil2(1) = 'MT1.r'
      cfil2(2) = 'MT1.z'
      cfil2(3) = 'MT1.t'
      cfil2(4) = 'MT1.n'
      cfil2(5) = 'MT1.e'     

      cfil2(6) = 'MT2.r'
      cfil2(7) = 'MT2.z'
      cfil2(8) = 'MT2.t'
      cfil2(9) = 'MT2.n'
      cfil2(10) = 'MT2.e'
      
      cfil2(11) = 'MT3.r'
      cfil2(12) = 'MT3.z'
      cfil2(13) = 'MT3.t'
      cfil2(14) = 'MT3.n'
      cfil2(15) = 'MT3.e'
      
      cfil2(16) = 'MT4.r'
      cfil2(17) = 'MT4.z'
      cfil2(18) = 'MT4.t'
      cfil2(19) = 'MT4.n'
      cfil2(20) = 'MT4.e'

      cfil2(21) = 'MT5.r'
      cfil2(22) = 'MT5.z'
      cfil2(23) = 'MT5.t'
      cfil2(24) = 'MT5.n'
      cfil2(25) = 'MT5.e'

      cfil2(26) = 'MT6.r'    
      cfil2(27) = 'MT6.z'
      cfil2(28) = 'MT6.t'
      cfil2(29) = 'MT6.n'
      cfil2(30) = 'MT6.e'

      cfil2(31) = 'SF1.r'
      cfil2(32) = 'SF1.z'
      cfil2(33) = 'SF1.t'
      cfil2(34) = 'SF1.n'
      cfil2(35) = 'SF1.e'

      cfil2(36) = 'SF2.r'
      cfil2(37) = 'SF2.z'
      cfil2(38) = 'SF2.t'
      cfil2(39) = 'SF2.n'
      cfil2(40) = 'SF2.e'

      cfil2(41) = 'SF3.r'     
      cfil2(42) = 'SF3.z'
      cfil2(43) = 'SF3.t'      
      cfil2(44) = 'SF3.n'
      cfil2(45) = 'SF3.e'            
                
      do j=1,nstat   !............ loop over stations
       ifenl = int(isamp(j)/4)
       if(istat(j).lt.10)then
        write(char3,'(a1,i1)')'.',istat(j)
       else if(istat(j).lt.100)then
        write(char3,'(a1,i2)')'.',istat(j)
       else
        write(char3,'(a1,i3)')'.',istat(j)
       end if

       igreen=9
       if(comp(j).eq.'r')then
        cfil2a(1) = cdir(1:ibl)//cfil2(1)//char3
        cfil2a(2) = cdir(1:ibl)//cfil2(6)//char3
        cfil2a(3) = cdir(1:ibl)//cfil2(11)//char3
        cfil2a(4) = cdir(1:ibl)//cfil2(16)//char3
        cfil2a(5) = cdir(1:ibl)//cfil2(21)//char3
        cfil2a(6) = cdir(1:ibl)//cfil2(26)//char3	
        cfil2a(7) = cdir(1:ibl)//cfil2(31)//char3
        cfil2a(8) = cdir(1:ibl)//cfil2(36)//char3
        cfil2a(9) = cdir(1:ibl)//cfil2(41)//char3
       else if(comp(j).eq.'z')then
        cfil2a(1) = cdir(1:ibl)//cfil2(2)//char3
        cfil2a(2) = cdir(1:ibl)//cfil2(7)//char3
        cfil2a(3) = cdir(1:ibl)//cfil2(12)//char3
        cfil2a(4) = cdir(1:ibl)//cfil2(17)//char3
        cfil2a(5) = cdir(1:ibl)//cfil2(22)//char3
        cfil2a(6) = cdir(1:ibl)//cfil2(27)//char3	
        cfil2a(7) = cdir(1:ibl)//cfil2(32)//char3
        cfil2a(8) = cdir(1:ibl)//cfil2(37)//char3
        cfil2a(9) = cdir(1:ibl)//cfil2(42)//char3
       else if(comp(j).eq.'t')then
        cfil2a(1) = cdir(1:ibl)//cfil2(3)//char3
        cfil2a(2) = cdir(1:ibl)//cfil2(8)//char3
        cfil2a(3) = cdir(1:ibl)//cfil2(13)//char3
        cfil2a(4) = cdir(1:ibl)//cfil2(18)//char3
        cfil2a(5) = cdir(1:ibl)//cfil2(23)//char3
        cfil2a(6) = cdir(1:ibl)//cfil2(28)//char3	
        cfil2a(7) = cdir(1:ibl)//cfil2(33)//char3
        cfil2a(8) = cdir(1:ibl)//cfil2(38)//char3
        cfil2a(9) = cdir(1:ibl)//cfil2(43)//char3
       else if(comp(j).eq.'n')then
        cfil2a(1) = cdir(1:ibl)//cfil2(4)//char3
        cfil2a(2) = cdir(1:ibl)//cfil2(9)//char3
        cfil2a(3) = cdir(1:ibl)//cfil2(14)//char3
        cfil2a(4) = cdir(1:ibl)//cfil2(19)//char3
        cfil2a(5) = cdir(1:ibl)//cfil2(24)//char3
        cfil2a(6) = cdir(1:ibl)//cfil2(29)//char3	
        cfil2a(7) = cdir(1:ibl)//cfil2(34)//char3
        cfil2a(8) = cdir(1:ibl)//cfil2(39)//char3
        cfil2a(9) = cdir(1:ibl)//cfil2(44)//char3
       else if(comp(j).eq.'e')then
        cfil2a(1) = cdir(1:ibl)//cfil2(5)//char3
        cfil2a(2) = cdir(1:ibl)//cfil2(10)//char3
        cfil2a(3) = cdir(1:ibl)//cfil2(15)//char3
        cfil2a(4) = cdir(1:ibl)//cfil2(20)//char3
        cfil2a(5) = cdir(1:ibl)//cfil2(25)//char3
        cfil2a(6) = cdir(1:ibl)//cfil2(30)//char3	
        cfil2a(7) = cdir(1:ibl)//cfil2(35)//char3
        cfil2a(8) = cdir(1:ibl)//cfil2(40)//char3
        cfil2a(9) = cdir(1:ibl)//cfil2(45)//char3
       end if

       do i=1,igreen

        if(igsein.eq.0)then   ! 13 bit ascii format
         open(1,file=cfil2a(i))
         read(1,'(a1)')cdummy
         read(1,110)nlay,mdeck,iso,idummy
         do k=1,nlay
          read(1,'(a1)')cdummy
         end do
         read(1,150) nstat1
         if(nstat1.ne.1)then
          write(*,*)'nstat1.ne.1 in ',cfil2a(i)
	  write(40,*)'nstat1.ne.1 in ',cfil2a(i)
          stop
         end if
         read(1,160)x1
         read(1,170)dtr
        else
          call rsac1(cfil2a(i),xxx,lx,abal,dtr,wmaxIBALLA,ierr)
          !call rsac1(cfil2a(i),xxx,lx,abal,dtr,wmax,ierr)
	  if (ierr.gt.0) then
	   write(*,"('reading error , ierr= ',i4)") ierr
           if(ierr.eq.108) write(*,*)'file ',cfil(i),' not existing'
           write(40,"('reading error , ierr= ',i4)") ierr
           if(ierr.eq.108) write(40,*)'file ',cfil(i),' not existing'        
	   stop
          else if (ierr.lt.0) then
           write(40,*)'too many points, file ',cfil(i),' cutted',ierr
          end if
          call getfhv('DIST',x1,ierr)
           r = x1
           !re(i) = x1
          call getfhv('AZ',x1,ierr)
           az = x1
          !call getnhv ('NEVID',komp,ierr)
          call getfhv ('USER1',x1,ierr)
           komp = int(x1)
          !call getfhv ('USER2',ABAL,ierr) ! bereits in rsac1 eingelesen
          call getkhv('KSTNM',cinput,ierr)
        end if
        if(abs(dt-dtr).ge.0.001) then
          write(40,*)'WARNING: dt.ne.dtr in ',cfil2a(i)
          !stop
        end if

        if(igsein.eq.0)then
         read(1,5001)r,ABAL,komp,lx,BALMAX,az
         if(iscal.gt.0)BALMAX = BALMAX*re(j)/1000.
         ! BALMAX = BALMAX*wt(j)
        else
         BALMAX = 1.
         if(iscal.gt.0)BALMAX = BALMAX*re(j)/1000.
        end if

        ibal = (abal/dtr)
        if(i.eq.1)then
         itmp=isampf
         i2 = iretardg(j) - ibal
         if(i2.lt.0)then
          write(*,*)cfil2a(i),j,i2
          write(*,*)'Error!!! i2 < 0, stop',j,iretardg(j),ibal
          write(40,*)cfil2a(i),j,i2
          write(40,*)'Error!!! i2 < 0, stop'
          stop
         end if
         kbeg = i2-1
         if(kbeg.lt.0)then
          write(*,*)cfil2a(i),j,i2,itmp,isamp(j),kbeg
          write(*,*)'Error!!! kbeg < 0, stop'
          write(40,*)cfil2a(i),j,i2,itmp,isamp(j),kbeg
          write(40,*)'Error!!! kbeg < 0, stop'
          stop
         end if
        end if

        if((abs(r-re(j))).gt.(r/100))then
         write(40,*)'WARNING: r=',r,'.ne.re(j)=',re(j),
     &                        ' in ',cfil2a(i),j
         !stop
        end if

        if(igsein.eq.0)then
         read(1,5002)(iballa(i1),i1=1,lx)
         do i1=kbeg+1,kbeg+isamp(j)
          x2=(iballa(i1)/10000. -5.0001)*balmax/4.999
          x2 = x2*ifen(j,i1-kbeg) / fscal
          cx(i1-kbeg)=cmplx(x2,0.e0)
         end do
        else
         do i1=kbeg+1,kbeg+isamp(j)
c          x2=xxx(k)*balmax
          x2=xxx(i1)*balmax
          x2 = x2*ifen(j,i1-kbeg) / fscal
          cx(i1-kbeg)=cmplx(x2,0.e0)
         end do
        end if
        do i1=kbeg+isamp(j)+1,kbeg+itmp
         x2=0.
         cx(i1-kbeg)=cmplx(x2,0.e0)
        end do
        ! already including taper in time domain

        call fork(itmp,cx,-1.0)
        do i1=1,lxx
         w=2.*pi*real(i1-1)/(real(isampf)*dt)
         !v6 add if
         if (v6i.ne.1) then
          if(fc.gt.0) cx(i1)=cx(i1)*cmplx(filt(i1),0.)
         end if
         gmatfre(j,i,i1)= cx(i1)
        end do
        ! already including BP filter in frequency dom.

        do i1=1,lxx
c         write(*,*)j,i,i1,gmatfre(j,i,i1)
	 ccx(i1)=gmatfre(j,i,i1)
	 ccx(isampf-i1+2)=conjg(ccx(i1))
	end do
        ccx(lxx)=cmplx(real(ccx(lxx)),0.)
        call fork(isampf,ccx,+1.0)

        do k=1,isampf
	 gmattim(j,i,k) =real( ccx(k) )
c         write(*,*)j,i,k,gmattim(j,i,k)
        end do
        ! filtered and taperes green functions in time dom.
	! are saved in gmattim(j,i,k)
 
        if(igsein.eq.0)then
         read(1,5003)r,sn,se,sz
         close(1)
        end if

       end do      
       ! end local loop over GF (1-9)

      end do      
      ! end local loop over traces

      do j=1,nstat
       re(j) = re(j)*1000.
       azi(j) = azi(j)/rad
       dazi = azi(j)
       w =  sa(j,1)
       sa(j,1) = w*dcos(dazi)
       sa(j,2) = w*dsin(dazi)
      end do


********************************************************************************
*     5. PREPARING INVERSION SECTION
********************************************************************************

      write(*,*)'Step 5: Preparing inversion'
      write(40,*)'Step 5: Preparing inversion'
      
      lx = isampf

      do ig = 1,20
       costlist(ig)=1.e+35
       cbest(ig)='-12345'
      end do
      iout = 9
      open(iout+1,file='result.dat')
      write(iout+1,*)'0, M0, M11, M12, M22, M13, M23, M33',
     & ' strike1 dip1 rake1 strike2 dip2 rake2 ISO%  DC% ',
     & ' P azi P plunge T azi T plunge SF0 SF1 SF2 SF3'
     
      !........ number of free source parameter depends on igame
      if(igame.eq.1)then
       k1 = 6
      else if(igame.eq.1)then
       k1 = 3
      else if(igame.eq.2)then
       k1 = 9
      else
       write(*,*)'Error : igame incorrect'
       write(40,*)'Error : igame incorrect'
       stop
      end if


C     INVERSION SWITCH
      if (v6i.eq.1) go to 740 


********************************************************************************
*     6. INVERSION IN TIME DOMAIN
********************************************************************************
      write(*,*)'Step 6: Inversion in time domain (not impl.)'
      write(40,*)'Step 6: Inversion in time domain (not impl.)'

      if (ifr.eq.0) then
       write(*,*)'Sorry, t-domain inv not yet implemented!'
       write(40,*)'Sorry, t-domain inv not yet implemented!'
       stop
      end if
      

********************************************************************************
*     7. INVERSION IN FREQ DOMAIN
********************************************************************************

      write(*,*)'Step 7: Inversion in frequency domain'      
      write(40,*)'Step 7: Inversion in frequency domain'      

      if (ifr.eq.1) then

cc     GF(w) are saved in: gmatfre(ntraces,ngreen(9),npt(lxx=isampf/2+1))
cc     U(w)  are saved in: fdat(ntraces,npt(lxx))

       do i=1,nstat
        do wi=1,lxx
	 uw(i)=fdat(i,wi)
        end do
       end do

       chinv=0

c      added for vers.9
       do wi=1,lxx
        wsing(wi)=0
       end do
c      <vers.9        

       do wi=1,lxx  ! loop over each sampled frequency 
c       do wi=2,lxx  ! loop over each sampled frequency 
 
c        write(*,*)'Frequency sample ',wi       
c       assign vector U(w_i), with dimension=number traces -> uw(i), i=1..nstat
        do i=1,nstat
	 uw(i)=fdat(i,wi)
        end do
        do i=nstat+1,nmax
	 uw(i)=(0.,0.)
        end do	

c       assign matrix G(w_i), with dimensions:
c                                       lines=number traces
c                                       columns=number free source parameters (9)
c                            -> gw(i,j), i=1..nstat, j=1..9        
        do i=1,nstat
	 do j=1,9
	  gw(i,j)=gmatfre(i,j,wi)
	 end do
	end do

c       calculate G^t, transposed of G ; dim G^t = (9,nstat)
        do i=1,nstat
	 do j=1,9
	  gwt(j,i)=gw(i,j)
	 end do
	end do        

c       calculate G^t G ; dim(9,9)
        do i=1,9
	 do j=1,9
	  hws(i,j)=(0.0)
	  do k=1,nstat
	   hws(i,j)=hws(i,j)+(gwt(i,k)*gw(k,j))
	  end do
	 end do
        end do

c       calculate (G^t G)^{-1}; dim(9,9)
c       the inverse matrix is saved in hwi
c       using external subroutines 
*       .... Gauss Elimination 

        do i=1,9
	 do j=1,9
	  hw(i,j)=hws(i,j)
	 end do
	end do
        dimmat=9
        nfpar=(dimmat*2)+1   !19

        call zgeco(hw,nfpar,dimmat,ipiv,deven,nonimp)
c        write(*,*)' condition of matrix = ',deven
        if(1.d0+deven .eq. 1.d0)then
         write(40,*)wi,' : acond singular to working precision'
         do j=1,9
          sourfre(j,wi)=(0.0)
          wsing(wi)=1
         end do
	 go to 700
c         stop
        end if
	chinv=chinv+1
        call zgedi(hw,nfpar,dimmat,ipiv,dete,bvect,01)
    
        do i=1,9
	 do j=1,9
	  hwi(i,j)=hw(i,j)
	 end do
	end do
    
c       end of calculation inverted complex matrix    

c       calculate (G^t G)^{-1}) G^t ; dim (gwi) = (9,nstat)
c       call matmulc (hwi,9,9,gwt,9,nstat,gwi,9,nstat)
        do i=1,9
	 do j=1,nstat
	  gwi(i,j)=(0.0)
	  do k=1,9
	   gwi(i,j)=gwi(i,j)+(hwi(i,k)*gwt(k,j))
	  end do
	 end do
        end do

c       added in version 9, S.C.
        do i=1,9
	 do j=1,nstat
          gwiw(wi,i,j)=gwi(i,j)
	 end do
        end do
c       <vers.9

c       calculate M = (G^t G)^{-1}) G^t U ; dim (mw) = (9,1)
c       call matmulc (gwi,9,nstat,uw,nstat,1,mw,9,1)
        do i=1,9
	 j=1
	 mw(i)=(0.0)
	 do k=1,nstat
	  mw(i)=mw(i)+(gwi(i,k)*uw(k))
	 end do
	end do

c       at this point M(w_i) have been determinated, with dimension=free sour. par. (9)                  
c                           -> mw(j), j=1..9

c       save spectra of components of source in sourfre(j,w_i). j=comp, w_i=freq
        do j=1,9
         sourfre(j,wi)=mw(j)
        end do

  700   continue
	
       end do !end loop over sampled frequencies

       if (chinv.eq.0) then
        write(*,*)'ERROR: Inversion incorrect'
	write(*,*)'Complex matrices at all frequencies can not be inverted'
        write(40,*)'ERROR: Inversion incorrect'
	write(40,*)'Complex matrices at all frequencies can not be inverted'
       end if

c      calculate synthetic spectra in syntfre(i,w_i), i=1..nstat, w_i=freq
       do wi=1,lxx
        do j=1,9
         mw(j)=sourfre(j,wi)
         do i=1,nstat
          gw(i,j)=gmatfre(i,j,wi)
         end do
        end do
c        call matmulc(gw,nstat,9,mw,9,1,sfre,nstat,1)
        do i=1,nstat
	 sfre(i)=(0.0)
	 do k=1,9
	  sfre(i)=sfre(i)+(gw(i,k)*mw(k))
	 end do
	end do
 
 	do i=1,nstat
         syntfre(i,wi)=sfre(i)
        end do
       end do
     
c      calculate ifft of M(w_i)
       do j=1,9
        do wi=1,lxx
         csourtim(wi)=sourfre(j,wi)
	 csourtim(isampf-wi+2)=conjg(csourtim(wi))	 
        end do
	call fork (isampf,csourtim,+1.0)
        do ti=1,isampf
	 sourtim(j,ti)=real(csourtim(ti))
	end do
       end do

c      calculate ifft of synthetic spectra
       do i=1,nstat
        do wi=1,lxx
         csourtim(wi)=syntfre(i,wi)
	 csourtim(isampf-wi+2)=conjg(csourtim(wi))	 
        end do
	call fork (isampf,csourtim,+1.0)
        do ti=1,isampf
	 synttim(i,ti)=real(csourtim(ti))
         synth(ti,i)=synttim(i,ti)
	end do
       end do

c      calculate ifft of observed spectra
       do i=1,nstat
        do wi=1,lxx
         dhelp(wi)=fdat(i,wi)
	 dhelp(isampf-wi+2)=conjg(dhelp(wi))	 
        end do
	call fork (isampf,dhelp,+1.0)
        do ti=1,isampf
	 fildat(i,ti)=real(dhelp(ti))
	end do
       end do

c      evaluate misfit between synthetic seismograms U and observed DAT       
c       call calmis(fildat,synttim,nstat,isampf,misf)
       call calmis(fildat,synttim,nmax,tmax,nstat,isamp,misf)

c****************************************
c      Error estimation (version 9)     *
c****************************************

c        for each frequency e=d-Gm=spectra(data)-spectra(synth)       
c        E = e^t e = sum (d-s)**2 (saved in var: ete=(eter,etei))

       do wi=1,lxx
        eter(wi)=0.
	etei(wi)=0.
       end do
       do wi=1,lxx
        if(wsing(wi).ne.1) then
         do i=1,nstat
	  pasval=fdat(i,wi)-syntfre(i,wi)
 	  eter(wi)=eter(wi)+(real(pasval))**2
	  pasval=pasval*(0.,-1.)
	  etei(wi)=etei(wi)+(real(pasval))**2
	 end do
         eter(wi)=eter(wi)/nstat
         etei(wi)=etei(wi)/nstat
        end if
       end do

c        Calculate m_j(w) and Dm_j(w)
c        they have a real and immaginary component
c        in var: rvaluem,ivaluem,rdeltam,ideltam
c        all are matrices (1..9,1..nfreq)
c
c        if m = (G^t G)^-1 G^t d = H d
c        then   
c        Dm_j = sum (H_jk^2 Dd_k^2) = Dd^2 sum (H_jk^2)

       do wi=1,lxx
        if(wsing(wi).ne.1) then
         do j=1,9
          deltam(j,wi)=0
 	  do i=1,nstat
           pasval=gwiw(wi,j,i)
	   deltam(j,wi)=deltam(j,wi)+(real(pasval))**2
           pasval=pasval*(0.,-1.)
	   deltam(j,wi)=deltam(j,wi)+(real(pasval))**2
	  end do
	  rdeltam(j,wi)=deltam(j,wi)*eter(wi)
	  ideltam(j,wi)=deltam(j,wi)*etei(wi)
	  rdeltam(j,wi)=sqrt(rdeltam(j,wi))
	  ideltam(j,wi)=sqrt(ideltam(j,wi))
          pasval=sourfre(j,wi)
          rvaluem(j,wi)=real(pasval)
	  pasval=(sourfre(j,wi))*(0,-1)
	  ivaluem(j,wi)=real(pasval)
	 end do
        end if
       end do 

c      BOOTSTRAP

c      Loop over tested solutions, e.g. 200
       do iboot=1,200

c       Loop over frequencies
        do wi=1,lxx  

         if(wsing(wi).ne.1) then
          do j=1,9
	  
c          Random number generator for real component 
c           tempmr=gennor(rvaluem(j,wi),rdeltam(j,wi))
c           jgas=-j*wi*iboot
c           tempmr=gasdev(jgas)
           tempmr=rgauss(0.,1.)
	   tempmr=rvaluem(j,wi)+(tempmr*rdeltam(j,wi))	   
c          Random number generator for immaginary component
c           tempir=gennor(ivaluem(j,wi),ivaluem(j,wi))
c           jgas=2*jgas
c 	    tempir=gasdev(jgas)
           tempir=rgauss(0.,1.)
	   tempir=ivaluem(j,wi)+(tempir*ideltam(j,wi))
	   
c           if (iboot.eq.10) then
c	    if (j.eq.1) then
c             write(*,*)'MR DR ',rvaluem(j,wi),rdeltam(j,wi)
c	     write(*,*)'mr ',tempmr
c             write(*,*)'MI DI ',ivaluem(j,wi),ideltam(j,wi)
c	     write(*,*)'mi ',tempir
c	    end if
c	   end if	  

c          Reconstruct new m(w)
           tempm(j,wi)=cmplx(tempmr,tempir)
          end do
         else
          tempm(j,wi)=sourfre(j,wi)
         end if
c       End loop freq.
        end do
	
c       Calculate ifft of modified M(w_i)
        do j=1,9
         do wi=1,lxx
          ttcsourtim(wi)=tempm(j,wi)
 	  ttcsourtim(isampf-wi+2)=conjg(ttcsourtim(wi))	 
         end do
	 call fork (isampf,ttcsourtim,+1.0)
         do ti=1,isampf
	  ttsourtim(iboot,j,ti)=real(ttcsourtim(ti))
	 end do
        end do

c      End loop over tested solutions
       end do
 
c      Loop over t
       do ti=1,isampf

c       Loop over source components
        do j=1,9
	 
         meanboot(j,ti)=0
	 sdevboot(j,ti)=0
	          
c        Loop over bootstrap solutions - calc media
         do iboot=1,200
          meanboot(j,ti)=meanboot(j,ti)+ttsourtim(iboot,j,ti)
         end do
         meanboot(j,ti)=meanboot(j,ti)/200
	
c        Loop over bootstrap solutions - calc s.dev.
         do iboot=1,200
          sdevboot(j,ti)=sdevboot(j,ti)+
     &           (ttsourtim(iboot,j,ti)-meanboot(j,ti))**2
         end do
         sdevboot(j,ti)=sdevboot(j,ti)/200
         sdevboot(j,ti)=sqrt(sdevboot(j,ti))
 
c       End loop components
        end do
 
c      End loop over t
       end do

c      Save variance, v(t)
       open(11,file='emt1.dat')
       open(12,file='emt2.dat')
       open(13,file='emt3.dat')
       open(14,file='emt4.dat')
       open(15,file='emt5.dat')
       open(16,file='emt6.dat')
       open(17,file='esf1.dat')
       open(18,file='esf2.dat')
       open(19,file='esf3.dat')            

       do i=1,9
       k=i+10
        do j=1,lx
         write(k,800)dt*(j-1),sdevboot(i,j)
	end do
	close(k)
       end do

c      Save mean+variance
       open(11,file='zmt1.dat')
       open(12,file='zmt2.dat')
       open(13,file='zmt3.dat')
       open(14,file='zmt4.dat')
       open(15,file='zmt5.dat')
       open(16,file='zmt6.dat')
       open(17,file='zsf1.dat')
       open(18,file='zsf2.dat')
       open(19,file='zsf3.dat')            

       do i=1,9
       k=i+10
        do j=1,lx
         write(k,800)dt*(j-1),sourtim(i,j)+sdevboot(i,j)
	end do
	close(k)
       end do

c      Save mean-variance
       open(11,file='ymt1.dat')
       open(12,file='ymt2.dat')
       open(13,file='ymt3.dat')
       open(14,file='ymt4.dat')
       open(15,file='ymt5.dat')
       open(16,file='ymt6.dat')
       open(17,file='ysf1.dat')
       open(18,file='ysf2.dat')
       open(19,file='ysf3.dat')            

       do i=1,9
       k=i+10
        do j=1,lx
         write(k,800)dt*(j-1),sourtim(i,j)-sdevboot(i,j)
	end do
	close(k)
       end do

      end if  !end options frequancy inversion   
      

 740  continue

      if (v6i.ne.1) go to 760


c     OPTION FORWARD MODELING

      write(*,*)'Step 7.5: Forward modeling'      
      write(40,*)'Step 7.5: Forward modeling'      

c      calculate ifft of observed spectra
       do i=1,nstat
        do wi=1,lxx
         dhelp(wi)=fdat(i,wi)
	 dhelp(isampf-wi+2)=conjg(dhelp(wi))	 
        end do
	call fork (isampf,dhelp,+1.0)
        do ti=1,isampf
	 fildat(i,ti)=real(dhelp(ti))
	end do
       end do
       
c      synthetic seismograms for constrain: f1=MT, f2=SF
       call sforw(lxx,gmatfre,nmax,wmax,tmax,nstat,
     &            sfun,isampf,syntforw,v6istf,v6i)

c      calculate best M0, and mul synthetics
       num(1)=0
       den(1)=0
       
       do i=1,nstat
c       V7: next line
        iwl=isamp(i)
        do ti=1,iwl
         num(1)=num(1)+(1000000*syntforw(i,ti)*
     &                  fildat(i,ti))
         den(1)=den(1)+(1000000*syntforw(i,ti)*
     &                  syntforw(i,ti))
        end do
       end do
       scalmom=num(1)/den(1)
       write(40,*)'iwl-scalmo',iwl,scalmom
       do i=1,nstat
c       V7: next line
        iwl=isamp(i)
        do ti=1,iwl
         syntforw(i,ti)= syntforw(i,ti)*scalmom
	 synth(ti,i)=syntforw(i,ti)
	end do
        do ti=iwl+1,isampf
         syntforw(i,ti)=0.
	 synth(ti,i)=0.
	end do
       end do     

c     additional output format
       do k=1,nstat
        if (istat(k).lt.10) then
         write(fileforw,'(a4,a1,a1,i1)')'syn.',comp(k),'.',istat(k)
        else if (istat(k).lt.100) then
         write(fileforw,'(a4,a1,a1,i2)')'syn.',comp(k),'.',istat(k)
        else 
         write(fileforw,'(a4,a1,a1,i3)')'syn.',comp(k),'.',istat(k)
        end if
        open(20,file=fileforw)
        do i=1,isamp(1)
         write(20,748)i*dt,synth(i,k)
        end do
        close(20)
        if (istat(k).lt.10) then
         write(fileforw,'(a4,a1,a1,i1)')'dis.',comp(k),'.',istat(k)
        else if (istat(k).lt.100) then
         write(fileforw,'(a4,a1,a1,i2)')'dis.',comp(k),'.',istat(k)
        else 
         write(fileforw,'(a4,a1,a1,i3)')'dis.',comp(k),'.',istat(k)
        end if
        open(20,file=fileforw)      
        do i=1,isamp(1)
         write(20,748)i*dt,fildat(k,i)
        end do
        close(20)
       end do
       
       
  748 format(F7.3,' ',E15.5)
          
c      evaluate misfit between synthetic seismograms U and observed DAT       
       call calmis(fildat,syntforw,nmax,tmax,nstat,isamp,misf)

c     write output file forward.str-dp.res
       
      istr=nint(v6str*180/PI)
      idp=nint(v6dip*180/PI)
      if (istr.lt.0) istr=istr+360
      if (istr.lt.10) then
       write(cstr,'(a2,i1)')'00',istr
      else if (istr.lt.100) then
       write(cstr,'(a1,i2)')'0',istr     
      else
       write(cstr,'(i3)')istr
      end if
      if (idp.lt.10) then
       write(cdp,'(a1,i1)')'0',idp     
      else
       write(cdp,'(i2)')idp     
      endif
      
      fileforw='forward.'//cstr//'-'//cdp//'.res'
      open(20,file=fileforw)
      write(20,750)istr,idp,misf,misf,scalmom
      close(20)
      fileforw='forward.'//cstr//'-'//cdp//'.mis'
      open(20,file=fileforw)
      write(20,751)istr,idp,misf
      close(20)
     
      
 750  format(I5,' ',I5,' ',F10.5,' ',F10.5,' ',E15.5)
 751  format(I5,' ',I5,' ',F10.5)
     
      go to 802

 760  continue

*****************************END INVERSION**************************************


c    >qui
c      additional output format
       do k=1,nstat
        if (istat(k).lt.10) then
         write(fileforw,'(a4,a1,a1,i1)')'syn.',comp(k),'.',istat(k)
        else if (istat(k).lt.100) then
         write(fileforw,'(a4,a1,a1,i2)')'syn.',comp(k),'.',istat(k)
        else 
         write(fileforw,'(a4,a1,a1,i3)')'syn.',comp(k),'.',istat(k)
        end if
        open(20,file=fileforw)
        do i=1,isamp(1)
         write(20,748)i*dt,synth(i,k)
        end do
        close(20)
        if (istat(k).lt.10) then
         write(fileforw,'(a4,a1,a1,i1)')'dis.',comp(k),'.',istat(k)
        else if (istat(k).lt.100) then
         write(fileforw,'(a4,a1,a1,i2)')'dis.',comp(k),'.',istat(k)
        else 
         write(fileforw,'(a4,a1,a1,i3)')'dis.',comp(k),'.',istat(k)
        end if
        open(20,file=fileforw)      
        do i=1,isamp(1)
         write(20,748)i*dt,fildat(k,i)
        end do
        close(20)
       end do
c      < add out format      

       write(*,*)'Step 8: Writing output files'       
       write(40,*)'Step 8: Writing output files'       

      !--------------------------------------------
      ! Source 
      !--------------------------------------------

c      independent source components

       open(11,file='imt1.dat')
       open(12,file='imt2.dat')
       open(13,file='imt3.dat')
       open(14,file='imt4.dat')
       open(15,file='imt5.dat')
       open(16,file='imt6.dat')
       open(17,file='isf1.dat')
       open(18,file='isf2.dat')
       open(19,file='isf3.dat')            

       somax=0
       do i=1,9
       sourekm=0
       sourekf=0
       sourek(i)=0
       k=i+10
c        do j=1,isampf
        do j=1,lx
         write(k,800)dt*(j-1),sourtim(i,j)
	 if (abs(sourtim(i,j)).ge.sourek(i)) then
	  sourek(i)=abs(sourtim(i,j))
	 end if
	end do
	close(k)
       end do
       do i=1,6
        if (sourek(i).ge.sourekm) sourekm=sourek(i)
       end do
       do i=7,9
        if (sourek(i).ge.sourekf) sourekf=sourek(i)
       end do

c      common mt 
       
       do i=1,6
        do j=1,lx
         amat(i,j)=sourtim(i,j)
        end do
       end do 
       m1=6
       n1=lx
       mda=lx
       do i=1,6
        do j=1,6
	 if (i.eq.j) bmat(i,j)=1
	 if (i.ne.j) bmat(i,j)=0
	end do
       end do
       mdb=6
       nb=6
       call svdrs(amat,mda,m1,n1,bmat,mdb,
     &               nb,svet,wsvd)
       do i=1,lx
        stfm(i)=amat(i,1)
       end do
       do i=1,6
        num(i)=0
	den(i)=0
	do j=1,lx
	 num(i)=num(i)+stfm(j)*sourtim(i,j)
	 den(i)=den(i)+stfm(j)*stfm(j)
	end do
       end do
       write(*,*)'coef ',num(1)/den(1)
       write(40,*)'COMMON MT'
       do i=1,6
        write(40,*)'SVET ',i,' ',svet(i)
       end do
       
       open(11,file='mmt1.dat')
       open(12,file='mmt2.dat')
       open(13,file='mmt3.dat')
       open(14,file='mmt4.dat')
       open(15,file='mmt5.dat')
       open(16,file='mmt6.dat')
       
       misffm=0
       do i=1,6
        mismax=0
        do j=1,lx
         write(i+10,800)dt*j,stfm(j)*num(i)/den(i)
         confunr(i,j)=stfm(j)*num(i)/den(i)
         confun(i,j)=cmplx(confunr(i,j),0.)
	end do         
	misfterm=0         
        do j=1,lx
	 misft2=abs((stfm(j)*num(i)/den(i))-sourtim(i,j))
	 misfterm=misfterm+misft2/sourekm
	end do
        write(40,*)',mt-sf ',i,' ',misfterm/lx
	misffm=misffm+misfterm/lx
       end do

       do i=11,16
        close(i)
       end do

  800 format(F10.5,' ',D20.10)


c      common sf 

       do i=1,3
        k=i+6
        do j=1,lx
         amat(i,j)=sourtim(k,j)
        end do
       end do 
       m1=3
       n1=lx
       mda=lx
       do i=1,3
        do j=1,3
	 if (i.eq.j) bmat(i,j)=1
	 if (i.ne.j) bmat(i,j)=0
	end do
       end do
       mdb=3
       nb=3
       call svdrs(amat,mda,m1,n1,bmat,mdb,
     &               nb,svet,wsvd)
       do i=1,lx
        stff(i)=amat(i,1)
       end do
       do i=1,3
        k=i+6
        num(k)=0
	den(k)=0
	do j=1,lx
	 num(k)=num(k)+stff(j)*sourtim(k,j)
	 den(k)=den(k)+stff(j)*stff(j)
	end do
       end do
       write(40,*)'COMMON SF'
       do i=1,3
        write(40,*)'SVET ',i,' ',svet(i)
       end do
         
       open(11,file='fsf1.dat')
       open(12,file='fsf2.dat')
       open(13,file='fsf3.dat')
       do i=1,3
        mismax=0
        k=i+6
        do j=1,lx
         write(i+10,800)dt*j,stff(j)*num(k)/den(k)
         confunr(k,j)=stff(j)*num(k)/den(k)
         confun(k,j)=cmplx(confunr(k,j),0.)
	end do         
	misfterm=0         
        do j=1,lx
	 misft2=abs((stff(j)*num(k)/den(k))-sourtim(k,j))
	 misfterm=misfterm+misft2/sourekf
	end do
        write(40,*)'mt-sf ',k,' ',misfterm/lx
	misffm=misffm+misfterm/lx       
       end do
       
       misffm=misffm/9

       do i=11,13
        close(i)
       end do
      
c      synthetic seismograms for constrain: f1=MT, f2=SF
       call sforw(lxx,gmatfre,nmax,wmax,tmax,nstat,
     &            confun,isampf,synmf,v6istf,v6i)

c      evaluate misfit between synthetic seismograms U and observed DAT       
       call calmis(fildat,synmf,nmax,tmax,nstat,isamp,misfitfm)


c      additional output format
       do k=1,nstat
        if (istat(k).lt.10) then
         write(fileforw,'(a4,a1,a1,i1)')'sy2.',comp(k),'.',istat(k)
        else if (istat(k).lt.100) then
         write(fileforw,'(a4,a1,a1,i2)')'sy2.',comp(k),'.',istat(k)
        else 
         write(fileforw,'(a4,a1,a1,i3)')'sy2.',comp(k),'.',istat(k)
        end if
        open(20,file=fileforw)
        do i=1,isamp(k)
         write(20,748)i*dt,synmf(k,i)
        end do
        close(20)
       end do
c       <qui


c      common source history (mt & sf)

       do i=1,9
        do j=1,lx
         amat(i,j)=sourtim(i,j)
        end do
       end do 
       m1=9
       n1=lx
       mda=lx
       do i=1,9
        do j=1,9
	 if (i.eq.j) bmat(i,j)=1
	 if (i.ne.j) bmat(i,j)=0
	end do
       end do
       mdb=9
       nb=9
       call svdrs(amat,mda,m1,n1,bmat,mdb,
     &               nb,svet,wsvd)
       do i=1,lx
        stfc(i)=amat(i,1)
       end do
       do i=1,9
        num(i)=0
	den(i)=0
	do j=1,lx
	 num(i)=num(i)+stfc(j)*sourtim(i,j)
	 den(i)=den(i)+stfc(j)*stfc(j)
	end do
       end do
       write(40,*)'COMMON MT&SF'
       do i=1,9
        write(40,*)'SVET ',i,' ',svet(i)
       end do
       
       open(11,file='cmt1.dat')
       open(12,file='cmt2.dat')
       open(13,file='cmt3.dat')
       open(14,file='cmt4.dat')
       open(15,file='cmt5.dat')
       open(16,file='cmt6.dat')
       open(17,file='csf1.dat')
       open(18,file='csf2.dat')
       open(19,file='csf3.dat')
       misfc=0
 
       do i=1,9
        mismax=0
        do j=1,lx
         write(i+10,800)dt*j,stfc(j)*num(i)/den(i)
         confunr(i,j)=stfc(j)*num(i)/den(i)
         confun(i,j)=cmplx(confunr(i,j),0.)
	end do
	misfterm=0         
        do j=1,lx
	 misft2=abs((stfc(j)*num(i)/den(i))-sourtim(i,j))
         if (i.le.6) then
	  misfterm=misfterm+(misft2/sourekm)
	 end if
         if (i.gt.6) then
	  misfterm=misfterm+(misft2/sourekf)
	 end if 
	end do
        write(40,*)'common ',i,' ',misfterm/lx
	misfc=misfc+misfterm/lx
       end do
       misfc=misfc/9
       
       do i=11,19
        close(i)
       end do


c      synthetic seismograms for constrain: f1=MT, f2=SF
       call sforw(lxx,gmatfre,nmax,wmax,tmax,nstat,
     &            confun,isampf,syncom,v6istf,v6i)


c      evaluate misfit between synthetic seismograms U and observed DAT       
c       call calmis(fildat,syncom,nstat,isampf,misfitcom)
       call calmis(fildat,syncom,nmax,tmax,nstat,isamp,misfitcom)

c      additional output format
       do k=1,nstat
        if (istat(k).lt.10) then
         write(fileforw,'(a4,a1,a1,i1)')'syc.',comp(k),'.',istat(k)
        else if (istat(k).lt.100) then
         write(fileforw,'(a4,a1,a1,i2)')'syc.',comp(k),'.',istat(k)
        else 
         write(fileforw,'(a4,a1,a1,i3)')'syc.',comp(k),'.',istat(k)
        end if
        open(20,file=fileforw)
        do i=1,isamp(k)
         write(20,748)i*dt,syncom(k,i)
        end do
        close(20)
       end do
c      <qui


      !--------------------------------------------
      ! Fit
      !--------------------------------------------
      
      open(11,file='result.fit')
      write(11,801)misf
      write(11,801)misfitfm
      write(11,801)misfitcom
      write(11,801)misffm
      write(11,801)misfc
      write(11,801)misffm*misf
      write(11,801)misfc*misf
   
  801 format(D20.10)
  802 continue
      close(11)

       
      !--------------------------------------------
      ! Waveforms 
      !--------------------------------------------

      if (v6i.eq.1) then
       iom=1
      else
       iom=3
      end if
      
      do io=1,iom  !loop over diff. constraints 

      ifr=0
 
      if(igseout.eq.1)then  ! ---------------- gse2 format --

       if (io.eq.1) then
        cfil0(1) = 'synth.r.gse'
        cfil0(2) = 'synth.z.gse'
        cfil0(3) = 'synth.t.gse'
        cfil0(7) = 'synth.n.gse'
        cfil0(8) = 'synth.e.gse'
       else if (io.eq.2) then
        cfil0(1) = 'synt2.r.gse'
        cfil0(2) = 'synt2.z.gse'
        cfil0(3) = 'synt2.t.gse'
        cfil0(7) = 'synt2.n.gse'
        cfil0(8) = 'synt2.e.gse'      
       else 
        cfil0(1) = 'syntc.r.gse'
        cfil0(2) = 'syntc.z.gse'
        cfil0(3) = 'syntc.t.gse'
        cfil0(7) = 'syntc.n.gse'
        cfil0(8) = 'syntc.e.gse'
       end if
        
       do j=1,3
        if(nstatc(j).gt.0)then

        open(ifil(j),file=cfil0(j))
        write(ifil(j),'(a50)')'Synthetisierte Raumwellen'
        write(ifil(j),'(a50)')'scale with: GSEcounts * calib *1.E-9 ym'
        WRITE(ifil(j),151) nstatc(j)
        WRITE(ifil(j),161)(re(1)/1000.,ie=1,nstatc(j))
        if(ifr.eq.0)then
         WRITE(ifil(j),171)0.,0.,dt,0.,0.
        else
         WRITE(ifil(j),171)0.,dt,1./(isampf*dt),0.,0.
        end if
      
        end if
       end do
       do j=7,8
        if(nstatc(j).gt.0)then

        open(ifil(j),file=cfil0(j))
        write(ifil(j),'(a50)')'Synthetisierte Raumwellen'
        write(ifil(j),'(a50)')'scale with: GSEcounts * calib *1.E-9 ym'
        WRITE(ifil(j),151) nstatc(j)
        WRITE(ifil(j),161)(re(1)/1000.,ie=1,nstatc(j))
        if(ifr.eq.0)then
         WRITE(ifil(j),171)0.,0.,dt,0.,0.
        else
         WRITE(ifil(j),171)0.,dt,1./(isampf*dt),0.,0.
        end if
      
        end if
       end do
     
      else      ! ------------------- 1 bit ascii files 

       if (io.eq.1) then
        cfil0(1) = 'synth.r'
        cfil0(2) = 'synth.z'
        cfil0(3) = 'synth.t'
        cfil0(7) = 'synth.n'
        cfil0(8) = 'synth.e'
       else if (io.eq.2) then
        cfil0(1) = 'synt2.r'
        cfil0(2) = 'synt2.z'
        cfil0(3) = 'synt2.t'
        cfil0(7) = 'synt2.n'
        cfil0(8) = 'synt2.e'      
       else 
        cfil0(1) = 'syntc.r'
        cfil0(2) = 'syntc.z'
        cfil0(3) = 'syntc.t'
        cfil0(7) = 'syntc.n'
        cfil0(8) = 'syntc.e'
       end if
       
       do j=1,3
        if(nstatc(j).gt.0)then

        open(ifil(j),file=cfil0(j))
        write(ifil(j),'(a50)')'Synthetisierte Raumwellen'
        write(ifil(j),111)1,0,0,j
        WRITE(ifil(j),141)vel(1)/1000.,vel(1)/1000.,vel(1)/1000.
        WRITE(ifil(j),151) nstatc(j)
        WRITE(ifil(j),161)(re(1)/1000.,ie=1,nstatc(j))
        if(ifr.eq.0)then           ! attention: assuming that ifreq is the same for all traces
         WRITE(ifil(j),171)0.,0.,dt,0.,0.
        else
         WRITE(ifil(j),171)0.,dt,1./(isampf*dt),0.,0.
        end if
      
        end if
       end do
       do j=7,8
        if(nstatc(j).gt.0)then

        open(ifil(j),file=cfil0(j))
        write(ifil(j),'(a50)')'Synthetisierte Raumwellen'
        write(ifil(j),111)1,0,0,j
        WRITE(ifil(j),141)vel(1)/1000.,vel(1)/1000.,vel(1)/1000.
        WRITE(ifil(j),151) nstatc(j)
        WRITE(ifil(j),161)(re(1)/1000.,ie=1,nstatc(j))
        if(ifr.eq.0)then           ! attention: assuming that ifreq is the same for all traces
         WRITE(ifil(j),171)0.,0.,dt,0.,0.
        else
         WRITE(ifil(j),171)0.,dt,1./(isampf*dt),0.,0.
        end if
      
        end if
       end do

      end if

      do ie=1,nstat
       if(ifr.eq.0)then
        lx0 = isamp(ie)
       else
        lx0 = lxx
       end if
       do j=1,lx0
        if (v6i.eq.1) then
	 wsynth(j,ie)=syntforw(ie,j)
	else
         if (io.eq.1) then
	  wsynth(j,ie)=synth(j,ie)
         else if (io.eq.2) then
	  wsynth(j,ie)=synmf(ie,j)
	 else 
	  wsynth(j,ie)=syncom(ie,j)
 	 end if
	end if
       end do
      end do 

      !...............  loop over traces to write
      do ie = 1,nstat

         abal = iretard(ie)*dt    ! immer der gleiche Spurbeginn

         if(ifr.eq.0)then
          lx0 = isamp(ie)
         else
          lx0 = lxx
         end if
         komp = ikomp(ie)

        !................  find real-valued maximum BAL of trace
        BALMAX = AMAX(wsynth(1,ie),lx0)
        BAL = abs(AMIN(wsynth(1,ie),lx0))
        IF(BAL.GT.BALMAX) BALMAX = BAL
        IF(BALMAX.EQ.0.) then
         bal = 1.
        else
         bal = BALMAX
        end if

        if(igseout.eq.1)then   ! ----------- gse2 format ---

         call init_hdrvars_mtinvers
         hdr_debug    = 4   !  muss wohl groesser als 3 sein ?
         hdr_nsamp    = lx0
          if(ifr.eq.0)then
           hdr_smprate = 1.0 / dt
          else
           hdr_smprate = isampf*dt
          end if

         call add_time_shift(abal) !-add t-shift abal, calculate new date

         hdr_station   = cphas(ie)
         hdr_chan      = '__'//comp(ie)
         if(komp.eq.1)then          ! rad
           hdr_hang    = azi(ie)
           hdr_vang    = 90.0
         else if(komp.eq.2)then     ! ver
           hdr_hang    = -1.
           hdr_vang    = 0.0
         else if(komp.eq.3)then     ! tra
           hdr_hang    = azi(ie)-90.
           hdr_vang    = 90.0
         else if(komp.eq.7)then     ! north
           hdr_hang    = 0.
           hdr_vang    = 90.0
         else if(komp.eq.8)then     ! east
           hdr_hang    = 90.
           hdr_vang    = 90.0
         end if

         hdr_calfac  = BAL*59.6  ! bal * 1E9/2**24
         !........conversion to 24 bit integer field based on micrometer
         do i=1,lx0
          wsynth(i,ie)=wsynth(i,ie)/hdr_calfac
          IBALLA(I)=1.E+9 *wsynth(i,ie)
         end do

         do i1=1,10  ! check free space in character variables
          if(hdr_station(i1:i1).eq.' ')hdr_station(i1:i1)='_'
         end do

         iout = ifil(ikomp(ie))
         call gseout(iout,cbuf,iballa,nchecksum)

        else                ! ----------- 13 bit ascii format ---

         do i=1,lx0
          wsynth(i,ie)=4.999*wsynth(i,ie)/BAL + 5.0001
          IBALLA(I)=10000.*wsynth(i,ie)
         end do
         write(ifil(ikomp(ie)),6001)re(ie)/1000.,ABAL,
     &    komp,lx0,BALMAX,rad*azi(ie)
         write(ifil(ikomp(ie)),6002)(IBALLA(i),i=1,lx0)
         write(ifil(ikomp(ie)),6003)re(ie)/1000.,
     &    sn,se,sz

        end if

      end do
      !..........................  close files
      do i=1,3
       if(nstatc(i).gt.0) close(ifil(i))
      end do

      end do       ! end loop diff. constraints
      
      
456   continue
      close(iout+1)
       
      write(*,*)'Ho finito!'
      write(40,*)'Ho finito!'
      close(40)

111   FORMAT(4I5,5X,'NLAY, MDECK,  ISO,  ISS(4)')
141   FORMAT(20X,3F10.4,5X,'REFL. ZONE')
151   FORMAT(I5,' ZAHL DER SEISMOGRAMME = ZAHL DER DISTANZEN')
161   FORMAT(7F10.3)
171   FORMAT(2F10.4,3F10.6,5X,'V-R - IN TIME - INC')
6001  FORMAT(2F15.5,I5/,I10,5X,E15.4,f15.5,e15.5)
6002  FORMAT(16I5)
6003  FORMAT(4F15.5)
      end

*****************************************************************
*     END PROGRAM VOLPIS
*****************************************************************



      subroutine introvolp  
******************************************************************************
*     Plot VOLPIS drawing
*     S. Cesca, June 2006                                    
****************************************************************************

      write(*,*)'                    VOLPIS'
      write(*,*)'VOlcanic Long Period Inversion for the Source'
      write(*,*)'    S. Cesca, University of Hamburg, 2005'
      write(*,*)''   
      write(*,*)'              ^ ^   VOLPIS   8 '
      write(*,*)'             <* *> -------~~~* '
      write(*,*)'               -      |  |     '
      write(*,*)'             !   !   !  !      '
      write(*,*)''
      write(*,*)''
      write(40,*)'                    VOLPIS'
      write(40,*)'VOlcanic Long Period Inversion for the Source'
      write(40,*)'    S. Cesca, University of Hamburg, 2005'
      write(40,*)''   
      write(40,*)'              ^ ^   VOLPIS   8 '
      write(40,*)'             <* *> -------~~~* '
      write(40,*)'               -      |  |     '
      write(40,*)'             !   !   !  !      '
      write(40,*)''
      write(40,*)''

      return
      end   


      subroutine calmis(ddo,ss,n1max,n2max,n1,calsam,mm)  
******************************************************************************
*     Calculate misfit (L2)
*     S. Cesca, June 2006                                    
****************************************************************************

      IMPLICIT NONE
      integer i,j,n1max,n2max,n1,n2
      integer calsam(n1max)
      real*8  mm,ddo(n1max,n2max),ss(n1max,n2max)
      real*8  dd(n1max,n2max)
      real*8  dmax,dabs,sumnum,sumden

      do i=1,n1
       n2=calsam(n1)
       do j=1,n2
        dd(i,j)=ddo(i,j)
       end do
      end do

      dmax=0
      do i=1,n1
       n2=calsam(n1)
       do j=1,n2
        dabs=abs(dd(i,j))
        if (dabs.ge.dmax) dmax=abs(dd(i,j))
       end do
      end do

      do i=1,n1
       n2=calsam(n1)
       do j=1,n2
        dd(i,j)=dd(i,j)/dmax
        ss(i,j)=ss(i,j)/dmax	
       end do
      end do

      write(*,*)'CALMIS'
      sumnum=0
      sumden=0     
      do i=1,n1
       n2=calsam(n1)
       do j=1,n2
        sumnum=sumnum+((dd(i,j)-ss(i,j))**2)
c        sumden=sumden+((abs(dd(i,j))+abs(ss(i,j)))**2)
	sumden=sumden+((dd(i,j))*(dd(i,j)))
       end do
      end do

      mm=sumnum/sumden
      
      return
      end



      subroutine sforw(sflxx,sfgmatfre,sfnmax,sfwmax,sftmax,sfnstat,
     &                 sfsfun,sfisampf,sfsynth,sfv6istf,sfv6i)  
******************************************************************************
*     Calculate synthetic seismograms for given source and GF
*     S. Cesca, June 2006                                    
****************************************************************************

      IMPLICIT NONE
      integer sflxx,sfnmax,sfwmax,sftmax,sfnstat,sfv6i
      complex*16 sfgmatfre(sfnmax,9,sfwmax/2+1)
      complex*8 sfsfun(9,sftmax)
      complex*16 sfmw(9),sfsfre(sfnmax),sfgw(sfnmax,9)
      complex*16 sfsyntfre(sfnstat,sfwmax)
      complex*16 sfsourfre(9,sfwmax)
      complex*8  sfcsourtim(sfwmax)
      real*8 sfsynttim(sfnmax,sfwmax)
      real*8   sfsynth(sfnmax,sftmax)
      integer sfisampf
      double precision sfv6istf
      integer wri,i,j,k,ti


c      calculate fft of M
c      assign M; save spectra of components of source in sourfre(j,w_i). j=comp, w_i=freq
       do j=1,9
        do wri=1,sfisampf       
         sfcsourtim(wri)=sfsfun(j,wri) 
        end do
	call fork(sfisampf,sfcsourtim,-1.0)
c        do wri=1,sflxx       
        do wri=1,sfisampf       
         sfsourfre(j,wri)=sfcsourtim(wri) 
        end do
       end do
   
c      calculate synthetic spectra in syntfre(i,w_i), i=1..nstat, w_i=freq
       do wri=1,sflxx
c       do wri=1,sfisampf       

        do j=1,9
         sfmw(j)=sfsourfre(j,wri)
         do i=1,sfnstat
          sfgw(i,j)=sfgmatfre(i,j,wri)
         end do
        end do

c        call matmulc(gw,nstat,9,mw,9,1,sfre,nstat,1)
        do i=1,sfnstat
	 sfsfre(i)=(0.0)
	 do k=1,9
	  if (sfv6i.eq.1) then
           if (sfv6istf.ne.0) then
            sfsfre(i)=sfsfre(i)+(sfgw(i,k)*sfmw(k))
           else
            sfsfre(i)=sfsfre(i)+(sfgw(i,k)*real(sfsfun(k,1)))	   
  	   end if
          else
	   sfsfre(i)=sfsfre(i)+(sfgw(i,k)*sfmw(k))
	  end if
	 end do
	end do

	do i=1,sfnstat
         sfsyntfre(i,wri)=sfsfre(i)
        end do
       end do
     
c      calculate ifft of synthetic spectra
       do i=1,sfnstat
        do wri=1,sflxx
         sfcsourtim(wri)=sfsyntfre(i,wri)
	 sfcsourtim(sfisampf-wri+2)=conjg(sfcsourtim(wri))	 
        end do
	call fork (sfisampf,sfcsourtim,+1.0)
        do ti=1,sfisampf
c	 sfsynttim(i,ti)=real(sfcsourtim(ti))
c         sfsynth(i,ti)=sfsynttim(i,ti)
         sfsynth(i,ti)=real(sfcsourtim(ti))
	end do
       end do

      return
      end



*******************************************************
      SUBROUTINE MATMULC(A,la,ca,B,lb,cb,C,lc,cc)
*******************************************************
c     VOLPIS subroutine, S.C.,  7.12.05
c     Product of complex matrices
c     C(a1,a3) = A(a1,a2) x B(a2,a3)
      integer    la,lb,lc,ca,cb,cc
      complex*16 A(la,ca),B(lb,cb),C(lc,cc),app
      integer    i,j,k
 
      write(*,*)'matrix a'
      do i=1,la
       write(*,*)'line ',la
       do j=1,ca
        write(*,*)i,j,A(i,j)
       end do
      end do
      
      write(*,*)'matrix b'
      do i=1,lb
       write(*,*)'line ',lb
       do j=1,cb
        write(*,*)i,j,B(i,j)
       end do
      end do
      
      do i=1,lc
       do j=1,cc
        write(*,*)'term ',i,j
	app=(0.,0.)
	write(*,*)'starting from ',app
        do k=1,ca
         app=app+(A(i,k)*B(k,j))
	 write(*,*)i,j,k,A(i,k),B(k,j),app
	end do
	C(i,j)=app
       end do
      end do
      
      return
      end



c
c     SUBROUTINES FROM MTINVERS:  
c  
*******************************************************
      SUBROUTINE MATMUL1(L,M,N,A,LA,B,MB,C,LC)
******************************************************      
*     CALCULATES C = A*B
      DOUBLE PRECISION A(LA,M), B(MB,N), C(LC,N)
      INTEGER I,J,K

      DO 1 I = 1,L
      DO 1 J = 1,N
	 C(I,J) = 0.D0
	 DO 2 K = 1,M
	    C(I,J) = C(I,J) + A(I,K)*B(K,J)
 2       CONTINUE
 1    CONTINUE
      RETURN
      END

*******************************************************
      SUBROUTINE MATMUL(L,M,N,A,LA,B,MB,C,LC)
******************************************************      
*     CALCULATES C = A*B
      DOUBLE PRECISION A(LA,M), B(MB,N), C(LC,N)
      INTEGER I,J,K

      DO 1 I = 1,L
      DO 1 J = 1,N
	 C(I,J) = 0.D0
	 DO 2 K = 1,M
	    C(I,J) = C(I,J) + A(I,K)*B(K,J)
 2       CONTINUE
 1    CONTINUE
      RETURN
      END

      FUNCTION AMAX(A,N)
      DIMENSION A(1)
      AMAX=A(1)
      DO 10 I=2,N
      IF(A(I)-AMAX) 10,10,9
9     AMAX=A(I)
10    CONTINUE
      RETURN
      END
 
      FUNCTION AMIN(A,N)
      DIMENSION A(1)
      AMIN=A(1)
      DO 10 I=2,N
      IF(A(I)-AMIN) 9,10,10
9     AMIN=A(I)
10    CONTINUE
      RETURN
      END

      SUBROUTINE FORK(LX,CX,SIGNI)
      COMPLEX*8 CX(LX),CARG,CW,CTEMP
      PI=3.14159265
      J=1
      SC=1./FLOAT(LX)
      SC=SQRT(SC)
      DO 5 I=1,LX
      IF(I.GT.J) GOTO 2
      CTEMP=CX(J)*SC
      CX(J)=CX(I)*SC
      CX(I)=CTEMP
2     M=LX/2
3     IF(J.LE.M) GOTO 5
      J=J-M
      M=M/2
      IF(M.GE.1) GOTO 3
5     J=J+M
      L=1
6     ISTEP=2*L
      DO 8 M=1,L
      CARG=CMPLX(0.,1.)*(PI*SIGNI*FLOAT(M-1))/FLOAT(L)
      CW=CEXP(CARG)
      DO 8 I=M,LX,ISTEP
      CTEMP=CW*CX(I+L)
      CX(I+L)=CX(I)-CTEMP
8     CX(I)=CX(I)+CTEMP
      L=ISTEP
      IF(L.LT.LX) GOTO 6
      RETURN
      END

      subroutine filtkaus(dt,f1,f2,ifenster,lx,fil,ifail)

      ! filterkoeff. fil(lx) (lx=pot von 2)
      ! eines kausalen tiefpasses  (einfach rechteck)
      ! %der mit
      ! %rechteck- (ifenster=1) oder cosinusfenster(ifenster=0) 
      ! %getapert wird (noch nicht eingebaut!!!!!!! ). 
      ! Die laenge des filters wird durch die
      ! laenge des fensters lx*dt bestimmt, dt ist die abtastfrequenz.
      ! f2 (in Hz) ist die Grenzferquenz
      ! fuer tiefpass ist f1<=0
      ! ifail = 0, wenn alles ok.

      implicit none

      integer lx,lxx,ifenster
      real dt,f1,f2,w
      real fil(lx) 
      integer ifail
      real pi
      real fnyq,w1,w2,xf,x1,x2,x,xk
      integer n,k

      lxx = lx/2 +1
c     f1 = 0.
      ifail = 0
      pi = 4.*atan(1.)
      fnyq = 1./(2.*dt)
      w2=f2*2.*pi
      w1=f1*2.*pi

      if(dt .le. 0 .or. f2.le.0)then
       write(*,*)'f oder  dt =',dt,' <= 0!'
       ifail=1
       ! end subroutine filtkaus
       return
      end if

      !-----> Tiefpassfilter im Frequenzbereich

       fil(1) =  w2/pi
      ! fil(lx+1) = fil(1)
       x2 = pi*w2/(2.*pi*fnyq) 

       do k = 1,lxx
        w=2.*pi*real(k-1)/(real(lx)*dt)
        ! if (ifenster .eq. 0) fft eines cos-fensters
        !  noch nicht implementiert                  
        ! else fft eines rechtecks
        if(w.ge.w1 .and. w.le.w2)then
         fil(k) = 1.
        else
         fil(k)=0.
        end if
        fil(lx-k+2) = fil(k)
       end do
       

      return
      end 

      subroutine slip(strike,dip,rake,sv,nv,bv)  
******************************************************************************
*
*     Program calculates slip vector, normal vector and b-axis vector
*     from Strike, dip and Rake of Double Couple
*                                         
****************************************************************************
      IMPLICIT NONE
      real*8 nv(3),sv(3),bv(3),strike,dip,rake
 
*     sv,nv und bv aus strike, dip und slip (aki+richards, 1980)

      sv(1) = dcos(rake)*dcos(strike)
     &        +dcos(dip)*dsin(rake)*dsin(strike)
      sv(2) = dcos(rake)*dsin(strike)
     &        -dcos(dip)*dsin(rake)*dcos(strike)
      sv(3) = -dsin(dip)*dsin(rake)

      nv(1) = -dsin(dip)*dsin(strike)
      nv(2) = +dsin(dip)*dcos(strike)
      nv(3) = -dcos(dip)
 
      bv(1) = sv(2)*nv(3) - sv(3)*nv(2) 
      bv(2) = sv(3)*nv(1) - sv(1)*nv(3)
      bv(3) = sv(1)*nv(2) - sv(2)*nv(1)
c     bv(1) = -bv(1) 
c     bv(2) = -bv(2)
c     bv(3) = -bv(3)

      return
      end

      subroutine searchbl1(cname,ilen,ibl)
*     sucht von vorne die erste .ne. ' ' Stelle (ibl)
*     vorandene blanks mit 0-en . 
      implicit none
      integer ilen,ibl,i
      character*50 cname
      i=0
 1    i=i+1
      if(i.gt.ilen) then
       ibl = ilen-1
       goto 2
      end if
      if(cname(i:i).eq.' ')then
       ibl=i-1
       goto 2
      end if
      goto 1
 2    return
      end
      subroutine comment(lll)
****************************************
*     spult input datei weiter auf erste zeile, die nicht pattern 
*     in 1 spalte enthaelt
****************************************
      integer lll 
      character*1 cpat,c1st
      cpat='%'
    1  read(lll,"(a1)",end=2)c1st
       if(cpat.eq.c1st) goto 1
    2  backspace(lll)
      return
      end



***********************************************************
      subroutine mech1(xm,m,xi,xi1,dc,rm,iso,pt)
***********************************************************
*     changes 19.04.93 iso wird zusaetlich uebergeben
*     xm(1) = 0.5*(M22-M11)
*     xm(2) = M12
*     xm(3) = M13
*     xm(4) = M23
*     xm(5) = -1/3 *(0.5(m(1)+m(3)) - m(6))
*     xm(6) = 1/3 *(m(1)+m(3)+m(6))
*     rm; totale seismische Moment nach Silver und Jordan 1982
*         rm = Tensornorm / sqrt(2)
*     iso relative isotrope Quellanteil nach Silver und Jordan 
*         iso = MI**2 / MT**2 
      implicit none
      integer iso
      real xm(6),m(6),xi(3),xi1(3),pt(6)
      real dc,rm

c     m(1)=-xm(1)+xm(5)+xm(6)
      m(1)=-xm(1)-xm(5)+xm(6)
      m(2)=xm(2)
c     m(3)=+xm(1)+xm(5)+xm(6)
      m(3)=+xm(1)-xm(5)+xm(6)
      m(4)=xm(3)
      m(5)=xm(4)
c     m(6)=xm(6)-2.*xm(5)
      m(6)=xm(6)+2.*xm(5)
      call moment1(m,xi,xi1,pt,dc,rm,iso)
      return
      end
      
      SUBROUTINE MOMENT1(M,XI,XI1,PT,DC,SMT,I1)
*
      IMPLICIT NONE
      INTEGER I,J,K,IACT,IPRT,ISWI,ISLI,IMYCRO,I1,I2,I3 
      INTEGER NMA,NMI,NNU,NROT 
      REAL A1(3),A2(3),A3(3),A11(3,3),A22(3,3),A33(3,3),EV(3)
      REAL Z(3,3),Z1(3,3),Z2(3,3),XI(3),ZETA(3),ETA(3),XI1(3),PT(6)
      REAL NV(3),SV(3),BV(3),P(3),T(3),SIG,DEL,GAM,SIG1,DEL1,GAM1
      REAL DELTAI(3),DELTAD(3),DELTAC(3),DELTAV(3),DELTA(3)
      REAL M(6),MCLVD(3,3),MDC(3,3),MTOT(3,3),V(3,3)
      REAL MDC1(3,3),MDC2(3,3),MDC3(3,3)
      REAL RM0,TOL,CODC,COCLVD,XF,FM,FMIN,FMAX,VORZ1,VORZ2
      REAL TRACE,FACT,EMIN,EMAX,DC,CLVD,CODC1,CODC2,CODC3,X1,X2,X3 
      REAL DDUN(3)
      REAL CNI,DELTD1,DELTD2,DELTD3,SMT,SMI,SMD,SMDC
      CHARACTER*1 IC(6),IC2(3)
      CHARACTER*6 LAB1(3)
      CHARACTER*12 LAB2(3)
*
      EQUIVALENCE (IC2,IC(4))
*
      DATA IC/'N','S','B','P','T','N'/
      DATA LAB1,LAB2/'STRIKE',' DIP  ',' RAKE ','COMPRESSION',
     & '  TENSION  ','   NULL    '/,TOL/1.E-3/
*
      DATA CNI/57.295780/
      DATA IPRT/0/
      DATA IACT/0/
      DATA ISWI/1/
      DATA ISLI/0/
      DATA IMYCRO/0/
*
      CALL MTOFI(MTOT,M)
*
* .....Berechnen der Eigenwerte des Momententesors
*
      CALL JACOBI(3,MTOT,DELTA,V,NROT)

*
* ......Berechnen der Eigenwerte des isoropen und des deviatorischen MT
*
      TRACE = DELTA(1) + DELTA(2) + DELTA(3)
      IF(ABS(TRACE).LE.TOL) TRACE = 0.
      DO 11 I=1,3
        DELTAI(I) = TRACE/3.
        DELTAD(I) = DELTA(I)-TRACE/3.
	IF(ABS(DELTAD(I)).LT.0.1*TOL) DELTAD(I) = 0.
	DDUN(I) = DELTAD(I)
 11   CONTINUE
      DELTD1=DELTAD(1)
      DELTD2=DELTAD(2)
      DELTD3=DELTAD(3)
* 
* ....Finden des max. und min deviatorischen Eigenvektors, umsortieren
*
      EMAX = ABS(DELTAD(1))
      EMIN = EMAX
      NMI = 1
      NMA = 1
      DO 2 K = 2,3
       IF(ABS(DELTAD(K)).GE.EMIN) GOTO 3
       EMIN = ABS(DELTAD(K))
       NMI = K
 3     IF(ABS(DELTAD(K)).LE.EMAX) GOTO 2
       EMAX = ABS(DELTAD(K))
       NMA = K
 2    CONTINUE
      NNU = 1
 4    IF(NNU.EQ.NMI .OR. NNU.EQ.NMA) THEN
	NNU = NNU + 1
	GO TO 4
      END IF
      EV(1) = DELTA(NMI)
      EV(2) = DELTA(NNU)
      EV(3) = DELTA(NMA)
      DO I=1,3
       DELTA(I)=EV(I)
      END DO
      EV(1) = DELTAD(NMI)
      EV(2) = DELTAD(NNU)
      EV(3) = DELTAD(NMA)
      DO  I=1,3
       DELTAD(I)=EV(I)
      END DO

*     Seismic Moment after Silver + Jordan

      SMT = (DELTA(1)**2+DELTA(2)**2+DELTA(3)**2)/2.
      SMI = TRACE**2/6.
      SMD = (DELTAD(1)**2+DELTAD(2)**2+DELTAD(3)**2)/2.
      FM = -DELTAD(1)/DELTAD(3)
      SMDC= ((1.-2.*FM)*DELTAD(3))**2
 
* ....Berechnen der Dyaden aus den Eigenvektoren
 
      DO 12 I = 1,3
	A1(I) = V(I,NMI)
	A2(I) = V(I,NNU)
 12     A3(I) = V(I,NMA)
      CALL DIADE(A1,A1,A11,3)
      CALL DIADE(A2,A2,A22,3)
      CALL DIADE(A3,A3,A33,3)
*
* .....3 elemtare Double Coubles
*
      CALL MATADD(1.,A11,-1.,A22,MDC1,3)
      CALL MATADD(1.,A22,-1.,A33,MDC2,3)
      CALL MATADD(1.,A33,-1.,A11,MDC3,3)
      CODC1 = (DELTA(1) - DELTA(2))/3.
      CODC2 = (DELTA(2) - DELTA(3))/3.
      CODC3 = (DELTA(3) - DELTA(1))/3.
      X1 = 1.
      X2 = 1.
      X3 = 1.
      IF(CODC1.LE.0.)THEN
	 CODC1=-1.*CODC1
	 X1 = -1.
      END IF
      IF(CODC2.LE.0.) THEN
	 CODC2 = -1.*CODC2
	 X2 = -1.
      END IF
      IF(CODC3.LE.0.) THEN
	 CODC3 = -1.*CODC3
	 X3 = -1.
      END IF
      DO 19 I=1,3
      DO 19 J = 1,3
        MDC1(I,J) = X1 * MDC1(I,J) 
        MDC2(I,J) = X2 * MDC2(I,J) 
 19     MDC3(I,J) = X3 * MDC3(I,J) 
*
* .....CLVD - ANTEIL + DOUBLE COUPLE
*
      CALL MATADD(2.,A33,-1.,A22,MCLVD,3)
      CALL MATADD(1.,MCLVD,-1.,A11,MCLVD,3)
      CALL MATADD(1.,A33,-1.,A22,MDC,3)
*
*     ........Prozentuale Anteil von DC%  und CLVD%
*
      FM=0.
      FMIN=AMIN1(ABS(DELTAD(1)),ABS(DELTAD(2)),ABS(DELTAD(3)))
      FMAX=AMAX1(ABS(DELTAD(1)),ABS(DELTAD(2)),ABS(DELTAD(3)))
      FM=FMIN/FMAX
      DC=100.*(1.-FM*2.)
c     CLVD=100.*FM*2.
*     ... den Wert dc durch 100 * -e2/max(|e1|,|e3|) belegen
*         sieher Definition Kawakatsu
      fmin = DELTAD(1)
      if(abs(deltad(2)).lt.abs(fmin))fmin = deltad(2)
      if(abs(deltad(3)).lt.abs(fmin))fmin = deltad(3)
      fmin = -fmin
*     ... fmax wie gehabt
c     dc=100.*fmin/fmax
*
* ....Berechnen der Eigenvectoren und seism.  Moment des major DC
*     Hier wird ausgenutzt, dass die Eigenvektoren des totalen 
*     Momententensors und des Major Double Anteils gleich sind
*     Im Grunde doppelt gemoppelt
*
      CALL FGMRTI(SV,NV,BV,RM0,MTOT)
*
*    Nochmal P und T aus SV und NV, da in FGMRTI nicht uebergeben (Z2=P)
*
      CALL FGPAF(SV,NV,P,T)
      DO 22 I=1,3
        Z2(I,1) = P(I)
        Z2(I,2) = T(I)
 22     Z2(I,3) = BV(I)
*
* ....Azimuth (ZETA) und Plunge (ETA) aus slip, normal und T-vector
*
      CALL PAI(ZETA,ETA,SV,NV,BV)
      PT(1)=ZETA(1)
      PT(2)=ETA(1)
      PT(3)=ZETA(2)
      PT(4)=ETA(2)
      PT(5)=ZETA(3)
      PT(6)=ETA(3)
*
* ... Strike, Dip und Rake (XI=SIG) (Z=NV) aus slip, normal und T
* 
      CALL FGI(SIG,DEL,GAM,SV,NV,BV)
      XI(1) = SIG
      XI(2) = DEL
      XI(3) = GAM
      DO 14 I = 1,3
        Z(I,1) = NV(I)
        Z(I,2) = SV(I)
 14     Z(I,3) = BV(I)
* ......auxiliary Plane: SV und NV vertauschen
      CALL FGI(SIG1,DEL1,GAM1,NV,SV,BV)
      XI1(1) = SIG1
      XI1(2) = DEL1
      XI1(3) = GAM1
      DO 15 I = 1,3
        Z1(I,1) = NV(I)
        Z1(I,2) = SV(I)
 15     Z1(I,3) = BV(I)
      
*
* ....................Berechnen des DOUBLE COUPLE MCD aus Eigenvectoren
*                     ist schon mit Dyaden berechnet worden
C     CALL FGMRTF(SV,NV,BV,1.,MDC)
*
*       .........Eigenwerte von DC
      CALL JACOBI(3,MDC,DELTAC,V,NROT)
*
*       .........Eigenwrte von CLVD
      CALL JACOBI(3,MCLVD,DELTAV,V,NROT)
*
* .........Normieren von MDC ind MCLVD mit RM0
*
      FACT = -1.*DELTAD(1)/( DELTAD(1) + DELTAD(3) )
      XF = -1.*DELTAD(1)/DELTAD(3)
      CODC = DELTAD(3)*(1.-2.*XF)
      COCLVD = -1.*DELTAD(1)
      VORZ1 = 1.
      VORZ2 = 1.
      IF(CODC.LT.0.) VORZ1 = -1.
      IF(COCLVD.LT.0.) VORZ2 = -1.
      CODC = VORZ1*CODC
      COCLVD = VORZ2*COCLVD
      IF(ABS(CODC).LE.TOL) VORZ1=0.0
      IF(ABS(COCLVD).LE.TOL) VORZ2=0.0
      DO 16 I=1,3
      DO 17 J=1,3
	MDC(I,J) = VORZ1 * MDC(I,J)
C       MCLVD(I,J) = VORZ2*FACT*MCLVD(I,J)
        MCLVD(I,J) = VORZ2*MCLVD(I,J)
        IF(ABS(MTOT(J,I)).LE.TOL) MTOT(J,I)=0.
        IF(ABS(MDC(J,I)).LE.TOL) MDC(J,I)=0.
        IF(ABS(MCLVD(J,I)).LE.TOL) MCLVD(J,I)=0.
 17   CONTINUE
 16   CONTINUE
 
*     ............Seismic Moment in  (dyne-cm) (if data in counts*s/ym , 
*                                               IMYCRO.eq.1)
      IF(SMI.ge.SMT)THEN
*       .......... Error ?
	GOTO 163
      ELSE IF(SMI.lt.0.001*SMT) THEN
	SMI=0.
      END IF
      IF(SMT.gt. 1.E-12) THEN
       I1=INT(100.*SMI/SMT)
c      I2=INT(100.*SMD/SMT)
c      I3=INT(100.*SMDC/SMT)
      ELSE
       I1=100*1.e+20*(SMI/(1.e+20*SMT))
c      I2=-1
c      I3=-1
      END IF
c     ISO durch TRACE
      I1 = sign(1.,TRACE)*I1
      
 163  CONTINUE
c     IF(IMYCRO.EQ.1) THEN
c       SMT=SQRT(SMT)*1.E14
c       SMI=SQRT(SMI)*1.E14
c       SMD=SQRT(SMD)*1.E14
c       SMDC=SQRT(SMDC)*1.E14
c       RM0=RM0*1.E14
c      ELSE
c       SMT=SQRT(SMT)*1.E20
c       SMI=SQRT(SMI)*1.E20
c       SMD=SQRT(SMD)*1.E20
c       SMDC=SQRT(SMDC)*1.E20
c        RM0=RM0*1.E20
c      END IF
      SMT=SQRT(SMT)
 
      RETURN
      END
      
      SUBROUTINE JACOBI(N,A,D,V,NROT)
C
C $$$$$ CALLS ONLY LIBRARY ROUTINES $$$$$
C
C   JACOBI FINDS ALL EIGENVALUES AND EIGENVECTORS OF THE N BY N REAL
C   SYMMETRIC MATRIX A BY JACOBI'S METHOD.  N MUST BE LESS THAN OR
C   EQUAL TO 50 (DETERMINED BY THE DIMENSIONS OF B AND Z BELOW).
C   THE UNORDERED EIGENVALUES ARE RETURNED IN REAL N VECTOR D AND
C   THE CORRESPONDING EIGENVECTORS ARE IN THE RESPECTIVE COLUMNS OF
C   N BY N REAL MATRIX V.  THE NUMBER OF JACOBI ROTATIONS NEEDED ARE
C   RETURNED IN NROT.  TAKEN FROM LINEAR ALGEBRA, VOLUME II, BY
C   WILKINSON AND REINSCH, 1971, SPRINGER-VERLAG.  THE ALGORITHM
C   IS PUBLISHED (IN ALGOL) BY H. RUTISHAUSER ON PAGES 202-211.
C   CONVERTED TO FORTRAN BY R. BULAND, 30 OCT. 1978.
C
      INTEGER P,Q,P1,P2,Q1,Q2
      DIMENSION A(N,N),D(N),V(N,N),B(50),Z(50)
      N1=N-1
      NN=N*N
C   INITIALIZE ARRAY STORAGE.
      DO 1 P=1,N
      DO 2 Q=1,N
 2    V(P,Q)=0.
      V(P,P)=1.
      B(P)=A(P,P)
      D(P)=B(P)
 1    Z(P)=0.
      NROT=0
C
C   MAKE UP TO 50 PASSES ROTATING EACH OFF DIAGONAL ELEMENT.
C
      DO 3 I=1,50
      SM=0.
      DO 4 P=1,N1
      P1=P+1
      DO 4 Q=P1,N
 4    SM=SM+ABS(A(P,Q))
C   EXIT IF ALL OFF DIAGONAL ELEMENTS HAVE UNDERFLOWED.
      IF(SM.EQ.0.) GO TO 13
      TRESH=0.
      IF(I.LT.4) TRESH=.2*SM/NN
C
C   LOOP OVER EACH OFF DIAGONAL ELEMENT.
C
      DO 5 P=1,N1
      P1=P+1
      P2=P-1
      DO 5 Q=P1,N
      Q1=Q+1
      Q2=Q-1
      G=100.*ABS(A(P,Q))
C   SKIP THIS ELEMENT IF IT HAS ALREADY UNDERFLOWED.
      IF((I.LE.4).OR.(ABS(D(P))+G.NE.ABS(D(P))).OR.
     1 (ABS(D(Q))+G.NE.ABS(D(Q)))) GO TO 6
      A(P,Q)=0.
      GO TO 5
 6    IF(ABS(A(P,Q)).LE.TRESH) GO TO 5
C   COMPUTE THE ROTATION.
      H=D(Q)-D(P)
      IF(ABS(H)+G.EQ.ABS(H)) GO TO 7
      THETA=.5*H/A(P,Q)
      T=1./(ABS(THETA)+SQRT(1.+THETA*THETA))
      IF(THETA.LT.0.) T=-T
      GO TO 14
 7    T=A(P,Q)/H
 14   C=1./SQRT(1.+T*T)
      S=T*C
      TAU=S/(1.+C)
C   ROTATE THE DIAGONAL.
      H=T*A(P,Q)
      Z(P)=Z(P)-H
      Z(Q)=Z(Q)+H
      D(P)=D(P)-H
      D(Q)=D(Q)+H
      A(P,Q)=0.
C   ROTATE THE OFF DIAGONAL.  NOTE THAT ONLY THE UPPER DIAGONAL
C   ELEMENTS ARE TOUCHED.  THIS ALLOWS THE RECOVERY OF MATRIX A LATER.
      IF(P2.LT.1) GO TO 15
      DO 8 J=1,P2
      G=A(J,P)
      H=A(J,Q)
      A(J,P)=G-S*(H+G*TAU)
 8    A(J,Q)=H+S*(G-H*TAU)
 15   IF(Q2.LT.P1) GO TO 16
      DO 9 J=P1,Q2
      G=A(P,J)
      H=A(J,Q)
      A(P,J)=G-S*(H+G*TAU)
 9    A(J,Q)=H+S*(G-H*TAU)
 16   IF(N.LT.Q1) GO TO 17
      DO 10 J=Q1,N
      G=A(P,J)
      H=A(Q,J)
      A(P,J)=G-S*(H+G*TAU)
 10   A(Q,J)=H+S*(G-H*TAU)
C   ROTATE THE EIGENVECTOR MATRIX.
 17   DO 11 J=1,N
      G=V(J,P)
      H=V(J,Q)
      V(J,P)=G-S*(H+G*TAU)
 11   V(J,Q)=H+S*(G-H*TAU)
      NROT=NROT+1
 5    CONTINUE
C   RESET THE TEMPORARY STORAGE FOR THE NEXT ROTATION PASS.
      DO 3 P=1,N
      D(P)=B(P)+Z(P)
      B(P)=D(P)
 3    Z(P)=0.
C
C   ALL FINISHED.  PREPARE FOR EXITING BY RECONSTRUCTING THE UPPER
C   TRIANGLE OF A FROM THE UNTOUCHED LOWER TRIANGLE.
C
 13   DO 12 P=1,N1
      P1=P+1
      DO 12 Q=P1,N
 12   A(P,Q)=A(Q,P)
      RETURN
      END
C    --- INCLUDE CROSS,FIXED,INSOURCE
      SUBROUTINE CROSS(A,B,C)
C
C $$$$$ CALLS NO OTHER ROUTINE $$$$$
C
C   CROSS COMPUTES THE VECTOR PRODUCT C = A X B, WHERE A, B, AND C ARE
C   ALL THREE VECTORS.  PROGRAMMED ON 30 OCT 78 BY R. BULAND.
C
      DIMENSION A(3),B(3),C(3)
      C(1)=A(2)*B(3)-B(2)*A(3)
      C(2)=A(3)*B(1)-B(3)*A(1)
      C(3)=A(1)*B(2)-B(1)*A(2)
      RETURN
      END
C    --- INCLUDE MTOFI,FIXED,INSOURCE
      SUBROUTINE MTOFI(M,G)
C
C $$$$$ CALLS NO OTHER ROUTINE $$$$$
C
      REAL M(3,3),G(3)
      K=0
      DO 2 J=1,3
      DO 2 L=1,J
      K=K+1
      M(L,J)=G(K)
 2    M(J,L)=G(K)
      RETURN
      END
C    --- INCLUDE FGMRTIXY,FIXED,INSOURCE
      SUBROUTINE FGMRTI(S,N,B,M0,M)
C     VERSION FOR X,Y,Z-(NORTH-EAST-DOWN)-COORDINATES
C
C $$$$$ CALLS CROSS, FGPAF, FGPAI, AND JACOBI $$$$$
C
      REAL S(3),N(3),B(3),M0,M(3,3),P(3),T(3),EV(3),VEC(3,3)
C ??? CALL FGPAF(S,N,P,T)
C   DECOMPOSE THE MOMENT TENSOR.
      CALL JACOBI(3,M,EV,VEC,NROT)
C   FIND THE MAXIMUM AND MINIMUM EIGENVALUES.
C   The Eigenvectors of the deviatoric Momenttensor are equivalent to the 
C   Eigenvectors of the composite Momenttensor
      EMAX=EV(1)
      EMIN=EMAX
      NP=1
      NT=1
      DO 2 K=2,3
      IF(EV(K).GE.EMIN) GO TO 3
      EMIN=EV(K)
      NP=K
 3    IF(EV(K).LE.EMAX) GO TO 2
      EMAX=EV(K)
      NT=K
 2    CONTINUE
C   COMPUTE S AND N.
      SP=SIGN(1.,VEC(3,NP))
      ST=SIGN(1.,VEC(3,NT))
      DO 4 K=1,3
      P(K)=-SP*VEC(K,NP)
 4    T(K)=-ST*VEC(K,NT)
      CALL CROSS(P,T,B)
* .......Vgl Jost + Herrmann (19)
C     M0=.5*(EV(NT)-EV(NP))
C     M0=.5*( ABS(EV(NT)) + ABS(EV(NP)) )
      M0=.5*( ABS(EMAX) + ABS(EMIN) )
C     PRINT*,' M0 = ',M0
      CALL FGPAI(S,N,P,T)
      RETURN
      END
C    --- INCLUDE PAIXY,FIXED,INSOURCE
      SUBROUTINE PAI(ZETA,ETA,S,N,B)
C     VERSION FOR X,Y,Z-COORDINATES
C
C $$$$$ CALLS FGPAF $$$$$
C
C $$$$$ THIS ROUTINE IS DEFINED IN A (X,Y,Z) COORDINATE SYSTEM WITH
C       ORIENTATION (NORTH,EAST,DOWN) OTHER THAN THE PREVIOUS COOR-
C       DINATES WITH (UP,SOUTH,EAST). SEE  AKI &RICHARDS P.119FF
C       EXTENDED ON JUNE 1, 1985 BY K.KOCH
C
      REAL ZETA(3),ETA(3),S(3),N(3),B(3),P(3),T(3)
      DATA CNI,TOL,TOL2/57.295780,1E-5,2E-5/
      CALL FGPAF(S,N,P,T)
      DO 2 K=1,3
 2    ZETA(K)=0.
      CE= SIGN(1.,P(3))
      IF(ABS(P(1))+ABS(P(2)).GT.TOL2) ZETA(1)=CNI*ATAN2(CE*P(2),
     1  CE*P(1))
      ETA(1)=CNI*ASIN(CE*P(3))
      CE= SIGN(1.,T(3))
      IF(ABS(T(1))+ABS(T(2)).GT.TOL2) ZETA(2)=CNI*ATAN2(CE*T(2),
     1  CE*T(1))
      ETA(2)=CNI*ASIN(CE*T(3))
      CE= SIGN(1.,B(3))
      IF(ABS(B(1))+ABS(B(2)).GT.TOL2) ZETA(3)=CNI*ATAN2(CE*B(2),
     1  CE*B(1))
      ETA(3)=CNI*ASIN(CE*B(3))
      DO 1 K=1,3
      IF(ZETA(K).LT.0.) ZETA(K)=ZETA(K)+360.
 1    CONTINUE
C   PUT ZETA INTO CANONICAL FORM IF ETA IS NEARLY ZERO.
      IF((ABS(P(3)).LE.TOL).AND.(ZETA(1).GT.180.)) ZETA(1)=ZETA(1)-180.
      IF((ABS(T(3)).LE.TOL).AND.(ZETA(2).GT.180.)) ZETA(2)=ZETA(2)-180.
      IF((ABS(B(3)).LE.TOL).AND.(ZETA(3).GT.180.)) ZETA(3)=ZETA(3)-180.
      RETURN
      END
C    --- INCLUDE FGIXY,FIXED,INSOURCE
      SUBROUTINE FGI(SIG,DEL,GAM,S,N,B)
C     VERSION FOR X,Y,Z-COORDINATES
C
C $$$$$ CALLS ONLY SYSTEM ROUTINES $$$$$
C
C $$$$$ THIS ROUTINE IS DEFINED IN A (X,Y,Z) COORDINATE SYSTEM WITH
C       ORIENTATION (NORTH,EAST,DOWN) OTHER THAN THE PREVIOUS COOR-
C       DINATES WITH (UP,SOUTH,EAST). SEE  AKI &RICHARDS P.115FF
C       EXTENDED ON JUNE 1, 1985 BY K.KOCH
C
      REAL S(3),N(3),B(3)
      DATA CNI/57.295780/
      DEL=ACOS(-N(3))
      SIG=ATAN2(-N(1),N(2))
      GAM=CNI*ATAN2(-S(3)/SIN(DEL),S(2)*SIN(SIG)+S(1)*COS(SIG))
      DEL=CNI*DEL
      SIG=CNI*SIG
      IF(DEL.GT.90.) THEN
        SIG=SIG-180.
        DEL=180.-DEL
        GAM=-GAM
* .....Hier muss rake gedreht werden,da man immer dieselbe Seite der
*      dipping fault betrachtet
        PRINT*,'RAKE UMGEDREHT'
      END IF
      IF(SIG.LT.0.) SIG=SIG+360.
C     IF(GAM.LT.0.) GAM=GAM+360.
      IF(GAM.GT.180.) GAM=GAM-360.
      IF(GAM.LE.-180.) GAM=GAM+360.
      RETURN
      END
C    --- INCLUDE FGPAFXY,FIXED,INSOURCE
      SUBROUTINE FGPAF(S,N,P,T)
C     VERSION FOR X,Y,Z-(NORTH-EAST-DOWN)-COORDINATES
C
C $$$$$ CALLS NO OTHER ROUTINE $$$$$
C
C   THERE ARE TWO FGPAXY ROUTINES.  FGPAFXY CONVERTS THE UNIT VECTORS
C   DEFINING THE FAULT GEOMETRY COORDINATE SYSTEM (SEE SUBROUTINE FGFXY)
C   TO THE UNIT VECTORS DEFINING THE PRINCIPAL AXIS COORDINATE SYSTEM:
C   P, T, AND B.  THESE UNIT VECTORS ARE COLINEAR WITH THE P, T, AND B
C   AXES SUCH THAT P(3)<=0 AND T = B X P (WHERE B = N X S).  B IS NOT
C   NEEDED AS IT IS COMMON TO BOTH SYSTEMS.  FGPAIXY PERFORMS THE IN-
C   VERSE TRANSFORMATION.  PROGRAMMED ON 30 OCT 78 BY R. BULAND.
C   UPDATED ON 11 SEP 1986 BY K.KOCH.
C
      REAL S(3),N(3),P(3),T(3)
      DATA CON/.70710678/
      CN=CON
      IF(N(3)-S(3).GT.0.) CN=-CN
      DO 1 K=1,3
      P(K)=CN*(N(K)-S(K))
 1    T(K)=CN*(N(K)+S(K))
      RETURN
      END
C    --- INCLUDE FGPAIXY,FIXED,INSOURCE
      SUBROUTINE FGPAI(S,N,P,T)
C     VERSION FOR X,Y,Z-(NORTH-EAST-DOWN)-COORDINATES
C
C $$$$$ CALLS NO OTHER ROUTINE $$$$$
C
      REAL S(3),N(3),P(3),T(3)
      DATA CON/.70710678/
      CN=CON
      IF(T(3)+P(3).GT.0.) CN=-CN
      DO 2 K=1,3
      N(K)=CN*(T(K)+P(K))
 2    S(K)=CN*(T(K)-P(K))
      RETURN
      END
C    --- INCLUDE FGMRTF,FIXED,INSOURCE
      SUBROUTINE FGMRTF(S,N,B,M0,M)
C
C $$$$$ CALLS NO OTHER ROUTINE $$$$$
C
C   THERE ARE TWO FGMRT ROUTINES.  FGMRTF CONVERTS THE UNIT VECTORS S
C   AND N (SEE SUBROUTINE FGF) AND THE MOMENT M0 TO A MOMENT TENSOR M.
C   FGMRTI PERFORMS THE INVERSE OPERATION.  PROGRAMMED ON 30 OCT 78 BY
C   R. BULAND.
C
      REAL S(3),N(3),B(3),M0,M(3,3)
      DO 1 K=1,3
      DO 1 J=1,3
 1    M(J,K)=M0*(S(K)*N(J)+S(J)*N(K))
      RETURN
      END
C
      SUBROUTINE MATADD(CA,A,CB,B,C,N)
C  CALCULATES C = aA + bB , C = C(N,N), B=B(N,N), A=A(N,N)
C
      INTEGER N,I,J
      REAL A(N,N), B(N,N), C(N,N),CA,CB
      DO 1 I = 1,N
      DO 1 J = 1,N
 1    C(I,J) = CA*A(I,J) + CB*B(I,J)
      RETURN 
      END
C
      SUBROUTINE DIADE(A,B,AB,N)
C  CALCULATES FROM VECTOR A(N) AND B(N) THE DYADIC AB 
C  LITERATURE: JOST + HERMANN SEISM RESEARCH LETTERS 60, 1989
C
      INTEGER I,J,N
      REAL A(N), B(N), AB(N,N)
      DO 1 I = 1,N
      DO 1 J = 1,N
 1    AB(I,J) = A(I) * B(J)
      RETURN
      END
       
      REAL FUNCTION RD(X)
      REAL X,STELLE
      INTEGER ISTELLE
      DATA ISTELLE / 3 /
*
      STELLE = 10.**ISTELLE
      RD = INT( STELLE * X) / STELLE
      RETURN
      END

      SUBROUTINE SSLIP(SIG,XFACT,XM)
 
*     .......SIG = Azimuth in degree
      IMPLICIT NONE
      INTEGER I,J
      REAL SIG,SIGR,XFACT,XM(3,3),CNI
      DATA CNI/57.295780/
      SIGR = SIG/CNI
      XM(1,1) = -1.*SIN(2*SIGR)
      XM(2,2) = SIN(2*SIGR)
      XM(3,3) = 0.
      XM(1,2) = COS(2.*SIGR)
      XM(2,1) = XM(1,2)
      XM(1,3) = 0.
      XM(3,1) = XM(1,3)
      XM(2,3) = 0.
      XM(3,2) = XM(2,3)
      DO 1 I = 1,3
      DO 1 J = 1,3
 1    XM(I,J) = XFACT*XM(I,J)
      RETURN
      END
       
      SUBROUTINE ISLIP(SIG,XFACT,XM)
 
*     .......SIG = Azimuth in degree
      IMPLICIT NONE
      INTEGER I,J
      REAL SIG,SIGR,XFACT,XM(3,3),CNI
      DATA CNI/57.295780/
      SIGR = SIG/CNI
      XM(1,1) = -1.*SIN(SIGR)**2
      XM(2,2) = -1.*COS(SIGR)**2
      XM(3,3) = 1.
      XM(1,2) = SIN(2.*SIGR)/2.
      XM(2,1) = XM(1,2)
      XM(1,3) = 0.
      XM(3,1) = XM(1,3)
      XM(2,3) = 0.
      XM(3,2) = XM(2,3)
      DO 1 I = 1,3
      DO 1 J = 1,3
 1    XM(I,J) = XFACT*XM(I,J)
      RETURN
      END
       
      SUBROUTINE DSLIP(SIG,XFACT,XM)
 
*     .......SIG = Azimuth in degree
      IMPLICIT NONE
      INTEGER I,J
      REAL SIG,SIGR,XFACT,XM(3,3),CNI
      DATA CNI/57.295780/
      SIGR = SIG/CNI
      XM(1,1) = 0.
      XM(2,2) = 0.
      XM(3,3) = 0.
      XM(1,2) = 0. 
      XM(2,1) = XM(1,2)
      XM(1,3) = SIN(SIGR)
      XM(3,1) = XM(1,3)
      XM(2,3) = -COS(SIGR)
      XM(3,2) = XM(2,3)
      DO 1 I = 1,3
      DO 1 J = 1,3
 1    XM(I,J) = XFACT*XM(I,J)
      RETURN
      END
       
      SUBROUTINE VECDIPO(XFACT,XM)
 
      IMPLICIT NONE
      INTEGER I,J
      REAL XFACT,XM(3,3)
      XM(1,1) = -1./3.
      XM(2,2) = -1./3.
      XM(3,3) = 2./3.
      XM(1,2) = 0. 
      XM(2,1) = XM(1,2)
      XM(1,3) = 0.
      XM(3,1) = XM(1,3)
      XM(2,3) = 0.
      XM(3,2) = XM(2,3)
      DO 1 I = 1,3
      DO 1 J = 1,3
 1    XM(I,J) = XFACT*XM(I,J)
      RETURN
      END
       
      subroutine mdekomp(rdp,rsl,re,rc,iswi)
*******************************************************************************
*     Function:  Program rotates a Tensor given in the mean
*                axis system (3 Eigenvalues) with  
*                STRIKE, DIP and SLIP defined for the  fault plane 
*                of a double couble  (in degree)
*     INPUT: 1) STRIKE, DIP and SLIP
*            2) Three eigenvalues of a Moment-tensor
*     OUTPUT:   Coefficients (RC(6) of rotatet elementary  Momenttensors
*     
*     Language : Fortran
*
*     Used Parameters:
*                       dip (f)   : dip
*                       rake (f)  : rake (all in degree)
*                       e(3)      : Eigenvalues
*                              abs(ev3) > abs(ev1) > abs(ev2) !!
*                       mij       : Moment-tensor elements
*                       c(6)      : Coefficients of elementary sources
*                       iswi      : switch
*
*     Author: Torsten Dahm  7.1991
*
* References: ..................
****************************************************************************** 
      implicit none
      integer  i,iswi
      real rdp,rsl,re(3),rc(6)
      double precision ev(3),c(6)
      double precision y1,y2
      double precision pi
      double precision dip,rake 
      double precision xdc,xclvd

      pi = 4.d0*atan(1.)
      do 11 i=1,3
 11     ev(i) = re(i)
      dip = rdp
      rake = rsl
      if(ev(3).gt.0) rake=rake-180.
      dip = pi*dip/180.
      rake = pi*rake/180.

*     Berechnen der Vorfaktoren

      xdc =  -ev(3) - 2.d0*ev(2)
      xclvd =  3.d0*ev(2)

*     festlegen der Coeffizienten

      y1=-dcos(2.*dip)*dsin(rake)
      y2=-0.5*dcos(2.*dip)*dsin(rake) 
     .                   -0.25*dsin(2.*dip)*dcos(rake)**2
      c(1) = xdc*y1 + xclvd*y2

      y1=-dcos(dip)*dcos(rake)
      y2=-0.5*dcos(dip)*dcos(rake)+0.25*dsin(dip)*dsin(2.*rake)
      c(2) = xdc*y1 + xclvd*y2

      y1=dsin(dip)*dcos(rake)
      y2=0.5*dsin(dip)*dcos(rake)+0.25*dcos(dip)*dsin(2.*rake)
      c(3) = xdc*y1 + xclvd*y2

      y1=dsin(2.*dip)*dsin(rake)
      y2=+0.5*(dsin(dip) + dcos(dip)*dsin(rake))**2
      c(4) = xdc*y1 + xclvd*y2

      c(5)=+0.5*xclvd*dcos(rake)**2

      c(6) = -1.d0*xclvd

      do 18 i=1,6
 18     rc(i)=real(c(i))
      
      return
      end

c------------------------------------------------------------
      subroutine init_hdrvars_mtinvers
c------------------------------------------------------------
      ! INTEGER 
      !  hdr_ifmtnum        ! input format number
      !  hdr_ofmtnum        ! output format number
      !  hdr_idiff          ! number of input differences
      !  hdr_odiff          ! number of output differences
      !  hdr_calunit        ! unit of calibration (0=displ, 1=vel, 2=acc)
      ! REAL    
      !  hdr_smprate        ! sample rate in Hz
      !  hdr_calfac         ! calibration factor
      !  hdr_calper         ! calibration period
      !  hdr_beamaz         ! beam azimuth
      !  hdr_beamslo        ! beam slowness
      !  hdr_hang           ! horizontal orientation
      !  hdr_vang           ! vertical orientation
      ! CHARACTER
      !  hdr_station*10     ! station name
      !  hdr_stadescr*10    ! station/chan description
      !  hdr_chan*3         ! channel name
      !  hdr_instr*6        ! instrument type
      !  hdr_network*9      ! network name
      !  hdr_coosys*12      ! name of coordinate system
      !  hdr_ifmtname*10    ! name of input format
      !  hdr_ofmtname*10    ! name of output format

      include 'codeco_common.f'

      !  hdr_debug = 4      ! debug level
      hdr_ifmtname = 'ascii'
      hdr_ifmtnum  = 0
      hdr_ofmtname = 'CM6' 
      hdr_ofmtnum  = 2   ! gse2 output format
      hdr_odiff    = 2   ! 2fache differenzen
      hdr_idiff    = 0

      hdr_year     = 2000
      hdr_month    = 1
      hdr_day      = 1
      hdr_jday     = 1
      hdr_hour     = 0
      hdr_min      = 0
      hdr_sec      = 0
      hdr_msec     = 0
      hdr_nsamp    = 0
      hdr_ifmtnum  = 0
      hdr_calunit  = 0
      hdr_smprate  = 0.0
      hdr_calfac   = 1.0
      hdr_calper   = 1.0
      hdr_stalat   = -999.0
      hdr_stalon   = -999.0
      hdr_staelev  = -999.0
      hdr_stadepth = -999.0
      hdr_beamaz   = -1.0
      hdr_beamslo  = -1.0
      hdr_hang     = -1.0
      hdr_vang     = -1.0
      hdr_station  = ' '
      hdr_stadescr = 'mtinvers'
      hdr_chan     = '   '
      hdr_instr    = 'unkown'
      hdr_network  = 'unknown'
      hdr_coosys   = 'cartesian'

      return
      end

c---------------------------------------------------
c adopted from H. Dufumier ? 
c---------------------------------------------------
      subroutine add_time_shift(sec) 
      implicit none
      include 'codeco_common.f'
      real sec

      hdr_sec=int(sec)
      hdr_msec=1000.*(sec-float(hdr_sec))
      if (hdr_msec .lt.0) then
           hdr_sec=hdr_sec -1
c --- just in case you worry about the 999: the float to integer
c     operation from above cuts for negative values towards 0
          hdr_msec=999+hdr_msec
      endif
      hdr_min=hdr_min+int(float(hdr_sec)/60.)
      hdr_sec=mod(hdr_sec,60)
      if (hdr_sec .lt.0) then
           hdr_min=hdr_min-1
           hdr_sec=60+hdr_sec
      endif
      hdr_hour=hdr_hour+int(float(hdr_min)/60.)
      hdr_min=mod(hdr_min,60)
      if (hdr_min .lt.0) then
           hdr_hour=hdr_hour-1
           hdr_min=60+hdr_min
      endif
      hdr_jday=hdr_jday+int(float(hdr_hour)/24.)
      hdr_hour=mod(hdr_hour,24)
      if (hdr_hour .lt.0) then
           hdr_jday=hdr_jday-1
           hdr_hour=24+hdr_hour
      endif
      if (hdr_jday.lt.1) then
       hdr_year=hdr_year-1
       call julien(31,12,hdr_year,hdr_jday,1)
      else
       call julien(hdr_day,hdr_month,hdr_year,hdr_jday,0)
      endif
      if  (hdr_debug .gt. 0)  write (*,*) '-- got start time ',
     &   hdr_year, hdr_jday, hdr_hour, hdr_min, hdr_sec,
     &   hdr_msec, sec

      return
      end

**************************************************************************
*
* HERE ARE SOME COMMENTS ON THE ORIGINAL VERSION OF MTINVERS USED AS
* STARTING SCRIPT FOR CODING VOLPIS.
* Some comments on the version used might be useful.
* Original version was coded by T. Dahm
* Some changes have been added by T. Dahm, F. Kruger and S. Cesca
*
* changes: 
* 17.08.05: - komp als real auf USER0 anstatt auf nevid 
* 30.09.00: -Zeilenumbruch in mtinvers.ball und mtinvers.ball1
*           -dispf.a dispf.b, dispf.l anstatt gross R, Z, T
*           -ebenso fuer synth.a, etc.
*           -itap eingefuehrt, um Taper zu steuern
*           - Zeitfenster beginnt jetzt bei iretard (bezogen auf origin),
*            (oder um Taperfenster voher wenn itap negativ)
*            und hat Laenge isamp (etwas laegr wenn itap < 0)
*            Hinter Fenster wird Spur mit Nullen aufgefuellt (vgl. isampf).
*    COSINUS FENSTER und Spuren numerisch testen
*           -Daten werden um iadd Punkte verlaengert
*
********************************************************************************
* Moment tensor inversion code (T. Dahm)
* This code is based on a program that has originally been used
* to estimate parameter of the extended earthquake source.
* In this version, only point source parameter can be estimated.
* No source time function is estimated (step function assumed)
* Also, acoustic emissions event have been analysed with a pre-version 
* of this code. Therefore, some parameter and statements are "commented out"
* or have no meaning for the present version.
*
* Input data are observed displacement seismograms (instrument deconvolved)
* and the corresponding Green functions (e.g. calculated with refgreen).
* Optionally, either waveforms or amplitude spectra are fitted.
* Different constraints can be used: 
*  (0) no constraint (6 moment tensor components)
*  (1) zero isotropic component constraint
*  (2) general dislocation source constraint (shear + tension dislocation)
*  (3) pure tension crack constraint
*  (4) pure double couple constraint
*
* An iterative inversion scheme is used, i.e. a starting model has to 
* be provided.
* For the time domain inversion and case (0) and (1) the inversion problem  
* is linear, and onlz then the starting model is uncritical; in all other
* cases different starting models should be tested to find the absolut 
* minimum misfit solution.
* Note also that the amplitude inversion can never resolve the "sign"
* of the moment tensor, i.e. P and T axes can be interchanged.
* 
* Input paramter are read from screen (or see mtinvers.cmd file).
* Input Files are 1. the controlling file (e.g. retard.dat)
*                 2. the waveforms (files termed  displ.n.k, where
*                    n is either z, n, e, and k is 1,2,3,....)
*                    The component n and numbers k are given in retard.dat
*                 3. the corresponding Green functions (displacement), 
*                    called PS_rSr.r.k, PS_rSz.k, etc. .
*                    Waveforms as well as Green functions are expected in 
*                    a special ascii format going back to Gerhard Muellers 
*                    reflectivity code.
*
* Output Files are:
*                 1. screen output
*                 2. result.best (list of the 10 best solutions)
*                 3. dispf.z, etc., containing filtered observations
*                 4. synth.z, etc, containing final fit solution traces 
*                    both may be plotted with e.g. Seismic Handler
*                 5 mtinvers.ball, psmeca moment tensor input file
*                 6 mtinvers.ball1, psmeca, best double couple input file
******************************************************************
*
*   some internals (not important for user)
*        xm(1) =0.5(M22-M11)      = sin(2sphi)
*        xm(2) =    M12           = cos(2sphi)
*        xm(3) =    M13           = sin(dphi)
*        xm(4) =    M23           = -cos(dphi))
*        xm(5) =1/3(0.5(M22+M11)) = VM
*        xm(6) =1/3(M11+M22+M33)  = Trace/3.
*
*     zu den Komponenten - Anteilen
*   | north |   | cos( azi )  -sin( azi )  |   | rad |
*   |       | = |                          | = |     |
*   | east  |   | sin( azi )   cos( azi )  |   | tra |
*
*   | rad   |   | cos( azi )   sin( azi )  |   | north |
*   |       | = |                          | = |       |
*   | tra   |   |-sin( azi )   cos( azi )  |   | east  |
*
*   | sv |   | cos(azi)*cos(toff)  sin( azi )*cos(toff) -sin(toff)|   |north|
*   | sh | = |-sin(azi)            cos(azi)              0        | = | east|
*   | p  |   | cos(azi)*sin(toff)  sin(azi)*sin(toff)    cos(toff)|   | down|
*
**************************************************************
