      program ircg

      use constants
      use settings
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
      use rm0
c      use rm
      use lstep
      use rm1
      use newmodc
      use arraycons

      implicit none      
c
        integer*4 n1,i,j,icd,ivarprd,k,iper,is,fa,fb,iter,x
        integer*4 iprd,iprd0,isite,ico,var
	character fexpdat*200,fmodel*200,fout*200
	real*4 prunit,period
        real::  g,g1,e,e1
        real::  zxy_re_e,zxy_im_e,zxy_re_h,zxy_im_h
        real::  sens_zxy_ree,sens_zxy_ime,sens_zxy_reh,sens_zxy_imh
        real::  ro_e,ro_h,phase_e,phase_h
        real::  sens_ro_e,sens_phase_e,sens_ro_h,sens_phase_h
        real::  ro_ed, phase_ed,ro_hd,phase_hd
        real::  dmro_e,dmphase_e,dmro_h,dmphase_h,koef,dmfnorm
        real,dimension(21) :: l,lt,l0,koef1,res1
        real,dimension(57) :: r
        real,dimension(57,3) :: xn      
        real :: beta,norm2,norm20,norm2a,norm2b,z1,z2
        real, dimension(80,21) :: sens,wfm,swd
        real, dimension(21,80) :: ist,ftw
        integer,dimension(80,80) :: wd,wd2
        real,dimension(80,1):: d,a1    
        real,dimension(1,21):: t1
        real,dimension(21,1):: p,q,lmat,ltmat,dlm,b1
        real,dimension(1,1):: z11
        integer,dimension(21,21) :: wm,wm2,wm3
               
c
c     Give name of the file with periods and sites of interest.
c     Remove the backslash symbol '\' from the input format if
c     non-Intel Fortran compilar is used!

	open(2,file='test_p7.dat',status='old')
c
c     Give name of the file with the structure of the 2D model.
c     Again, the backslash '\' may cause an error with some
c     compilers.

	open(3,file='test_sh7.dat',status='old')
c
c     Give name for the output file. Again, take care with the
c     backslash '\' in the format!
c     I have changed,introduce 51,52,53
        open(4,file='test_res0.dat',status='old')
        open(5,file='test_res1.dat',status='old')
        open(10,file='e.dat',status='old')
        open(11,file='esens.dat',status='old')
        open(12,file='h.dat',status='old')
        open(13,file='hsens.dat',status='old')
        open(30,file='e1.dat',status='unknown')
        open(31,file='esens1.dat',status='unknown')
        open(32,file='h1.dat',status='unknown')
        open(41,file='data1.dat',status='unknown')
        open(42,file='data2.dat',status='unknown')
        open(43,file='data3.dat',status='unknown')
        open(44,file='data4.dat',status='unknown')
        open(33,file='hsens1.dat',status='unknown')
        open(14,file='sensresults.dat',status='unknown')
        open(114,file='sensresults1.dat',status='unknown')
        open(15,file='data.dat',status='old')
        open(16,file='misfit.dat',status='unknown')
        open(116,file='misfit1.dat',status='unknown')
        open(17,file='wd.dat',status='old',action='read')
        open(18,file='wm.dat',status='old',action='read')
        open(20,file='aprphase.dat',status='unknown')
        open(120,file='aprphase1.dat',status='unknown')
        open(220,file='aprphase2.dat',status='unknown')
        open(22,file='deltam0.dat',status='old',action='read')
        open(122,file='deltam1.dat',status='unknown')
        open(222,file='deltam.dat',status='unknown')
        open(100,file='koef.dat',status='unknown')
        open(101,file='koef1.dat',status='unknown')
        open(102,file='res.dat',status='unknown')
        open(103,file='res1.dat',status='unknown')
        open(104,file='resnew.dat',status='unknown')
        open(105,file='apriori.dat',status='old',action='read')
        open(128,file='z11.dat',status='unknown')
        open(200,file='newmodel.dat',status='unknown')
c
      call readmodelzero
c     Solution of the direct 2D problem and sensitivity evaluation
      call solve_mt2d_directsens 
      rewind(14)
      rewind(16)
      call calc0(norm20,l0)

c     do  iter=1,itermax
c     if  iter.eq.1
       iter=1
       rewind(3)
       rewind(2)
      call readmodelone
c     else
c      call redmodel
      call version2_solve_mt2d_directsens
      rewind(114)
      rewind(116)
  
         do j=1,21
         do i=1,80
            read(114,*)sens(i,j)
               end do 
              end do

        ist=transpose(sens)
        write(88,1001)ist
1001    format(80e12.5)

        rewind(17)
        do j=1,80
         do i=1,80
            read(17,*)wd(i,j)
         end do
       end do

       wd2=matmul(wd,wd)
       ftw=matmul(ist,wd2)           
    
       read(116,*)d
1002   format(e12.5)
       p=matmul(ftw,d)

c     calculate α*Wm^2*(m-mapr)
c     ovde alfa mora da se racuna drugacije 
c     po formuli iz zhdanova 
        x=iter-1   
        alpha=300*(0.92)**x

         do i=1,21
            read(18,*)(wm(i,j),j=1,21)
               end do 
       wm2=matmul(wm,wm)
       wm3=alpha*wm2
 
       read(122,1006)dlm
1006   format(e10.3) 
       q=matmul(wm3,dlm)
c      lstep is calculated l=ftranspose*wd2*r+alpha*wm2(m-mapr)
       lmat=p+q
       l=reshape(lmat,shape=shape(l))
       lbig(:,iter)=l(:)
       write(*,*)l

      norm2= (sum(l(1:21)**2))
      norm2big(iter)=norm2
      select case(iter)
      case(1)
      beta=norm2big(1)/norm20
      case default
      beta=norm2big(iter+1)/norm2big(iter)
      end select
      betabig(iter)=beta

      select case(iter)
      case(1)
      lt=lbig(:,1)+betabig(1)*l0
      write(*,*)lt
      case default
      lt=lbig(:,iter)+betabig(iter)*ltbig(:,iter-1)
      end select
      ltbig(:,iter)=lt(:)    

c     calculate koefficient k
      ltmat=reshape(lt,shape=shape(ltmat))  
      t1=transpose(ltmat)
      z11=matmul(t1,lmat)
      write(128,*)z11    
      
      rewind(114)
      do j=1,21
        do i=1,80
        read(114,*)sens(i,j)
        end do 
       end do

      rewind(17)
        do j=1,80
         do i=1,80
            read(17,*)wd(i,j)
         end do
       end do

       swd=matmul(wd,sens)
       a1=matmul(swd,ltmat)
       norm2a=(sum(a1(1:80,1:1)**2))
       write(*,*)norm2a

       rewind(18)
       do i=1,21
         read(18,*)(wm(i,j),j=1,21)
       end do 
     
       b1=matmul(wm,ltmat)
       norm2b=(sum(b1(1:21,1:1)**2))
       z2=norm2a+alpha*norm2b
       rewind(128)
       read(128,*)z1
       koef=z1/z2
       write(100,*)koef
c      calculate -k*l
       koef1=-koef*lt
       rewind(102)
       read(102,*)res1 
       res1=res1+koef1 
       write(103,501)res1
501    format(f8.3)
       rewind(103)
       rewind(5)
       call newmodel(r)  
       write(*,*)r
       rewind(5)
       call arrayc(xn)  
       xn(:,2)=r(:)
       write(*,*)xn(:,2) 
c       e nd do      
         end
c            
c
c
      subroutine solve_mt2d_directsens
c     ================================
c     Subroutine for the solution of a 2D MT model (isotropic) and for the
c     sensitivity evaluations w.r.t. the conductivities of selected homogeneous
c     domains of the model
c
      use constants
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c
      implicit none      
c
      integer*4 iper,isite,i,ico,iprd
      real*4 period,prunit
	complex*8 hxsurf
c
c     Compute boundary models for the 1D margins at the left and right hand
c     sides of the 2D model. The boundary models are the same for all periods
c     if the 2D model does not change
      call gebolr
c
c     Loop over all periods considered      
      do iper=1,np
        period=per(iper)
c
c     Solve the 2D model for the E-mode case for one period
        call d2emod(iper)
c
c     Evaluate the sensitivities of the 2D solution w.r.t. all the required
c     conductivities for the E-mode case for one period
        call sensd2emod(iper)
c
c     Solve the 2D for the H-mode case for one period
        call d2hmod(iper,hxsurf)
c        
c     Evaluate the sensitivities of the 2D solution w.r.t. all the required
c     conductivities for the H-mode case for one period        
        call sensd2hmod(iper,hxsurf)
c
c     Store the MT functions for the particular period in a larger array of
c     all the MT functions for all periods. Here, the impedances are output
c     at points of interest only
	  do i=1,nsites
	    isite=isit(i)
	    zmd(iper,i,1)=zmd1p(i,1)
	    zmd(iper,i,2)=zmd1p(i,2)
	    zmd(iper,i,3)=zmd1p(i,3)
	    zmd(iper,i,4)=zmd1p(i,4)

c
c     Auxiliary output, file 81: E and H-mode impedances at sites of interest
c     for one period
c	    write(81,'(2i5,5x,4e15.6)')iper,isite,
c     &      zmd(iper,i,1),zmd(iper,i,2),	        
c     &      zmd(iper,i,3),zmd(iper,i,4)
c     
        enddo    
