global_main.f90 22.3 KB
Newer Older
capdevil's avatar
capdevil committed
1
2
3
4
5
6
7
!------------------------------------------------------
module global_main_param
!------------------------------------------------------
  use def_gparam
  implicit none
  public:: recepteur_dat  &
         ,nms_dat,NBT,dt,t0,f1,f2,f3,f4,RA,f0    &
capdevil's avatar
   
capdevil committed
8
         ,eigenfileT,eigenfileS,amp_source   &
capdevil's avatar
capdevil committed
9
10
11
         ,get_NBE,init_global_main,GEOC,source_type &
         ,geocentric_correction,nbcomp,NBR,NBE,stations_name      &
         ,get_source_dat,sources_name,logunit,euler&
yann's avatar
yann committed
12
         ,coord_stations,get_receivers,load_source,czero,cone     &
capdevil's avatar
capdevil committed
13
         ,set_logunit_in_global_main,t_delay_sources &
capdevil's avatar
   
capdevil committed
14
         ,SEXP,M0,source_delay,comp,compr,half_duration    &
15
         ,local_backup_dir,reset_sources_delay,get_Mnt,local_grad &
capdevil's avatar
   
capdevil committed
16
         ,load_receivers,get_NBR,rotation ,stacking_sources,nmax_in,nmin &
capdevil's avatar
   
capdevil committed
17
         ,source_force,station_response,response_constant,zeros,poles    &
18
         ,NBZ,NBP,channel,DEP,VEL,ACC,GRA,GRA_COR,homogenization,s2,is2,II,NBF2 
capdevil's avatar
capdevil committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
  private
!configuration files:
  character(len=30)            :: source_dat         ='sources.dat'
  character(len=30)            :: recepteur_dat      ='receivers.dat'
  character(len=30), parameter :: nms_dat            ='nms.dat'
  character(len=1), dimension(3), parameter  :: compr=(/'Z','R','T'/)
  character(len=1), dimension(3), parameter  :: comp=(/'Z','N','E'/)
!sources, stations, 
  integer :: NBE
  integer :: NBR
  integer, parameter :: nbcomp=3!nombre de composante (1== veritcal, 3 == all)
  integer :: iNNmax
  real(DP), dimension(:,:), allocatable :: coord_stations,Mtmp
  real(DP), dimension(:,:), allocatable :: coord_sources
capdevil's avatar
   
capdevil committed
33
  real(DP), dimension(:  ), allocatable :: t_delay_sources,half_duration
capdevil's avatar
capdevil committed
34
  logical :: geoc_corr,rotation
35
  logical :: local_grad
capdevil's avatar
capdevil committed
36
  character(len=8), dimension(: ), allocatable :: sources_name
capdevil's avatar
   
capdevil committed
37
  logical :: sources_available=.false.,stacking_sources,source_delay,source_force
capdevil's avatar
capdevil committed
38
39
  doubleprecision, parameter :: SEXP=1.d18
  character(len=4), dimension(: ), allocatable :: stations_name
40
!character(len=6), dimension(: ), allocatable :: stations_name
capdevil's avatar
capdevil committed
41
!temps:
capdevil's avatar
   
capdevil committed
42
43
  integer :: NBT
  real(DP):: dt,t0,amp_source
capdevil's avatar
capdevil committed
44
45
46
47
48
49
50
51
52
  real(DP):: f1,f2,f3,f4,f0 
  character(len=6), parameter :: source_type='heavis'
  logical :: local_backup
  character(len=100)  :: local_backup_dir
!frequence
  integer  :: NBF,NBF2 
   real(DP):: dtf,duree,deltaf,dtf2
  character(len=100) :: eigenfileT,eigenfileS
!
capdevil's avatar
   
capdevil committed
53
54
55
  real(DP), parameter ::  GEOC=0.993277d0
!RA will be reinitialized by util_fctp:
  real(DP) :: RA=-1.d0
capdevil's avatar
capdevil committed
56
57
58
59
60
61
  integer :: nmax_in,nmin
!
!logical unit for log file:
  integer :: logunit
