readanigemini.f 3.82 KB
Newer Older
1
2
c this is <readanigemini.f>
c ----------------------------------------------------------------------------
thomas.forbriger's avatar
thomas.forbriger committed
3
c   ($Id: readanigemini.f,v 1.4 2005-06-16 12:57:47 tforb Exp $)
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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
thomas.forbriger's avatar
thomas.forbriger committed
19
      character filename*(*), text*(*)
20
21
22
23
24
25
      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)
thomas.forbriger's avatar
thomas.forbriger committed
26
      integer iflso(maxlayer), nco(maxlayer)
27
28
29
30
cE
      character cnlay*2,form*11
      double precision fref, rboc
      integer i,j,ifanis
thomas.forbriger's avatar
thomas.forbriger committed
31
      character junk*80
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

thomas.forbriger's avatar
thomas.forbriger committed
65
66
67
68
      text='-'
  1   read(lu,'(a72)') junk
      if (text.eq.'-') text=junk
      if (junk(1:1).eq.'#') then
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
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)
thomas.forbriger's avatar
thomas.forbriger committed
94
95
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)
96
97
98
99
100
101
102
103
104
105
106

        if (qm(i).le.0.) then          !
            qm(i) = 0.                 !
            iflso(i)=1                 !
            nloc = i                   ! layer number of outer core
        else                           ! 
            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)
thomas.forbriger's avatar
thomas.forbriger committed
107
c       write(88,*) rho(i,j),vpv(i,j),vph(i,j),vsv(i,j),vsh(i,j),eta(i,j)
108
109
        enddo
        read(lu,'(1x)')
thomas.forbriger's avatar
thomas.forbriger committed
110
c       write(88,'(1x)')
111
112
113
114

      enddo


thomas.forbriger's avatar
thomas.forbriger committed
115
      read(lu,*) rb(nlayer)           ! Get Earth radius
116
117

c  Check for ocean.
thomas.forbriger's avatar
thomas.forbriger committed
118
      if(iflso(nlayer).eq.1) stop '~~~~ I don`t like oceans. ~~~~'
119
120
121
122
123
124
125
126
      rboc = rb(nloc)                  ! Radius of outer core

      close(lu)
      return
   99 stop 'ERROR (gemini_getani): opening file'
      end
c
c ----- END OF readanigemini.f -----