c	    
c     Complete array of all required sensitivities w.r.t. flagged domains at
c     points of interest for all periods. For larger models, this array may
c     become prohibitely large!!!
        do i=1,nsites
          isite=isit(i)
	    do ico=1,nc
 	      if(ivarp(ico).ne.0)then
	        iprd=ivarp(ico)
	        zmdsens(iper,i,1,iprd)=zmd1psens(i,1,iprd)
	        zmdsens(iper,i,2,iprd)=zmd1psens(i,2,iprd)
	        zmdsens(iper,i,3,iprd)=zmd1psens(i,3,iprd)
	        zmdsens(iper,i,4,iprd)=zmd1psens(i,4,iprd)
c
c     Auxiliary output, file 82: E and H-mode impedance sensitivities at sites
c     of interest for one period	        
c	        write(82,'(3i5,4e15.6)')iper,isite,iprd,
c     &          zmdsens(iper,i,1,iprd),zmdsens(iper,i,2,iprd),	        
c     &          zmdsens(iper,i,3,iprd),zmdsens(iper,i,4,iprd)
c     
	      endif
	    enddo
	  enddo    
c
c     Go on to the next period        
      enddo      
c
      return
      end
c
c
c
      subroutine d2emod(iper)
c     =======================
c     Direct 2D modeling, E-mode
c
      use constants
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
c
      implicit none
c
      integer*4 iper,ieh,n1,m2,n2,kj,kj1,ki,ki1
      real*4 period,cof      
      complex*8 x1,x2,x3,x4,x5
c
      period=per(iper)
c
c     Define the polarization, IEH=1 means E-mode
      ieh=1
c
c     Compute boundary conditions (E-fields) for the 2D model 
      call pobo(ieh,period)
c
c     Form the discretized linear equations at all mesh nodes of the
c     2D model. Use the boundary conditions whenever the nodes on the
c     margins of the model are referred to. Remember that the multiplicative
c     factor COF takes into consideration that the mesh steps are given in km!
      cof=0.2*pi*pi/period
      m1=m-1
      n1=n-1
      m2=m-2
      n2=n-2
      a=0.
      b=0.
      nk=0
      do 2 kj=1,n2
      kj1=kj+1
      do 3 ki=1,m2
      ki1=ki+1
      nk=nk+1
      x1=0.5*(sz(ki)+sz(ki1))/sy(kj)
      x2=0.5*(sy(kj)+sy(kj1))/sz(ki)
      x3=0.5*(sy(kj)+sy(kj1))/sz(ki1)
      x4=0.5*(sz(ki)+sz(ki1))/sy(kj1)
      x5=-(x1+x2+x3+x4)+imc*cof*(cond(ic(ki,kj))*sy(kj)*sz(ki)+
     +                           cond(ic(ki,kj1))*sy(kj1)*sz(ki)+
     +                           cond(ic(ki1,kj))*sy(kj)*sz(ki1)+
     +                           cond(ic(ki1,kj1))*sy(kj1)*sz(ki1))
      a(nk,1)=x5
      if(ki.lt.m2)then
        a(nk,2)=x3
      else if(ki.eq.m2)then
        a(nk,2)=0.
        b(nk)=b(nk)-x3*bob(ieh,kj1)
      endif
      if(kj.lt.n2)then
        a(nk,m1)=x4
      else if(kj.eq.n2)then
        a(nk,m1)=0.
        b(nk)=b(nk)-x4*bor(ieh,ki1)
      endif
      if(kj.eq.1)b(nk)=b(nk)-x1*bol(ieh,ki1)
      if(ki.eq.1)b(nk)=b(nk)-x2*bot(ieh,kj1)
    3 continue
    2 continue
c
c     Solve the normal system by Gaussian elimination. The eliminated
c     form of the matrix A is saved to speed up the sensitivity computations!    
      call gsres
c
c     Solution the the 2D direct problem, E-mode, for the selected period
c     is finished here.
c
      return
      end
c
c
c
      subroutine d2hmod(iper,b0)
c     ==========================
c     Direct 2D modeling, H-mode
c
      use constants
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
c      
      implicit none
c
      integer*4 iper,ieh,n1,m2,n2,kj,kj1,kii,ki,ki1
      real*4 period,cof
      complex*8 x1,x2,x3,x4,x5
      complex*8 b0
c
      period=per(iper)
c
c     Define the polarization, IEH=2 means H-mode
      ieh=2
c      
c     Compute boundary conditions (H-fields) for the 2D model 
      call pobo(ieh,period)
c
c     Form the discretized linear equations at all mesh nodes in the conductive
c     domain of the 2D model (no field is computed in the air in the H-mode case, 
c     as H is constant in the air). Use the boundary conditions whenever the nodes
c     on the margins of the model are referred to. Remember that the multiplicative
c     factor COF takes into account that the mesh steps are given in km!
      b0=bol(ieh,ma)
      cof=0.2*pi*pi/period
      m1=m-ma
      n1=n-1
      m2=m1-1
      n2=n-2
      a=0.
      b=0.
      nk=0
      do 2 kj=1,n2
      kj1=kj+1
      do 3 kii=1,m2
      ki=kii+ma-1
      ki1=ki+1
      nk=nk+1
      x1=0.5*(sz(ki)*res(ic(ki,kj))+sz(ki1)*res(ic(ki1,kj)))/sy(kj)
      x2=0.5*(sy(kj)*res(ic(ki,kj))+sy(kj1)*res(ic(ki,kj1)))/sz(ki)
      x3=0.5*(sy(kj)*res(ic(ki1,kj))+sy(kj1)*res(ic(ki1,kj1)))/sz(ki1)
      x4=0.5*(sz(ki)*res(ic(ki,kj1))+sz(ki1)*res(ic(ki1,kj1)))/sy(kj1)
      x5=-(x1+x2+x3+x4)+imc*cof*(sy(kj)*sz(ki)+sy(kj1)*sz(ki)+
     +                           sy(kj)*sz(ki1)+sy(kj1)*sz(ki1))
      a(nk,1)=x5
      if(kii.lt.m2)then
        a(nk,2)=x3
      else if(kii.eq.m2)then
        a(nk,2)=0.
        b(nk)=b(nk)-x3*bob(ieh,kj1)
      endif
      if(kj.lt.n2)then
        a(nk,m1)=x4
      else if(kj.eq.n2)then
        a(nk,m1)=0.
        b(nk)=b(nk)-x4*bor(ieh,ki1)
      endif
      if(kj.eq.1)b(nk)=b(nk)-x1*bol(ieh,ki1)
      if(kii.eq.1)b(nk)=b(nk)-x2*bot(ieh,kj1)
    3 continue
    2 continue
c    
c     Solve the normal system by Gaussian elimination. The eliminated
c     form of the matrix A is saved to speed up the sensitivity computations!    
      call gsres
c
c     Solution the the 2D direct problem, H-mode, for the selected period
c     is finished here.
c
      return
      end
c
c
c
      subroutine gsres
c     ================
c     Full Gaussian elimination for the symmetric, banded system of equations
c     obtained by FV-discretizing the partial differential equations for E/H fields
c     in the 2-D model
c
      use constants
      use params
      use fdsystem
c      
      implicit none
c
      integer*4 nk1,nkm,mg,i,i1,me,k,ik,j,l,ne
      real*4 aik,ank
      complex*8 c
c
      nk1=nk-1
      nkm=nk-m1+2
      mg=m1
c
c     Forward elimination step
      do i=1,nk1
        i1=i-1
        if(i.ge.nkm)mg=mg-1
        me=2
        do k=2,mg
          aik=abs(real(a(i,k)))
          if(aik.ge.mach_real4)then
            c=a(i,k)/a(i,1)
            ik=i1+k
            j=0
            do l=me,mg
              j=j+1
              a(ik,j)=a(ik,j)-c*a(i,l)
            enddo
            b(ik)=b(ik)-c*b(i)
          endif  
          me=me+1
        enddo
      enddo
c
c     Backward substitution step
      ne=nk+1
  310 continue
      ne=ne-1
      if(ne.le.0)goto 350
      l=ne
      do k=2,m1
        l=l+1
        if(l.gt.nk)exit
        ank=abs(real(a(ne,k)))
        if(ank.gt.mach_real4)then
          b(ne)=b(ne)-a(ne,k)*b(l)
        endif  
      enddo
      b(ne)=b(ne)/a(ne,1)
      goto 310
  350 continue
c
c     Solution to the linear system is stored in B, the eliminated
c     form of the matrix A is saved for subsequent sensitivity
c     calculations        
c
      return
      end
c
c
c
      subroutine boem(ieh,period)
c     ===========================
c     Calculation of the boundary conditions at the left and right hand side
c     margins of the 2D model for one period and for the polarization of
c     the field given by IEH (1 for E-mode, 2 for H-mode)
c
      use constants
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c      
      implicit none
c
      integer*4 i1,i,ieh
      real*4 period
      complex*8 elf(mmax),mgf(mmax),z1(mmax)
      complex*8 z1sens(mmax,mmax)
