No,it doesn't work.I calculate these for arrays in 2 different subroutines.I will try to attach the whole code.
program ircg
use definitions
use constants
use settings
use params
use mt2dmod
use fdsystem
use mt2ddat
use mt2dsens
use rm0
use lstep
use rm1
use newmodc
use arraycons
implicit none
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='a11.dat',status='unknown')
open(200,file='newmodel.dat',status='unknown')
open(320,file='aproba.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
iter=1
if (iter.eq.1) then
rewind(3)
rewind(2)
call readmodelone
else
rewind(3)
rewind(2)
call readmodel
end if
call version2_solve_mt2d_directsens
c if (dmfnorm.le.10) stop
c else
c continue
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
rewind(18)
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)
call vfm
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
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)
rewind(5)
call arrayc(xn)
xn
,2)=r
)
xnbig
,:,iter+1)=xn
c end 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,j,f
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)
call transenpas
c Calculates Frobenius norm square of datamisfit
call norm
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
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
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))
write(320,*) aprph1d(i,1), aprph1d(i,2)
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,iprd)=2*e*(real(sens_zxy)*real(zxy)+
& aimag(sens_zxy)*aimag(zxy))
write(14,*) apr1dsens(j,iprd)
apr2dsens(j,iprd)=e1*(aimag(sens_zxy)*real(zxy)-
& aimag(zxy)*real(sens_zxy))
write(14,*) apr2dsens(j,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
apr3dsens(j,iprd)=2*g*(real(sens_zyx)*real(zyx)
& +aimag(sens_zyx)*aimag(zyx))
write(14,*) apr3dsens(j,iprd)
apr4dsens(j,iprd)=g1*(aimag(sens_zyx)*real(zyx)
& -aimag(zyx)*real(sens_zyx))
write(14,*) apr4dsens(j,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
if(ieh.eq.1)then
botsens(i)=bolsens(1)+syy(i)*(borsens(1)-bolsens(1))/syy
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,iprd)=2*e*(real(sens_zxy)*real(zxy)+
& aimag(sens_zxy)*aimag(zxy))
apr2dsens(j,iprd)=e1*(aimag(sens_zxy)*real(zxy)-
& aimag(zxy)*real(sens_zxy))
write(114,*)apr1dsens(j,iprd)
write(114,*)apr2dsens(j,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
apr3dsens(j,iprd)=2*g*(real(sens_zyx)*real(zyx)
& +aimag(sens_zyx)*aimag(zyx))
apr4dsens(j,iprd)=g1*(aimag(sens_zyx)*real(zyx)
& -aimag(zyx)*real(sens_zyx))
write(114,*)apr3dsens(j,iprd)
write(114,*)apr4dsens(j,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 rm5
c ==============================
real :: icd,res,ivar
real,dimension(57,3) :: xn
real,dimension(57) :: resar1,ivarar1,icdar1
icdar1= xn
,1)
resar1= xn
,2)
ivarar1=xn
,3)
do i=1,57
icd=icdar1(i)
res=resar1(i)
ivar=ivarar1(i)
end do
end subroutine
c ==============================
subroutine vfm
c ==============================
real,dimension(1) :: z1
real,dimension(1,1) :: z11
z1=reshape(z11,shape(z1))
end subroutine
c ===============================
subroutine norm
c ===============================
real:: dmfnorm
real :: dmfnorm1,dmfnorm2,dmfnorm3,dmfnorm4
real,dimension(20,1) :: dmf
dmfnorm1= (sum(dmf(1:20,1:1)**2))
dmfnorm2= (sum(dmf(1:20,2:2)**2))
dmfnorm3= (sum(dmf(1:20,3:3)**2))
dmfnorm4= (sum(dmf(1:20,4:4)**2))
dmfnorm=dmfnorm1+dmfnorm2+dmfnorm3+dmfnorm4
write(*,*)dmfnorm
end subroutine
c ================================
subroutine readmodel
c ================================
use constants
use settings
use params
use mt2dmod
use fdsystem
use mt2ddat
use mt2dsens
read(3,*)n,m,ma
n1=n-1
m1=m-1
c
c Read the horizontal and vertical mesh steps, in km
read(3,*)(sy(i),i=1,n1)
read(3,*)(sz(i),i=1,m1)
syy(1)=0.
do i=1,n1
syy(i+1)=syy(i)+sy(i)
enddo
szz(ma)=0.
do i=ma,m1
szz(i+1)=szz(i)+sz(i)
enddo
do i=ma-1,1,-1
szz(i)=szz(i+1)-sz(i)
enddo
do i=1,ma-1
do j=1,n1
ic(i,j)=1
enddo
enddo
do i=ma,m1
read(3,*)(ic(i,j),j=1,n1)
enddo
icd=1
res(icd)=-1.
cond(icd)=0.
ivar(icd)=0
read(3,*)nc
nc=nc+1
ivarprd=0
close(3)
c read array from 5
xnbig
,:,iter)=xn
call rm5
cond(icd)=1./res(icd)
if(ivar(icd).ne.0)then
ivarprd=ivarprd+1
ivarp(ivarprd)=icd
endif
c
c
c Inputs from the period/sites file.
c Actual number of periods is read and a list of periods in seconds
read(2,*)np
read(2,*)(per(i),i=1,np)
c
read(2,*)nsites
read(2,*)(isit(i),i=1,nsites)
c
c Flag whether the reciprocity approach is to be used for the sensitivity
c calculations (IRCPR=1) or not (IRCPR=0)
read(2,*)ircpr
c
c That is all the input from the period/structure file
c close(2)
c
c Conversion of units (SI or practical) for the impedances
if(lsiunits)then
prev_zunit=prev_siunit
else
prev_zunit=prev_pracunit
endif
c
end subroutine
c ========================================
subroutine transenpas
c ========================================
use params
real :: group(n4mat,nfmax,ncmax)
real :: apr1dsens(nfmax,ncmax)
real :: apr2dsens(nfmax,ncmax)
real :: apr3dsens(nfmax,ncmax)
real :: apr4dsens(nfmax,ncmax)
group(1,:,
= apr1dsens
group(2,:,
= apr2dsens
group(3,:,
= apr3dsens
group(4,:,
= apr4dsens
write(*,*)group(1,:,
c
end subroutine
