merge rank 3 arrays

[milenko76]

I have problem like this:real*4 apr1dsens(nmax,nfmax,ncmax)
apr1dsens(j,1,iprd),apr1dsens(j,2,iprd),apr1dsens(j,3,iprd),apr1dsens(j,4,iprd)
All have 420 numerical values.How to join them into single array?Can merge intinsic function be used for 4 arrays?I have tried like this but got some apsurd results:
do i=1,4
apr4sens(j,i,iprd)=apr1dsens(j,1,iprd)
apr4sens(j,i,iprd)=apr1dsens(j,2,iprd)
apr4sens(j,i,iprd)=apr1dsens(j,3,iprd)
apr4sens(j,i,iprd)=apr1dsens(j,4,iprd)
end do

 

[xwb]

Don't quite understand what you're trying to do. You're assigning the same value in apr1dsens to the same variable in apr4sens 4 times.

When you say merge, what exactly do you really want to do. The merge function does the following. Say you have 2 arrays

xxx = reshape((/ 1, 2, 3, 4, 5, 6 /), (/ 2, 3 /))
yyy = reshape((/ 7, 8, 9, 0, 1, 2 /), (/ 2, 3 /))
mask= reshape((/ .true., .false., .false., .true., .true., .false. /), (/ 2, 3 /))

result = merge (xxx, yyy, mask)

The mask will items from xxx if it is true, and yyy if it is false. So result will be 1, 8, 9, 4, 5, 2. Is that what you want?

 

[milenko76]

No,english is not my first language.I want to put values from 4 arrays into one and then later to play with these.I have 1680 numerical value in every array.I tried this but it doesn't work.
subroutine transenpas
c ========================================
use params
use mt2dsens
integer :: k