c
c     1D MT impedances at all nodes on the left hand side of 
c     the model and elementary sensitivities of these impedances
c     w.r.t. the conductivities of the elementary (mesh) layers
      call z1idus(period,hbl,sbl,z1,z1sens,nbl)
c
c     Electric and magnetic fields and their elementary sensitivities
c     at all nodes on the left hand side of the model      
	mgf(1)=1.d0
	call h1iuds(period,hbl,sbl,z1,z1sens,
     &            mgf,elf,mgflsens,elflsens,nbl)
c
c     Boundary conditions for the EM fields on the left hand side 
c     of the model
c     ... in the conducting Earth     
	i1=0
	do i=ma,m
	  i1=i1+1
        if(ieh.eq.1)then
          bol(ieh,i)=elf(i1)
        else if(ieh.eq.2)then
          bol(ieh,i)=mgf(i1)
        endif
	enddo
c
c     ... and in the insulating air layer	
	do i=ma-1,1,-1
	  if(ieh.eq.1)then
	    bol(ieh,i)=8.d-4*pi*pi*imc*szz(i)*mgf(1)/period+elf(1)
	  else if(ieh.eq.2)then
	    bol(ieh,i)=mgf(1)	  
	  endif
	enddo
c
c     1D MT impedances at all nodes on the right hand side of 
c     the model and elementary sensitivities of these impedances
c     w.r.t. the conductivities of the elementary (mesh) layers
      call z1idus(period,hbr,sbr,z1,z1sens,nbr)
c
c     Electric and magnetic fields and their elementary sensitivities
c     at all nodes on the right hand side of the model      
	mgf(1)=1.d0
	call h1iuds(period,hbr,sbr,z1,z1sens,
     &            mgf,elf,mgfrsens,elfrsens,nbr)
c
c     Boundary conditions for the EM fields on the right hand side 
c     of the model.
c     ... in the conducting Earth     
	i1=0
	do i=ma,m
	  i1=i1+1
        if(ieh.eq.1)then
          bor(ieh,i)=elf(i1)
        else if(ieh.eq.2)then
          bor(ieh,i)=mgf(i1)
        endif
	enddo
c
c     ... and in the insulating air layer	
	do i=ma-1,1,-1
	  if(ieh.eq.1)then
	    bor(ieh,i)=8.d-4*pi*pi*imc*szz(i)*mgf(1)/period+elf(1)
	  else if(ieh.eq.2)then
	    bor(ieh,i)=mgf(1)	  
	  endif
	enddo
c
c     Boundary conditions on the left and right hand sides of the 2D model
c     are computed here completely, and they are used in forming the normal
c     FD linear systems for the direct 2D modeling. For the sensitivities,
c     only the elementary field sensitivities are evaluated here (i.e. 
c     sensitivities at marginal mesh nodes with respect to the conductivities
c     of the elementary (mesh) layers) and stored. The full bounadary conditions
c     for the sensitivities are constructed only in the routine BOEMSENS
c     
c
      return
      end
c
c
c
      subroutine gebolr
c     =================
c     Defines the 1D layered models on the left and right hand side
c     of the 2D model
c
      use constants
      use params
      use mt2dmod
c      
      implicit none
c
      integer*4 m1,n1,j,jma
c
      m1=m-1
      n1=n-1
c
c     Form the elementary 1D layered model on the left hand side margin
c     of the model. The model assumes that each mesh cell in the first (left-most)
c     column of the conductivity map continues as a layer towards plus
c     infinity
	nbl=m1-ma+2
	do j=1,nbl
	  jma=j+ma-1
	  if(j.eq.nbl)then
	    hbl(j)=0.
	    sbl(j)=cond(ic(jma-1,1))
	    isbl(j)=ic(jma-1,1)
	  else  
	    hbl(j)=1000.*sz(jma)
	    sbl(j)=cond(ic(jma,1))
	    isbl(j)=ic(jma,1)
	  endif  
	enddo
c
c     Form the elementary 1D layered model on the right hand side margin
c     of the model. The model assumes that each mesh cell in the last (right-most)
c     column of the conductivity map continues as a layer towards minus
c     infinity
	nbr=m1-ma+2
	do j=1,nbr
	  jma=j+ma-1
	  if(j.eq.nbr)then
	    hbr(j)=0.
	    sbr(j)=cond(ic(jma-1,n1))
	    isbr(j)=ic(jma-1,n1)
	  else  
	    hbr(j)=1000.*sz(jma)
	    sbr(j)=cond(ic(jma,n1))
	    isbr(j)=ic(jma,n1)
	  endif  
	enddo
c
      return
      end
c
c
c
      subroutine pobo(ieh,period)
c     ===========================
c     Boundary conditions for the 2-D model at all margins (left, right, top, bottom)
c     (here, the linear interpolation of values from the left and right hand
c     side margins of the model is used)
c
      use constants
      use params
      use mt2dmod
c
      implicit none
c
      integer*4 ieh,i
      real*4 period            
c
c     Compute the boundary conditions at left and right hand side margins
c     of the 2D model
      call boem(ieh,period)
c
c     Boundary conditions on the top and bottom of the 2D model      
      do 6 i=1,n
c
c     The bottom field is computed by linearly interpolating the fields from the
c     left bottom and right bottom corners of the model      
      bob(ieh,i)=bol(ieh,m)+syy(i)*(bor(ieh,m)-bol(ieh,m))/syy(n)
c
c     For the top boundary conditions, the linear interpolation is used in case
c     of the E-mode only; for the H-mode, the top (surface) boundary conditions
c     is constant
      if(ieh.eq.1)then
        bot(ieh,i)=bol(ieh,1)+syy(i)*(bor(ieh,1)-bol(ieh,1))/syy(n)
      else if(ieh.eq.2)then
        bot(ieh,i)=bol(ieh,ma)
      endif
    6 continue
c
      return
      end
c
c
c
	real*4 function phase(z)
c	------------------------
c     Computes the phase (argument) of a complex number Z in RAD's
c
      use constants
c
      implicit none      
c
      complex*8 z
      real*4 zr,zi
c
      zr=real(z)
      zi=aimag(z)
c                  
	if(abs(zr).lt.mach_real4)then
	  if(zi.gt.0.)then
	    phase=0.5*pi
	  else if(zi.lt.0.)then
	    phase=-0.5*pi
	  else if(zi.eq.0.)then
	    phase=0.d0
	  endif
	else
	  phase=atan(zi/zr)
	  if(zr.lt.0)phase=phase+pi
        if(phase.gt.pi)phase=phase-2.*pi
	endif
c
	return
	end
c
c
c
c     ---------------------------
      subroutine sensd2emod(iper)
c     ---------------------------
c     Sensitivities of the E-mode impedances with respect to the
c     conductivities of homogeneous subdomains of the 2D model. Here,
c     geomagnetic transfer functions (induction arrows) are not considered!
c
      use constants
      use settings
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c      
      implicit none
c
      integer*4 iper,m2,ma1,iprd0,i,i1,i2,ik1,ik0,ik2,iprd,j
      real*4 period,cof,spom,dw
      complex*8 sensv(nmax,(nmax-2)*(mmax-2)),sensw_t,sensw_c,sensw_b
      complex*8 v((nmax-2)*(mmax-2)),rp((nmax-2)*(mmax-2))
      complex*8 q0,q1,q2,hy,zxy,sens_zxy
      complex*8 zpom(nmax)
      real :: k,e,e1
c
      period=per(iper)
c
c     Pre-compute the site and MT function-specific vectors for sensitivity 
c     computations first. Only the "sites of interest" are treated here which
c     were specified on the input
      cof=1250.*period/(pi*pi)
      m2=m-2
      ma1=ma-1
      sensv=0.
      iprd0=0
c      write(50,'("---- E-mode impedance, period",f10.4)')period
      do i=1,nsites
c
c     Only sites from within the internal nodes of the mesh are considered.
c     Sites on the margins (1, N) are omitted.      
        if(isit(i).gt.1.and.isit(i).lt.n)then
          i1=isit(i)-1
          i2=i1+1
          ik1=(isit(i)-2)*m2+ma1      
          ik0=ik1-1
          ik2=ik1+1
          q0=imc*cof*sz(ma)/(sz(ma1)*(sz(ma1)+sz(ma)))
          q2=-imc*cof*sz(ma1)/(sz(ma)*(sz(ma1)+sz(ma)))
          spom=(sy(i1)*cond(ic(ma,i1))+sy(i2)*cond(ic(ma,i2)))/
     &         (sy(i1)+sy(i2))
          dw=500.*sz(ma1)*sz(ma)*spom/(sz(ma1)+sz(ma))
          q1=-(q0+q2)+dw
          hy=q0*b(ik0)+q1*b(ik1)+q2*b(ik2)
c
c     E-mode impedance          
          zxy=b(ik1)/hy
          zmd1p(i,1)=prev_zunit*real(zxy)
          zmd1p(i,2)=prev_zunit*aimag(zxy)
           k=period/(2*mu0)
          aprph1d(i,1)=k*(real(zxy)**2+aimag(zxy)**2)
          aprph1d(i,2)=atan(aimag(zxy)/real(zxy))
          read(41,*)data(i,1)
          read(42,*)data(i,2)
          dmf(i,1)=data(i,1)-aprph1d(i,1)
          dmf(i,2)=data(i,2)-aprph1d(i,2)
          write(16,*)dmf(i,1)
          write(16,*)dmf(i,2)
          e=(real(zxy)**2+aimag(zxy)**2)
          e1=1/e
