All Classes Files Functions Variables Pages
module_solver.f90
Go to the documentation of this file.
1 !=====================================================================================================================================
2 ! EFISPEC3D !
3 ! (Elements FInis SPECtraux 3D) !
4 ! !
5 ! http://efispec.free.fr !
6 ! !
7 ! !
8 ! This file is part of EFISPEC3D !
9 ! Please refer to http://efispec.free.fr if you use it or part of it !
10 ! !
11 ! !
12 !1 ---> French License: CeCILL V2 !
13 ! !
14 ! Copyright BRGM 2009 contributeurs : Florent DE MARTIN !
15 ! David MICHEA !
16 ! Philippe THIERRY !
17 ! !
18 ! Contact: f.demartin at brgm.fr !
19 ! !
20 ! Ce logiciel est un programme informatique servant a resoudre l'equation du !
21 ! mouvement en trois dimensions via une methode des elements finis spectraux. !
22 ! !
23 ! Ce logiciel est regi par la licence CeCILL soumise au droit francais et !
24 ! respectant les principes de diffusion des logiciels libres. Vous pouvez !
25 ! utiliser, modifier et/ou redistribuer ce programme sous les conditions de la !
26 ! licence CeCILL telle que diffusee par le CEA, le CNRS et l'INRIA sur le site !
27 ! "http://www.cecill.info". !
28 ! !
29 ! En contrepartie de l'accessibilite au code source et des droits de copie, de !
30 ! modification et de redistribution accordes par cette licence, il n'est offert !
31 ! aux utilisateurs qu'une garantie limitee. Pour les memes raisons, seule une !
32 ! responsabilite restreinte pese sur l'auteur du programme, le titulaire des !
33 ! droits patrimoniaux et les concedants successifs. !
34 ! !
35 ! A cet egard l'attention de l'utilisateur est attiree sur les risques associes !
36 ! au chargement, a l'utilisation, a la modification et/ou au developpement et a !
37 ! la reproduction du logiciel par l'utilisateur etant donne sa specificite de !
38 ! logiciel libre, qui peut le rendre complexe a manipuler et qui le reserve donc !
39 ! a des developpeurs et des professionnels avertis possedant des connaissances !
40 ! informatiques approfondies. Les utilisateurs sont donc invites a charger et !
41 ! tester l'adequation du logiciel a leurs besoins dans des conditions permettant !
42 ! d'assurer la securite de leurs systemes et ou de leurs donnees et, plus !
43 ! generalement, a l'utiliser et l'exploiter dans les memes conditions de !
44 ! securite. !
45 ! !
46 ! Le fait que vous puissiez acceder a cet en-tete signifie que vous avez pris !
47 ! connaissance de la licence CeCILL et que vous en avez accepte les termes. !
48 ! !
49 ! !
50 !2 ---> International license: GNU GPL V3 !
51 ! !
52 ! EFISPEC3D is a computer program that solves the three-dimensional equations of !
53 ! motion using a finite spectral-element method. !
54 ! !
55 ! Copyright (C) 2009 Florent DE MARTIN !
56 ! !
57 ! Contact: f.demartin at brgm.fr !
58 ! !
59 ! This program is free software: you can redistribute it and/or modify it under !
60 ! the terms of the GNU General Public License as published by the Free Software !
61 ! Foundation, either version 3 of the License, or (at your option) any later !
62 ! version. !
63 ! !
64 ! This program is distributed in the hope that it will be useful, but WITHOUT ANY !
65 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A !
66 ! PARTICULAR PURPOSE. See the GNU General Public License for more details. !
67 ! !
68 ! You should have received a copy of the GNU General Public License along with !
69 ! this program. If not, see http://www.gnu.org/licenses/. !
70 ! !
71 ! !
72 !3 ---> Third party libraries !
73 ! !
74 ! EFISPEC3D uses the following of third party libraries: !
75 ! !
76 ! --> METIS 5.1.0 !
77 ! see http://glaros.dtc.umn.edu/gkhome/metis/metis/overview !
78 ! !
79 ! --> Lib_VTK_IO !
80 ! see S. Zaghi's website: https://github.com/szaghi/Lib_VTK_IO !
81 ! !
82 ! --> INTERP_LINEAR !
83 ! see J. Burkardt website: http://people.sc.fsu.edu/~jburkardt/ !
84 ! !
85 ! --> SAC !
86 ! http://ds.iris.edu/files/sac-manual/ !
87 ! !
88 ! --> EXODUS II !
89 ! http://sourceforge.net/projects/exodusii/ !
90 ! !
91 ! --> NETCDF !
92 ! http://www.unidata.ucar.edu/software/netcdf/ !
93 ! !
94 ! --> HDF5 !
95 ! http://www.hdfgroup.org/HDF5/ !
96 ! !
97 ! Some of these libraries are located in directory lib and pre-compiled !
98 ! with intel compiler for x86_64 platform. !
99 ! !
100 ! !
101 !4 ---> Related Articles !
102 ! !
103 ! De Martin, F., Matsushima, M., Kawase, H. (BSSA, 2013) !
104 ! Impact of geometric effects on near-surface Green's functions !
105 ! doi:10.1785/0120130039 !
106 ! !
107 ! Aochi, H., Ducellier, A., Dupros, F., Delatre, M., Ulrich, T., De Martin, F., !
108 ! Yoshimi, M., (Pure Appl. Geophys. 2013) !
109 ! Finite Difference Simulations of Seismic Wave Propagation for the 2007 Mw 6.6 !
110 ! Niigata-ken Chuetsu-Oki Earthquake: Validity of Models and Reliable !
111 ! Input Ground Motion in the Near-Field !
112 ! doi:10.1007/s00024-011-0429-5 !
113 ! !
114 ! De Martin, F. (BSSA, 2011) !
115 ! Verification of a Spectral-Element Method Code for the Southern California !
116 ! Earthquake Center LOH.3 Viscoelastic Case !
117 ! doi:10.1785/0120100305 !
118 ! !
119 !=====================================================================================================================================
120 
123 
126 
128 
129  use mpi
130 
131  implicit none
132 
133  public :: newmark_ini
134  public :: newmark_end
138  public :: compute_absorption_forces
139  public :: compute_external_force
140 
141  contains
142 
143 !
144 !
149 !***********************************************************************************************************************************************************************************
150  subroutine newmark_ini()
151 !***********************************************************************************************************************************************************************************
152 
153  use mpi
154 
155  use mod_global_variables, only :&
156  ig_ngll_total&
157  ,rg_gll_displacement&
158  ,rg_gll_velocity&
159  ,rg_gll_acceleration&
160  ,rg_dt&
161  ,rg_dt2&
162  ,rg_newmark_gamma
163 
164  implicit none
165 
166  integer :: igll
167 
168 !
169 !------->compute displacement at step n+1
170  do igll = 1,ig_ngll_total
171  rg_gll_displacement(1,igll) = rg_gll_displacement(1,igll) + rg_dt*rg_gll_velocity(1,igll) + rg_dt2*0.5*rg_gll_acceleration(1,igll) !displacement x
172  rg_gll_displacement(2,igll) = rg_gll_displacement(2,igll) + rg_dt*rg_gll_velocity(2,igll) + rg_dt2*0.5*rg_gll_acceleration(2,igll) !dispalcement y
173  rg_gll_displacement(3,igll) = rg_gll_displacement(3,igll) + rg_dt*rg_gll_velocity(3,igll) + rg_dt2*0.5*rg_gll_acceleration(3,igll) !displacement z
174  enddo
175 
176 !
177 !------->compute velocity at step n+1/2
178  do igll = 1,ig_ngll_total
179  rg_gll_velocity(1,igll) = rg_gll_velocity(1,igll) + rg_dt*(1.0-rg_newmark_gamma)*rg_gll_acceleration(1,igll) !velocity x
180  rg_gll_velocity(2,igll) = rg_gll_velocity(2,igll) + rg_dt*(1.0-rg_newmark_gamma)*rg_gll_acceleration(2,igll) !velocity y
181  rg_gll_velocity(3,igll) = rg_gll_velocity(3,igll) + rg_dt*(1.0-rg_newmark_gamma)*rg_gll_acceleration(3,igll) !velocity z
182  enddo
183 
184 !
185 !------->flush array rg_gll_acceleration to zero. This array is used as external/internal forces array in subroutines compute_external_force/compute_internal_forces, respectively
186  do igll = 1,ig_ngll_total
187  rg_gll_acceleration(1,igll) = 0.0
188  rg_gll_acceleration(2,igll) = 0.0
189  rg_gll_acceleration(3,igll) = 0.0
190  enddo
191 
192  return
193 !***********************************************************************************************************************************************************************************
194  end subroutine newmark_ini
195 !***********************************************************************************************************************************************************************************
196 
197 !
198 !
203 !***********************************************************************************************************************************************************************************
204  subroutine newmark_end()
205 !***********************************************************************************************************************************************************************************
206 
207  use mpi
208 
209  use mod_global_variables, only :&
210  ig_ngll_total&
211  ,rg_gll_velocity&
212  ,rg_gll_acceleration&
213  ,rg_dt&
214  ,rg_newmark_gamma&
215  ,rg_gll_mass_matrix
216  implicit none
217 
218  integer :: igll
219 
220  do igll = 1,ig_ngll_total
221  rg_gll_acceleration(1,igll) = rg_gll_acceleration(1,igll)*rg_gll_mass_matrix(igll) !acceleration x step n+1
222  rg_gll_acceleration(2,igll) = rg_gll_acceleration(2,igll)*rg_gll_mass_matrix(igll) !acceleration y step n+1
223  rg_gll_acceleration(3,igll) = rg_gll_acceleration(3,igll)*rg_gll_mass_matrix(igll) !acceleration z step n+1
224  enddo
225 
226  do igll = 1,ig_ngll_total
227  rg_gll_velocity(1,igll) = rg_gll_velocity(1,igll) + rg_dt*rg_newmark_gamma*rg_gll_acceleration(1,igll) !velocity x step n+1
228  rg_gll_velocity(2,igll) = rg_gll_velocity(2,igll) + rg_dt*rg_newmark_gamma*rg_gll_acceleration(2,igll) !velocity y step n+1
229  rg_gll_velocity(3,igll) = rg_gll_velocity(3,igll) + rg_dt*rg_newmark_gamma*rg_gll_acceleration(3,igll) !velocity z step n+1
230  enddo
231 
232  return
233 !***********************************************************************************************************************************************************************************
234  end subroutine newmark_end
235 !***********************************************************************************************************************************************************************************
236 
237 !
238 !
244 !***********************************************************************************************************************************************************************************
245  subroutine compute_internal_forces_order4(elt_start,elt_end)
246 !***********************************************************************************************************************************************************************************
247 
248  use mpi
249 
250  use mod_global_variables, only :&
251  ig_ngll&
252  ,ig_ndof&
253  ,rg_gll_lagrange_deriv&
254  ,rg_gll_displacement&
255  ,rg_gll_acceleration&
256  ,rg_gll_weight&
257  ,rg_gll_abscissa&
258  ,ig_nreceiver_hexa&
259  ,ig_idt&
260  ,ig_receiver_saving_incr&
261  ,ig_nhexa&
262  ,ig_hexa_gll_glonum&
263  ,rg_hexa_gll_dxidx&
264  ,rg_hexa_gll_dxidy&
265  ,rg_hexa_gll_dxidz&
266  ,rg_hexa_gll_detdx&
267  ,rg_hexa_gll_detdy&
268  ,rg_hexa_gll_detdz&
269  ,rg_hexa_gll_dzedx&
270  ,rg_hexa_gll_dzedy&
271  ,rg_hexa_gll_dzedz&
272  ,rg_hexa_gll_jacobian_det&
273  ,ig_myrank&
274  ,ig_nhexa_outer&
275  ,rg_hexa_gll_rho&
276  ,rg_hexa_gll_rhovs2&
277  ,rg_hexa_gll_rhovp2&
278  ,rg_hexa_gll_rhovs2&
279  ,rg_hexa_gll_rhovp2&
280  ,rg_hexa_gll_wkqs&
281  ,rg_hexa_gll_wkqp&
282  ,rg_hexa_gll_ksixx&
283  ,rg_hexa_gll_ksiyy&
284  ,rg_hexa_gll_ksizz&
285  ,rg_hexa_gll_ksixy&
286  ,rg_hexa_gll_ksixz&
287  ,rg_hexa_gll_ksiyz&
288  ,rg_relax_coeff&
289  ,rg_mem_var_exp&
290  ,ig_nrelax&
291  ,rg_dt&
292  ,lg_visco&
293  ,ig_nhexa
294 
295  implicit none
296 
297  integer, intent(in) :: elt_start
298  integer, intent(in) :: elt_end
299 
300  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx1
301  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx2
302  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx3
303  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy1
304  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy2
305  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy3
306  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz1
307  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz2
308  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz3
309  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: rl_displacement_gll
310  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: rl_acceleration_gll
311 
312  real :: duxdxi
313  real :: duxdet
314  real :: duxdze
315  real :: duydxi
316  real :: duydet
317  real :: duydze
318  real :: duzdxi
319  real :: duzdet
320  real :: duzdze
321  real :: duxdx
322  real :: duydy
323  real :: duzdz
324  real :: duxdy
325  real :: duxdz
326  real :: duydx
327  real :: duydz
328  real :: duzdx
329  real :: duzdy
330  real :: dxidx
331  real :: dxidy
332  real :: dxidz
333  real :: detdx
334  real :: detdy
335  real :: detdz
336  real :: dzedx
337  real :: dzedy
338  real :: dzedz
339  real :: tauxx
340  real :: tauyy
341  real :: tauzz
342  real :: tauxy
343  real :: tauxz
344  real :: tauyz
345  real :: tauxx_n12
346  real :: tauyy_n12
347  real :: tauzz_n12
348  real :: tauxy_n12
349  real :: tauxz_n12
350  real :: tauyz_n12
351  real :: trace_tau
352  real :: tmpx1
353  real :: tmpx2
354  real :: tmpx3
355  real :: tmpx4
356  real :: tmpy1
357  real :: tmpy2
358  real :: tmpy3
359  real :: tmpz1
360  real :: tmpz2
361  real :: tmpz3
362  real :: fac1
363  real :: fac2
364  real :: fac3
365 
366  integer :: iel
367  integer :: k
368  integer :: l
369  integer :: m
370  integer :: igll
371  integer :: imem_var
372 
373 
374  do iel = elt_start,elt_end
375 
376  !
377  !------->fill local displacement
378  do k = 1,ig_ngll !zeta
379  do l = 1,ig_ngll !eta
380  do m = 1,ig_ngll !xi
381 
382  igll = ig_hexa_gll_glonum(m,l,k,iel)
383 
384  rl_displacement_gll(1,m,l,k) = rg_gll_displacement(1,igll)
385  rl_displacement_gll(2,m,l,k) = rg_gll_displacement(2,igll)
386  rl_displacement_gll(3,m,l,k) = rg_gll_displacement(3,igll)
387 
388  enddo
389  enddo
390  enddo
391  !
392  !
393  !******************************************************************************
394  !->compute integrale at GLL nodes + assemble forces in global gll grid for hexa
395  !******************************************************************************
396  do k = 1,ig_ngll !zeta
397  do l = 1,ig_ngll !eta
398  do m = 1,ig_ngll !xi
399 
400  !
401  !---------->derivative of displacement with respect to local coordinate xi, eta and zeta at the gll node klm
402 
403  duxdxi = rl_displacement_gll(1,1,l,k)*rg_gll_lagrange_deriv(1,m) &
404  + rl_displacement_gll(1,2,l,k)*rg_gll_lagrange_deriv(2,m) &
405  + rl_displacement_gll(1,3,l,k)*rg_gll_lagrange_deriv(3,m) &
406  + rl_displacement_gll(1,4,l,k)*rg_gll_lagrange_deriv(4,m) &
407  + rl_displacement_gll(1,5,l,k)*rg_gll_lagrange_deriv(5,m)
408 
409  duydxi = rl_displacement_gll(2,1,l,k)*rg_gll_lagrange_deriv(1,m) &
410  + rl_displacement_gll(2,2,l,k)*rg_gll_lagrange_deriv(2,m) &
411  + rl_displacement_gll(2,3,l,k)*rg_gll_lagrange_deriv(3,m) &
412  + rl_displacement_gll(2,4,l,k)*rg_gll_lagrange_deriv(4,m) &
413  + rl_displacement_gll(2,5,l,k)*rg_gll_lagrange_deriv(5,m)
414 
415  duzdxi = rl_displacement_gll(3,1,l,k)*rg_gll_lagrange_deriv(1,m) &
416  + rl_displacement_gll(3,2,l,k)*rg_gll_lagrange_deriv(2,m) &
417  + rl_displacement_gll(3,3,l,k)*rg_gll_lagrange_deriv(3,m) &
418  + rl_displacement_gll(3,4,l,k)*rg_gll_lagrange_deriv(4,m) &
419  + rl_displacement_gll(3,5,l,k)*rg_gll_lagrange_deriv(5,m)
420 
421  duxdet = rl_displacement_gll(1,m,1,k)*rg_gll_lagrange_deriv(1,l) &
422  + rl_displacement_gll(1,m,2,k)*rg_gll_lagrange_deriv(2,l) &
423  + rl_displacement_gll(1,m,3,k)*rg_gll_lagrange_deriv(3,l) &
424  + rl_displacement_gll(1,m,4,k)*rg_gll_lagrange_deriv(4,l) &
425  + rl_displacement_gll(1,m,5,k)*rg_gll_lagrange_deriv(5,l)
426 
427  duydet = rl_displacement_gll(2,m,1,k)*rg_gll_lagrange_deriv(1,l) &
428  + rl_displacement_gll(2,m,2,k)*rg_gll_lagrange_deriv(2,l) &
429  + rl_displacement_gll(2,m,3,k)*rg_gll_lagrange_deriv(3,l) &
430  + rl_displacement_gll(2,m,4,k)*rg_gll_lagrange_deriv(4,l) &
431  + rl_displacement_gll(2,m,5,k)*rg_gll_lagrange_deriv(5,l)
432 
433  duzdet = rl_displacement_gll(3,m,1,k)*rg_gll_lagrange_deriv(1,l) &
434  + rl_displacement_gll(3,m,2,k)*rg_gll_lagrange_deriv(2,l) &
435  + rl_displacement_gll(3,m,3,k)*rg_gll_lagrange_deriv(3,l) &
436  + rl_displacement_gll(3,m,4,k)*rg_gll_lagrange_deriv(4,l) &
437  + rl_displacement_gll(3,m,5,k)*rg_gll_lagrange_deriv(5,l)
438 
439  duxdze = rl_displacement_gll(1,m,l,1)*rg_gll_lagrange_deriv(1,k) &
440  + rl_displacement_gll(1,m,l,2)*rg_gll_lagrange_deriv(2,k) &
441  + rl_displacement_gll(1,m,l,3)*rg_gll_lagrange_deriv(3,k) &
442  + rl_displacement_gll(1,m,l,4)*rg_gll_lagrange_deriv(4,k) &
443  + rl_displacement_gll(1,m,l,5)*rg_gll_lagrange_deriv(5,k)
444 
445  duydze = rl_displacement_gll(2,m,l,1)*rg_gll_lagrange_deriv(1,k) &
446  + rl_displacement_gll(2,m,l,2)*rg_gll_lagrange_deriv(2,k) &
447  + rl_displacement_gll(2,m,l,3)*rg_gll_lagrange_deriv(3,k) &
448  + rl_displacement_gll(2,m,l,4)*rg_gll_lagrange_deriv(4,k) &
449  + rl_displacement_gll(2,m,l,5)*rg_gll_lagrange_deriv(5,k)
450 
451  duzdze = rl_displacement_gll(3,m,l,1)*rg_gll_lagrange_deriv(1,k) &
452  + rl_displacement_gll(3,m,l,2)*rg_gll_lagrange_deriv(2,k) &
453  + rl_displacement_gll(3,m,l,3)*rg_gll_lagrange_deriv(3,k) &
454  + rl_displacement_gll(3,m,l,4)*rg_gll_lagrange_deriv(4,k) &
455  + rl_displacement_gll(3,m,l,5)*rg_gll_lagrange_deriv(5,k)
456 
457  !
458  !---------->derivative of displacement at step n+1 with respect to global coordinate x, y and z at the gll node klm
459  dxidx = rg_hexa_gll_dxidx(m,l,k,iel)
460  dxidy = rg_hexa_gll_dxidy(m,l,k,iel)
461  dxidz = rg_hexa_gll_dxidz(m,l,k,iel)
462  detdx = rg_hexa_gll_detdx(m,l,k,iel)
463  detdy = rg_hexa_gll_detdy(m,l,k,iel)
464  detdz = rg_hexa_gll_detdz(m,l,k,iel)
465  dzedx = rg_hexa_gll_dzedx(m,l,k,iel)
466  dzedy = rg_hexa_gll_dzedy(m,l,k,iel)
467  dzedz = rg_hexa_gll_dzedz(m,l,k,iel)
468 
469  duxdx = duxdxi*dxidx + duxdet*detdx + duxdze*dzedx
470  duxdy = duxdxi*dxidy + duxdet*detdy + duxdze*dzedy
471  duxdz = duxdxi*dxidz + duxdet*detdz + duxdze*dzedz
472  duydx = duydxi*dxidx + duydet*detdx + duydze*dzedx
473  duydy = duydxi*dxidy + duydet*detdy + duydze*dzedy
474  duydz = duydxi*dxidz + duydet*detdz + duydze*dzedz
475  duzdx = duzdxi*dxidx + duzdet*detdx + duzdze*dzedx
476  duzdy = duzdxi*dxidy + duzdet*detdy + duzdze*dzedy
477  duzdz = duzdxi*dxidz + duzdet*detdz + duzdze*dzedz
478  !
479  !---------->compute elastic stress (elastic simulation) or unrelaxed elastic stress (viscoelastic simulation)
480  trace_tau = (rg_hexa_gll_rhovp2(m,l,k,iel) - 2.0*rg_hexa_gll_rhovs2(m,l,k,iel))*(duxdx+duydy+duzdz)
481  tauxx = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duxdx
482  tauyy = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duydy
483  tauzz = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duzdz
484  tauxy = rg_hexa_gll_rhovs2(m,l,k,iel)*(duxdy+duydx)
485  tauxz = rg_hexa_gll_rhovs2(m,l,k,iel)*(duxdz+duzdx)
486  tauyz = rg_hexa_gll_rhovs2(m,l,k,iel)*(duydz+duzdy)
487  !
488  !----------->compute viscoelastic stress
489 
490  if (lg_visco) then
491 
492  do imem_var = 1,ig_nrelax
493 
494  tmpx1 = rg_mem_var_exp(imem_var)
495 
496  tmpx2 = rg_hexa_gll_rhovp2(m,l,k,iel)*rg_hexa_gll_wkqp(imem_var,m,l,k,iel)
497  tmpx3 = 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*rg_hexa_gll_wkqs(imem_var,m,l,k,iel)
498 
499  tmpx4 = (duxdx+duydy+duzdz)*(tmpx2 - tmpx3)
500 
501  !
502  !----------------->anelastic stress at step n+1/2 following s. ma and p. liu (2006) using epsnp1 and unrelaxed material modulus
503  tauxx_n12 = tmpx1*rg_hexa_gll_ksixx(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duxdx + tmpx4)
504  tauyy_n12 = tmpx1*rg_hexa_gll_ksiyy(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duydy + tmpx4)
505  tauzz_n12 = tmpx1*rg_hexa_gll_ksizz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duzdz + tmpx4)
506  tauxy_n12 = tmpx1*rg_hexa_gll_ksixy(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duxdy+duydx))
507  tauxz_n12 = tmpx1*rg_hexa_gll_ksixz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duxdz+duzdx))
508  tauyz_n12 = tmpx1*rg_hexa_gll_ksiyz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duydz+duzdy))
509  !
510  !----------------->compute final stress at step n+1 according to day and minster (1984) + ma and liu (2006) : tauxx - SUM anelastic stress at step n+1
511  tauxx = tauxx - 0.5*(tauxx_n12 + rg_hexa_gll_ksixx(imem_var,m,l,k,iel))
512  tauyy = tauyy - 0.5*(tauyy_n12 + rg_hexa_gll_ksiyy(imem_var,m,l,k,iel))
513  tauzz = tauzz - 0.5*(tauzz_n12 + rg_hexa_gll_ksizz(imem_var,m,l,k,iel))
514  tauxy = tauxy - 0.5*(tauxy_n12 + rg_hexa_gll_ksixy(imem_var,m,l,k,iel))
515  tauxz = tauxz - 0.5*(tauxz_n12 + rg_hexa_gll_ksixz(imem_var,m,l,k,iel))
516  tauyz = tauyz - 0.5*(tauyz_n12 + rg_hexa_gll_ksiyz(imem_var,m,l,k,iel))
517  !
518  !----------------->update memory for stress (step n+1/2)
519  rg_hexa_gll_ksixx(imem_var,m,l,k,iel) = tauxx_n12
520  rg_hexa_gll_ksiyy(imem_var,m,l,k,iel) = tauyy_n12
521  rg_hexa_gll_ksizz(imem_var,m,l,k,iel) = tauzz_n12
522  rg_hexa_gll_ksixy(imem_var,m,l,k,iel) = tauxy_n12
523  rg_hexa_gll_ksixz(imem_var,m,l,k,iel) = tauxz_n12
524  rg_hexa_gll_ksiyz(imem_var,m,l,k,iel) = tauyz_n12
525 
526  enddo
527 
528  endif
529  !
530  !---------->store members of integration of the gll node klm
531  intpx1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*dxidx+tauxy*dxidy+tauxz*dxidz)
532  intpx2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*detdx+tauxy*detdy+tauxz*detdz)
533  intpx3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*dzedx+tauxy*dzedy+tauxz*dzedz)
534 
535  intpy1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*dxidx+tauyy*dxidy+tauyz*dxidz)
536  intpy2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*detdx+tauyy*detdy+tauyz*detdz)
537  intpy3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*dzedx+tauyy*dzedy+tauyz*dzedz)
538 
539  intpz1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*dxidx+tauyz*dxidy+tauzz*dxidz)
540  intpz2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*detdx+tauyz*detdy+tauzz*detdz)
541  intpz3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*dzedx+tauyz*dzedy+tauzz*dzedz)
542 
543  enddo !xi
544  enddo !eta
545  enddo !zeta
546 
547  !
548  !->finish integration for hexa (internal forces at step n+1)
549  do k = 1,ig_ngll
550  do l = 1,ig_ngll
551  do m = 1,ig_ngll
552 
553  tmpx1 = intpx1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
554  + intpx1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
555  + intpx1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
556  + intpx1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
557  + intpx1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5)
558 
559  tmpy1 = intpy1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
560  + intpy1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
561  + intpy1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
562  + intpy1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
563  + intpy1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5)
564 
565  tmpz1 = intpz1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
566  + intpz1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
567  + intpz1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
568  + intpz1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
569  + intpz1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5)
570 
571  tmpx2 = intpx2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
572  + intpx2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
573  + intpx2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
574  + intpx2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
575  + intpx2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5)
576 
577  tmpy2 = intpy2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
578  + intpy2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
579  + intpy2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
580  + intpy2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
581  + intpy2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5)
582 
583  tmpz2 = intpz2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
584  + intpz2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
585  + intpz2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
586  + intpz2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
587  + intpz2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5)
588 
589  tmpx3 = intpx3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
590  + intpx3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
591  + intpx3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
592  + intpx3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
593  + intpx3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5)
594 
595  tmpy3 = intpy3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
596  + intpy3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
597  + intpy3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
598  + intpy3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
599  + intpy3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5)
600 
601  tmpz3 = intpz3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
602  + intpz3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
603  + intpz3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
604  + intpz3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
605  + intpz3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5)
606 
607  fac1 = rg_gll_weight(l)*rg_gll_weight(k)
608  fac2 = rg_gll_weight(m)*rg_gll_weight(k)
609  fac3 = rg_gll_weight(m)*rg_gll_weight(l)
610 
611  rl_acceleration_gll(1,m,l,k) = (fac1*tmpx1 + fac2*tmpx2 + fac3*tmpx3)
612  rl_acceleration_gll(2,m,l,k) = (fac1*tmpy1 + fac2*tmpy2 + fac3*tmpy3)
613  rl_acceleration_gll(3,m,l,k) = (fac1*tmpz1 + fac2*tmpz2 + fac3*tmpz3)
614 
615  enddo
616  enddo
617  enddo
618 
619  do k = 1,ig_ngll !zeta
620  do l = 1,ig_ngll !eta
621  do m = 1,ig_ngll !xi
622 
623  igll = ig_hexa_gll_glonum(m,l,k,iel)
624 
625  rg_gll_acceleration(1,igll) = rg_gll_acceleration(1,igll) - rl_acceleration_gll(1,m,l,k)
626  rg_gll_acceleration(2,igll) = rg_gll_acceleration(2,igll) - rl_acceleration_gll(2,m,l,k)
627  rg_gll_acceleration(3,igll) = rg_gll_acceleration(3,igll) - rl_acceleration_gll(3,m,l,k)
628 
629  enddo
630  enddo
631  enddo
632 
633  enddo !loop on hexahedron elements
634 
635  return
636 !***********************************************************************************************************************************************************************************
637  end subroutine compute_internal_forces_order4
638 !***********************************************************************************************************************************************************************************
639 
640 !
641 !
647 !***********************************************************************************************************************************************************************************
648  subroutine compute_internal_forces_order5(elt_start,elt_end)
649 !***********************************************************************************************************************************************************************************
650 
651  use mpi
652 
653  use mod_global_variables, only :&
654  ig_ngll&
655  ,ig_ndof&
656  ,rg_gll_lagrange_deriv&
657  ,rg_gll_displacement&
658  ,rg_gll_acceleration&
659  ,rg_gll_weight&
660  ,rg_gll_abscissa&
661  ,ig_nreceiver_hexa&
662  ,ig_idt&
663  ,ig_receiver_saving_incr&
664  ,ig_nhexa&
665  ,ig_hexa_gll_glonum&
666  ,rg_hexa_gll_dxidx&
667  ,rg_hexa_gll_dxidy&
668  ,rg_hexa_gll_dxidz&
669  ,rg_hexa_gll_detdx&
670  ,rg_hexa_gll_detdy&
671  ,rg_hexa_gll_detdz&
672  ,rg_hexa_gll_dzedx&
673  ,rg_hexa_gll_dzedy&
674  ,rg_hexa_gll_dzedz&
675  ,rg_hexa_gll_jacobian_det&
676  ,ig_myrank&
677  ,ig_nhexa_outer&
678  ,rg_hexa_gll_rho&
679  ,rg_hexa_gll_rhovs2&
680  ,rg_hexa_gll_rhovp2&
681  ,rg_hexa_gll_rhovs2&
682  ,rg_hexa_gll_rhovp2&
683  ,rg_hexa_gll_wkqs&
684  ,rg_hexa_gll_wkqp&
685  ,rg_hexa_gll_ksixx&
686  ,rg_hexa_gll_ksiyy&
687  ,rg_hexa_gll_ksizz&
688  ,rg_hexa_gll_ksixy&
689  ,rg_hexa_gll_ksixz&
690  ,rg_hexa_gll_ksiyz&
691  ,rg_relax_coeff&
692  ,rg_mem_var_exp&
693  ,ig_nrelax&
694  ,rg_dt&
695  ,lg_visco&
696  ,ig_nhexa
697 
698  implicit none
699 
700  integer, intent(in) :: elt_start
701  integer, intent(in) :: elt_end
702 
703  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx1
704  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx2
705  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx3
706  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy1
707  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy2
708  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy3
709  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz1
710  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz2
711  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz3
712  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: rl_displacement_gll
713  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: rl_acceleration_gll
714 
715  real :: duxdxi
716  real :: duxdet
717  real :: duxdze
718  real :: duydxi
719  real :: duydet
720  real :: duydze
721  real :: duzdxi
722  real :: duzdet
723  real :: duzdze
724  real :: duxdx
725  real :: duydy
726  real :: duzdz
727  real :: duxdy
728  real :: duxdz
729  real :: duydx
730  real :: duydz
731  real :: duzdx
732  real :: duzdy
733  real :: dxidx
734  real :: dxidy
735  real :: dxidz
736  real :: detdx
737  real :: detdy
738  real :: detdz
739  real :: dzedx
740  real :: dzedy
741  real :: dzedz
742  real :: tauxx
743  real :: tauyy
744  real :: tauzz
745  real :: tauxy
746  real :: tauxz
747  real :: tauyz
748  real :: tauxx_n12
749  real :: tauyy_n12
750  real :: tauzz_n12
751  real :: tauxy_n12
752  real :: tauxz_n12
753  real :: tauyz_n12
754  real :: trace_tau
755  real :: tmpx1
756  real :: tmpx2
757  real :: tmpx3
758  real :: tmpx4
759  real :: tmpy1
760  real :: tmpy2
761  real :: tmpy3
762  real :: tmpz1
763  real :: tmpz2
764  real :: tmpz3
765  real :: fac1
766  real :: fac2
767  real :: fac3
768 
769  integer :: iel
770  integer :: k
771  integer :: l
772  integer :: m
773  integer :: igll
774  integer :: imem_var
775 
776 
777  do iel = elt_start,elt_end
778 
779  !
780  !------->flush to zero local acceleration
781  do k = 1,ig_ngll
782  do l = 1,ig_ngll
783  do m = 1,ig_ngll
784 
785  rl_acceleration_gll(1,m,l,k) = 0.0
786  rl_acceleration_gll(2,m,l,k) = 0.0
787  rl_acceleration_gll(3,m,l,k) = 0.0
788 
789  enddo
790  enddo
791  enddo
792 
793  !
794  !------->fill local displacement
795  do k = 1,ig_ngll !zeta
796  do l = 1,ig_ngll !eta
797  do m = 1,ig_ngll !xi
798 
799  igll = ig_hexa_gll_glonum(m,l,k,iel)
800 
801  rl_displacement_gll(1,m,l,k) = rg_gll_displacement(1,igll)
802  rl_displacement_gll(2,m,l,k) = rg_gll_displacement(2,igll)
803  rl_displacement_gll(3,m,l,k) = rg_gll_displacement(3,igll)
804 
805  enddo
806  enddo
807  enddo
808  !
809  !
810  !******************************************************************************
811  !->compute integrale at gll nodes + assemble force in global gll grid for hexa
812  !******************************************************************************
813  do k = 1,ig_ngll !zeta
814  do l = 1,ig_ngll !eta
815  do m = 1,ig_ngll !xi
816  !
817  !---------->derivative of displacement with respect to local coordinate xi, eta and zeta at the gll node klm
818 
819 
820  duxdxi = rl_displacement_gll(1,1,l,k)*rg_gll_lagrange_deriv(1,m) &
821  + rl_displacement_gll(1,2,l,k)*rg_gll_lagrange_deriv(2,m) &
822  + rl_displacement_gll(1,3,l,k)*rg_gll_lagrange_deriv(3,m) &
823  + rl_displacement_gll(1,4,l,k)*rg_gll_lagrange_deriv(4,m) &
824  + rl_displacement_gll(1,5,l,k)*rg_gll_lagrange_deriv(5,m) &
825  + rl_displacement_gll(1,6,l,k)*rg_gll_lagrange_deriv(6,m)
826 
827  duydxi = rl_displacement_gll(2,1,l,k)*rg_gll_lagrange_deriv(1,m) &
828  + rl_displacement_gll(2,2,l,k)*rg_gll_lagrange_deriv(2,m) &
829  + rl_displacement_gll(2,3,l,k)*rg_gll_lagrange_deriv(3,m) &
830  + rl_displacement_gll(2,4,l,k)*rg_gll_lagrange_deriv(4,m) &
831  + rl_displacement_gll(2,5,l,k)*rg_gll_lagrange_deriv(5,m) &
832  + rl_displacement_gll(2,6,l,k)*rg_gll_lagrange_deriv(6,m)
833 
834  duzdxi = rl_displacement_gll(3,1,l,k)*rg_gll_lagrange_deriv(1,m) &
835  + rl_displacement_gll(3,2,l,k)*rg_gll_lagrange_deriv(2,m) &
836  + rl_displacement_gll(3,3,l,k)*rg_gll_lagrange_deriv(3,m) &
837  + rl_displacement_gll(3,4,l,k)*rg_gll_lagrange_deriv(4,m) &
838  + rl_displacement_gll(3,5,l,k)*rg_gll_lagrange_deriv(5,m) &
839  + rl_displacement_gll(3,6,l,k)*rg_gll_lagrange_deriv(6,m)
840 
841  duxdet = rl_displacement_gll(1,m,1,k)*rg_gll_lagrange_deriv(1,l) &
842  + rl_displacement_gll(1,m,2,k)*rg_gll_lagrange_deriv(2,l) &
843  + rl_displacement_gll(1,m,3,k)*rg_gll_lagrange_deriv(3,l) &
844  + rl_displacement_gll(1,m,4,k)*rg_gll_lagrange_deriv(4,l) &
845  + rl_displacement_gll(1,m,5,k)*rg_gll_lagrange_deriv(5,l) &
846  + rl_displacement_gll(1,m,6,k)*rg_gll_lagrange_deriv(6,l)
847 
848  duydet = rl_displacement_gll(2,m,1,k)*rg_gll_lagrange_deriv(1,l) &
849  + rl_displacement_gll(2,m,2,k)*rg_gll_lagrange_deriv(2,l) &
850  + rl_displacement_gll(2,m,3,k)*rg_gll_lagrange_deriv(3,l) &
851  + rl_displacement_gll(2,m,4,k)*rg_gll_lagrange_deriv(4,l) &
852  + rl_displacement_gll(2,m,5,k)*rg_gll_lagrange_deriv(5,l) &
853  + rl_displacement_gll(2,m,6,k)*rg_gll_lagrange_deriv(6,l)
854 
855  duzdet = rl_displacement_gll(3,m,1,k)*rg_gll_lagrange_deriv(1,l) &
856  + rl_displacement_gll(3,m,2,k)*rg_gll_lagrange_deriv(2,l) &
857  + rl_displacement_gll(3,m,3,k)*rg_gll_lagrange_deriv(3,l) &
858  + rl_displacement_gll(3,m,4,k)*rg_gll_lagrange_deriv(4,l) &
859  + rl_displacement_gll(3,m,5,k)*rg_gll_lagrange_deriv(5,l) &
860  + rl_displacement_gll(3,m,6,k)*rg_gll_lagrange_deriv(6,l)
861 
862  duxdze = rl_displacement_gll(1,m,l,1)*rg_gll_lagrange_deriv(1,k) &
863  + rl_displacement_gll(1,m,l,2)*rg_gll_lagrange_deriv(2,k) &
864  + rl_displacement_gll(1,m,l,3)*rg_gll_lagrange_deriv(3,k) &
865  + rl_displacement_gll(1,m,l,4)*rg_gll_lagrange_deriv(4,k) &
866  + rl_displacement_gll(1,m,l,5)*rg_gll_lagrange_deriv(5,k) &
867  + rl_displacement_gll(1,m,l,6)*rg_gll_lagrange_deriv(6,k)
868 
869  duydze = rl_displacement_gll(2,m,l,1)*rg_gll_lagrange_deriv(1,k) &
870  + rl_displacement_gll(2,m,l,2)*rg_gll_lagrange_deriv(2,k) &
871  + rl_displacement_gll(2,m,l,3)*rg_gll_lagrange_deriv(3,k) &
872  + rl_displacement_gll(2,m,l,4)*rg_gll_lagrange_deriv(4,k) &
873  + rl_displacement_gll(2,m,l,5)*rg_gll_lagrange_deriv(5,k) &
874  + rl_displacement_gll(2,m,l,6)*rg_gll_lagrange_deriv(6,k)
875 
876  duzdze = rl_displacement_gll(3,m,l,1)*rg_gll_lagrange_deriv(1,k) &
877  + rl_displacement_gll(3,m,l,2)*rg_gll_lagrange_deriv(2,k) &
878  + rl_displacement_gll(3,m,l,3)*rg_gll_lagrange_deriv(3,k) &
879  + rl_displacement_gll(3,m,l,4)*rg_gll_lagrange_deriv(4,k) &
880  + rl_displacement_gll(3,m,l,5)*rg_gll_lagrange_deriv(5,k) &
881  + rl_displacement_gll(3,m,l,6)*rg_gll_lagrange_deriv(6,k)
882 
883  !
884  !---------->derivative of displacement at step n+1 with respect to global coordinate x, y and z at the gll node klm
885  dxidx = rg_hexa_gll_dxidx(m,l,k,iel)
886  dxidy = rg_hexa_gll_dxidy(m,l,k,iel)
887  dxidz = rg_hexa_gll_dxidz(m,l,k,iel)
888  detdx = rg_hexa_gll_detdx(m,l,k,iel)
889  detdy = rg_hexa_gll_detdy(m,l,k,iel)
890  detdz = rg_hexa_gll_detdz(m,l,k,iel)
891  dzedx = rg_hexa_gll_dzedx(m,l,k,iel)
892  dzedy = rg_hexa_gll_dzedy(m,l,k,iel)
893  dzedz = rg_hexa_gll_dzedz(m,l,k,iel)
894 
895  duxdx = duxdxi*dxidx + duxdet*detdx + duxdze*dzedx
896  duxdy = duxdxi*dxidy + duxdet*detdy + duxdze*dzedy
897  duxdz = duxdxi*dxidz + duxdet*detdz + duxdze*dzedz
898  duydx = duydxi*dxidx + duydet*detdx + duydze*dzedx
899  duydy = duydxi*dxidy + duydet*detdy + duydze*dzedy
900  duydz = duydxi*dxidz + duydet*detdz + duydze*dzedz
901  duzdx = duzdxi*dxidx + duzdet*detdx + duzdze*dzedx
902  duzdy = duzdxi*dxidy + duzdet*detdy + duzdze*dzedy
903  duzdz = duzdxi*dxidz + duzdet*detdz + duzdze*dzedz
904  !
905  !---------->compute elastic stress (elastic simulation) or unrelaxed elastic stress (viscoelastic simulation)
906  trace_tau = (rg_hexa_gll_rhovp2(m,l,k,iel) - 2.0*rg_hexa_gll_rhovs2(m,l,k,iel))*(duxdx+duydy+duzdz)
907  tauxx = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duxdx
908  tauyy = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duydy
909  tauzz = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duzdz
910  tauxy = rg_hexa_gll_rhovs2(m,l,k,iel)*(duxdy+duydx)
911  tauxz = rg_hexa_gll_rhovs2(m,l,k,iel)*(duxdz+duzdx)
912  tauyz = rg_hexa_gll_rhovs2(m,l,k,iel)*(duydz+duzdy)
913  !
914  !---------->compute viscoelastic stress
915 
916  if(lg_visco) then
917 
918  do imem_var = 1,ig_nrelax
919 
920  tmpx1 = rg_mem_var_exp(imem_var)
921 
922  tmpx2 = rg_hexa_gll_rhovp2(m,l,k,iel)*rg_hexa_gll_wkqp(imem_var,m,l,k,iel)
923  tmpx3 = 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*rg_hexa_gll_wkqs(imem_var,m,l,k,iel)
924 
925  tmpx4 = (duxdx+duydy+duzdz)*(tmpx2 - tmpx3)
926 
927  !
928  !------------->anelastic stress at step n+1/2 following s. ma and p. liu (2006) using epsnp1 and unrelaxed material modulus
929  tauxx_n12 = tmpx1*rg_hexa_gll_ksixx(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duxdx + tmpx4)
930  tauyy_n12 = tmpx1*rg_hexa_gll_ksiyy(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duydy + tmpx4)
931  tauzz_n12 = tmpx1*rg_hexa_gll_ksizz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duzdz + tmpx4)
932  tauxy_n12 = tmpx1*rg_hexa_gll_ksixy(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duxdy+duydx))
933  tauxz_n12 = tmpx1*rg_hexa_gll_ksixz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duxdz+duzdx))
934  tauyz_n12 = tmpx1*rg_hexa_gll_ksiyz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duydz+duzdy))
935  !
936  !------------->compute final stress at step n+1 according to day and minster (1984) + ma and liu (2006) : tauxx - SUM anelastic stress at step n+1
937  tauxx = tauxx - 0.5*(tauxx_n12 + rg_hexa_gll_ksixx(imem_var,m,l,k,iel))
938  tauyy = tauyy - 0.5*(tauyy_n12 + rg_hexa_gll_ksiyy(imem_var,m,l,k,iel))
939  tauzz = tauzz - 0.5*(tauzz_n12 + rg_hexa_gll_ksizz(imem_var,m,l,k,iel))
940  tauxy = tauxy - 0.5*(tauxy_n12 + rg_hexa_gll_ksixy(imem_var,m,l,k,iel))
941  tauxz = tauxz - 0.5*(tauxz_n12 + rg_hexa_gll_ksixz(imem_var,m,l,k,iel))
942  tauyz = tauyz - 0.5*(tauyz_n12 + rg_hexa_gll_ksiyz(imem_var,m,l,k,iel))
943  !
944  !------------->update memory for stress (step n+1/2)
945  rg_hexa_gll_ksixx(imem_var,m,l,k,iel) = tauxx_n12
946  rg_hexa_gll_ksiyy(imem_var,m,l,k,iel) = tauyy_n12
947  rg_hexa_gll_ksizz(imem_var,m,l,k,iel) = tauzz_n12
948  rg_hexa_gll_ksixy(imem_var,m,l,k,iel) = tauxy_n12
949  rg_hexa_gll_ksixz(imem_var,m,l,k,iel) = tauxz_n12
950  rg_hexa_gll_ksiyz(imem_var,m,l,k,iel) = tauyz_n12
951 
952  enddo
953 
954  endif
955  !
956  !---------->store members of integration of the gll node klm
957  intpx1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*dxidx+tauxy*dxidy+tauxz*dxidz)
958  intpx2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*detdx+tauxy*detdy+tauxz*detdz)
959  intpx3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*dzedx+tauxy*dzedy+tauxz*dzedz)
960 
961  intpy1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*dxidx+tauyy*dxidy+tauyz*dxidz)
962  intpy2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*detdx+tauyy*detdy+tauyz*detdz)
963  intpy3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*dzedx+tauyy*dzedy+tauyz*dzedz)
964 
965  intpz1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*dxidx+tauyz*dxidy+tauzz*dxidz)
966  intpz2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*detdx+tauyz*detdy+tauzz*detdz)
967  intpz3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*dzedx+tauyz*dzedy+tauzz*dzedz)
968  enddo !xi
969  enddo !eta
970  enddo !zeta
971 
972  !
973  !->finish integration for hexa (internal forces at step n+1)
974  do k = 1,ig_ngll
975  do l = 1,ig_ngll
976  do m = 1,ig_ngll
977 
978  tmpx1 = intpx1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
979  + intpx1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
980  + intpx1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
981  + intpx1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
982  + intpx1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5) &
983  + intpx1(6,l,k)*rg_gll_lagrange_deriv(m,6)*rg_gll_weight(6)
984 
985  tmpy1 = intpy1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
986  + intpy1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
987  + intpy1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
988  + intpy1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
989  + intpy1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5) &
990  + intpy1(6,l,k)*rg_gll_lagrange_deriv(m,6)*rg_gll_weight(6)
991 
992  tmpz1 = intpz1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
993  + intpz1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
994  + intpz1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
995  + intpz1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
996  + intpz1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5) &
997  + intpz1(6,l,k)*rg_gll_lagrange_deriv(m,6)*rg_gll_weight(6)
998 
999  tmpx2 = intpx2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
1000  + intpx2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
1001  + intpx2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
1002  + intpx2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
1003  + intpx2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5) &
1004  + intpx2(m,6,k)*rg_gll_lagrange_deriv(l,6)*rg_gll_weight(6)
1005 
1006  tmpy2 = intpy2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
1007  + intpy2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
1008  + intpy2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
1009  + intpy2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
1010  + intpy2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5) &
1011  + intpy2(m,6,k)*rg_gll_lagrange_deriv(l,6)*rg_gll_weight(6)
1012 
1013  tmpz2 = intpz2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
1014  + intpz2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
1015  + intpz2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
1016  + intpz2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
1017  + intpz2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5) &
1018  + intpz2(m,6,k)*rg_gll_lagrange_deriv(l,6)*rg_gll_weight(6)
1019 
1020  tmpx3 = intpx3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
1021  + intpx3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
1022  + intpx3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
1023  + intpx3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
1024  + intpx3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5) &
1025  + intpx3(m,l,6)*rg_gll_lagrange_deriv(k,6)*rg_gll_weight(6)
1026 
1027  tmpy3 = intpy3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
1028  + intpy3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
1029  + intpy3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
1030  + intpy3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
1031  + intpy3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5) &
1032  + intpy3(m,l,6)*rg_gll_lagrange_deriv(k,6)*rg_gll_weight(6)
1033 
1034  tmpz3 = intpz3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
1035  + intpz3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
1036  + intpz3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
1037  + intpz3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
1038  + intpz3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5) &
1039  + intpz3(m,l,6)*rg_gll_lagrange_deriv(k,6)*rg_gll_weight(6)
1040 
1041  fac1 = rg_gll_weight(l)*rg_gll_weight(k)
1042  fac2 = rg_gll_weight(m)*rg_gll_weight(k)
1043  fac3 = rg_gll_weight(m)*rg_gll_weight(l)
1044 
1045  rl_acceleration_gll(1,m,l,k) = rl_acceleration_gll(1,m,l,k) + (fac1*tmpx1 + fac2*tmpx2 + fac3*tmpx3)
1046  rl_acceleration_gll(2,m,l,k) = rl_acceleration_gll(2,m,l,k) + (fac1*tmpy1 + fac2*tmpy2 + fac3*tmpy3)
1047  rl_acceleration_gll(3,m,l,k) = rl_acceleration_gll(3,m,l,k) + (fac1*tmpz1 + fac2*tmpz2 + fac3*tmpz3)
1048 
1049  enddo
1050  enddo
1051  enddo
1052 
1053  do k = 1,ig_ngll !zeta
1054  do l = 1,ig_ngll !eta
1055  do m = 1,ig_ngll !xi
1056 
1057  igll = ig_hexa_gll_glonum(m,l,k,iel)
1058 
1059  rg_gll_acceleration(1,igll) = rg_gll_acceleration(1,igll) - rl_acceleration_gll(1,m,l,k)
1060  rg_gll_acceleration(2,igll) = rg_gll_acceleration(2,igll) - rl_acceleration_gll(2,m,l,k)
1061  rg_gll_acceleration(3,igll) = rg_gll_acceleration(3,igll) - rl_acceleration_gll(3,m,l,k)
1062 
1063  enddo
1064  enddo
1065  enddo
1066 
1067  enddo !loop on hexahedron elements
1068 
1069  return
1070 !***********************************************************************************************************************************************************************************
1071  end subroutine compute_internal_forces_order5
1072 !***********************************************************************************************************************************************************************************
1073 
1074 !
1075 !
1081 !***********************************************************************************************************************************************************************************
1082  subroutine compute_internal_forces_order6(elt_start,elt_end)
1083 !***********************************************************************************************************************************************************************************
1084 
1085  use mpi
1086 
1087  use mod_global_variables, only :&
1088  ig_ngll&
1089  ,ig_ndof&
1090  ,rg_gll_lagrange_deriv&
1091  ,rg_gll_displacement&
1092  ,rg_gll_acceleration&
1093  ,rg_gll_weight&
1094  ,rg_gll_abscissa&
1095  ,ig_nreceiver_hexa&
1096  ,ig_idt&
1097  ,ig_receiver_saving_incr&
1098  ,ig_nhexa&
1099  ,ig_hexa_gll_glonum&
1100  ,rg_hexa_gll_dxidx&
1101  ,rg_hexa_gll_dxidy&
1102  ,rg_hexa_gll_dxidz&
1103  ,rg_hexa_gll_detdx&
1104  ,rg_hexa_gll_detdy&
1105  ,rg_hexa_gll_detdz&
1106  ,rg_hexa_gll_dzedx&
1107  ,rg_hexa_gll_dzedy&
1108  ,rg_hexa_gll_dzedz&
1109  ,rg_hexa_gll_jacobian_det&
1110  ,ig_myrank&
1111  ,ig_nhexa_outer&
1112  ,rg_hexa_gll_rho&
1113  ,rg_hexa_gll_rhovs2&
1114  ,rg_hexa_gll_rhovp2&
1115  ,rg_hexa_gll_rhovs2&
1116  ,rg_hexa_gll_rhovp2&
1117  ,rg_hexa_gll_wkqs&
1118  ,rg_hexa_gll_wkqp&
1119  ,rg_hexa_gll_ksixx&
1120  ,rg_hexa_gll_ksiyy&
1121  ,rg_hexa_gll_ksizz&
1122  ,rg_hexa_gll_ksixy&
1123  ,rg_hexa_gll_ksixz&
1124  ,rg_hexa_gll_ksiyz&
1125  ,rg_relax_coeff&
1126  ,rg_mem_var_exp&
1127  ,ig_nrelax&
1128  ,rg_dt&
1129  ,lg_visco&
1130  ,ig_nhexa
1131 
1132  implicit none
1133 
1134  integer, intent(in) :: elt_start
1135  integer, intent(in) :: elt_end
1136 
1137  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx1
1138  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx2
1139  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpx3
1140  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy1
1141  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy2
1142  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpy3
1143  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz1
1144  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz2
1145  real, dimension(IG_NGLL,IG_NGLL,IG_NGLL) :: intpz3
1146  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: rl_displacement_gll
1147  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: rl_acceleration_gll
1148 
1149  real :: duxdxi
1150  real :: duxdet
1151  real :: duxdze
1152  real :: duydxi
1153  real :: duydet
1154  real :: duydze
1155  real :: duzdxi
1156  real :: duzdet
1157  real :: duzdze
1158  real :: duxdx
1159  real :: duydy
1160  real :: duzdz
1161  real :: duxdy
1162  real :: duxdz
1163  real :: duydx
1164  real :: duydz
1165  real :: duzdx
1166  real :: duzdy
1167  real :: dxidx
1168  real :: dxidy
1169  real :: dxidz
1170  real :: detdx
1171  real :: detdy
1172  real :: detdz
1173  real :: dzedx
1174  real :: dzedy
1175  real :: dzedz
1176  real :: tauxx
1177  real :: tauyy
1178  real :: tauzz
1179  real :: tauxy
1180  real :: tauxz
1181  real :: tauyz
1182  real :: tauxx_n12
1183  real :: tauyy_n12
1184  real :: tauzz_n12
1185  real :: tauxy_n12
1186  real :: tauxz_n12
1187  real :: tauyz_n12
1188  real :: trace_tau
1189  real :: tmpx1
1190  real :: tmpx2
1191  real :: tmpx3
1192  real :: tmpx4
1193  real :: tmpy1
1194  real :: tmpy2
1195  real :: tmpy3
1196  real :: tmpz1
1197  real :: tmpz2
1198  real :: tmpz3
1199  real :: fac1
1200  real :: fac2
1201  real :: fac3
1202 
1203  integer :: iel
1204  integer :: k
1205  integer :: l
1206  integer :: m
1207  integer :: igll
1208  integer :: imem_var
1209 
1210  do iel = elt_start,elt_end
1211 
1212  !
1213  !------->flush to zero local acceleration
1214  do k = 1,ig_ngll
1215  do l = 1,ig_ngll
1216  do m = 1,ig_ngll
1217 
1218  rl_acceleration_gll(1,m,l,k) = 0.0
1219  rl_acceleration_gll(2,m,l,k) = 0.0
1220  rl_acceleration_gll(3,m,l,k) = 0.0
1221 
1222  enddo
1223  enddo
1224  enddo
1225 
1226  !
1227  !------->fill local displacement
1228  do k = 1,ig_ngll !zeta
1229  do l = 1,ig_ngll !eta
1230  do m = 1,ig_ngll !xi
1231 
1232  igll = ig_hexa_gll_glonum(m,l,k,iel)
1233 
1234  rl_displacement_gll(1,m,l,k) = rg_gll_displacement(1,igll)
1235  rl_displacement_gll(2,m,l,k) = rg_gll_displacement(2,igll)
1236  rl_displacement_gll(3,m,l,k) = rg_gll_displacement(3,igll)
1237 
1238  enddo
1239  enddo
1240  enddo
1241  !
1242  !
1243  !******************************************************************************
1244  !->compute integrale at gll nodes + assemble force in global gll grid for hexa
1245  !******************************************************************************
1246  do k = 1,ig_ngll !zeta
1247  do l = 1,ig_ngll !eta
1248  do m = 1,ig_ngll !xi
1249  !
1250  !---------->derivative of displacement with respect to local coordinate xi, eta and zeta at the gll node klm
1251 
1252 
1253  duxdxi = rl_displacement_gll(1,1,l,k)*rg_gll_lagrange_deriv(1,m) &
1254  + rl_displacement_gll(1,2,l,k)*rg_gll_lagrange_deriv(2,m) &
1255  + rl_displacement_gll(1,3,l,k)*rg_gll_lagrange_deriv(3,m) &
1256  + rl_displacement_gll(1,4,l,k)*rg_gll_lagrange_deriv(4,m) &
1257  + rl_displacement_gll(1,5,l,k)*rg_gll_lagrange_deriv(5,m) &
1258  + rl_displacement_gll(1,6,l,k)*rg_gll_lagrange_deriv(6,m) &
1259  + rl_displacement_gll(1,7,l,k)*rg_gll_lagrange_deriv(7,m)
1260 
1261  duydxi = rl_displacement_gll(2,1,l,k)*rg_gll_lagrange_deriv(1,m) &
1262  + rl_displacement_gll(2,2,l,k)*rg_gll_lagrange_deriv(2,m) &
1263  + rl_displacement_gll(2,3,l,k)*rg_gll_lagrange_deriv(3,m) &
1264  + rl_displacement_gll(2,4,l,k)*rg_gll_lagrange_deriv(4,m) &
1265  + rl_displacement_gll(2,5,l,k)*rg_gll_lagrange_deriv(5,m) &
1266  + rl_displacement_gll(2,6,l,k)*rg_gll_lagrange_deriv(6,m) &
1267  + rl_displacement_gll(2,7,l,k)*rg_gll_lagrange_deriv(7,m)
1268 
1269  duzdxi = rl_displacement_gll(3,1,l,k)*rg_gll_lagrange_deriv(1,m) &
1270  + rl_displacement_gll(3,2,l,k)*rg_gll_lagrange_deriv(2,m) &
1271  + rl_displacement_gll(3,3,l,k)*rg_gll_lagrange_deriv(3,m) &
1272  + rl_displacement_gll(3,4,l,k)*rg_gll_lagrange_deriv(4,m) &
1273  + rl_displacement_gll(3,5,l,k)*rg_gll_lagrange_deriv(5,m) &
1274  + rl_displacement_gll(3,6,l,k)*rg_gll_lagrange_deriv(6,m) &
1275  + rl_displacement_gll(3,7,l,k)*rg_gll_lagrange_deriv(7,m)
1276 
1277  duxdet = rl_displacement_gll(1,m,1,k)*rg_gll_lagrange_deriv(1,l) &
1278  + rl_displacement_gll(1,m,2,k)*rg_gll_lagrange_deriv(2,l) &
1279  + rl_displacement_gll(1,m,3,k)*rg_gll_lagrange_deriv(3,l) &
1280  + rl_displacement_gll(1,m,4,k)*rg_gll_lagrange_deriv(4,l) &
1281  + rl_displacement_gll(1,m,5,k)*rg_gll_lagrange_deriv(5,l) &
1282  + rl_displacement_gll(1,m,6,k)*rg_gll_lagrange_deriv(6,l) &
1283  + rl_displacement_gll(1,m,7,k)*rg_gll_lagrange_deriv(7,l)
1284 
1285  duydet = rl_displacement_gll(2,m,1,k)*rg_gll_lagrange_deriv(1,l) &
1286  + rl_displacement_gll(2,m,2,k)*rg_gll_lagrange_deriv(2,l) &
1287  + rl_displacement_gll(2,m,3,k)*rg_gll_lagrange_deriv(3,l) &
1288  + rl_displacement_gll(2,m,4,k)*rg_gll_lagrange_deriv(4,l) &
1289  + rl_displacement_gll(2,m,5,k)*rg_gll_lagrange_deriv(5,l) &
1290  + rl_displacement_gll(2,m,6,k)*rg_gll_lagrange_deriv(6,l) &
1291  + rl_displacement_gll(2,m,7,k)*rg_gll_lagrange_deriv(7,l)
1292 
1293  duzdet = rl_displacement_gll(3,m,1,k)*rg_gll_lagrange_deriv(1,l) &
1294  + rl_displacement_gll(3,m,2,k)*rg_gll_lagrange_deriv(2,l) &
1295  + rl_displacement_gll(3,m,3,k)*rg_gll_lagrange_deriv(3,l) &
1296  + rl_displacement_gll(3,m,4,k)*rg_gll_lagrange_deriv(4,l) &
1297  + rl_displacement_gll(3,m,5,k)*rg_gll_lagrange_deriv(5,l) &
1298  + rl_displacement_gll(3,m,6,k)*rg_gll_lagrange_deriv(6,l) &
1299  + rl_displacement_gll(3,m,7,k)*rg_gll_lagrange_deriv(7,l)
1300 
1301  duxdze = rl_displacement_gll(1,m,l,1)*rg_gll_lagrange_deriv(1,k) &
1302  + rl_displacement_gll(1,m,l,2)*rg_gll_lagrange_deriv(2,k) &
1303  + rl_displacement_gll(1,m,l,3)*rg_gll_lagrange_deriv(3,k) &
1304  + rl_displacement_gll(1,m,l,4)*rg_gll_lagrange_deriv(4,k) &
1305  + rl_displacement_gll(1,m,l,5)*rg_gll_lagrange_deriv(5,k) &
1306  + rl_displacement_gll(1,m,l,6)*rg_gll_lagrange_deriv(6,k) &
1307  + rl_displacement_gll(1,m,l,7)*rg_gll_lagrange_deriv(7,k)
1308 
1309  duydze = rl_displacement_gll(2,m,l,1)*rg_gll_lagrange_deriv(1,k) &
1310  + rl_displacement_gll(2,m,l,2)*rg_gll_lagrange_deriv(2,k) &
1311  + rl_displacement_gll(2,m,l,3)*rg_gll_lagrange_deriv(3,k) &
1312  + rl_displacement_gll(2,m,l,4)*rg_gll_lagrange_deriv(4,k) &
1313  + rl_displacement_gll(2,m,l,5)*rg_gll_lagrange_deriv(5,k) &
1314  + rl_displacement_gll(2,m,l,6)*rg_gll_lagrange_deriv(6,k) &
1315  + rl_displacement_gll(2,m,l,7)*rg_gll_lagrange_deriv(7,k)
1316 
1317  duzdze = rl_displacement_gll(3,m,l,1)*rg_gll_lagrange_deriv(1,k) &
1318  + rl_displacement_gll(3,m,l,2)*rg_gll_lagrange_deriv(2,k) &
1319  + rl_displacement_gll(3,m,l,3)*rg_gll_lagrange_deriv(3,k) &
1320  + rl_displacement_gll(3,m,l,4)*rg_gll_lagrange_deriv(4,k) &
1321  + rl_displacement_gll(3,m,l,5)*rg_gll_lagrange_deriv(5,k) &
1322  + rl_displacement_gll(3,m,l,6)*rg_gll_lagrange_deriv(6,k) &
1323  + rl_displacement_gll(3,m,l,7)*rg_gll_lagrange_deriv(7,k)
1324 
1325  !
1326  !---------->derivative of displacement at step n+1 with respect to global coordinate x, y and z at the gll node klm
1327  dxidx = rg_hexa_gll_dxidx(m,l,k,iel)
1328  dxidy = rg_hexa_gll_dxidy(m,l,k,iel)
1329  dxidz = rg_hexa_gll_dxidz(m,l,k,iel)
1330  detdx = rg_hexa_gll_detdx(m,l,k,iel)
1331  detdy = rg_hexa_gll_detdy(m,l,k,iel)
1332  detdz = rg_hexa_gll_detdz(m,l,k,iel)
1333  dzedx = rg_hexa_gll_dzedx(m,l,k,iel)
1334  dzedy = rg_hexa_gll_dzedy(m,l,k,iel)
1335  dzedz = rg_hexa_gll_dzedz(m,l,k,iel)
1336 
1337  duxdx = duxdxi*dxidx + duxdet*detdx + duxdze*dzedx
1338  duxdy = duxdxi*dxidy + duxdet*detdy + duxdze*dzedy
1339  duxdz = duxdxi*dxidz + duxdet*detdz + duxdze*dzedz
1340  duydx = duydxi*dxidx + duydet*detdx + duydze*dzedx
1341  duydy = duydxi*dxidy + duydet*detdy + duydze*dzedy
1342  duydz = duydxi*dxidz + duydet*detdz + duydze*dzedz
1343  duzdx = duzdxi*dxidx + duzdet*detdx + duzdze*dzedx
1344  duzdy = duzdxi*dxidy + duzdet*detdy + duzdze*dzedy
1345  duzdz = duzdxi*dxidz + duzdet*detdz + duzdze*dzedz
1346  !
1347  !---------->compute elastic stress (elastic simulation) or unrelaxed elastic stress (viscoelastic simulation)
1348  trace_tau = (rg_hexa_gll_rhovp2(m,l,k,iel) - 2.0*rg_hexa_gll_rhovs2(m,l,k,iel))*(duxdx+duydy+duzdz)
1349  tauxx = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duxdx
1350  tauyy = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duydy
1351  tauzz = trace_tau + 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*duzdz
1352  tauxy = rg_hexa_gll_rhovs2(m,l,k,iel)*(duxdy+duydx)
1353  tauxz = rg_hexa_gll_rhovs2(m,l,k,iel)*(duxdz+duzdx)
1354  tauyz = rg_hexa_gll_rhovs2(m,l,k,iel)*(duydz+duzdy)
1355  !
1356  !---------->compute viscoelastic stress
1357 
1358  if(lg_visco) then
1359 
1360  do imem_var = 1,ig_nrelax
1361 
1362  tmpx1 = rg_mem_var_exp(imem_var)
1363 
1364  tmpx2 = rg_hexa_gll_rhovp2(m,l,k,iel)*rg_hexa_gll_wkqp(imem_var,m,l,k,iel)
1365  tmpx3 = 2.0*rg_hexa_gll_rhovs2(m,l,k,iel)*rg_hexa_gll_wkqs(imem_var,m,l,k,iel)
1366 
1367  tmpx4 = (duxdx+duydy+duzdz)*(tmpx2 - tmpx3)
1368 
1369  !
1370  !------------->anelastic stress at step n+1/2 following s. ma and p. liu (2006) using epsnp1 and unrelaxed material modulus
1371  tauxx_n12 = tmpx1*rg_hexa_gll_ksixx(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duxdx + tmpx4)
1372  tauyy_n12 = tmpx1*rg_hexa_gll_ksiyy(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duydy + tmpx4)
1373  tauzz_n12 = tmpx1*rg_hexa_gll_ksizz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*duzdz + tmpx4)
1374  tauxy_n12 = tmpx1*rg_hexa_gll_ksixy(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duxdy+duydx))
1375  tauxz_n12 = tmpx1*rg_hexa_gll_ksixz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duxdz+duzdx))
1376  tauyz_n12 = tmpx1*rg_hexa_gll_ksiyz(imem_var,m,l,k,iel) + (1.0 - tmpx1) * (tmpx3*0.5*(duydz+duzdy))
1377  !
1378  !------------->compute final stress at step n+1 according to day and minster (1984) + ma and liu (2006) : tauxx - SUM anelastic stress at step n+1
1379  tauxx = tauxx - 0.5*(tauxx_n12 + rg_hexa_gll_ksixx(imem_var,m,l,k,iel))
1380  tauyy = tauyy - 0.5*(tauyy_n12 + rg_hexa_gll_ksiyy(imem_var,m,l,k,iel))
1381  tauzz = tauzz - 0.5*(tauzz_n12 + rg_hexa_gll_ksizz(imem_var,m,l,k,iel))
1382  tauxy = tauxy - 0.5*(tauxy_n12 + rg_hexa_gll_ksixy(imem_var,m,l,k,iel))
1383  tauxz = tauxz - 0.5*(tauxz_n12 + rg_hexa_gll_ksixz(imem_var,m,l,k,iel))
1384  tauyz = tauyz - 0.5*(tauyz_n12 + rg_hexa_gll_ksiyz(imem_var,m,l,k,iel))
1385  !
1386  !------------->update memory for stress (step n+1/2)
1387  rg_hexa_gll_ksixx(imem_var,m,l,k,iel) = tauxx_n12
1388  rg_hexa_gll_ksiyy(imem_var,m,l,k,iel) = tauyy_n12
1389  rg_hexa_gll_ksizz(imem_var,m,l,k,iel) = tauzz_n12
1390  rg_hexa_gll_ksixy(imem_var,m,l,k,iel) = tauxy_n12
1391  rg_hexa_gll_ksixz(imem_var,m,l,k,iel) = tauxz_n12
1392  rg_hexa_gll_ksiyz(imem_var,m,l,k,iel) = tauyz_n12
1393 
1394  enddo
1395 
1396  endif
1397  !
1398  !---------->store members of integration of the gll node klm
1399  intpx1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*dxidx+tauxy*dxidy+tauxz*dxidz)
1400  intpx2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*detdx+tauxy*detdy+tauxz*detdz)
1401  intpx3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxx*dzedx+tauxy*dzedy+tauxz*dzedz)
1402 
1403  intpy1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*dxidx+tauyy*dxidy+tauyz*dxidz)
1404  intpy2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*detdx+tauyy*detdy+tauyz*detdz)
1405  intpy3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxy*dzedx+tauyy*dzedy+tauyz*dzedz)
1406 
1407  intpz1(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*dxidx+tauyz*dxidy+tauzz*dxidz)
1408  intpz2(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*detdx+tauyz*detdy+tauzz*detdz)
1409  intpz3(m,l,k) = rg_hexa_gll_jacobian_det(m,l,k,iel)*(tauxz*dzedx+tauyz*dzedy+tauzz*dzedz)
1410  enddo !xi
1411  enddo !eta
1412  enddo !zeta
1413 
1414  !
1415  !->finish integration for hexa (internal forces at step n+1)
1416  do k = 1,ig_ngll
1417  do l = 1,ig_ngll
1418  do m = 1,ig_ngll
1419 
1420  tmpx1 = intpx1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
1421  + intpx1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
1422  + intpx1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
1423  + intpx1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
1424  + intpx1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5) &
1425  + intpx1(6,l,k)*rg_gll_lagrange_deriv(m,6)*rg_gll_weight(6) &
1426  + intpx1(7,l,k)*rg_gll_lagrange_deriv(m,7)*rg_gll_weight(7)
1427 
1428  tmpy1 = intpy1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
1429  + intpy1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
1430  + intpy1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
1431  + intpy1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
1432  + intpy1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5) &
1433  + intpy1(6,l,k)*rg_gll_lagrange_deriv(m,6)*rg_gll_weight(6) &
1434  + intpy1(7,l,k)*rg_gll_lagrange_deriv(m,7)*rg_gll_weight(7)
1435 
1436  tmpz1 = intpz1(1,l,k)*rg_gll_lagrange_deriv(m,1)*rg_gll_weight(1) &
1437  + intpz1(2,l,k)*rg_gll_lagrange_deriv(m,2)*rg_gll_weight(2) &
1438  + intpz1(3,l,k)*rg_gll_lagrange_deriv(m,3)*rg_gll_weight(3) &
1439  + intpz1(4,l,k)*rg_gll_lagrange_deriv(m,4)*rg_gll_weight(4) &
1440  + intpz1(5,l,k)*rg_gll_lagrange_deriv(m,5)*rg_gll_weight(5) &
1441  + intpz1(6,l,k)*rg_gll_lagrange_deriv(m,6)*rg_gll_weight(6) &
1442  + intpz1(7,l,k)*rg_gll_lagrange_deriv(m,7)*rg_gll_weight(7)
1443 
1444  tmpx2 = intpx2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
1445  + intpx2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
1446  + intpx2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
1447  + intpx2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
1448  + intpx2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5) &
1449  + intpx2(m,6,k)*rg_gll_lagrange_deriv(l,6)*rg_gll_weight(6) &
1450  + intpx2(m,7,k)*rg_gll_lagrange_deriv(l,7)*rg_gll_weight(7)
1451 
1452  tmpy2 = intpy2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
1453  + intpy2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
1454  + intpy2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
1455  + intpy2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
1456  + intpy2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5) &
1457  + intpy2(m,6,k)*rg_gll_lagrange_deriv(l,6)*rg_gll_weight(6) &
1458  + intpy2(m,7,k)*rg_gll_lagrange_deriv(l,7)*rg_gll_weight(7)
1459 
1460  tmpz2 = intpz2(m,1,k)*rg_gll_lagrange_deriv(l,1)*rg_gll_weight(1) &
1461  + intpz2(m,2,k)*rg_gll_lagrange_deriv(l,2)*rg_gll_weight(2) &
1462  + intpz2(m,3,k)*rg_gll_lagrange_deriv(l,3)*rg_gll_weight(3) &
1463  + intpz2(m,4,k)*rg_gll_lagrange_deriv(l,4)*rg_gll_weight(4) &
1464  + intpz2(m,5,k)*rg_gll_lagrange_deriv(l,5)*rg_gll_weight(5) &
1465  + intpz2(m,6,k)*rg_gll_lagrange_deriv(l,6)*rg_gll_weight(6) &
1466  + intpz2(m,7,k)*rg_gll_lagrange_deriv(l,7)*rg_gll_weight(7)
1467 
1468  tmpx3 = intpx3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
1469  + intpx3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
1470  + intpx3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
1471  + intpx3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
1472  + intpx3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5) &
1473  + intpx3(m,l,6)*rg_gll_lagrange_deriv(k,6)*rg_gll_weight(6) &
1474  + intpx3(m,l,7)*rg_gll_lagrange_deriv(k,7)*rg_gll_weight(7)
1475 
1476  tmpy3 = intpy3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
1477  + intpy3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
1478  + intpy3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
1479  + intpy3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
1480  + intpy3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5) &
1481  + intpy3(m,l,6)*rg_gll_lagrange_deriv(k,6)*rg_gll_weight(6) &
1482  + intpy3(m,l,7)*rg_gll_lagrange_deriv(k,7)*rg_gll_weight(7)
1483 
1484  tmpz3 = intpz3(m,l,1)*rg_gll_lagrange_deriv(k,1)*rg_gll_weight(1) &
1485  + intpz3(m,l,2)*rg_gll_lagrange_deriv(k,2)*rg_gll_weight(2) &
1486  + intpz3(m,l,3)*rg_gll_lagrange_deriv(k,3)*rg_gll_weight(3) &
1487  + intpz3(m,l,4)*rg_gll_lagrange_deriv(k,4)*rg_gll_weight(4) &
1488  + intpz3(m,l,5)*rg_gll_lagrange_deriv(k,5)*rg_gll_weight(5) &
1489  + intpz3(m,l,6)*rg_gll_lagrange_deriv(k,6)*rg_gll_weight(6) &
1490  + intpz3(m,l,7)*rg_gll_lagrange_deriv(k,7)*rg_gll_weight(7)
1491 
1492  fac1 = rg_gll_weight(l)*rg_gll_weight(k)
1493  fac2 = rg_gll_weight(m)*rg_gll_weight(k)
1494  fac3 = rg_gll_weight(m)*rg_gll_weight(l)
1495 
1496  rl_acceleration_gll(1,m,l,k) = rl_acceleration_gll(1,m,l,k) + (fac1*tmpx1 + fac2*tmpx2 + fac3*tmpx3)
1497  rl_acceleration_gll(2,m,l,k) = rl_acceleration_gll(2,m,l,k) + (fac1*tmpy1 + fac2*tmpy2 + fac3*tmpy3)
1498  rl_acceleration_gll(3,m,l,k) = rl_acceleration_gll(3,m,l,k) + (fac1*tmpz1 + fac2*tmpz2 + fac3*tmpz3)
1499 
1500  enddo
1501  enddo
1502  enddo
1503 
1504  do k = 1,ig_ngll !zeta
1505  do l = 1,ig_ngll !eta
1506  do m = 1,ig_ngll !xi
1507 
1508  igll = ig_hexa_gll_glonum(m,l,k,iel)
1509 
1510  rg_gll_acceleration(1,igll) = rg_gll_acceleration(1,igll) - rl_acceleration_gll(1,m,l,k)
1511  rg_gll_acceleration(2,igll) = rg_gll_acceleration(2,igll) - rl_acceleration_gll(2,m,l,k)
1512  rg_gll_acceleration(3,igll) = rg_gll_acceleration(3,igll) - rl_acceleration_gll(3,m,l,k)
1513 
1514  enddo
1515  enddo
1516  enddo
1517 
1518  enddo !loop on hexahedron elements
1519 
1520  return
1521 !***********************************************************************************************************************************************************************************
1522  end subroutine compute_internal_forces_order6
1523 !***********************************************************************************************************************************************************************************
1524 
1525 !
1526 !
1530 !***********************************************************************************************************************************************************************************
1532 !***********************************************************************************************************************************************************************************
1533 
1534  use mpi
1535 
1536  use mod_global_variables, only : &
1537  rg_gll_velocity&
1538  ,rg_gll_acceleration&
1539  ,rg_dt&
1540  ,rg_newmark_gamma&
1541  ,ig_ngll&
1542  ,rg_gll_weight&
1543  ,ig_nquad_parax&
1544  ,rg_quadp_gll_jaco_det&
1545  ,rg_quadp_gll_normal&
1546  ,ig_quadp_gll_glonum&
1547  ,rg_quadp_gll_rhovp&
1548  ,rg_quadp_gll_rhovs
1549 
1550  implicit none
1551 
1552  real :: vx
1553  real :: vy
1554  real :: vz
1555  real :: jaco
1556  real :: nx
1557  real :: ny
1558  real :: nz
1559  real :: tx
1560  real :: ty
1561  real :: tz
1562  real :: vn
1563 
1564  integer :: iquad
1565  integer :: k
1566  integer :: l
1567  integer :: igll
1568 
1569  do iquad = 1,ig_nquad_parax
1570 
1571  do k = 1,ig_ngll
1572  do l = 1,ig_ngll
1573 
1574  jaco = rg_quadp_gll_jaco_det(l,k,iquad)
1575  nx = rg_quadp_gll_normal(1,l,k,iquad)
1576  ny = rg_quadp_gll_normal(2,l,k,iquad)
1577  nz = rg_quadp_gll_normal(3,l,k,iquad)
1578 
1579  igll = ig_quadp_gll_glonum(l,k,iquad)
1580 
1581  vx = rg_gll_velocity(1,igll)
1582  vy = rg_gll_velocity(2,igll)
1583  vz = rg_gll_velocity(3,igll)
1584 
1585  vn = vx*nx+vy*ny+vz*nz
1586 
1587  tx = rg_quadp_gll_rhovp(l,k,iquad)*vn*nx + rg_quadp_gll_rhovs(l,k,iquad)*(vx-vn*nx)
1588  ty = rg_quadp_gll_rhovp(l,k,iquad)*vn*ny + rg_quadp_gll_rhovs(l,k,iquad)*(vy-vn*ny)
1589  tz = rg_quadp_gll_rhovp(l,k,iquad)*vn*nz + rg_quadp_gll_rhovs(l,k,iquad)*(vz-vn*nz)
1590 
1591  rg_gll_acceleration(1,igll) = rg_gll_acceleration(1,igll) - jaco*rg_gll_weight(l)*rg_gll_weight(k)*tx
1592  rg_gll_acceleration(2,igll) = rg_gll_acceleration(2,igll) - jaco*rg_gll_weight(l)*rg_gll_weight(k)*ty
1593  rg_gll_acceleration(3,igll) = rg_gll_acceleration(3,igll) - jaco*rg_gll_weight(l)*rg_gll_weight(k)*tz
1594 
1595  enddo
1596  enddo
1597 
1598  enddo
1599 
1600  return
1601 !***********************************************************************************************************************************************************************************
1602  end subroutine compute_absorption_forces
1603 !***********************************************************************************************************************************************************************************
1604 
1609 
1610  use mpi
1611 
1612  use mod_global_variables, only :&
1613  ig_ngll&
1614  ,rg_simu_current_time&
1615  ,rg_gll_acceleration&
1616  ,ig_ndcsource&
1617  ,ig_nsfsource&
1618  ,tg_dcsource&
1619  ,tg_sfsource&
1620  ,rg_dcsource_gll_force&
1621  ,ig_hexa_gll_glonum&
1622  ,ig_idt&
1623  ,rg_dt&
1624  ,rg_dcsource_user_func&
1625  ,rg_sfsource_user_func
1626 
1628 
1629  implicit none
1630 
1631  real :: s
1632  real :: fac
1633  real :: val
1634  real, save :: val_dc
1635  real, save :: val_pf
1636 
1637  integer :: iso
1638  integer :: ipf
1639  integer :: k
1640  integer :: l
1641  integer :: m
1642  integer :: igll
1643  integer :: idir
1644 
1645 
1646  !
1647  !
1648  !*****************************************************************************************************
1649  !double couple point source located at any location
1650  !*****************************************************************************************************
1651  do iso = 1,ig_ndcsource
1652 
1653  if (tg_dcsource(iso)%icur == 3) then
1654 
1655  val = gabor(rg_simu_current_time,tg_dcsource(iso)%shift_time,tg_dcsource(iso)%rise_time,1.0)
1656 
1657  elseif (tg_dcsource(iso)%icur == 4) then
1658 
1659  val = expcos(rg_simu_current_time,tg_dcsource(iso)%shift_time,tg_dcsource(iso)%rise_time)
1660 
1661  elseif (tg_dcsource(iso)%icur == 5) then
1662 
1663  if (ig_idt == 1) val_dc = 0.0
1664  call ispli3(rg_simu_current_time,tg_dcsource(iso)%rise_time,s)
1665  val_dc = val_dc + s*rg_dt
1666  val = val_dc
1667 
1668  elseif (tg_dcsource(iso)%icur == 6) then
1669 
1670  val = ricker(rg_simu_current_time,tg_dcsource(iso)%shift_time,tg_dcsource(iso)%rise_time,1.0)
1671 
1672  elseif (tg_dcsource(iso)%icur == 7) then
1673 
1674  val = spiexp(rg_simu_current_time,tg_dcsource(iso)%rise_time,1.0)
1675 
1676  elseif (tg_dcsource(iso)%icur == 8) then
1677 
1678  val = fctanh(rg_simu_current_time,tg_dcsource(iso)%shift_time,tg_dcsource(iso)%rise_time,1.0)
1679 
1680  elseif (tg_dcsource(iso)%icur == 9) then
1681 
1682  val = fctanh_dt(rg_simu_current_time,tg_dcsource(iso)%shift_time,tg_dcsource(iso)%rise_time,1.0)
1683 
1684  elseif (tg_dcsource(iso)%icur == 10) then
1685 
1686  val = rg_dcsource_user_func(ig_idt)
1687 
1688  endif
1689 
1690  do k = 1,ig_ngll
1691  do l = 1,ig_ngll
1692  do m = 1,ig_ngll
1693  igll = ig_hexa_gll_glonum(m,l,k,tg_dcsource(iso)%iel)
1694  rg_gll_acceleration(1,igll) = rg_gll_acceleration(1,igll) + rg_dcsource_gll_force(1,m,l,k,iso)*val
1695  rg_gll_acceleration(2,igll) = rg_gll_acceleration(2,igll) + rg_dcsource_gll_force(2,m,l,k,iso)*val
1696  rg_gll_acceleration(3,igll) = rg_gll_acceleration(3,igll) + rg_dcsource_gll_force(3,m,l,k,iso)*val
1697  enddo
1698  enddo
1699  enddo
1700 
1701  enddo
1702 
1703 
1704  !
1705  !
1706  !*****************************************************************************************************
1707  !single force point source located at gll nodes
1708  !*****************************************************************************************************
1709  do ipf = 1,ig_nsfsource
1710 
1711  if (tg_sfsource(ipf)%icur == 3) then
1712 
1713  val = gabor(rg_simu_current_time,tg_sfsource(ipf)%shift_time,tg_sfsource(ipf)%rise_time,1.0)
1714 
1715  elseif (tg_sfsource(ipf)%icur == 4) then
1716 
1717  val = expcos(rg_simu_current_time,tg_sfsource(ipf)%shift_time,tg_sfsource(ipf)%rise_time)
1718 
1719  elseif (tg_sfsource(ipf)%icur == 5) then
1720 
1721  if (ig_idt == 1) val_pf = 0.0
1722  call ispli3(rg_simu_current_time,tg_sfsource(ipf)%rise_time,s)
1723  val_pf = val_pf + s*rg_dt
1724  val = val_pf
1725 
1726  elseif (tg_sfsource(ipf)%icur == 6) then
1727 
1728  val = ricker(rg_simu_current_time,tg_sfsource(ipf)%shift_time,tg_sfsource(ipf)%rise_time,1.0)
1729 
1730  elseif (tg_sfsource(ipf)%icur == 7) then
1731 
1732  val = spiexp(rg_simu_current_time,tg_sfsource(ipf)%rise_time,1.0)
1733 
1734  elseif (tg_sfsource(ipf)%icur == 8) then
1735 
1736  val = fctanh(rg_simu_current_time,tg_sfsource(ipf)%shift_time,tg_sfsource(ipf)%rise_time,1.0)
1737 
1738  elseif (tg_sfsource(ipf)%icur == 9) then
1739 
1740  val = fctanh_dt(rg_simu_current_time,tg_sfsource(ipf)%shift_time,tg_sfsource(ipf)%rise_time,1.0)
1741 
1742  elseif (tg_sfsource(ipf)%icur == 10) then
1743 
1744  val = rg_sfsource_user_func(ig_idt)
1745 
1746  endif
1747 
1748  igll = tg_sfsource(ipf)%iequ
1749  idir = tg_sfsource(ipf)%idir
1750  fac = tg_sfsource(ipf)%fac
1751  rg_gll_acceleration(idir,igll) = rg_gll_acceleration(idir,igll) + val*fac
1752 
1753  enddo
1754 
1755  return
1756 
1757 end subroutine compute_external_force
1758 
1759 end module mod_solver
subroutine, public compute_external_force()
This subroutine sets external forces of the system for double couple and single force point sources...
subroutine, public compute_internal_forces_order5(elt_start, elt_end)
This subroutine computes internal forces for spectral-elements of order 5. Stress-strain relationshi...
This module contains subroutines to compute Newmark explicit time marching scheme, external forces , internal forces and boundary traction forces of the system .
real function, public spiexp(t, tp, a)
function to compute spice project exponential source function :
real function, public gabor(t, ts, l, a)
function to compute real part of Gabor wavelet :
This module defines all global variables of EFISPEC3D. Scalar variables are initialized directly in t...
subroutine, public compute_internal_forces_order4(elt_start, elt_end)
This subroutine computes internal forces for spectral-elements of order 4. Stress-strain relationshi...
real function, public ricker(t, ts, tp, a)
function to compute order 2 Ricker wavelet :
real function, public expcos(t, ts, fp)
function to compute euroseistest project source function for case can1 :
This module contains functions and subroutines to compute tsource functions's time history...
subroutine, public compute_absorption_forces()
This subroutine computes absorption forces for any spectral-elements order. A so-called 'P1' explici...
subroutine, public ispli3(t, du, s)
function to compute order 3 spline : see Wikipedia
real function, public fctanh_dt(t, ts, tp, a)
function to compute first derivative of hyperbolic tangent function : with
real function, public fctanh(t, ts, tp, a)
function to compute hyperbolic tangent function (also called 'Bouchon pulse') : ...
subroutine, public compute_internal_forces_order6(elt_start, elt_end)
This subroutine computes internal forces for spectral-elements of order 6. Stress-strain relationshi...
subroutine, public newmark_end()
This subroutine finalizes Newmark time marching scheme at step n+1.
subroutine, public newmark_ini()
This subroutine initializes Newmark time marching scheme at step n+1.