!pour les FFT
  integer :: coef_nbt
capdevil's avatar
   
capdevil committed
62
  integer :: homogenization=-1 ! pas : -1, order 0 : 0, order 1 : 1
capdevil's avatar
capdevil committed
63
64
65
66
67
!station_response
  logical :: station_response
  integer, dimension(:), allocatable :: NBZ,NBP
  complex*16, dimension(:,:), allocatable :: zeros,poles
  doubleprecision, dimension(:), allocatable :: response_constant
capdevil's avatar
   
capdevil committed
68
!
69
  integer, parameter :: DEP=1,VEL=2,ACC=3,GRA=5,GRA_COR=4
capdevil's avatar
   
capdevil committed
70
  integer :: channel
yann's avatar
yann committed
71
72
73
74
  real(DP), parameter :: s2=sqrt(2._DP),is2=1._DP/s2
  complex(DP), parameter :: II=cmplx(0._DP,1._DP)
  complex(DP), parameter :: czero=cmplx(0._DP,0._DP)
  complex(DP), parameter :: cone =cmplx(1._DP,0._DP)
capdevil's avatar
capdevil committed
75
76
77
78
79
80
81
82
83
84
85
86
!
contains 
!---------------------------------------------------
  subroutine init_global_main(logunit_)
!---------------------------------------------------
    implicit none
    integer, intent(in) :: logunit_
!
    logunit=logunit_
!
    call load_nms()
    call init_rcv_src()
capdevil's avatar
capdevil committed
87
    if (station_response) call load_station_response()
capdevil's avatar
capdevil committed
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
!---------------------------------------------------
  end subroutine init_global_main
!---------------------------------------------------
!---------------------------------------------------
  subroutine set_logunit_in_global_main(log)
!---------------------------------------------------
    implicit none
    integer, intent(in) :: log
    logunit=log
!---------------------------------------------------
  end subroutine set_logunit_in_global_main
!---------------------------------------------------
!-------------------------------------------------------------------------
    logical function geocentric_correction()
!-------------------------------------------------------------------------
      implicit none
      geocentric_correction=geoc_corr
!-------------------------------------------------------------------------
    end function geocentric_correction
!-------------------------------------------------------------------------
!---------------------------------------------------
  subroutine load_nms()
!---------------------------------------------------
    implicit none
    integer, parameter:: unit=13
    integer :: i
    
!
!    nbcomp=1
    open(unit,file=nms_dat,status='old')
    read(unit,*)
    read(unit,*) dt
    read(unit,*)
    read(unit,*) NBT
    read(unit,*)
    read(unit,*) t0
    read(unit,*)
capdevil's avatar
   
capdevil committed
125
126
    read(unit,*) amp_source
    read(unit,*) 
capdevil's avatar
capdevil committed
127
128
129
130
131
132
133
134
135
136
137
138
139
    read(unit,*) f1,f2,f3,f4
    read(unit,*)
    read(unit,*) nmin,nmax_in
    if (nmin<0) nmin=0
    if (nmax_in<0) nmax_in=-1
    read(unit,*)
    read(unit,*) geoc_corr
    read(unit,*)
    read(unit,*) rotation
    read(unit,*)
    read(unit,'(a)') eigenfileS
    read(unit,*)
    read(unit,'(a)') eigenfileT
capdevil's avatar
   
capdevil committed
140
141
    read(unit,*) 
    read(unit,*) source_force
capdevil's avatar
capdevil committed
142
    read(unit,*)     
capdevil's avatar
   
capdevil committed
143
144
    read(unit,*) channel
    read(unit,*)     
capdevil's avatar
capdevil committed
145
    read(unit,*) station_response
capdevil's avatar
   
capdevil committed
146
147
    read(unit,*)! #tag homogenization: -1: none; 0: 0 order , 1: first order
    read(unit,*) homogenization