c          write(*,*)e1
          zpom(i)=zxy/hy
         
          write(10,'(f12.4,2i5,2e15.6)')period,isit(i),iprd0,
     &                                  prev_zunit*zxy
          sensv(i,ik0)=-zpom(i)*q0
          sensv(i,ik1)=1./hy-zpom(i)*q1
          sensv(i,ik2)=-zpom(i)*q2
          if(ircpr.ne.0)then
            v=0.
            v(ik0)=sensv(i,ik0)
            v(ik1)=sensv(i,ik1)
            v(ik2)=sensv(i,ik2)
            call gsres0el(v)
            sensv(i,:)=v
          endif
        else
          write(*,*)' Site positioned at a margin',isit(i)
        endif
      enddo      
c      write(51,'("---- E-mode sensits, period",f10.4)')period
c
c     Loop over all conductivity domains flagged for sensitivity 
c     evaluations
	do i=1,nc
	  if(ivarp(i).ne.0)then
	    iprd=ivarp(i)
	    call rp2emod(iprd,period,rp)
	    v=rp
	    if(ircpr.eq.0)call gsres0el(v)	    
	    do j=1,nsites
            if(isit(j).gt.1.and.isit(j).lt.n)then
              i1=isit(j)-1
              i2=i1+1
              ik1=(isit(j)-2)*m2+ma1      
              ik0=ik1-1
              ik2=ik1+1
              if(ircpr.eq.0)then
                sens_zxy=sensv(j,ik0)*v(ik0)+
     &                   sensv(j,ik1)*v(ik1)+
     &                   sensv(j,ik2)*v(ik2)
              else
                sens_zxy=dot_product(conjg(sensv(j,1:nk)),v(1:nk))
              endif    
              sensw_t=0.
              sensw_c=0.
              if(ic(ma,i1).eq.iprd)then
                sensw_c=sensw_c+500.*sz(ma1)*sz(ma)*sy(i1)/
     &                  ((sz(ma1)+sz(ma))*(sy(i1)+sy(i2)))
              endif
              if(ic(ma,i2).eq.iprd)then
                sensw_c=sensw_c+500.*sz(ma1)*sz(ma)*sy(i2)/
     &                  ((sz(ma1)+sz(ma))*(sy(i1)+sy(i2)))
              endif
              sensw_b=0.      
c
c     Sensitivity of the E-mode impedance with respect to the CONDUCTIVITY
c     of the domain IPRD
              sens_zxy=sens_zxy-zpom(j)*
     &                 (sensw_t*b(ik0)+
     &                 sensw_c*b(ik1)+
     &                 sensw_b*b(ik2))
c
c     If the sensitivity of the E-mode impedance with respect to the RESISTIVITY
c     of the domain IPRD is required
             if(.not.lcondder)then
 	         sens_zxy=-sens_zxy/res(iprd)**2.
 	       endif  
c
c     Sensitivity of the E-mode impedance with respect to IPRD on the output
c     and then stored in ZMD1PSENS (indices 1,2)
              write(11,'(f12.4,2i5,2e15.6)')period,isit(j),iprd,
     &                                      prev_zunit*sens_zxy    
              zmd1psens(j,1,iprd)=prev_zunit*real(sens_zxy)
              zmd1psens(j,2,iprd)=prev_zunit*aimag(sens_zxy)
          e=(real(zxy)**2+aimag(zxy)**2)
          e1=1/e
          apr1dsens(j,1,iprd)=2*e*(real(sens_zxy)*real(zxy)+
     &    aimag(sens_zxy)*aimag(zxy))
          write(14,*) apr1dsens(j,1,iprd)
          apr1dsens(j,2,iprd)=e1*(aimag(sens_zxy)*real(zxy)-
     &    aimag(zxy)*real(sens_zxy))
          write(14,*) apr1dsens(j,2,iprd)
          

            else
c
c     Sites at any of the margins are not considered here
              write(*,*)' Site positioned at a margin',isit(j)
c              
            endif
	    enddo
	  endif  	  
	enddo
c
      return
      end
c
c
c
c     ----------------------------------
      subroutine sensd2hmod(iper,hxsurf)
c     ----------------------------------
c     Sensitivities of the H-mode impedances with respect to the
c     conductivities of homogeneous subdomains of the 2D model. 
c
      use constants
      use settings
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c
      implicit none      
c
      integer*4 iper,n4,m2,ma1,iprd0,i,i1,i2,ik2,iprd,j
      real*4 period,cof,cpom,zrh,zimh
      complex*8 hxsurf
      complex*8 sensv(nmax,(nmax-2)*(mmax-2)),sensw_b,sensgam
      complex*8 v((nmax-2)*(mmax-2)),rp((nmax-2)*(mmax-2))
      complex*8 q1,q2,ey,hx,zyx,sens_zyx
      real :: k,g,g1
c
      period=per(iper)
c
c     Pre-compute the site and MT function-specific vectors for sensitivity 
c     computations first. Only the "sites of interest" are treated here which
c     were specified on the input
      cof=0.0004*pi*pi/period
      n4=n-4
      m2=m-ma-1
      ma1=ma-1
      iprd0=0
c      write(52,'("---- H-mode impedance, period",f10.4)')period
      do i=1,nsites
c      
c     Only sites from within the internal nodes of the mesh are considered.
c     Sites on the margins (1, N) are omitted.      
        if(isit(i).gt.1.and.isit(i).lt.n)then
          i1=isit(i)-1
          i2=i1+1
          ik2=(isit(i)-2)*m2+1      
          hx=hxsurf
          q2=0.001*(sy(i1)*res(ic(ma,i1))+sy(i2)*res(ic(ma,i2)))/
     /             (sz(ma)*(sy(i1)+sy(i2)))      
          q1=imc*cof*sz(ma)-q2
          ey=q1*hxsurf+q2*b(ik2)
c
c     H-mode impedance, reversed sign used to make it comparable with 
c     the E-mode impedance
          zyx=-ey/hx
          if(.not.lneghsign)zyx=-zyx
          zmd1p(i,3)=prev_zunit*real(zyx)
          zmd1p(i,4)=prev_zunit*aimag(zyx)
          zrh=zmd1p(i,3)
          zimh=zmd1p(i,4)
          k=period/(2*mu0)
          aprph1d(i,3)=k*(real(zyx)**2+aimag(zyx)**2)
          aprph1d(i,4)=atan(aimag(zyx)/real(zyx))
          g=(real(zyx)**2+aimag(zyx)**2)
          g1=1/g
c          write(*,*)g1
          read(43,*)data(i,3)
          read(44,*)data(i,4)
          dmf(i,3)=data(i,3)-aprph1d(i,3)
          dmf(i,4)=data(i,4)-aprph1d(i,4)
          write(16,*)dmf(i,3)
          write(16,*)dmf(i,4)
          write(12,'(f12.4,2i5,2e15.6)')period,isit(i),iprd0,
     &                                  zyx
          sensv(i,ik2)=-q2/hxsurf
          if(ircpr.ne.0)then
            v=0.
            v(ik2)=sensv(i,ik2)
            call gsres0el(v)
            sensv(i,:)=v
          endif
        else
          write(*,*)' Site positioned at a margin',isit(i)
        endif
      enddo
c      write(53,'("---- H-mode sensits, period",f10.4)')period
c
c     Loop over all conductivity domains flagged for sensitivity 
c     evaluations
	do i=1,nc
	  if(ivarp(i).ne.0)then
	    iprd=ivarp(i)
	    call rp2hmod(iprd,period,rp)
	    v=rp
	    if(ircpr.eq.0)then
	      call gsres0el(v)	    
	    endif    
	    do j=1,nsites
            if(isit(j).gt.1.and.isit(j).lt.n)then
              i1=isit(j)-1
              i2=i1+1
              ik2=(isit(j)-2)*m2+1      
              if(ircpr.eq.0)then
                sens_zyx=sensv(j,ik2)*v(ik2)
              else
                sens_zyx=dot_product(conjg(sensv(j,1:nk)),v(1:nk))
              endif    
              sensw_b=0.
              sensgam=0.
              if(ic(ma,i1).eq.iprd)then
                cpom=0.001*sy(i1)/(sz(ma)*(sy(i1)+sy(i2)))
                sensw_b=sensw_b-cpom/hxsurf
                sensgam=sensgam+cpom
              endif
              if(ic(ma,i2).eq.iprd)then
                cpom=0.001*sy(i2)/(sz(ma)*(sy(i1)+sy(i2)))
                sensw_b=sensw_b-cpom/hxsurf
                sensgam=sensgam+cpom
              endif
c
c     Sensitivity of the H-mode impedance with respect to the RESISTIVITY
c     of the domain IPRD
              sens_zyx=sens_zyx+sensw_b*b(ik2)+sensgam
c
c     Sensitivity of the H-mode impedance with respect to the CONDUCTIVITY
c     of the domain IPRD
              if(lcondder)then
                sens_zyx=-sens_zyx*res(iprd)**2.    
              endif  
