I have already posted regarding this,but I will try to explain.I calculate some functions,store them in arrays and then want to send them to main program but there I get some wrong values.My module:
module s2em
contains
c ---------------------------
subroutine sensd2emod_mod1(iper,as1,as2)
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)
real,dimension(20,21) :: apr1dsens,apr2dsens
real,dimension(20,21),intent(out) :: as1,as2
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
as1=apr1dsens
as2=apr2dsens
return
end subroutine
end module
module s2em
contains
c ---------------------------
subroutine sensd2emod_mod1(iper,as1,as2)
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)
real,dimension(20,21) :: apr1dsens,apr2dsens
real,dimension(20,21),intent(out) :: as1,as2
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
as1=apr1dsens
as2=apr2dsens
return
end subroutine
end module