capdevil's avatar
capdevil committed
148
149
150
151
152
153
154
155
156
    write(logunit,*) "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
    write(logunit,*) "time step:",dt
    write(logunit,*) "Number of time step:",NBT
    write(logunit,*) "Source t0 time:",t0
    write(logunit,*) "source frequency band:",sngl(f1),sngl(f2),sngl(f3),sngl(f4)
    write(logunit,*) "Geocentric correction?:",geoc_corr 
    write(logunit,*) "Egeinfunctions files:"
    write(logunit,'(a)') eigenfileS
    write(logunit,'(a)') eigenfileT
capdevil's avatar
   
capdevil committed
157
158
    select case (channel)
    case(DEP)
capdevil's avatar
   
capdevil committed
159
        write(logunit,*) 'Output channel is displacement (for a step source time function) '
160
        print*, 'Output channel is displacement (for a step source time function) '
capdevil's avatar
   
capdevil committed
161
    case(VEL)
capdevil's avatar
   
capdevil committed
162
        write(logunit,*) 'Output channel is velocity (for a step source time function) or displacement (for a Dirac source time function)'
163
        print*,'Output channel is velocity (for a step source time function) or displacement (for a Dirac source time function)'
capdevil's avatar
   
capdevil committed
164
    case(ACC)
capdevil's avatar
   
capdevil committed
165
        write(logunit,*) 'Output channel is acceleration (for step source time function) or velocity (for a Dirac source time function) '
166
        print*,'Output channel is acceleration (for step source time function) or velocity (for a Dirac source time function) '
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
167
    case(GRA)
168
169
170
171
172
        write(logunit,*) 'Output channel is gradient + rotational (for step source time function) or  rate (for a Dirac source time function) '
        print*,'Output channel is gradient + rotational (for step source time function) or  rate (for a Dirac source time function) '
    case(GRA_COR)
        write(logunit,*) 'Output channel is PARTIAL (no sphericity term) gradient + rotational (for step source time function) or rate (for a Dirac source time function) '
        print*,'Output channel is PARTIAL (no sphericity term) gradient + rotational (for step source time function) or rate (for a Dirac source time function) '
capdevil's avatar
   
capdevil committed
173
    case default
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
174
       STOP 'Channel must 1, 2, 3, 4 or 5'
capdevil's avatar
   
capdevil committed
175
176
    end select
    write(logunit,*) 'Corrector flag:',homogenization
177
178
179
180
181
182
183
    if (channel==GRA_COR) then
       channel=GRA
       local_grad=.true.
    else
       local_grad=.false.
    endif
!       
capdevil's avatar
capdevil committed
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
!---------------------------------------------------
  end subroutine load_nms
!---------------------------------------------------

!---------------------------------------------------
  subroutine get_NBE(NBE_)
!---------------------------------------------------
    implicit none
    integer, intent(out) :: NBE_
    NBE_=NBE
!---------------------------------------------------
  end subroutine get_NBE
!---------------------------------------------------
!---------------------------------------------------
  subroutine get_NBR(NBR_)
!---------------------------------------------------
    implicit none
    integer, intent(out) :: NBR_
    NBR_=NBR
!---------------------------------------------------
  end subroutine get_NBR
!---------------------------------------------------
!============SOURCES:
!-------------------------------------------------------------------------
  subroutine init_rcv_src()
!-------------------------------------------------------------------------
    implicit none
    logical :: err
    call load_source()
    call load_receivers(err)
    if (err) STOP 'Error when loading the receivers file'
!-------------------------------------------------------------------------
  end subroutine init_rcv_src
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

!-------------------------------------------------------------------------
  subroutine load_source()
!-------------------------------------------------------------------------
    implicit none
    character(len=8) :: name
    integer     :: i
!pour les corrections geocentrique 
    character(len=100) :: source_file
    integer :: units=527,j,ibid
!
!lecture du fichier de source:
    open(units,file=source_dat,status='old')
    read(units,*) 
    read(units,*) NBE
    read(units,*) 
    read(units,*) stacking_sources
    read(units,*) 
capdevil's avatar
   
capdevil committed
237
238
    allocate(coord_sources(3,NBE),Mtmp(6,NBE),sources_name(NBE),half_duration(NBE), &
            t_delay_sources(NBE)) 
capdevil's avatar
capdevil committed
239
240
    source_delay=.false.
    do i=1,NBE