c
c     Sensitivity of the H-mode impedance with respect to the RES or COND
c     of the domain IPRD with reversed sign for comparison with the E-mode
c     sensitivity
              if(.not.lneghsign)then
                sens_zyx=-sens_zyx          
              endif  
c              
c     Sensitivity of the H-mode impedance with respect to IPRD on the output
c     and then stored in ZMD1PSENS (indices 3,4)
              write(13,'(f12.4,2i5,2e15.6)')period,isit(j),iprd,
     &                                      prev_zunit*sens_zyx    
              zmd1psens(j,3,iprd)=prev_zunit*real(sens_zyx)
              zmd1psens(j,4,iprd)=prev_zunit*aimag(sens_zyx)
c              write(*,*)g1
              apr1dsens(j,3,iprd)=2*g*(real(sens_zyx)*real(zyx)
     &        +aimag(sens_zyx)*aimag(zyx))
              write(14,*) apr1dsens(j,3,iprd)
              apr1dsens(j,4,iprd)=g1*(aimag(sens_zyx)*real(zyx)
     &        -aimag(zyx)*real(sens_zyx))
              write(14,*) apr1dsens(j,4,iprd)
c              
            else
c            
c     Sites at any of the margins are not considered here
              write(*,*)' Site positioned at a margin',isit(j)
c              
            endif
	    enddo
	  endif  	  
	enddo
c
      return
      end
c
c
c
c     ----------------------------------
      subroutine rp2emod(iprd,period,rp)
c     ----------------------------------
c     Right hand side of the normal system for sensitivity evaluations,
c     E-mode case, single period, single homogeneous domain
c
      use constants
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c
      implicit none      
c
      integer*4 iprd,ieh,n1,m2,n2,kj,kj1,ki,ki1
      real*4 cof,period
      complex*8 x1,x2,x3,x4,x5
      complex*8 rp((nmax-2)*(mmax-2))
c
c     Define the E-mode
      ieh=1
c
c     Boundary conditions for the sensitivity evaluation, mode and domain
c     specifications are given as routine parameters      
      call boemsens(ieh,iprd)
c
c     Form the right hand side of the linear system for the sensitivity
c     calculation from derivatives of the matrix A and from boundary
c     conditions
      cof=0.2*pi*pi/period
      m1=m-1
      n1=n-1
      m2=m-2
      n2=n-2
      rp=0.
      nk=0
      do 2 kj=1,n2
      kj1=kj+1
      do 3 ki=1,m2
      ki1=ki+1
      nk=nk+1
      x1=0.5*(sz(ki)+sz(ki1))/sy(kj)
      x2=0.5*(sy(kj)+sy(kj1))/sz(ki)
      x3=0.5*(sy(kj)+sy(kj1))/sz(ki1)
      x4=0.5*(sz(ki)+sz(ki1))/sy(kj1)
      x5=-(x1+x2+x3+x4)+imc*cof*(cond(ic(ki,kj))*sy(kj)*sz(ki)+
     +                           cond(ic(ki,kj1))*sy(kj1)*sz(ki)+
     +                           cond(ic(ki1,kj))*sy(kj)*sz(ki1)+
     +                           cond(ic(ki1,kj1))*sy(kj1)*sz(ki1))
      if(ki.ge.(ma-1))then
        if(ic(ki,kj).eq.iprd)rp(nk)=rp(nk)
     &                      -imc*cof*sy(kj)*sz(ki)*b(nk)
        if(ic(ki,kj1).eq.iprd)rp(nk)=rp(nk)
     &                        -imc*cof*sy(kj1)*sz(ki)*b(nk)
        if(ic(ki1,kj).eq.iprd)rp(nk)=rp(nk)
     &                        -imc*cof*sy(kj)*sz(ki1)*b(nk)
        if(ic(ki1,kj1).eq.iprd)rp(nk)=rp(nk)
     &                         -imc*cof*sy(kj1)*sz(ki1)*b(nk)
      endif
c
      if(ki.eq.m2)rp(nk)=rp(nk)-x3*bobsens(kj1)
      if(kj.eq.n2)rp(nk)=rp(nk)-x4*borsens(ki1)
      if(kj.eq.1)rp(nk)=rp(nk)-x1*bolsens(ki1)
      if(ki.eq.1)rp(nk)=rp(nk)-x2*botsens(kj1)
c      
    3 continue
    2 continue
c
      return
      end
c
c
c
c     ----------------------------------
      subroutine rp2hmod(iprd,period,rp)
c     ----------------------------------
c     Right hand side of the normal system for sensitivity evaluations,
c     H-mode case, single period, single homogeneous domain
c
      use constants
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c
      implicit none      
c
      integer*4 iprd,ieh,n1,m2,n2,kj,kj1,kii,ki,ki1
      real*4 cof,period
      complex*8 x1,x2,x3,x4,x5
      complex*8 rp((nmax-2)*(mmax-2)),b0
c
c     Define the E-mode
      ieh=2
c      
c     Boundary conditions for the sensitivity evaluation, mode and domain
c     specifications are given as routine parameters      
      call boemsens(ieh,iprd)
c
c     Form the right hand side of the linear system for the sensitivity
c     calculation from derivatives of the matrix A and from boundary
c     conditions
      b0=bol(ieh,ma)
      cof=0.2*pi*pi/period
      m1=m-ma
      n1=n-1
      m2=m1-1
      n2=n-2
      rp=0.
      nk=0
      do 2 kj=1,n2
      kj1=kj+1
      do 3 kii=1,m2
      ki=kii+ma-1
      ki1=ki+1
      nk=nk+1
      x1=0.5*(sz(ki)*res(ic(ki,kj))+sz(ki1)*res(ic(ki1,kj)))/sy(kj)
      x2=0.5*(sy(kj)*res(ic(ki,kj))+sy(kj1)*res(ic(ki,kj1)))/sz(ki)
      x3=0.5*(sy(kj)*res(ic(ki1,kj))+sy(kj1)*res(ic(ki1,kj1)))/sz(ki1)
      x4=0.5*(sz(ki)*res(ic(ki,kj1))+sz(ki1)*res(ic(ki1,kj1)))/sy(kj1)
      x5=-(x1+x2+x3+x4)+imc*cof*(sy(kj)*sz(ki)+sy(kj1)*sz(ki)+
     +                           sy(kj)*sz(ki1)+sy(kj1)*sz(ki1))
      if(ic(ki,kj).eq.iprd)then
        rp(nk)=rp(nk)+0.5*(sz(ki)/sy(kj)+sy(kj)/sz(ki))*b(nk)
        if(kii.gt.1.and.kj.gt.1)then
          rp(nk)=rp(nk)
     &    -0.5*sz(ki)*b(nk-m2)/sy(kj)
     &    -0.5*sy(kj)*b(nk-1)/sz(ki)
        else if(kii.eq.1.and.kj.gt.1)then
          rp(nk)=rp(nk)
     &    -0.5*sz(ki)*b(nk-m2)/sy(kj)
     &    -0.5*sy(kj)*b0/sz(ki)
        else if(kii.gt.1.and.kj.eq.1)then
          rp(nk)=rp(nk)
     &    -0.5*sz(ki)*bol(ieh,ki1)/sy(kj)
     &    -0.5*sy(kj)*b(nk-1)/sz(ki)
        else if(kii.eq.1.and.kj.eq.1)then
          rp(nk)=rp(nk)
     &    -0.5*sz(ki)*bol(ieh,ki1)/sy(kj)
     &    -0.5*sy(kj)*b0/sz(ki)               
        endif
      endif
      if(ic(ki,kj1).eq.iprd)then
        rp(nk)=rp(nk)+0.5*(sy(kj1)/sz(ki)+sz(ki)/sy(kj1))*b(nk)
        if(kii.gt.1.and.kj.lt.n2)then
          rp(nk)=rp(nk)
     &    -0.5*sy(kj1)*b(nk-1)/sz(ki)
     &    -0.5*sz(ki)*b(nk+m2)/sy(kj1)
        else if(kii.eq.1.and.kj.lt.n2)then
          rp(nk)=rp(nk)
     &    -0.5*sy(kj1)*b0/sz(ki)
     &    -0.5*sz(ki)*b(nk+m2)/sy(kj1)
        else if(kii.gt.1.and.kj.eq.n2)then
          rp(nk)=rp(nk)
     &    -0.5*sy(kj1)*b(nk-1)/sz(ki)
     &    -0.5*sz(ki)*bor(ieh,ki1)/sy(kj1)
        else if(kii.eq.1.and.kj.eq.n2)then
          rp(nk)=rp(nk)
     &    -0.5*sy(kj1)*b0/sz(ki)
     &    -0.5*sz(ki)*bor(ieh,ki1)/sy(kj1)
        endif                    
      endif
      if(ic(ki1,kj).eq.iprd)then
        rp(nk)=rp(nk)+0.5*(sz(ki1)/sy(kj)+sy(kj)/sz(ki1))*b(nk)
        if(kii.lt.m2.and.kj.gt.1)then
          rp(nk)=rp(nk)
     &    -0.5*sz(ki1)*b(nk-m2)/sy(kj)
     &    -0.5*sy(kj)*b(nk+1)/sz(ki1)
        else if(kii.eq.m2.and.kj.gt.1)then
          rp(nk)=rp(nk)
     &    -0.5*sz(ki1)*b(nk-m2)/sy(kj)
     &    -0.5*sy(kj)*bob(ieh,kj1)/sz(ki1)
        else if(kii.lt.m2.and.kj.eq.1)then
          rp(nk)=rp(nk)
     &    -0.5*sz(ki1)*bol(ieh,ki1)/sy(kj)
     &    -0.5*sy(kj)*b(nk+1)/sz(ki1)
        else if(kii.eq.m2.and.kj.eq.1)then               
          rp(nk)=rp(nk)
     &    -0.5*sz(ki1)*bol(ieh,ki1)/sy(kj)
     &    -0.5*sy(kj)*bob(ieh,kj1)/sz(ki1)
        endif
      endif
      if(ic(ki1,kj1).eq.iprd)then
        rp(nk)=rp(nk)+0.5*(sy(kj1)/sz(ki1)+sz(ki1)/sy(kj1))*b(nk)
        if(kii.lt.m2.and.kj.lt.n2)then
          rp(nk)=rp(nk)
     &    -0.5*sy(kj1)*b(nk+1)/sz(ki1)
     &    -0.5*sz(ki1)*b(nk+m2)/sy(kj1)
        else if(kii.eq.m2.and.kj.lt.n2)then
          rp(nk)=rp(nk)
     &    -0.5*sy(kj1)*bob(ieh,kj1)/sz(ki1)
     &    -0.5*sz(ki1)*b(nk+m2)/sy(kj1)
        else if(kii.lt.m2.and.kj.eq.n2)then
          rp(nk)=rp(nk)
     &    -0.5*sy(kj1)*b(nk+1)/sz(ki1)
     &    -0.5*sz(ki1)*bor(ieh,ki1)/sy(kj1)               
        else if(kii.eq.m2.and.kj.eq.n2)then
          rp(nk)=rp(nk)
     &    -0.5*sy(kj1)*bob(ieh,kj1)/sz(ki1)
     &    -0.5*sz(ki1)*bor(ieh,ki1)/sy(kj1)
        endif
      endif   
