I have a problem,actually I want to join two FORTRAN codes.One of them is main and the other calculates some issues(geophysics stuff).I am not sure that the second code can be represented as subroutine because it has 2000+ lines,with 10 subroutines inside itself.How to solve this?Is it possible that one executable calls the other?I am working on Linux,with Intel Fortran compiler.
Tek-Tips is the largest IT community on the Internet today!
Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!
-
Congratulations SkipVought on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!
FORTRAN 90 subroutine
- Thread starter milenko76
- Start date
NickFort
Technical User
- Jun 10, 2010
- 113
- 0
- 0
Why would the size of the subroutine matter? You just need an explicit interface from the subroutine(s) in the driver program. Alternatively, you can just put all of the subroutines into a module, and then "use" that module in the main program.
However, I'm not entirely certain I understand your problem. Maybe give us a bit more info?
--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.
However, I'm not entirely certain I understand your problem. Maybe give us a bit more info?
--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.
GerritGroot
Technical User
I suspect that the 2000+ code hasn't ANY separate file with subroutines or modules, right?
There are indeed two ways, both already mentioned by NickFort. The quick and dirty way is to place the word "CONTAINS" before the first (in the main program included) subroutine and replace the END PROGRAM of the main to the bottom of it all calling it END SUBROUTINE. The header of the main program must be changed of course to a subroutine with the proper INTENT(IN)'s and OUT's to give you the right stuff.
For the quick and dirty way, I like Nick's "dump it in a module" solution, which is possibly the one that will be working inmediately.
The clean way is to take each subroutine out of the main code and place them in a separate file (if there are a lot of them, you can also automise that... ...in fortran reading and writing lines, recognising the word "SUBROUTINE" and "END SUBROUYINE")
There are indeed two ways, both already mentioned by NickFort. The quick and dirty way is to place the word "CONTAINS" before the first (in the main program included) subroutine and replace the END PROGRAM of the main to the bottom of it all calling it END SUBROUTINE. The header of the main program must be changed of course to a subroutine with the proper INTENT(IN)'s and OUT's to give you the right stuff.
For the quick and dirty way, I like Nick's "dump it in a module" solution, which is possibly the one that will be working inmediately.
The clean way is to take each subroutine out of the main code and place them in a separate file (if there are a lot of them, you can also automise that... ...in fortran reading and writing lines, recognising the word "SUBROUTINE" and "END SUBROUYINE")
- Thread starter
- #4
NickFort
Technical User
- Jun 10, 2010
- 113
- 0
- 0
Yes -- what you can do is take the file with all of the subroutines in it, put "module <your module name>" at the very beginning, then "end module <your module name>" at the very end. In your driver program, BEFORE any declarations, just add "use <your module name>". Now you can call the subroutines from the main program. One thing to note: you have to compile the module BEFORE the driver program. So, e.g.
If you specify the .f90 files in the reverse order, you'll get a compiling error. (Disclaimer: I've never used ifort, but I imagine it's the same with all compilers)
--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.
Code:
ifort your_module_name.f90 your_driver.f90 -o my_progIf you specify the .f90 files in the reverse order, you'll get a compiling error. (Disclaimer: I've never used ifort, but I imagine it's the same with all compilers)
--------------------------------------
Background: Chemical engineer, familiar mostly with MATLAB, but now branching out into real programming.
- Thread starter
- #6
- Thread starter
- #7
I got problems:
milenko@milenkons:~/ircg$ ifort -c mt2ddib1.for
mt2ddib1.for(23): error #6274: This statement must not appear in the specification part of a module
open(2,file='test_p7.dat',status='old')
--------^
mt2ddib1.for(29): error #6274: This statement must not appear in the specification part of a module
open(3,file='test_sh7.dat',status='old')
--------^
mt2ddib1.for(30): error #6274: This statement must not appear in the specification part of a module
open(4,file='test_res.dat',status='old')
--------^
mt2ddib1.for(35): error #6274: This statement must not appear in the specification part of a module
open(50,file='e.dat',status='unknown')
--------^
mt2ddib1.for(36): error #6274: This statement must not appear in the specification part of a module
open(51,file='esens.dat',status='unknown')
--------^
mt2ddib1.for(37): error #6274: This statement must not appear in the specification part of a module
open(52,file='h.dat',status='unknown')
--------^
milenko@milenkons:~/ircg$ ifort -c mt2ddib1.for
mt2ddib1.for(23): error #6274: This statement must not appear in the specification part of a module
open(2,file='test_p7.dat',status='old')
--------^
mt2ddib1.for(29): error #6274: This statement must not appear in the specification part of a module
open(3,file='test_sh7.dat',status='old')
--------^
mt2ddib1.for(30): error #6274: This statement must not appear in the specification part of a module
open(4,file='test_res.dat',status='old')
--------^
mt2ddib1.for(35): error #6274: This statement must not appear in the specification part of a module
open(50,file='e.dat',status='unknown')
--------^
mt2ddib1.for(36): error #6274: This statement must not appear in the specification part of a module
open(51,file='esens.dat',status='unknown')
--------^
mt2ddib1.for(37): error #6274: This statement must not appear in the specification part of a module
open(52,file='h.dat',status='unknown')
--------^
- Thread starter
- #9
I have changed but it doesn't work!Example:
module resph
use constants
implicit none
integer:: is,iprd0,k
real:: zxy_re_e,zxy_im_e,zxy_re_h,zxy_im_h
real:: ro_e,ro_h,phase_e,phase_h,period
do n=1,20
read(10,150)period,is,iprd0,zxy_re_e,zxy_im_e
150 format(f12.4,2i5,2e15.6)
read(12,152)period,is,iprd0,zxy_re_h,zxy_im_h
152 format(f12.4,2i5,2e15.6)
k=period/(2*mu0)
ro_e=k*(zxy_re_e**2+zxy_im_e**2)
phase_e=atan(zxy_im_e/zxy_re_e)
ro_h=k*(zxy_re_h**2+zxy_im_h**2)
phase_h=atan(zxy_im_h/zxy_re_h)
write(20,175)ro_e,phase_e,ro_h,phase_h
175 format(f11.6,2x,f11.8,2x,f11.6,2x,f11.8)
end do
end module
resph.f90(10): error #6274: This statement must not appear in the specification part of a module
do n=1,20
-------^
resph.f90(11): error #6274: This statement must not appear in the specification part of a module
read(10,150)period,is,iprd0,zxy_re_e,zxy_im_e
-------^
resph.f90(13): error #6274: This statement must not appear in the specification part of a module
read(12,152)period,is,iprd0,zxy_re_h,zxy_im_h
-------^
resph.f90(15): error #6274: This statement must not appear in the specification part of a module
k=period/(2*mu0)
module resph
use constants
implicit none
integer:: is,iprd0,k
real:: zxy_re_e,zxy_im_e,zxy_re_h,zxy_im_h
real:: ro_e,ro_h,phase_e,phase_h,period
do n=1,20
read(10,150)period,is,iprd0,zxy_re_e,zxy_im_e
150 format(f12.4,2i5,2e15.6)
read(12,152)period,is,iprd0,zxy_re_h,zxy_im_h
152 format(f12.4,2i5,2e15.6)
k=period/(2*mu0)
ro_e=k*(zxy_re_e**2+zxy_im_e**2)
phase_e=atan(zxy_im_e/zxy_re_e)
ro_h=k*(zxy_re_h**2+zxy_im_h**2)
phase_h=atan(zxy_im_h/zxy_re_h)
write(20,175)ro_e,phase_e,ro_h,phase_h
175 format(f11.6,2x,f11.8,2x,f11.6,2x,f11.8)
end do
end module
resph.f90(10): error #6274: This statement must not appear in the specification part of a module
do n=1,20
-------^
resph.f90(11): error #6274: This statement must not appear in the specification part of a module
read(10,150)period,is,iprd0,zxy_re_e,zxy_im_e
-------^
resph.f90(13): error #6274: This statement must not appear in the specification part of a module
read(12,152)period,is,iprd0,zxy_re_h,zxy_im_h
-------^
resph.f90(15): error #6274: This statement must not appear in the specification part of a module
k=period/(2*mu0)
- Thread starter
- #10
Of course, the compiler is right with the message This statement must not appear in the specification part of a module
You forgot CONTAINS.
Read first something about modules - for example here:
You forgot CONTAINS.
Read first something about modules - for example here:
- Thread starter
- #12
It works fine for one subroutine but when I try something like this,I get :m0.for(13): error #6274: This statement must not appear in the specification part of a module
call readmodel
------^
m0.for(14): error #6274: This statement must not appear in the specification part of a module
call solve_mt2d_directsens
module m0
use constants
use settings
use params
use mt2dmod
use fdsystem
use mt2ddat
use mt2dsens
implicit none
call readmodel
call solve_mt2d_directsens
contains
subroutine readmodel
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
c
do i=ma,m1
read(3,*)(ic(i,j),j=1,n1)
enddo
icd=1
res(icd)=-1.
cond(icd)=0.
ivar(icd)=0
c
c Read the actual number of different conductivities used in the conductivity
c map. This number is for conductivities in the Earth only, the air domain is
c NOT included here!
read(3,*)nc
nc=nc+1
ivarprd=0
c
c Read the resistivities of the uniform domains in the 2D model. The index of
c the particular domain is read first (corresponding to the index in the
c conductivity map above, IC), then the resistivity of the domain, in ohm*m,
c and finally a flag (0/1) which indicates whether the sensitivity w.r.t. this
c domain's conductivity is to be computed.
do i=2,nc
read(4,*)icd,res(icd),ivar(icd)
cond(icd)=1./res(icd)
if(ivar(icd).ne.0)then
ivarprd=ivarprd+1
ivarp(ivarprd)=icd
endif
enddo
close (4)
c Conversion of units (SI or practical) for the impedances
if(lsiunits)then
prev_zunit=prev_siunit
else
prev_zunit=prev_pracunit
endif
end subroutine
c
c Solution of the direct 2D problem and sensitivity evaluation
c call solve_mt2d_directsens
subroutine solve_mt2d_directsens
c ================================
c Subroutine for the solution of a 2D MT model (isotropic) and for the
c sensitivity evaluations w.r.t. the conductivities of selected homogeneous
c domains of the model
c
use constants
use params
use mt2dmod
use fdsystem
use mt2ddat
use mt2dsens
c
implicit none
c
integer*4 iper,isite,i,ico,iprd
real*4 period,prunit
complex*8 hxsurf
c
c Compute boundary models for the 1D margins at the left and right hand
c sides of the 2D model. The boundary models are the same for all periods
c if the 2D model does not change
call gebolr
c
c Loop over all periods considered
do iper=1,np
period=per(iper)
c
c Solve the 2D model for the E-mode case for one period
call d2emod(iper)
c
c Evaluate the sensitivities of the 2D solution w.r.t. all the required
c conductivities for the E-mode case for one period
call sensd2emod(iper)
c
c Solve the 2D for the H-mode case for one period
call d2hmod(iper,hxsurf)
c
c Evaluate the sensitivities of the 2D solution w.r.t. all the required
c conductivities for the H-mode case for one period
call sensd2hmod(iper,hxsurf)
c
c Store the MT functions for the particular period in a larger array of
c all the MT functions for all periods. Here, the impedances are output
c at points of interest only
do i=1,nsites
isite=isit(i)
zmd(iper,i,1)=zmd1p(i,1)
zmd(iper,i,2)=zmd1p(i,2)
zmd(iper,i,3)=zmd1p(i,3)
zmd(iper,i,4)=zmd1p(i,4)
c
c Auxiliary output, file 81: E and H-mode impedances at sites of interest
c for one period
c write(81,'(2i5,5x,4e15.6)')iper,isite,
c & zmd(iper,i,1),zmd(iper,i,2),
c & zmd(iper,i,3),zmd(iper,i,4)
c
enddo
c
c Complete array of all required sensitivities w.r.t. flagged domains at
c points of interest for all periods. For larger models, this array may
c become prohibitely large!!!
do i=1,nsites
isite=isit(i)
do ico=1,nc
if(ivarp(ico).ne.0)then
iprd=ivarp(ico)
zmdsens(iper,i,1,iprd)=zmd1psens(i,1,iprd)
zmdsens(iper,i,2,iprd)=zmd1psens(i,2,iprd)
zmdsens(iper,i,3,iprd)=zmd1psens(i,3,iprd)
zmdsens(iper,i,4,iprd)=zmd1psens(i,4,iprd)
c
c Auxiliary output, file 82: E and H-mode impedance sensitivities at sites
c of interest for one period
c write(82,'(3i5,4e15.6)')iper,isite,iprd,
c & zmdsens(iper,i,1,iprd),zmdsens(iper,i,2,iprd),
c & zmdsens(iper,i,3,iprd),zmdsens(iper,i,4,iprd)
c
endif
enddo
enddo
c
c Go on to the next period
enddo
c
return
end subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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)
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)
c write(*,*)zmd1p(i,1),zmd1p(i,2)
c
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)
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 subroutine
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
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)
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
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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
end module m0
call readmodel
------^
m0.for(14): error #6274: This statement must not appear in the specification part of a module
call solve_mt2d_directsens
module m0
use constants
use settings
use params
use mt2dmod
use fdsystem
use mt2ddat
use mt2dsens
implicit none
call readmodel
call solve_mt2d_directsens
contains
subroutine readmodel
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
c
do i=ma,m1
read(3,*)(ic(i,j),j=1,n1)
enddo
icd=1
res(icd)=-1.
cond(icd)=0.
ivar(icd)=0
c
c Read the actual number of different conductivities used in the conductivity
c map. This number is for conductivities in the Earth only, the air domain is
c NOT included here!
read(3,*)nc
nc=nc+1
ivarprd=0
c
c Read the resistivities of the uniform domains in the 2D model. The index of
c the particular domain is read first (corresponding to the index in the
c conductivity map above, IC), then the resistivity of the domain, in ohm*m,
c and finally a flag (0/1) which indicates whether the sensitivity w.r.t. this
c domain's conductivity is to be computed.
do i=2,nc
read(4,*)icd,res(icd),ivar(icd)
cond(icd)=1./res(icd)
if(ivar(icd).ne.0)then
ivarprd=ivarprd+1
ivarp(ivarprd)=icd
endif
enddo
close (4)
c Conversion of units (SI or practical) for the impedances
if(lsiunits)then
prev_zunit=prev_siunit
else
prev_zunit=prev_pracunit
endif
end subroutine
c
c Solution of the direct 2D problem and sensitivity evaluation
c call solve_mt2d_directsens
subroutine solve_mt2d_directsens
c ================================
c Subroutine for the solution of a 2D MT model (isotropic) and for the
c sensitivity evaluations w.r.t. the conductivities of selected homogeneous
c domains of the model
c
use constants
use params
use mt2dmod
use fdsystem
use mt2ddat
use mt2dsens
c
implicit none
c
integer*4 iper,isite,i,ico,iprd
real*4 period,prunit
complex*8 hxsurf
c
c Compute boundary models for the 1D margins at the left and right hand
c sides of the 2D model. The boundary models are the same for all periods
c if the 2D model does not change
call gebolr
c
c Loop over all periods considered
do iper=1,np
period=per(iper)
c
c Solve the 2D model for the E-mode case for one period
call d2emod(iper)
c
c Evaluate the sensitivities of the 2D solution w.r.t. all the required
c conductivities for the E-mode case for one period
call sensd2emod(iper)
c
c Solve the 2D for the H-mode case for one period
call d2hmod(iper,hxsurf)
c
c Evaluate the sensitivities of the 2D solution w.r.t. all the required
c conductivities for the H-mode case for one period
call sensd2hmod(iper,hxsurf)
c
c Store the MT functions for the particular period in a larger array of
c all the MT functions for all periods. Here, the impedances are output
c at points of interest only
do i=1,nsites
isite=isit(i)
zmd(iper,i,1)=zmd1p(i,1)
zmd(iper,i,2)=zmd1p(i,2)
zmd(iper,i,3)=zmd1p(i,3)
zmd(iper,i,4)=zmd1p(i,4)
c
c Auxiliary output, file 81: E and H-mode impedances at sites of interest
c for one period
c write(81,'(2i5,5x,4e15.6)')iper,isite,
c & zmd(iper,i,1),zmd(iper,i,2),
c & zmd(iper,i,3),zmd(iper,i,4)
c
enddo
c
c Complete array of all required sensitivities w.r.t. flagged domains at
c points of interest for all periods. For larger models, this array may
c become prohibitely large!!!
do i=1,nsites
isite=isit(i)
do ico=1,nc
if(ivarp(ico).ne.0)then
iprd=ivarp(ico)
zmdsens(iper,i,1,iprd)=zmd1psens(i,1,iprd)
zmdsens(iper,i,2,iprd)=zmd1psens(i,2,iprd)
zmdsens(iper,i,3,iprd)=zmd1psens(i,3,iprd)
zmdsens(iper,i,4,iprd)=zmd1psens(i,4,iprd)
c
c Auxiliary output, file 82: E and H-mode impedance sensitivities at sites
c of interest for one period
c write(82,'(3i5,4e15.6)')iper,isite,iprd,
c & zmdsens(iper,i,1,iprd),zmdsens(iper,i,2,iprd),
c & zmdsens(iper,i,3,iprd),zmdsens(iper,i,4,iprd)
c
endif
enddo
enddo
c
c Go on to the next period
enddo
c
return
end subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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)
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)
c write(*,*)zmd1p(i,1),zmd1p(i,2)
c
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,
=vendif
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)
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 subroutine
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
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)
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,
=vendif
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
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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 subroutine
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
end module m0
It's because you cannot call in module a subroutine. A module can only contain declarations of variables, functions and subroutines.milenko76 said:
And btw:
[ol]
[li]Your code seems to be buggy - I cannot see end module in the code you posted[/li]
[li]Next time post your source code betwwen [ignore]
Code:
[/ignore][/b] and [b][ignore][/ol]
- Status