capdevil's avatar
   
capdevil committed
241
242
243
!        read(units,'(i2.2,1x,a8,1x,6(e9.2,1x),3(f6.2,1x),f6.2,x,f6.2)') ibid,name, &
!           (Mtmp (j,i),j=1,6),(coord_sources(j,i),j=1,3),half_duration(i),t_delay_sources(i)
        read(units,*) ibid,name, &
capdevil's avatar
   
capdevil committed
244
           (Mtmp (j,i),j=1,6),(coord_sources(j,i),j=1,3),half_duration(i),t_delay_sources(i)
capdevil's avatar
capdevil committed
245
246
        sources_name(i)=name
        if (abs(t_delay_sources(i)) > 1.E-4_DP)  source_delay=.true.
capdevil's avatar
   
capdevil committed
247
!conversion latitude->colatitude      
capdevil's avatar
   
capdevil committed
248
       coord_sources(2,i)=90.d0-coord_sources(2,i) 
capdevil's avatar
capdevil committed
249
250
251
252
    enddo
    close(units)
    write(logunit,*) "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
    write(logunit,*) "Number of sources to used:",NBE
capdevil's avatar
   
capdevil committed
253
254
    do i=1,NBE
        write(logunit,'(i2.2,1x,a8,1x,6(e10.3,1x),3(f6.2,1x),f6.2,x,f6.2)') i,sources_name(i), &
capdevil's avatar
   
capdevil committed
255
                (Mtmp (j,i),j=1,6),coord_sources(1,i),90.d0-coord_sources(2,i),coord_sources(3,i),half_duration(i),t_delay_sources(i)
capdevil's avatar
   
capdevil committed
256
    enddo
capdevil's avatar
capdevil committed
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
    if (source_delay)  then
       write(logunit,*) " Sources delay mode (not all sources will triggered at t=0)"
    else
       write(logunit,*) " All triggered at t=0"
    endif
    if (rotation.and.stacking_sources) then
       print*,'Warning, can''t rotate component when stacking sources'
       rotation=.false.
    endif
!
!-------------------------------------------------------------------------
  end subroutine load_source
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
  subroutine get_Mnt(Mout,is)
!-------------------------------------------------------------------------
    implicit none
    real, dimension(:), intent(out):: Mout
    integer, intent(in) :: is
    Mout(:)=real(Mtmp(:,is))
!-------------------------------------------------------------------------
  end subroutine get_Mnt
!-------------------------------------------------------------------------


!-------------------------------------------------------------------------
  subroutine load_receivers(error)
!-------------------------------------------------------------------------
    implicit none
    logical, intent(out) :: error
    integer :: irec,unitrec
!    
!=========================================
!Lecture des donnees dans recepteur.dat
!=========================================
    unitrec=135
    open(unitrec, file=recepteur_dat,status='old',err=100)
    read(unitrec,*)      
    read(unitrec,*) NBR
    read(unitrec,*)      
    allocate(stations_name(NBR),coord_stations(2,NBR))
    do irec= 1,NBR
capdevil's avatar
   
capdevil committed
299
300
301
!       read(unitrec,'(a4,1x,f8.3,1x,f8.3)') stations_name(irec) &
!                ,coord_stations(1,irec),coord_stations(2,irec)
       read(unitrec,*) stations_name(irec) &
capdevil's avatar
capdevil committed
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
                ,coord_stations(1,irec),coord_stations(2,irec)
!latitude-> colatitude
       coord_stations(1,irec)=90.d0-coord_stations(1,irec)
    enddo

!Fermeture du fichier recepteur.dat
    close(unitrec)
    error=.false.
    write(logunit,*) "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
    write(logunit,*) "Number of receivers:",NBR
    return
100 print*,'No recepteur.dat file found! ... skeeping receivers management..'     
    error=.true.
!-------------------------------------------------------------------------
  end subroutine load_receivers
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
  subroutine write_receivers(weight)
!-------------------------------------------------------------------------
    implicit none
    integer :: irec,unitrec
    real(DP), dimension(:), optional :: weight