c
      if(kii.eq.m2)rp(nk)=rp(nk)+x3*bobsens(kj1)/res(iprd)**2.
      if(kj.eq.n2)rp(nk)=rp(nk)+x4*borsens(ki1)/res(iprd)**2.
      if(kj.eq.1)rp(nk)=rp(nk)+x1*bolsens(ki1)/res(iprd)**2.
      if(kii.eq.1)rp(nk)=rp(nk)+x2*botsens(kj1)/res(iprd)**2.
c          
    3 continue
    2 continue
c
      return
      end
c
c
c
c     ----------------------
      subroutine gsres0el(v)
c     ----------------------
c     Gaussian elimination with eliminated form of a symmetric banded matrix used
c     for repeated solutions with a number of different right hand sides
c
      use constants
      use params
      use fdsystem
c      
      implicit none
c
      integer*4 nk1,nkm,mg,i,i1,me,k,ik,ne,l
      real*4 aik,ank
      complex*8 c,v((nmax-2)*(mmax-2))
c
      nk1=nk-1
      nkm=nk-m1+2
      mg=m1
c
c     Reduced form of the forward elimination, applied to the
c     right hand side of the system only
      do i=1,nk1
        i1=i-1
        if(i.ge.nkm)mg=mg-1
        me=2
        do k=2,mg
          aik=real(a(i,k))
          if(abs(aik).ge.tiny(aik))then
            c=a(i,k)/a(i,1)
            ik=i1+k
            v(ik)=v(ik)-c*v(i)
          endif  
          me=me+1
        enddo
      enddo
c
c     Backward substitution step identical with that for the full
c     Gaussian elimination
      ne=nk+1
  310 continue
      ne=ne-1
      if(ne.le.0)goto 350
      l=ne
      do k=2,m1
        l=l+1
        if(l.gt.nk)exit
        ank=abs(real(a(ne,k)))
        if(ank.ge.mach_real4)then
          v(ne)=v(ne)-a(ne,k)*v(l)
        endif  
      enddo
      v(ne)=v(ne)/a(ne,1)
      goto 310
  350 continue
c
      return
      end
c
c
c
      subroutine boemsens(ieh,iprd)
c     -----------------------------
c     Boundary conditions at the margins of the 2D model for the sensitivity
c     calculations, mode specification IEH, integer IPRD defines the domain
c     w.r.t. whose conductivity is differentiated
c
      use constants
      use params
      use mt2dmod
      use mt2ddat
      use mt2dsens
c
      implicit none      
c
      integer*4 i,i1,j,j1,ieh,iprd
c
      bolsens=0.
      borsens=0.
      do i=ma,m-1
        i1=i-ma+1
        if(isbl(i1).eq.iprd)then
          if(ieh.eq.1)then
            do j=1,m-ma+1
              j1=j+ma-1
              bolsens(j1)=bolsens(j1)+elflsens(j,i1)
            enddo
          else if(ieh.eq.2)then
            do j=1,m-ma+1
              j1=j+ma-1
              bolsens(j1)=bolsens(j1)+mgflsens(j,i1)
            enddo
          endif
        endif
        if(isbr(i1).eq.iprd)then
          if(ieh.eq.1)then
            do j=1,m-ma+1
              j1=j+ma-1
              borsens(j1)=borsens(j1)+elfrsens(j,i1)
            enddo
          else if(ieh.eq.2)then
            do j=1,m-ma+1
              j1=j+ma-1
              borsens(j1)=borsens(j1)+mgfrsens(j,i1)
            enddo
          endif
        endif
      enddo 
c
      if(ieh.eq.1)then
        do j=ma-1,1,-1
          bolsens(j)=bolsens(ma)
          borsens(j)=borsens(ma)
        enddo
      endif
c                 
      do i=1,n
        bobsens(i)=bolsens(m)+syy(i)*(borsens(m)-bolsens(m))/syy(n)
        if(ieh.eq.1)then
          botsens(i)=bolsens(1)+syy(i)*(borsens(1)-bolsens(1))/syy(n)
        else if(ieh.eq.2)then
          botsens(i)=0.
        endif
      enddo
c
      return
      end
c
c
c
      subroutine z1idus(period,h,s,z1,z1sens,nl)
c     ------------------------------------------
c     Impedances and sensitivities for a 1D layered model. Both the impedances
c     and sensitivities are computed at all layer boundaries, the sensitivities
c     are computed with respect to the conductivities of all the layers present.
c     A stable bottom-to-top impedance propagation routine is used.
c
      use constants
      use params
c
      implicit none	
c
      real*4 period
      real*4 h(mmax),s(mmax)
      complex*8 z1(mmax),z1sens(mmax,mmax)
      complex*8 k,commi,cimmo,zet,th,vm1,ph,sech,kh
      real*4 omega
	integer*4 nl,l,l1,ll
c
      complex*8 dth,dsech,cx
      dth(cx)=(1.-cexp(-2.*cx))/(1.+cexp(-2.*cx))
      dsech(cx)=2.*cexp(-cx)/(1.+cexp(-2.*cx))
c
      omega=2.*pi/period
      cimmo=imc*omega*mu0
      commi=(1.-imc)*sqrt(0.5*omega*mu0)
c
      z1sens=0.
      do l=nl,1,-1
        k=commi*sqrt(s(l))
	  zet=cimmo/k
        if(l.eq.nl)then
          z1(l)=k/s(l)
          z1sens(l,l)=-0.5*z1(l)/s(l)
        else
          l1=l+1
          vm1=z1(l1)/zet
          kh=k*h(l)
          th=dth(kh)
          sech=dsech(kh)
          ph=sech/(1.-vm1*th)
	    z1(l)=zet*(vm1-th)/(1.-vm1*th)
	    do ll=nl,l1,-1
	      z1sens(l,ll)=z1sens(l1,ll)*ph**2.
	    enddo
	    z1sens(l,l)=-(0.5/s(l))*
     &	              (z1(l)-zet*(vm1-(1.-vm1**2.)*kh)*ph**2.)
        endif
      enddo
c
      return
      end
c
c
c
      subroutine h1iuds(period,h,s,z1,z1sens,
     &                  mgf,elf,mgfsens,elfsens,nl)
c     ---------------------------------------------
c     Fields and sensitivities for a 1D layered model. Both the H and E fields
c     and sensitivities are computed at all layer boundaries, the sensitivities
c     are computed with respect to the conductivities of all the layers present.
c     A stable top-to-bottom field propagation routine is applied and impedances and
c     their sensitivities from the routine Z1IDUS are used.
c
      use constants
      use params
c
      implicit none	
c
      real*4 period
      real*4 h(mmax),s(mmax)
      complex*8 z1(mmax),z1sens(mmax,mmax)
      complex*8 mgf(mmax),elf(mmax)
      complex*8 mgfsens(mmax,mmax),elfsens(mmax,mmax)
      complex*8 k,commi,cimmo,zet,th,ph,sech,vm1
      complex*8 kh,termlleqlm1,termllgtlm1
      real*4 omega
	integer*4 nl,l,l1,ll