do k=1,4
apr4sens(k,:,= apr1dsens(j,iprd)
apr4sens(k,:,= apr2dsens(j,iprd)
apr4sens(k,:,= apr3dsens(j,iprd)
apr4sens(k,:,= apr4dsens(j,iprd)
end do
module mt2dsens
c ===============
c Parameters used in sensitivity calculations
use params
c
implicit none
c
c Compressed version if the sensitivity flag IVAR
integer*4 ivarp(ncmax)
c
c Sensitivities of the MT functions for one period of the EM filed
c at all surface nodes and for all domains which have been flagged
c for the sensitivity evaluation. In the present version, only MT
c impedances are considered, not the induction arrows!
real*4 zmd1psens(nmax,nfmax,ncmax)
real*4 zmd1psens2(nmax,nfmax,ncmax)
real*4 apr1dsens(nfmax,ncmax)
real*4 apr2dsens(nfmax,ncmax)
real*4 apr3dsens(nfmax,ncmax)
real*4 apr4dsens(nfmax,ncmax)
c
c Aggregate version of the above array for all periods. May become
c pretty huge!!!
real*4 zmdsens(npmax,nmax,nfmax,ncmax)
real*4 zmdsens2(npmax,nmax,nfmax,ncmax)
c Aggregate version od sensitivity matrix which need to be reshaped
real*4 apr4sens(n4mat,nfmax,ncmax)

 

[FJacq]

       do k=1,4
       apr4sens(k,:,:)= apr1dsens(j,iprd)
       apr4sens(k,:,:)= apr2dsens(j,iprd)
       apr4sens(k,:,:)= apr3dsens(j,iprd)
       apr4sens(k,:,:)= apr4dsens(j,iprd)
       end do

Do you really understand what you have programmed ? From a logical point of view, these instructions are just non sense !

Try instead, without any loop :

       real :: group(nmax,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

I have change the name "apr4sens" into "group" because of the possible confusion with "apr4dsens" which obfuscates your code.

François Jacq

 

[milenko76]

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

 

[milenko76]

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

 

[milenko76]

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

 

[milenko76]

Works better now but still not perfect.I have changed nfmax and ncmax to 20 and 21.Then as a result at he screen I got 420 numerical values,but some of them are just 0:
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
1.2908742E-18 1.1608330E-17 1.1927672E-16 1.2435267E-14 4.0228403E-14
2.1946319E-13 -4.8073974E-13 -1.6162873E-12 -1.1688290E-11 -1.6727184E-11
-2.0666685E-11 -4.4495192E-11 -7.6163804E-11 -3.5470935E-11 -8.4115935E-12
-1.5911483E-13 -1.1818656E-14 -3.7578939E-15 8.9163326E-16 1.1130369E-16
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
1.5810671E-18 1.3264045E-17 2.5806894E-16 -1.5020346E-15 -2.5819510E-13
-3.5149103E-13 -7.3843286E-13 -1.1032785E-12 -4.5275910E-12 -5.9744037E-12
-6.4341674E-12 -1.0753063E-11 -6.5364615E-12 -3.0312400E-12 1.3167024E-10
-5.8838056E-12 -1.6515093E-13 -5.4458662E-15 -7.3791884E-15 -4.9404335E-16
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
1.2103064E-18 8.6202477E-18 1.3681297E-16 -2.1310601E-15 -7.6861714E-14
-3.6190448E-14 -4.7557376E-14 -5.8863850E-14 -1.1567471E-13 -1.0792758E-13
-7.2387829E-14 -1.5461306E-13 -9.1826740E-14 -3.9773290E-13 -7.3846550E-13
-3.0486087E-12 5.6226128E-11 -5.2306640E-12 -5.9632421E-14 -1.5679719E-15
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
9.4598716E-18 1.4132135E-16 2.0056516E-13 -7.6410225E-11 -5.7890356E-11
-3.8615608E-12 6.7739965E-13 6.3318740E-13 9.4493976E-13 6.5950132E-13
3.9332822E-13 2.1042835E-13 5.1335279E-14 1.9013698E-13 2.2994478E-14
3.6856335E-15 -6.4241552E-16 -3.4615460E-15 1.2640444E-15 1.7410016E-16
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
1.2260472E-17 1.4296709E-16 4.8845658E-15 2.7806158E-15 -1.0331626E-10
-3.6651970E-10 -1.0284996E-10 -1.8013107E-11 -3.4645127E-11 -2.2542404E-11
-9.7404671E-12 -8.5855377E-12 4.3947284E-13 4.5772766E-13 -5.0359288E-13
-3.0335299E-14 -7.9230700E-15 -4.5027090E-15 2.0727538E-16 1.0118936E-16
0.0000000E+00

 

[FJacq]

Sorry but what do you mean by "it does not work" ?

Of course, if you tried the routine tansenpas as shown at the end of your program, it probably cannot work !

        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

You did not give the contents of the module "params" but I assume that the arrays apr1dsens, apr2dsens, apr3dsens and apr4dsens are not defined inside. so these arrays must be passed by arguments !

Something like :

        subroutine transenpas(group,apr1dsens,apr2dsens
      *                      ,apr3dsens,apr2dsens)
c      ========================================
       use params
       real,intent(out) :: group(n4mat,nfmax,ncmax)
       real,intent(in)  :: apr1dsens(nfmax,ncmax)
       real,intent(in)  :: apr2dsens(nfmax,ncmax)
       real,intent(in)  :: apr3dsens(nfmax,ncmax)
       real,intent(in)  :: apr4dsens(nfmax,ncmax)
       group(1,:,:)= apr1dsens
       group(2,:,:)= apr2dsens
       group(3,:,:)= apr3dsens
       group(4,:,:)= apr4dsens
       write(*,*)group(1,:,:)
c
        end subroutine

Please, before saying "it does not work", try to program with care and logic !

Here it was not difficult to see that the arrays apr*dsens were never defined.


François Jacq

 

[milenko76]

Well actualy there were define at mt2dsens module.
module mt2dsens
c ===============
c Parameters used in sensitivity calculations
use params
c
implicit none
c
c Compressed version if the sensitivity flag IVAR
integer*4 ivarp(ncmax)
c
c Sensitivities of the MT functions for one period of the EM filed
c at all surface nodes and for all domains which have been flagged
c for the sensitivity evaluation. In the present version, only MT
c impedances are considered, not the induction arrows!
real*4 zmd1psens(nmax,nfmax,ncmax)
real*4 zmd1psens2(nmax,nfmax,ncmax)
real*4 apr1dsens(npmax,ncs)
real*4 apr2dsens(npmax,ncs)
real*4 apr3dsens(npmax,ncs)
real*4 apr4dsens(npmax,ncs)
c
c Aggregate version of the above array for all periods. May become
c pretty huge!!!
real*4 zmdsens(npmax,nmax,nfmax,ncmax)
real*4 zmdsens2(npmax,nmax,nfmax,ncmax)
c Aggregate version od sensitivity matrix which need to be reshaped
real*4 apr4sens(n4mat,nfmax,ncmax)
c
c Elementary sensitivities of the 1D marginal models at the left and
c right hand side of the 2D model. By elementary sensitivities, we
c mean sensitivities at each vertical node w.r.t. each elementary
c layer of the model, defined as a layer between two successive
c horizontal mesh lines
complex*8 elflsens(mmax,mmax),mgflsens(mmax,mmax)
complex*8 elfrsens(mmax,mmax),mgfrsens(mmax,mmax)
c
c Boundary conditions for the sensitivity evaluations at the Left,
c Right, Top, and Bottom margins of the 2D model. These boundary
c conditions are computed on the fly for each polarization and for
c each derivative domain. After they have been used in the particular
c sensitivity evaluation, they are overwritten.
complex*8 bolsens(mmax),borsens(mmax)
complex*8 bobsens(nmax),botsens(nmax)
end module

 

[milenko76]

Now I have actually pasted your code,but I have problems like this:
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
ircg 08050FBF transenpas_ 2113 ircg.for
ircg 0804D5F9 solve_mt2d_direct 258 ircg.for
ircg 0804AA31 MAIN__ 73 ircg.for

 

[FJacq]

Sorry but I don't see neither the call statement to transenpas nor the declaration of the arguments in the calling routine.

In particular, have you declare the array "group" in solve_mt2d_direct according to its definition in transenpas ?

François Jacq

 

[milenko76]

Well I am using module mt2sens so group is defined in solve_mt2d_direct.I have problem in compiling the subroutine.
ifort -c ircg.for
ircg.for(2127): error #6366: The shapes of the array expressions do not conform. [GROUP]
group(1,:,= apr1dsens
-------^
compilation aborted for ircg.for (code 1)
subroutine transenpas(group,apr1dsens,apr2dsens,
& apr3dsens,apr4dsens)
c ========================================
use params
real,intent(out) :: group(n4mat,npmax,ncs)
real,intent(in) :: apr1dsens(npmax,ncs)
real,intent(in) :: apr2dsens(nfmax,ncs)
real,intent(in) :: apr3dsens(nfmax,ncs)
real,intent(in) :: apr4dsens(nfmax,ncs)
group(1,:,= apr1dsens
group(2,:,= apr2dsens
group(3,:,= apr3dsens
group(4,:,= apr4dsens
write(*,*)group(1,:,
c
end subroutine

 

[milenko76]

ok ,its npmax instead of nfmax.