!    
!=========================================
!Ecriture des donnees dans recepteur_out.dat
!=========================================
    unitrec=135
    open(unitrec, file='recepteur_out.dat')
    write(unitrec,'(a)')'#nombre de recepteurs:'
    write(unitrec,*) NBR
    write(unitrec,'(a)') '#nom (4 characters), latitude, longitude en degres &
                     &, 1/poids de la station'
    if (.not.present(weight) ) then
       do irec= 1,NBR
capdevil's avatar
   
capdevil committed
336
          write(unitrec,'(a4,1x,f8.3,1x,f8.3)') stations_name(irec) &
capdevil's avatar
capdevil committed
337
338
339
340
               ,90.d0-coord_stations(1,irec),coord_stations(2,irec)
       enddo
    else
       do irec= 1,NBR
capdevil's avatar
   
capdevil committed
341
          write(unitrec,'(a4,1x,f8.3,1x,f8.3,1x,f5.2)') stations_name(irec) &
capdevil's avatar
capdevil committed
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
               ,90.d0-coord_stations(1,irec),coord_stations(2,irec),weight(irec)
       enddo
    endif

!Fermeture du fichier recepteur_out.dat
    close(unitrec)
!cleaning
!-------------------------------------------------------------------------
  end subroutine write_receivers
!-------------------------------------------------------------------------
!-------------------------------------------------
  subroutine get_receivers(rec_coord)
!-------------------------------------------------
    use def_gparam
    implicit none
    real(DP), dimension(:,:), intent(out)  :: rec_coord
!
    rec_coord(:,:)=coord_stations(:,:)
    if (geoc_corr)                             &
         rec_coord(1,:)=rad2deg*(PI/2.d0-atan(GEOC*tan((PI/2.d0-  &
         rec_coord(1,:)/rad2deg))))
!on transforme les degres en radian
    rec_coord(:,:)=rec_coord(:,:)*deg2rad
!------------------------------------------------------------------------
  end subroutine get_receivers
!------------------------------------------------------------------------


!-------------------------------------------------------------------------
  doubleprecision function M0(M)
!calcul la norme de la matrice 3x3 symetrique
!| M(1) M(4) M(5) |
!|      M(2) M(6) |
!|           M(3) |
!-------------------------------------------------------------------------
    implicit none
    doubleprecision, dimension(6) ::M
    M0=sqrt(M(1)**2+M(2)**2+M(3)**2+2.d0*(M(4)**2+M(5)**2+M(6)**2))
!-------------------------------------------------------------------------
  end function M0
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
    subroutine get_source_dat(coord_out,Mtmp_out,name_out)
!-------------------------------------------------------------------------
      implicit none
      doubleprecision, dimension(3,NBE), intent(out)  :: coord_out
      doubleprecision, dimension(6,NBE), intent(out)  :: Mtmp_out
      character(len=8),dimension(NBE), optional, intent(out) :: name_out
      integer :: i,j
!
capdevil's avatar
   
capdevil committed
392
    if (RA<0.d0) STOP 'load_source: RA is not initialized!'
capdevil's avatar
capdevil committed
393
!conversion profondeur (en km)->rayon (en m)
capdevil's avatar
   
capdevil committed
394
395
396
397
398
399
    coord_out(1,:)=RA-coord_sources(1,:)*1000.d0
    coord_out(2:,:)=coord_sources(2:,:)
    Mtmp_out =Mtmp
    if (geoc_corr)                             &
         coord_out(2,:)=rad2deg*(PI/2.d0-atan(GEOC*tan((PI/2.d0-  &
         coord_out(2,:)/rad2deg))))      
capdevil's avatar
capdevil committed
400
!passage en radian
capdevil's avatar
   
capdevil committed
401
402
    coord_out(2:3,:)=coord_out(2:3,:)*deg2rad
    if (present(name_out)) name_out(:)=sources_name(:)
capdevil's avatar
capdevil committed
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
!
!-------------------------------------------------------------------------
    end subroutine get_source_dat
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
    subroutine reset_sources_delay
