Commit 1096960c authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

read transversal isotropic parameters

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 1754
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 046af724
# this is <Makefile>
# ----------------------------------------------------------------------------
# $Id: Makefile,v 1.2 2005-06-16 11:18:12 tforb Exp $
#
# Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach)
#
# compile earth model library
#
# REVISIONS and CHANGES
# 16/06/2005 V1.0 Thomas Forbriger
#
# ============================================================================
#
all:
flist: Makefile $(wildcard *.f)
echo $^ | tr ' ' '\n' | sort > $@
.PHONY: edit
edit: flist; vim $<
.PHONY: clean
clean: ;
-find . -name \*.bak | xargs --no-run-if-empty /bin/rm -v
-/bin/rm -vf flist *.o
#
# Makefile for prog/lib/libemod
#
......@@ -9,7 +36,8 @@ F2CFLAGS=-f -u
#CC=gcc
CFLAGS=-O2 -I${SERVERINCLUDEDIR} -I${LOCINCLUDEDIR}
LIBEMODSUB=psexdat.o readgemini.o efa.o readrefmet.o earthmod.o
LIBEMODSUB=psexdat.o readgemini.o efa.o readrefmet.o earthmod.o \
readanigemini.o
docs: $(DOCS)
......@@ -18,7 +46,7 @@ earthmod.o: earthmod.f
$(CC) $(CFLAGS) $(<:.f=.c) -c
@rm $(<:.f=.c)
readgemini.o: readgemini.f
readgemini.o readanigemini.o: %.o: %.f
f2c -f $<
$(CC) $(CFLAGS) $(<:.f=.c) -c
@rm $(<:.f=.c)
......@@ -33,9 +61,6 @@ psexdat.o: psexdat.f
$(CC) $(CFLAGS) $(<:.f=.c) -c
@rm $(<:.f=.c)
clean:
-/bin/rm *.o *.bak
libemod.a: $(LIBEMODSUB)
ar rcv libemod.a $(LIBEMODSUB)
ranlib libemod.a
......@@ -43,3 +68,5 @@ libemod.a: $(LIBEMODSUB)
efa.doc: efa.f
extractdoc.tcl efa.f > efa.doc
# ----- END OF Makefile -----
c this is <readanigemini.f>
c ----------------------------------------------------------------------------
c ($Id: readanigemini.f,v 1.1 2005-06-16 11:18:12 tforb Exp $)
c
c Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach)
c
c read transversal isotropic gemini model
c
c REVISIONS and CHANGES
c 16/06/2005 V1.0 Thomas Forbriger
c
c ============================================================================
c
subroutine gemini_getani(filename, lu, maxlayer, eta,
& rb, qm, qk, rho, vpv, vph, vsv, vsh, nlayer, iflso, nco, text)
c
c this routines reads a gemini model file
c
character filename*(*), text*72
integer nlayer, maxlayer, lu
double precision qm(maxlayer), qk(maxlayer), rho(maxlayer, 4)
double precision vpv(maxlayer, 4), vsv(maxlayer, 4)
double precision vph(maxlayer, 4), vsh(maxlayer, 4)
double precision eta(maxlayer, 4)
double precision rb(0:maxlayer)
integer iflso(maxlayer), nco(maxlayer), nlay
cE
character cnlay*2,form*11
double precision fref, rboc
integer i,j,ifanis
c
open(lu,file=filename,status='old',err=99)
c
c------------------------------------------------------
c code was taken from:
c
c G T M A N I S
c
c Reads tranversal isotropic Earth model parameters described by polynomials
c from file.
c See comments at the read-statements for the meaning of the parameter.
c
c------------------------------------------------------
do i=1,maxlayer
qm(i)=0.d0
qk(i)=0.d0
rb(i)=0.d0
iflso(i)=0
nco(i)=0
do j=1,4
rho(i,j)=0.d0
vpv(i,j)=0.d0
vph(i,j)=0.d0
vsv(i,j)=0.d0
vsh(i,j)=0.d0
eta(i,j)=0.d0
enddo
enddo
c Read header with comments and print them on stdout
c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1 read(lu,'(a72)') text
if (text(1:1).eq.'#') then
c print '(a72)', text
goto 1
endif
backspace lu
read(lu,'(i2)') nlayer ! Number of layers
if (nlayer.gt.maxlayer) stop 'ERROR: too many layers!'
write(cnlay,'(i2)') nlayer !
form='('//cnlay//'i2)' ! Number of polynomal
read(lu,form) (nco(i),i=1,nlayer) ! coefficients for each layer
read(lu,*) fref ! reference frequency of Qs in Hertz
read(lu,*) ifanis ! Transversal isotropic? 1=y, else=n
read(lu,'(1x/1x/)')
c Reading radii, density, P-velo-vert, P-velo-hori, S-velo-vert,
c S-velo-hori, Qmu, Qkappa, eta
do i = 1, nlayer
read(lu,*) rb(i-1),rho(i,1),vpv(i,1),vph(i,1),vsv(i,1),vsh(i,1),
& qm(i),qk(i),eta(i,1)
c write(88,*) rb(i-1),rho(i,1),vpv(i,1),vph(i,1),vsv(i,1),vsh(i,1),
c & qm(i),qk(i),eta(i,1)
qk(i)=1./qk(i) ! reciprocal quality factors
if (qm(i).le.0.) then !
qm(i) = 0. !
iflso(i)=1 !
nloc = i ! layer number of outer core
else !
qm(i)=1./qm(i) !
iflso(i)=0 !
endif !
do j = 2, nco(i) ! polynomal coefficients
read(lu,*) rho(i,j),vpv(i,j),vph(i,j),vsv(i,j),vsh(i,j),eta(i,j)
c write(88,*) rho(i,j),vpv(i,j),vph(i,j),vsv(i,j),vsh(i,j),eta(i,j)
enddo
read(lu,'(1x)')
c write(88,'(1x)')
enddo
read(lu,*) rb(nlay) ! Get Earth radius
c Check for ocean.
if(iflso(nlay).eq.1) stop '~~~~ I don`t like oceans. ~~~~'
rboc = rb(nloc) ! Radius of outer core
close(lu)
return
99 stop 'ERROR (gemini_getani): opening file'
end
c
c ----- END OF readanigemini.f -----
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment