Commit 94f5a99b authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

debugging

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: 1755
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 1096960c
c this is <readanigemini.f>
c ----------------------------------------------------------------------------
c ($Id: readanigemini.f,v 1.1 2005-06-16 11:18:12 tforb Exp $)
c ($Id: readanigemini.f,v 1.2 2005-06-16 11:43:30 tforb Exp $)
c
c Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach)
c
......@@ -88,25 +88,23 @@ c S-velo-hori, Qmu, Qkappa, eta
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)
write(88,*) 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)
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)
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)')
write(88,'(1x)')
enddo
......
c this is <gemmodpg.f>
c
c $Id: gemmodpg.f,v 1.2 2000-04-14 15:28:08 thof Exp $
c $Id: gemmodpg.f,v 1.3 2005-06-16 11:43:26 tforb Exp $
c
c plot GEMINI earth model
c
c V1.0 10/01/97 Thomas Forbriger
c V1.1 15/01/97 define radius range
c V1.2 14/04/00 use up-to-date commandline function
c V1.3 16/07/2005 plot files defining transverse isotropy
c
c----------------------------------------------------------------------
character*70 version
parameter(version='GEMMODPG V1.1 plot gemini earth model')
parameter(version='GEMMODPG V1.3 plot gemini earth model')
c gemini
integer maxlayer, nlayer
parameter(maxlayer=30)
double precision alpha(maxlayer, 4), beta(maxlayer, 4)
double precision alphah(maxlayer, 4), betah(maxlayer, 4)
double precision rho(maxlayer, 4), qm(maxlayer), qk(maxlayer)
double precision eta(maxlayer, 4)
double precision rb(0:maxlayer)
integer iflso(maxlayer), nco(maxlayer)
character text*72
c
character*80 filename, device, title
logical debug
logical debug, optani
integer fin, iargc, lu
parameter(lu=10)
real rmin, rmax
integer i,j
c commandline
integer maxopt, lastarg, iargc
parameter(maxopt=3)
parameter(maxopt=4)
character*2 optid(maxopt)
character*40 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
c here are the keys to our commandline options
data optid/2h-d,2h-D,2h-r/
data opthasarg/.TRUE.,.FALSE.,.TRUE./
data optarg/3hx11,1h-,7h0.,1.e4/
data optid/2h-d,2h-D,2h-r,2h-t/
data opthasarg/.TRUE.,.FALSE.,.TRUE.,.false./
data optarg/3hx11,1h-,7h0.,1.e4,1h-/
c----------------------------------------------------------------------
print *,version
print *,'Usage: gemmodpg [ -d device ] [ -r rmin,rmax ] file'
print *,'Usage: gemmodpg [-t] [ -d device ] [ -r rmin,rmax ] file'
print *,' or: gemmodpg -help'
if (iargc().lt.1) stop 'ERROR: missing parameters'
call getarg(1, filename)
......@@ -46,6 +50,7 @@ c----------------------------------------------------------------------
print *,' '
print *,'-d device any pgplot device'
print *,'-r rmin,rmax plot range from radius rmin to radius rmax'
print *,'-t earth model has transverse isotropy'
print *,' '
call pgp_showdevices
stop
......@@ -57,31 +62,72 @@ c set options
device=optarg(1)
debug=optset(2)
read(optarg(3), *) rmin,rmax
optani=optset(4)
if (iargc().eq.(lastarg)) stop 'ERROR: no file'
call getarg(iargc(), filename)
fin=index(filename,' ')
c----------------------------------------------------------------------
c read file
if (debug) print *,'DEBUG: filename ',filename(1:fin)
call gemini_getmod(filename, lu, maxlayer,
& rb, qm, qk, rho, alpha, beta, nlayer, iflso, nco, text)
if (optani) then
call gemini_getani(filename, lu, maxlayer, eta,
& rb, qm, qk, rho, alpha, alphah, beta, betah,
& nlayer, iflso, nco, text)
else
call gemini_getmod(filename, lu, maxlayer,
& rb, qm, qk, rho, alpha, beta, nlayer, iflso, nco, text)
endif
c----------------------------------------------------------------------
if (debug) then
do i=1,nlayer
print *,alpha(i,1),alphah(i,1),beta(i,1),
& betah(i,1),eta(i,1),rho(i,1),qm(i),qk(i)
do j=2,nco(i)
print *,alpha(i,j),alphah(i,j),beta(i,j),
& betah(i,j),eta(i,j),rho(i,j)
enddo
print *,' '
enddo
endif
c----------------------------------------------------------------------
c check ranges
rmin=max(rmin,sngl(rb(0)))
rmax=min(rmax,sngl(rb(nlayer)))
c----------------------------------------------------------------------
c start plot
call pgp_setdevice(device, 5, 0)
if (.false.) then
call pgp_setdevice(device, 4, 2)
else
call pgp_setdevice(device, 5, 0)
endif
call pgsch(2.3)
c plot velocities
title='P-velocity [km/s]'
call plotcurve(maxlayer, nlayer, rb, alpha, nco, iflso, rmin, rmax,
& title, debug)
call pgmtxt('T',1.,0.,0.,text)
title='S-velocity [km/s]'
call plotcurve(maxlayer, nlayer, rb, beta, nco, iflso, rmin, rmax,
& title, debug)
title='density [g/cm^3]'
if (.false.) then
title='Pv-velocity (km/s)'
call plotcurve(maxlayer, nlayer, rb, alpha, nco, iflso, rmin, rmax,
& title, debug)
title='Ph-velocity (km/s)'
call plotcurve(maxlayer, nlayer, rb, alphah, nco, iflso, rmin, rmax,
& title, debug)
title='Sv-velocity (km/s)'
call plotcurve(maxlayer, nlayer, rb, beta, nco, iflso, rmin, rmax,
& title, debug)
title='Sh-velocity (km/s)'
call plotcurve(maxlayer, nlayer, rb, betah, nco, iflso, rmin, rmax,
& title, debug)
title='eta'
call plotcurve(maxlayer, nlayer, rb, eta, nco, iflso, rmin, rmax,
& title, debug)
else
title='P-velocity (km/s)'
call plotcurve(maxlayer, nlayer, rb, alpha, nco, iflso, rmin, rmax,
& title, debug)
title='S-velocity (km/s)'
call plotcurve(maxlayer, nlayer, rb, beta, nco, iflso, rmin, rmax,
& title, debug)
endif
title='density (g/cm^3)'
call plotcurve(maxlayer, nlayer, rb, rho, nco, iflso, rmin, rmax,
& title, debug)
title='Qmu'
......@@ -128,7 +174,7 @@ c
if (debug) print *,'DEBUG: go plot xmin/xmax',xmin,xmax
call pgslw(1)
call pgenv(xmin, xmax, ymin, ymax, 0, 2)
call pglab(title, 'radius [km]', ' ')
call pglab(title, 'radius (km)', ' ')
call pgslw(4)
call pgline(nsteps, x, y)
call pgslw(1)
......@@ -169,7 +215,7 @@ c
if (debug) print *,'DEBUG: go plot xmin/xmax',xmin,xmax
call pgslw(1)
call pgenv(xmin, xmax, ymin, ymax, 0, 2)
call pglab(title, 'radius [km]', ' ')
call pglab(title, 'radius (km)', ' ')
call pgslw(4)
call pgline(nsteps, x, y)
call pgslw(1)
......
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