!-------------------------------------------------------------------------
      use def_gparam
      implicit none
      write(logunit,*)'Warning, forcing source delays to zero !!!!'
      t_delay_sources(:)=0.0_DP
!-------------------------------------------------------------------------
    end subroutine reset_sources_delay
!-------------------------------------------------------------------------

!-------------------------------------------------------------------------
  real(DP) function dist3D(pt1,pt2)
!pseudo distance dans le volume en radian entre pt1 et pt2
!-------------------------------------------------------------------------
    use def_gparam
    implicit none
    real(DP), dimension(3), intent(in) :: pt1,pt2
!
    real(DP) :: a,b,g
!
capdevil's avatar
capdevil committed
428
    if (RA<0.d0) print*,'load_source: RA is not initialized!'
capdevil's avatar
capdevil committed
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
    call euler(pt1(2),pt1(3),pt2(2),pt2(3),a,b,g)
    dist3D=sqrt(b**2+((pt1(1)-pt2(1))/RA)**2)
!-------------------------------------------------------------------------
  end function dist3D
!-------------------------------------------------------------------------

      subroutine euler(t1,p1,t2,p2,a,b,g)
!_____________________________________________________________________
!   
!     calcule les 3 angles d'euler de la rotation definie par:
!     D(a b g)=D(-p1 -t1 0)oD(0 t2 p2)
!     (notation cf Edmonds)
!
!     Entrees:
!        angles p1 p2 dans [-pi,pi] (en fait peu importe mais en
!                                    radians)
!        angles t1 t2 DANS [ 0 ,pi]! (pour etre avec utilise avec 
!        la regle  de sommation d'harmoniques spheriques)
!     Sorties:
!        angles a et g dans [-pi,pi]
!        angle  b      dans [ 0 ,pi]
!_____________________________________________________________________
!
      IMPLICIT NONE
!
      real*8 t1,t2,p1,p2,a,b,g,ca,sa,cg,sg,sb,cb,pi
      real*8 rt1(3,3),rt2(3,3),rp1(3,3),rp2(3,3),prod1(3,3), &
           prod2(3,3),rot(3,3)
      integer i,j,k
!
      pi = 3.141592653589793d0
      do i=1,3
         do j=1,3
            rt1(i,j)  =0.0d0
            rt2(i,j)  =0.0d0
            rp1(i,j)  =0.0d0
            rp2(i,j)  =0.0d0
            prod1(i,j)=0.0d0
            prod2(i,j)=0.0d0
            rot(i,j)  =0.0d0
         enddo
      enddo
!      
      rt1(1,1)=dcos(t1)
      rt1(2,2)=1.0d0
      rt1(3,3)=rt1(1,1)
      rt1(1,3)=dsin(t1)
      rt1(3,1)=-rt1(1,3)
!
      rt2(1,1)=dcos(t2)
      rt2(2,2)=1.0d0
      rt2(3,3)=rt2(1,1)
      rt2(1,3)=-dsin(t2)
      rt2(3,1)=-rt2(1,3)
!
      rp1(1,1)=dcos(p1)
      rp1(2,2)=rp1(1,1)
      rp1(3,3)=1.0d0
      rp1(1,2)=-dsin(p1)
      rp1(2,1)=-rp1(1,2)
!
      rp2(1,1)=dcos(p2)
      rp2(2,2)=rp2(1,1)
      rp2(3,3)=1.0d0
      rp2(1,2)=dsin(p2)
      rp2(2,1)=-rp2(1,2)
!
      do i=1,3
         do j=1,3
            do k=1,3
               prod1(i,j)=prod1(i,j)+rp2(i,k)*rt2(k,j)
            enddo
         enddo
      enddo
      do i=1,3
         do j=1,3
            do k=1,3
               prod2(i,j)=prod2(i,j)+rp1(i,k)*prod1(k,j)
            enddo
         enddo
      enddo
      do i=1,3
         do j=1,3
            do k=1,3
               rot(i,j)=rot(i,j)+rt1(i,k)*prod2(k,j)
            enddo
         enddo
      enddo
!
      cb=rot(3,3)
      if (cb>=1.d0) then
         b=0.d0
      else if (cb<=-1.d0) then
         b=pi
      else
         b=dacos(cb)
      endif
      sb=dsin(b)
      if (abs(sb).le.1.0d-15) then
         a=p2-p1
         g=0.
      else 
         ca=rot(3,1)/sb
         sa=rot(3,2)/sb
         if (abs(ca-1.0d0).lt.1.0d-8) then
            a=0.0d0
         else 
            if (abs(ca+1.0d0).lt.1.0d-8) then
               a=pi
            else
               a=dacos(ca)
            endif
         endif
         if (sa.lt.(0.0d0)) a=-1.0d0*a
         cg=-rot(1,3)/sb
         sg=rot(2,3)/sb
         if (abs(cg-1.0d0).lt.1.0d-8) then
            g=0.0d0
         else
            if (abs(cg+1.0d0).lt.1.0d-8) then
               g=pi
            else
               g=dacos(cg)
            endif
         endif
         if (sg.lt.(0.0d0)) g=-1.0d0*g
      endif
!
      end subroutine euler
!-------------------------------------------------------------------------------
capdevil's avatar
capdevil committed
559
560
561
      subroutine load_station_response()
        implicit none
        character(len=12) :: name,ctmp
capdevil's avatar
minor    
capdevil committed
562
        character(len=8) :: ctmp8
capdevil's avatar
capdevil committed
563
564
        integer ::         NBZmax,NBPmax
        complex*16 :: bid
capdevil's avatar
capdevil committed
565
        real :: a,b
capdevil's avatar
capdevil committed
566
567
568
569
570
571
572
573
        integer :: i,ii,j
        
!looking for the maximum of ZEROz and Poles
        NBZmax=0;NBPmax=0
        allocate(response_constant(NBR),NBZ(NBR),NBP(NBR))
        do i=1,NBR
           write(name,'(a4,"response")') stations_name(i)
           open(11,file=name,status='old')
capdevil's avatar
   
capdevil committed
574
575
!           read(11,'(a5,x,i)') ctmp,ii
           read(11,*) ctmp,ii
capdevil's avatar
capdevil committed
576
577
578
           NBZ(i)=ii
           NBZmax=max(ii,NBZmax)
           do j=1,ii
capdevil's avatar
   
capdevil committed
579
580
581
!              read(11,*) bid
              read(11,*) a,b
              bid=cmplx(a,b)
capdevil's avatar
capdevil committed
582
           enddo
capdevil's avatar
   
capdevil committed
583
584
!           read(11,'(a5,x,i)') ctmp,ii
           read(11,*) ctmp,ii
capdevil's avatar
capdevil committed
585
586
587
588
589
590
591
592
           NBP(i)=ii
           NBPmax=max(ii,NBPmax)
           close(11)
        enddo
        allocate(zeros(NBZmax,NBR),poles(NBPmax,NBR))
        do i=1,NBR
           write(name,'(a4,"response")') stations_name(i)
           open(11,file=name)
capdevil's avatar
   
capdevil committed
593
594
!           read(11,'(a5,x,i)') ctmp,ii
           read(11,*) ctmp,ii
capdevil's avatar
capdevil committed
595
596
597
598
           do j=1,ii
              read(11,*) a,b
              zeros(j,i)=cmplx(a,b)
           enddo
capdevil's avatar
   
capdevil committed
599
600
!           read(11,'(a5,x,i)') ctmp,ii
           read(11,*) ctmp,ii
capdevil's avatar
capdevil committed
601
           do j=1,ii
capdevil's avatar
capdevil committed
602
603
              read(11,*) a,b
              poles(j,i)=cmplx(a,b)
capdevil's avatar
capdevil committed
604
           enddo
capdevil's avatar
   
capdevil committed
605
606
!           read(11,'(a8,x,g)') ctmp,response_constant(i)
           read(11,*) ctmp,response_constant(i)
capdevil's avatar
capdevil committed
607
608
609
           close(11)
        enddo
      end subroutine load_station_response
capdevil's avatar
capdevil committed
610
611
612
!------------------------------------------------------
end module global_main_param
!------------------------------------------------------