c
      complex*8 dth,dsech,cx
      dth(cx)=(1.-cexp(-2.*cx))/(1.+cexp(-2.*cx))
      dsech(cx)=2.*cexp(-cx)/(1.+cexp(-2.*cx))
c
      omega=2.*pi/period
      cimmo=imc*omega*mu0
      commi=(1.-imc)*sqrt(0.5*omega*mu0)
c
      mgfsens=0.
      mgf(1)=1.
      elf(1)=z1(1)*mgf(1)
      do ll=1,nl
        elfsens(1,ll)=z1sens(1,ll)*mgf(1)
      enddo  
      do l=2,nl
        k=commi*sqrt(s(l-1))
	  zet=cimmo/k
        vm1=z1(l)/zet
        kh=k*h(l-1)
        th=dth(kh)
        sech=dsech(kh)
        ph=sech/(1.-vm1*th)
        mgf(l)=ph*mgf(l-1)
        elf(l)=z1(l)*mgf(l)
        do ll=1,nl
          mgfsens(l,ll)=ph*mgfsens(l-1,ll)
          if(ll.eq.(l-1))then
            termlleqlm1=(0.5/s(ll))*((vm1-kh)*th+vm1*kh)/(1.-vm1*th)
            mgfsens(l,ll)=mgfsens(l,ll)+termlleqlm1*mgf(l)
          else if(ll.gt.(l-1))then
            termllgtlm1=z1sens(l,ll)*th/(zet*(1.-vm1*th))
            mgfsens(l,ll)=mgfsens(l,ll)+termllgtlm1*mgf(l)  
          endif
          elfsens(l,ll)=z1sens(l,ll)*mgf(l)+z1(l)*mgfsens(l,ll)
        enddo
      enddo
c
      return
      end

      subroutine version2_solve_mt2d_directsens
c     ================================
c     Subroutine for the solution of a 2D MT model (isotropic) and for the
c     sensitivity evaluations w.r.t. the conductivities of selected homogeneous
c     domains of the model
c
      use constants
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c
      implicit none      
c
      integer*4 iper,isite,i,ico,iprd
      real*4 period,prunit
	complex*8 hxsurf
c
c     Compute boundary models for the 1D margins at the left and right hand
c     sides of the 2D model. The boundary models are the same for all periods
c     if the 2D model does not change
      call gebolr
c
c     Loop over all periods considered      
      do iper=1,np
        period=per(iper)
c
c     Solve the 2D model for the E-mode case for one period
        call d2emod(iper)
       
c     Evaluate the sensitivities of the 2D solution w.r.t. all the required
c     conductivities for the E-mode case for one period
        call sensd2emod_mod1(iper)
c
c     Solve the 2D for the H-mode case for one period
        call d2hmod(iper,hxsurf)
c        
c     Evaluate the sensitivities of the 2D solution w.r.t. all the required
c     conductivities for the H-mode case for one period        
        call sensd2hmod_mod1(iper,hxsurf)
c
c     Store the MT functions for the particular period in a larger array of
c     all the MT functions for all periods. Here, the impedances are output
c     at points of interest only
	  do i=1,nsites
	    isite=isit(i)
	    zmd2(iper,i,1)=zmd2p(i,1)
            zmd2(iper,i,2)=zmd2p(i,2)
	    zmd2(iper,i,3)=zmd2p(i,3)
	    zmd2(iper,i,4)=zmd2p(i,4)

c
c     Auxiliary output, file 81: E and H-mode impedances at sites of interest
c     for one period
c	    write(81,'(2i5,5x,4e15.6)')iper,isite,
c     &      zmd(iper,i,1),zmd(iper,i,2),	        
c     &      zmd(iper,i,3),zmd(iper,i,4)
c     
        enddo    
c	    
c     Complete array of all required sensitivities w.r.t. flagged domains at
c     points of interest for all periods. For larger models, this array may
c     become prohibitely large!!!
        do i=1,nsites
          isite=isit(i)
	    do ico=1,nc
 	      if(ivarp(ico).ne.0)then
	        iprd=ivarp(ico)
	        zmdsens2(iper,i,1,iprd)=zmd1psens2(i,1,iprd)
	        zmdsens2(iper,i,2,iprd)=zmd1psens2(i,2,iprd)
	        zmdsens2(iper,i,3,iprd)=zmd1psens2(i,3,iprd)
	        zmdsens2(iper,i,4,iprd)=zmd1psens2(i,4,iprd)
c
c     Auxiliary output, file 82: E and H-mode impedance sensitivities at sites
c     of interest for one period	        
c	        write(82,'(3i5,4e15.6)')iper,isite,iprd,
c     &          zmdsens(iper,i,1,iprd),zmdsens(iper,i,2,iprd),	        
c     &          zmdsens(iper,i,3,iprd),zmdsens(iper,i,4,iprd)
c     
	      endif
	    enddo
	  enddo    
c
c     Go on to the next period        
      enddo      
c
      return
      end

c     ---------------------------
      subroutine sensd2emod_mod1(iper)
c     ---------------------------
c     Sensitivities of the E-mode impedances with respect to the
c     conductivities of homogeneous subdomains of the 2D model. Here,
c     geomagnetic transfer functions (induction arrows) are not considered!
c
      use constants
      use settings
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c      
      implicit none
c
      integer*4 iper,m2,ma1,iprd0,i,i1,i2,ik1,ik0,ik2,iprd,j
      real*4 period,cof,spom,dw,k,e,e1
      complex*8 sensv(nmax,(nmax-2)*(mmax-2)),sensw_t,sensw_c,sensw_b
      complex*8 v((nmax-2)*(mmax-2)),rp((nmax-2)*(mmax-2))
      complex*8 q0,q1,q2,hy,zxy,sens_zxy
      complex*8 zpom(nmax)
c
      period=per(iper)
c
c     Pre-compute the site and MT function-specific vectors for sensitivity 
c     computations first. Only the "sites of interest" are treated here which
c     were specified on the input
      cof=1250.*period/(pi*pi)
      m2=m-2
      ma1=ma-1
      sensv=0.
      iprd0=0
c      write(50,'("---- E-mode impedance, period",f10.4)')period
      do i=1,nsites
c
c     Only sites from within the internal nodes of the mesh are considered.
c     Sites on the margins (1, N) are omitted.      
        if(isit(i).gt.1.and.isit(i).lt.n)then
          i1=isit(i)-1
          i2=i1+1
          ik1=(isit(i)-2)*m2+ma1      
          ik0=ik1-1
          ik2=ik1+1
          q0=imc*cof*sz(ma)/(sz(ma1)*(sz(ma1)+sz(ma)))
          q2=-imc*cof*sz(ma1)/(sz(ma)*(sz(ma1)+sz(ma)))
          spom=(sy(i1)*cond(ic(ma,i1))+sy(i2)*cond(ic(ma,i2)))/
     &         (sy(i1)+sy(i2))
          dw=500.*sz(ma1)*sz(ma)*spom/(sz(ma1)+sz(ma))
          q1=-(q0+q2)+dw
          hy=q0*b(ik0)+q1*b(ik1)+q2*b(ik2)
c
c     E-mode impedance          
          zxy=b(ik1)/hy
          zmd2p(i,1)=prev_zunit*real(zxy)
          zmd2p(i,2)=prev_zunit*aimag(zxy)
           k=period/(2*mu0)
          aprph1d(i,1)=k*(real(zxy)**2+aimag(zxy)**2)
          aprph1d(i,2)=atan(aimag(zxy)/real(zxy))
          rewind(41)
          rewind(42)
          read(41,*)data(i,1)
          read(42,*)data(i,2)
          dmf(i,1)=data(i,1)-aprph1d(i,1)
          dmf(i,2)=data(i,2)-aprph1d(i,2)
          write(116,*)dmf(i,1)
          write(116,*)dmf(i,2)
          e=(real(zxy)**2+aimag(zxy)**2)
          e1=1/e
c         
          zpom(i)=zxy/hy
          write(30,'(f12.4,2i5,2e15.6)')period,isit(i),iprd0,
     &                                  prev_zunit*zxy
          sensv(i,ik0)=-zpom(i)*q0
          sensv(i,ik1)=1./hy-zpom(i)*q1
          sensv(i,ik2)=-zpom(i)*q2
          if(ircpr.ne.0)then
            v=0.
            v(ik0)=sensv(i,ik0)
            v(ik1)=sensv(i,ik1)
            v(ik2)=sensv(i,ik2)
            call gsres0el(v)
            sensv(i,:)=v
          endif
        else
          write(*,*)' Site positioned at a margin',isit(i)
        endif
      enddo      
c      write(51,'("---- E-mode sensits, period",f10.4)')period
c
c     Loop over all conductivity domains flagged for sensitivity 
c     evaluations
	do i=1,nc
	  if(ivarp(i).ne.0)then
	    iprd=ivarp(i)
	    call rp2emod(iprd,period,rp)
	    v=rp
	    if(ircpr.eq.0)call gsres0el(v)	    
	    do j=1,nsites
            if(isit(j).gt.1.and.isit(j).lt.n)then
              i1=isit(j)-1
              i2=i1+1
              ik1=(isit(j)-2)*m2+ma1      
              ik0=ik1-1
              ik2=ik1+1
              if(ircpr.eq.0)then
                sens_zxy=sensv(j,ik0)*v(ik0)+
     &                   sensv(j,ik1)*v(ik1)+
     &                   sensv(j,ik2)*v(ik2)
              else
                sens_zxy=dot_product(conjg(sensv(j,1:nk)),v(1:nk))
              endif    
              sensw_t=0.
              sensw_c=0.
              if(ic(ma,i1).eq.iprd)then
                sensw_c=sensw_c+500.*sz(ma1)*sz(ma)*sy(i1)/
     &                  ((sz(ma1)+sz(ma))*(sy(i1)+sy(i2)))
              endif
              if(ic(ma,i2).eq.iprd)then
                sensw_c=sensw_c+500.*sz(ma1)*sz(ma)*sy(i2)/
     &                  ((sz(ma1)+sz(ma))*(sy(i1)+sy(i2)))
              endif
              sensw_b=0.      
c
c     Sensitivity of the E-mode impedance with respect to the CONDUCTIVITY
c     of the domain IPRD
              sens_zxy=sens_zxy-zpom(j)*
     &                 (sensw_t*b(ik0)+
     &                 sensw_c*b(ik1)+
     &                 sensw_b*b(ik2))
c
c     If the sensitivity of the E-mode impedance with respect to the RESISTIVITY
c     of the domain IPRD is required
             if(.not.lcondder)then
 	         sens_zxy=-sens_zxy/res(iprd)**2.
 	       endif  
c
c     Sensitivity of the E-mode impedance with respect to IPRD on the output
c     and then stored in ZMD1PSENS (indices 1,2)
              write(31,'(f12.4,2i5,2e15.6)')period,isit(j),iprd,
     &                                      prev_zunit*sens_zxy    
              zmd1psens2(j,1,iprd)=prev_zunit*real(sens_zxy)
              zmd1psens2(j,2,iprd)=prev_zunit*aimag(sens_zxy)
                    e=(real(zxy)**2+aimag(zxy)**2)
                    e1=1/e
          apr1dsens(j,1,iprd)=2*e*(real(sens_zxy)*real(zxy)+
     &    aimag(sens_zxy)*aimag(zxy))
          apr1dsens(j,2,iprd)=e1*(aimag(sens_zxy)*real(zxy)-
     &    aimag(zxy)*real(sens_zxy))
         write(114,*)apr1dsens(j,1,iprd)
         write(114,*)apr1dsens(j,2,iprd)
c              
            else
c
c     Sites at any of the margins are not considered here
              write(*,*)' Site positioned at a margin',isit(j)
c              
            endif
	    enddo
	  endif  	  
	enddo
c
      return
      end
c     ----------------------------------
      subroutine sensd2hmod_mod1(iper,hxsurf)
c     ----------------------------------
c     Sensitivities of the H-mode impedances with respect to the
c     conductivities of homogeneous subdomains of the 2D model. 
c
      use constants
      use settings
      use params
      use mt2dmod
      use fdsystem
      use mt2ddat
      use mt2dsens
c
      implicit none      
c
      integer*4 iper,n4,m2,ma1,iprd0,i,i1,i2,ik2,iprd,j
      real*4 period,cof,cpom,zrh,zimh,k,g,g1
      complex*8 hxsurf
      complex*8 sensv(nmax,(nmax-2)*(mmax-2)),sensw_b,sensgam
      complex*8 v((nmax-2)*(mmax-2)),rp((nmax-2)*(mmax-2))
      complex*8 q1,q2,ey,hx,zyx,sens_zyx
c
      period=per(iper)
c
c     Pre-compute the site and MT function-specific vectors for sensitivity 
c     computations first. Only the "sites of interest" are treated here which
c     were specified on the input
      cof=0.0004*pi*pi/period
      n4=n-4
      m2=m-ma-1
      ma1=ma-1
      iprd0=0
c      write(52,'("---- H-mode impedance, period",f10.4)')period
      do i=1,nsites
c      
c     Only sites from within the internal nodes of the mesh are considered.
c     Sites on the margins (1, N) are omitted.      
        if(isit(i).gt.1.and.isit(i).lt.n)then
          i1=isit(i)-1
          i2=i1+1
          ik2=(isit(i)-2)*m2+1      
          hx=hxsurf
          q2=0.001*(sy(i1)*res(ic(ma,i1))+sy(i2)*res(ic(ma,i2)))/
     /             (sz(ma)*(sy(i1)+sy(i2)))      
          q1=imc*cof*sz(ma)-q2
          ey=q1*hxsurf+q2*b(ik2)
c
c     H-mode impedance, reversed sign used to make it comparable with 
c     the E-mode impedance
          zyx=-ey/hx
          if(.not.lneghsign)zyx=-zyx
          zmd2p(i,3)=prev_zunit*real(zyx)
          zmd2p(i,4)=prev_zunit*aimag(zyx)
          zrh=zmd2p(i,3)
          zimh=zmd2p(i,4)
          k=period/(2*mu0)
          aprph1d(i,3)=k*(real(zyx)**2+aimag(zyx)**2)
          aprph1d(i,4)=atan(aimag(zyx)/real(zyx))
          rewind(43)
          rewind(44)
          read(43,*)data(i,3)
          read(44,*)data(i,4)
          dmf(i,3)=data(i,3)-aprph1d(i,3)
          dmf(i,4)=data(i,4)-aprph1d(i,4)
          write(116,*)dmf(i,3)
          write(116,*)dmf(i,4)
          write(32,'(f12.4,2i5,2e15.6)')period,isit(i),iprd0,
     &                                  zyx
          sensv(i,ik2)=-q2/hxsurf
          if(ircpr.ne.0)then
            v=0.
            v(ik2)=sensv(i,ik2)
            call gsres0el(v)
            sensv(i,:)=v
          endif
        else
          write(*,*)' Site positioned at a margin',isit(i)
        endif
      enddo
c      write(53,'("---- H-mode sensits, period",f10.4)')period
c
c     Loop over all conductivity domains flagged for sensitivity 
c     evaluations
	do i=1,nc
	  if(ivarp(i).ne.0)then
	    iprd=ivarp(i)
	    call rp2hmod(iprd,period,rp)
	    v=rp
	    if(ircpr.eq.0)then
	      call gsres0el(v)	    
	    endif    
	    do j=1,nsites
            if(isit(j).gt.1.and.isit(j).lt.n)then
              i1=isit(j)-1
              i2=i1+1
              ik2=(isit(j)-2)*m2+1      
              if(ircpr.eq.0)then
                sens_zyx=sensv(j,ik2)*v(ik2)
              else
                sens_zyx=dot_product(conjg(sensv(j,1:nk)),v(1:nk))
              endif    
              sensw_b=0.
              sensgam=0.
              if(ic(ma,i1).eq.iprd)then
                cpom=0.001*sy(i1)/(sz(ma)*(sy(i1)+sy(i2)))
                sensw_b=sensw_b-cpom/hxsurf
                sensgam=sensgam+cpom
              endif
              if(ic(ma,i2).eq.iprd)then
                cpom=0.001*sy(i2)/(sz(ma)*(sy(i1)+sy(i2)))
                sensw_b=sensw_b-cpom/hxsurf
                sensgam=sensgam+cpom
              endif
c
c     Sensitivity of the H-mode impedance with respect to the RESISTIVITY
c     of the domain IPRD
              sens_zyx=sens_zyx+sensw_b*b(ik2)+sensgam
c
c     Sensitivity of the H-mode impedance with respect to the CONDUCTIVITY
c     of the domain IPRD
              if(lcondder)then
                sens_zyx=-sens_zyx*res(iprd)**2.    
              endif  
c
c     Sensitivity of the H-mode impedance with respect to the RES or COND
c     of the domain IPRD with reversed sign for comparison with the E-mode
c     sensitivity
              if(.not.lneghsign)then
                sens_zyx=-sens_zyx          
              endif  
c              
c     Sensitivity of the H-mode impedance with respect to IPRD on the output
c     and then stored in ZMD1PSENS (indices 3,4)
              write(33,'(f12.4,2i5,2e15.6)')period,isit(j),iprd,
     &                                      prev_zunit*sens_zyx    
              zmd1psens2(j,3,iprd)=prev_zunit*real(sens_zyx)
              zmd1psens2(j,4,iprd)=prev_zunit*aimag(sens_zyx)
              g=real(zyx)**2+aimag(zyx)**2
              g1=1/g
              apr1dsens(j,3,iprd)=2*g*(real(sens_zyx)*real(zyx)
     &        +aimag(sens_zyx)*aimag(zyx))
              apr1dsens(j,4,iprd)=g1*(aimag(sens_zyx)*real(zyx)
     &        -aimag(zyx)*real(sens_zyx))
               write(114,*)apr1dsens(j,3,iprd)
               write(114,*)apr1dsens(j,4,iprd)
c              
            else
c            
c     Sites at any of the margins are not considered here
              write(*,*)' Site positioned at a margin',isit(j)
c              
            endif
	    enddo
	  endif  	  
	enddo
c
      return
      end



c
        
c
     
