[cig-commits] r22447 - in seismo/3D/SPECFEM3D_GLOBE/trunk: . setup src/specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Fri Jun 28 14:43:25 PDT 2013
Author: dkomati1
Date: 2013-06-28 14:43:25 -0700 (Fri, 28 Jun 2013)
New Revision: 22447
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/part1_classical.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/part1_undo_att.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/part2_classical.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/part2_undo_att.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_Dev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
Log:
merged the high-order time scheme LDDRK4-6 into the trunk
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/part1_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/part1_classical.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/part1_classical.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -93,7 +93,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,deltat,two_omega_earth, &
@@ -115,7 +116,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
endif
! Stacey absorbing boundaries
@@ -217,7 +219,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,deltat,two_omega_earth, &
@@ -238,7 +241,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
endif
do while (iphase <= 7) ! make sure the last communications are finished and processed
@@ -316,7 +320,8 @@
R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -356,7 +361,8 @@
R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
endif
! Deville routine
@@ -392,7 +398,8 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -425,7 +432,8 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
endif
! Stacey
@@ -461,7 +469,7 @@
accel_crust_mantle,sourcearrays, &
DT,t0,tshift_cmt,hdur_gaussian,ibool_crust_mantle, &
islice_selected_source,ispec_selected_source,it, &
- hdur,xi_source,eta_source,gamma_source,nu_source)
+ hdur,xi_source,eta_source,gamma_source,nu_source,istage)
! add adjoint sources only if adjoint simulation is performed for source inversion only
!! DK DK UNDO_ATTENUATION this must remain here even when SIMULATION_TYPE == 3 because it applies to array
@@ -616,7 +624,8 @@
R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -656,7 +665,8 @@
R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
endif
! Deville routine
@@ -692,7 +702,8 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -725,7 +736,8 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
endif
! assemble all the contributions between slices using MPI
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/part1_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/part1_undo_att.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/part1_undo_att.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -6,67 +6,87 @@
! Newmark time scheme update
- ! mantle
- do i=1,NGLOB_CRUST_MANTLE
- displ_crust_mantle(:,i) = displ_crust_mantle(:,i) &
- + deltat*veloc_crust_mantle(:,i) + deltatsqover2*accel_crust_mantle(:,i)
- veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) &
- + deltatover2*accel_crust_mantle(:,i)
- accel_crust_mantle(:,i) = 0._CUSTOM_REAL
- enddo
- ! outer core
- do i=1,NGLOB_OUTER_CORE
- displ_outer_core(i) = displ_outer_core(i) &
- + deltat*veloc_outer_core(i) + deltatsqover2*accel_outer_core(i)
- veloc_outer_core(i) = veloc_outer_core(i) &
- + deltatover2*accel_outer_core(i)
- accel_outer_core(i) = 0._CUSTOM_REAL
- enddo
- ! inner core
- do i=1,NGLOB_INNER_CORE
- displ_inner_core(:,i) = displ_inner_core(:,i) &
- + deltat*veloc_inner_core(:,i) + deltatsqover2*accel_inner_core(:,i)
- veloc_inner_core(:,i) = veloc_inner_core(:,i) &
- + deltatover2*accel_inner_core(:,i)
- accel_inner_core(:,i) = 0._CUSTOM_REAL
- enddo
- ! integral of strain for adjoint movie volume
- if(MOVIE_VOLUME .and. (MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3) ) then
-! do *NOT* use array syntax for that loop, otherwise you will get a compiler error when MOVIE_VOLUME is off
-! because the shape of the arrays will not match (due to some arrays purposely declared with a dummy size of 1)
- do ispec = 1,NSPEC_CRUST_MANTLE
- call compute_element_strain_undo_att_noDev(ispec,nglob_crust_mantle,nspec_crust_mantle,&
+ do istage = 1, NSTAGE_TIME_SCHEME !ZN begin loop of istage
+ if(USE_LDDRK)then
+ ! mantle
+ accel_crust_mantle(:,:) = 0._CUSTOM_REAL
+ ! outer core
+ accel_outer_core(:) = 0._CUSTOM_REAL
+ ! inner core
+ accel_inner_core(:,:) = 0._CUSTOM_REAL
+ else
+ ! mantle
+ do i=1,NGLOB_CRUST_MANTLE
+ displ_crust_mantle(:,i) = displ_crust_mantle(:,i) &
+ + deltat*veloc_crust_mantle(:,i) + deltatsqover2*accel_crust_mantle(:,i)
+ veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) &
+ + deltatover2*accel_crust_mantle(:,i)
+ accel_crust_mantle(:,i) = 0._CUSTOM_REAL
+ enddo
+ ! outer core
+ do i=1,NGLOB_OUTER_CORE
+ displ_outer_core(i) = displ_outer_core(i) &
+ + deltat*veloc_outer_core(i) + deltatsqover2*accel_outer_core(i)
+ veloc_outer_core(i) = veloc_outer_core(i) &
+ + deltatover2*accel_outer_core(i)
+ accel_outer_core(i) = 0._CUSTOM_REAL
+ enddo
+ ! inner core
+ do i=1,NGLOB_INNER_CORE
+ displ_inner_core(:,i) = displ_inner_core(:,i) &
+ + deltat*veloc_inner_core(:,i) + deltatsqover2*accel_inner_core(:,i)
+ veloc_inner_core(:,i) = veloc_inner_core(:,i) &
+ + deltatover2*accel_inner_core(:,i)
+ accel_inner_core(:,i) = 0._CUSTOM_REAL
+ enddo
+ endif
+
+ if(istage == 1)then
+ ! integral of strain for adjoint movie volume
+ if(MOVIE_VOLUME .and. (MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3) ) then
+ ! do *NOT* use array syntax for that loop, otherwise you will get a compiler error when MOVIE_VOLUME is off
+ ! because the shape of the arrays will not match (due to some arrays purposely declared with a dummy size of 1)
+ do ispec = 1,NSPEC_CRUST_MANTLE
+ call compute_element_strain_undo_att_noDev(ispec,nglob_crust_mantle,nspec_crust_mantle,&
displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,&
epsilondev_loc_crust_mantle,eps_trace_over_3_loc_crust_mantle)
- Iepsilondev_crust_mantle(:,:,:,:,ispec) = Iepsilondev_crust_mantle(:,:,:,:,ispec) &
+ Iepsilondev_crust_mantle(:,:,:,:,ispec) = Iepsilondev_crust_mantle(:,:,:,:,ispec) &
+ deltat*epsilondev_loc_crust_mantle(:,:,:,:)
- Ieps_trace_over_3_crust_mantle(:,:,:,ispec) = Ieps_trace_over_3_crust_mantle(:,:,:,ispec) &
+ Ieps_trace_over_3_crust_mantle(:,:,:,ispec) = Ieps_trace_over_3_crust_mantle(:,:,:,ispec) &
+ deltat*eps_trace_over_3_loc_crust_mantle(:,:,:)
- enddo
- endif
+ enddo
+ endif
- ! compute the maximum of the norm of the displacement
- ! in all the slices using an MPI reduction
- ! and output timestamp file to check that simulation is running fine
- if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin+4 .or. it == it_end) then
- call check_simulation_stability(it,displ_crust_mantle,displ_inner_core,displ_outer_core, &
- 1,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
- it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
+ ! compute the maximum of the norm of the displacement
+ ! in all the slices using an MPI reduction
+ ! and output timestamp file to check that simulation is running fine
+ if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin+4 .or. it == it_end) then
+ call check_simulation_stability(it,displ_crust_mantle,displ_inner_core,displ_outer_core, &
+ 1,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
+ it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
+ endif
endif
-
! ****************************************************
! big loop over all spectral elements in the fluid
! ****************************************************
-
- ! compute internal forces in the fluid region
- if(CUSTOM_REAL == SIZE_REAL) then
- time = sngl((dble(it-1)*DT-t0)*scale_t_inv)
+ if(USE_LDDRK)then
+ ! compute internal forces in the fluid region
+ if(CUSTOM_REAL == SIZE_REAL) then
+ time = sngl((dble(it-1)*DT+dble(C_LDDRK(istage))*DT-t0)*scale_t_inv)
+ else
+ time = (dble(it-1)*DT+dble(C_LDDRK(istage))*DT-t0)*scale_t_inv
+ endif
else
- time = (dble(it-1)*DT-t0)*scale_t_inv
+ ! compute internal forces in the fluid region
+ if(CUSTOM_REAL == SIZE_REAL) then
+ time = sngl((dble(it-1)*DT-t0)*scale_t_inv)
+ else
+ time = (dble(it-1)*DT-t0)*scale_t_inv
+ endif
endif
iphase = 0 ! do not start any non blocking communications at this stage
@@ -93,7 +113,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,deltat,two_omega_earth, &
@@ -115,7 +136,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
endif
! Stacey absorbing boundaries
@@ -217,7 +239,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,deltat,two_omega_earth, &
@@ -238,7 +261,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar,iphase,icall, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
endif
do while (iphase <= 7) ! make sure the last communications are finished and processed
@@ -258,10 +282,20 @@
enddo
! multiply by the inverse of the mass matrix and update velocity
- do i=1,NGLOB_OUTER_CORE
- accel_outer_core(i) = accel_outer_core(i)*rmass_outer_core(i)
- veloc_outer_core(i) = veloc_outer_core(i) + deltatover2*accel_outer_core(i)
- enddo
+ if(USE_LDDRK)then
+ do i=1,NGLOB_OUTER_CORE
+ accel_outer_core(i) = accel_outer_core(i)*rmass_outer_core(i)
+ veloc_outer_core_lddrk(i) = ALPHA_LDDRK(istage) * veloc_outer_core_lddrk(i) + deltat * accel_outer_core(i)
+ displ_outer_core_lddrk(i) = ALPHA_LDDRK(istage) * displ_outer_core_lddrk(i) + deltat * veloc_outer_core(i)
+ veloc_outer_core(i) = veloc_outer_core(i) + BETA_LDDRK(istage) * veloc_outer_core_lddrk(i)
+ displ_outer_core(i) = displ_outer_core(i) + BETA_LDDRK(istage) * displ_outer_core_lddrk(i)
+ enddo
+ else
+ do i=1,NGLOB_OUTER_CORE
+ accel_outer_core(i) = accel_outer_core(i)*rmass_outer_core(i)
+ veloc_outer_core(i) = veloc_outer_core(i) + deltatover2*accel_outer_core(i)
+ enddo
+ endif
! ------------------- new non blocking implementation -------------------
! ****************************************************
@@ -315,7 +349,8 @@
R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -355,7 +390,8 @@
R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
endif
! Deville routine
@@ -391,7 +427,8 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -424,7 +461,8 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
endif
! Stacey
@@ -460,7 +498,7 @@
accel_crust_mantle,sourcearrays, &
DT,t0,tshift_cmt,hdur_gaussian,ibool_crust_mantle, &
islice_selected_source,ispec_selected_source,it, &
- hdur,xi_source,eta_source,gamma_source,nu_source)
+ hdur,xi_source,eta_source,gamma_source,nu_source,istage)
! add adjoint sources only if adjoint simulation is performed for source inversion only
!! DK DK UNDO_ATTENUATION this must remain here even when SIMULATION_TYPE == 3 because it applies to array
@@ -603,7 +641,8 @@
R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -643,7 +682,8 @@
R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
endif
! Deville routine
@@ -679,7 +719,8 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -712,7 +753,8 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
endif
! assemble all the contributions between slices using MPI
@@ -793,9 +835,20 @@
! Newmark time scheme - corrector for elastic parts
! mantle
- do i=1,NGLOB_CRUST_MANTLE
- veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + deltatover2*accel_crust_mantle(:,i)
- enddo
+ if(USE_LDDRK)then
+ do i=1,NGLOB_CRUST_MANTLE
+ veloc_crust_mantle_lddrk(:,i) = ALPHA_LDDRK(istage) * veloc_crust_mantle_lddrk(:,i) &
+ + deltat * accel_crust_mantle(:,i)
+ displ_crust_mantle_lddrk(:,i) = ALPHA_LDDRK(istage) * displ_crust_mantle_lddrk(:,i) &
+ + deltat * veloc_crust_mantle(:,i)
+ veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + BETA_LDDRK(istage) * veloc_crust_mantle_lddrk(:,i)
+ displ_crust_mantle(:,i) = displ_crust_mantle(:,i) + BETA_LDDRK(istage) * displ_crust_mantle_lddrk(:,i)
+ enddo
+ else
+ do i=1,NGLOB_CRUST_MANTLE
+ veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + deltatover2*accel_crust_mantle(:,i)
+ enddo
+ endif
! inner core
do i=1,NGLOB_INNER_CORE
accel_inner_core(1,i) = accel_inner_core(1,i)*rmass_inner_core(i) &
@@ -804,8 +857,18 @@
- two_omega_earth*veloc_inner_core(1,i)
accel_inner_core(3,i) = accel_inner_core(3,i)*rmass_inner_core(i)
- veloc_inner_core(:,i) = veloc_inner_core(:,i) + deltatover2*accel_inner_core(:,i)
+ if(USE_LDDRK)then
+ veloc_inner_core_lddrk(:,i) = ALPHA_LDDRK(istage) * veloc_inner_core_lddrk(:,i) &
+ + deltat * accel_inner_core(:,i)
+ displ_inner_core_lddrk(:,i) = ALPHA_LDDRK(istage) * displ_inner_core_lddrk(:,i) &
+ + deltat * veloc_inner_core(:,i)
+ veloc_inner_core(:,i) = veloc_inner_core(:,i) + BETA_LDDRK(istage) * veloc_inner_core_lddrk(:,i)
+ displ_inner_core(:,i) = displ_inner_core(:,i) + BETA_LDDRK(istage) * displ_inner_core_lddrk(:,i)
+ else
+ veloc_inner_core(:,i) = veloc_inner_core(:,i) + deltatover2*accel_inner_core(:,i)
+ endif
enddo
+ enddo !ZN end loop of istage
! write the seismograms with time shift
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/part2_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/part2_classical.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/part2_classical.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -89,7 +89,8 @@
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
else
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
@@ -110,7 +111,8 @@
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
endif
endif
@@ -220,7 +222,8 @@
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
@@ -242,7 +245,8 @@
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
endif
do while (b_iphase <= 7) ! make sure the last communications are finished and processed
@@ -325,7 +329,8 @@
b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -365,7 +370,8 @@
b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
endif
endif
@@ -402,7 +408,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -435,7 +442,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
endif
endif
@@ -601,7 +609,8 @@
b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -641,7 +650,8 @@
b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
endif
! Deville routine
@@ -677,7 +687,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -710,7 +721,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
endif
! assemble all the contributions between slices using MPI
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/part2_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/part2_undo_att.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/part2_undo_att.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -93,7 +93,8 @@
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
else
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
@@ -114,7 +115,8 @@
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
endif
endif
@@ -251,7 +253,8 @@
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
else
! div_displ_outer_core is initialized to zero in the following subroutine.
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
@@ -273,7 +276,8 @@
b_buffer_send_chunkcorn_scalar,b_buffer_recv_chunkcorn_scalar,b_iphase,b_icall, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool_outer_core,MOVIE_VOLUME)
+ ibool_outer_core,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
endif
do while (b_iphase <= 7) ! make sure the last communications are finished and processed
@@ -356,7 +360,8 @@
b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -396,7 +401,8 @@
b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
endif
endif
@@ -433,7 +439,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -466,7 +473,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
endif
endif
@@ -666,7 +674,8 @@
b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -706,7 +715,8 @@
b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL)
endif
! Deville routine
@@ -742,7 +752,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -775,7 +786,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL)
endif
! assemble all the contributions between slices using MPI
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2013-06-28 21:43:25 UTC (rev 22447)
@@ -696,3 +696,19 @@
integer, parameter :: NSPEC_DOUBLING_BASICBRICK = 8
integer, parameter :: NGLOB_DOUBLING_BASICBRICK = 27
+! this for LDDRK
+ integer, parameter :: NSTAGE = 6
+ real(kind=CUSTOM_REAL), parameter, dimension(NSTAGE) :: ALPHA_LDDRK = &
+ (/0.0_CUSTOM_REAL,-0.737101392796_CUSTOM_REAL, -1.634740794341_CUSTOM_REAL,&
+ -0.744739003780_CUSTOM_REAL,-1.469897351522_CUSTOM_REAL,-2.813971388035_CUSTOM_REAL/)
+
+ real(kind=CUSTOM_REAL), parameter, dimension(NSTAGE) :: BETA_LDDRK = &
+ (/0.032918605146_CUSTOM_REAL,0.823256998200_CUSTOM_REAL,0.381530948900_CUSTOM_REAL,&
+ 0.200092213184_CUSTOM_REAL,1.718581042715_CUSTOM_REAL,0.27_CUSTOM_REAL/)
+
+ real(kind=CUSTOM_REAL), parameter, dimension(NSTAGE) :: C_LDDRK = &
+ (/0.0_CUSTOM_REAL,0.032918605146_CUSTOM_REAL,0.249351723343_CUSTOM_REAL,&
+ 0.466911705055_CUSTOM_REAL,0.582030414044_CUSTOM_REAL,0.847252983783_CUSTOM_REAL/)
+
+ logical, parameter :: USE_LDDRK = .false. ! .true.
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -29,7 +29,7 @@
accel_crust_mantle,sourcearrays, &
DT,t0,tshift_cmt,hdur_gaussian,ibool_crust_mantle, &
islice_selected_source,ispec_selected_source,it, &
- hdur,xi_source,eta_source,gamma_source,nu_source)
+ hdur,xi_source,eta_source,gamma_source,nu_source,istage)
implicit none
@@ -65,6 +65,9 @@
double precision :: f0
double precision, external :: comp_source_time_function_rickr
+!for LDDRK
+ integer :: istage
+
do isource = 1,NSOURCES
@@ -88,7 +91,12 @@
!endif
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ if(USE_LDDRK)then
+ stf_used = FACTOR_FORCE_SOURCE * &
+ comp_source_time_function_rickr(dble(it-1)*DT + dble(C_LDDRK(istage))*DT-t0-tshift_cmt(isource),f0)
+ else
+ stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ endif
! we use a force in a single direction along one of the components:
! x/y/z or E/N/Z-direction would correspond to 1/2/3 = COMPONENT_FORCE_SOURCE
@@ -97,9 +105,13 @@
+ sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
else
+ if(USE_LDDRK)then
+ stf = comp_source_time_function(dble(it-1)*DT + &
+ dble(C_LDDRK(istage))*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ else
+ stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ endif
- stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
-
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
stf_used = sngl(stf)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -157,6 +157,11 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
@@ -486,6 +491,11 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
@@ -999,6 +1009,11 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
@@ -1257,7 +1272,8 @@
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
c44store,muvstore, &
- epsilondev_loc_nplus1,epsilondev_loc)
+ epsilondev_loc_nplus1,epsilondev_loc,&
+ istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,deltat)
! crust mantle
! update memory variables based upon the Runge-Kutta scheme
@@ -1304,6 +1320,12 @@
integer :: i_memory
+! for LDDRK
+ integer :: istage
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_lddrk
+ real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+ real(kind=CUSTOM_REAL) :: deltat
+
! use Runge-Kutta scheme to march in time
! get coefficients for that standard linear solid
@@ -1320,11 +1342,21 @@
factor_common_c44_muv(:,:,:) = factor_common_c44_muv(:,:,:) * muvstore(:,:,:,ispec)
endif
- do i_memory = 1,5
- R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
- + factor_common_c44_muv(:,:,:) &
- * (betaval(i_SLS) * epsilondev_loc(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:))
- enddo
+ if(USE_LDDRK)then
+ do i_memory = 1,5
+ R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec) = ALPHA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec) + &
+ deltat * (factor_common_c44_muv(:,:,:)*epsilondev_loc(i_memory,:,:,:) - &
+ R_memory(i_memory,i_SLS,:,:,:,ispec)*(1._CUSTOM_REAL/tau_sigma_CUSTOM_REAL(i_SLS)))
+ R_memory(i_memory,i_SLS,:,:,:,ispec) = R_memory(i_memory,i_SLS,:,:,:,ispec) + &
+ BETA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec)
+ enddo
+ else
+ do i_memory = 1,5
+ R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+ + factor_common_c44_muv(:,:,:) &
+ * (betaval(i_SLS) * epsilondev_loc(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:))
+ enddo
+ endif
enddo ! i_SLS
end subroutine compute_element_att_memory_cr
@@ -1337,7 +1369,8 @@
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
muvstore, &
- epsilondev_loc_nplus1,epsilondev_loc)
+ epsilondev_loc_nplus1,epsilondev_loc,&
+ istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,deltat)
! inner core
! update memory variables based upon the Runge-Kutta scheme
@@ -1384,6 +1417,12 @@
integer :: i_memory
+! for LDDRK
+ integer :: istage
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_lddrk
+ real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+ real(kind=CUSTOM_REAL) :: deltat
+
! use Runge-Kutta scheme to march in time
! get coefficients for that standard linear solid
@@ -1392,11 +1431,28 @@
do i_SLS = 1,N_SLS
factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec)
- do i_memory = 1,5
- R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
- + muvstore(:,:,:,ispec) * factor_common_use(:,:,:) * &
- (betaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
- enddo
+! do i_memory = 1,5
+! R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+! + muvstore(:,:,:,ispec) * factor_common_use(:,:,:) * &
+! (betaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+! enddo
+
+ if(USE_LDDRK)then
+ do i_memory = 1,5
+ R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec) = ALPHA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec) + &
+ deltat * (muvstore(:,:,:,ispec) * factor_common_use(:,:,:)*epsilondev_loc(i_memory,:,:,:) - &
+ R_memory(i_memory,i_SLS,:,:,:,ispec)*(1._CUSTOM_REAL/tau_sigma_CUSTOM_REAL(i_SLS)))
+ R_memory(i_memory,i_SLS,:,:,:,ispec) = R_memory(i_memory,i_SLS,:,:,:,ispec) + &
+ BETA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec)
+ enddo
+ else
+ do i_memory = 1,5
+ R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+ + muvstore(:,:,:,ispec) * factor_common_use(:,:,:) * &
+ (betaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+ enddo
+ endif
+
enddo
end subroutine compute_element_att_memory_ic
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -53,7 +53,8 @@
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool,ispec_is_tiso, &
R_memory,one_minus_sum_beta,deltat,veloc_crust_mantle, &
- alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,PARTIAL_PHYS_DISPERSION_ONLY)
+ alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -228,6 +229,12 @@
integer NSPEC2D_BOTTOM_INNER_CORE,iend,ispec_glob
integer, dimension(NSPEC2D_BOTTOM_INNER_CORE) :: ibelm_bottom_inner_core
+!for LDDRK
+ integer :: istage
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_lddrk
+ real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+
+
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -571,7 +578,8 @@
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
c44store,muvstore, &
- epsilondev_loc_nplus1,epsilondev_loc)
+ epsilondev_loc_nplus1,epsilondev_loc,&
+ istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,deltat)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -53,7 +53,8 @@
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool,ispec_is_tiso, &
R_memory,one_minus_sum_beta,deltat,veloc_crust_mantle, &
- alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,PARTIAL_PHYS_DISPERSION_ONLY)
+ alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL)
implicit none
@@ -255,6 +256,11 @@
integer NSPEC2D_BOTTOM_INNER_CORE
integer, dimension(NSPEC2D_BOTTOM_INNER_CORE) :: ibelm_bottom_inner_core
+!for LDDRK
+ integer :: istage
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_lddrk
+ real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -382,6 +388,11 @@
ispec_strain = ispec
endif
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
@@ -933,30 +944,29 @@
factor_common_c44_muv = factor_common_c44_muv * muvstore(:,:,:,ispec)
endif
- R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * &
- R_memory(i_memory,i_SLS,:,:,:,ispec) + &
- factor_common_c44_muv * &
-!ZN (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + &
-!ZN gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
- (betaval(i_SLS) * epsilondev_loc(i_memory,:,:,:) + &
- gammaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:))
+! R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * &
+! R_memory(i_memory,i_SLS,:,:,:,ispec) + &
+! factor_common_c44_muv * &
+! (betaval(i_SLS) * epsilondev_loc(i_memory,:,:,:) + &
+! gammaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:))
+
+ if(USE_LDDRK)then
+ R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec) = ALPHA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec) + &
+ deltat * (factor_common_c44_muv(:,:,:)*epsilondev_loc(i_memory,:,:,:) - &
+ R_memory(i_memory,i_SLS,:,:,:,ispec)*(1._CUSTOM_REAL/tau_sigma_CUSTOM_REAL(i_SLS)))
+ R_memory(i_memory,i_SLS,:,:,:,ispec) = R_memory(i_memory,i_SLS,:,:,:,ispec) + &
+ BETA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec)
+ else
+ R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+ + factor_common_c44_muv(:,:,:) &
+ * (betaval(i_SLS) * epsilondev_loc(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:))
+ endif
+
enddo
enddo
endif
-! save deviatoric strain for Runge-Kutta scheme
-!ZN if(COMPUTE_AND_STORE_STRAIN) then
-!ZN !epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
-!ZN do k=1,NGLLZ
-!ZN do j=1,NGLLY
-!ZN do i=1,NGLLX
-!ZN epsilondev(:,i,j,k,ispec) = epsilondev_loc(:,i,j,k)
-!ZN enddo
-!ZN enddo
-!ZN enddo
-!ZN endif
-
enddo ! spectral element loop NSPEC_CRUST_MANTLE
end subroutine compute_forces_crust_mantle
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -49,7 +49,8 @@
kappavstore,muvstore,ibool,idoubling, &
c11store,c33store,c12store,c13store,c44store,R_memory,one_minus_sum_beta,deltat,veloc_inner_core,&
alphaval,betaval,gammaval,factor_common, &
- vx,vy,vz,vnspec,PARTIAL_PHYS_DISPERSION_ONLY)
+ vx,vy,vz,vnspec,PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -239,6 +240,11 @@
real(kind=CUSTOM_REAL) templ
+!for LDDRK
+ integer :: istage
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_lddrk
+ real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -752,7 +758,8 @@
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
muvstore, &
- epsilondev_loc_nplus1,epsilondev_loc)
+ epsilondev_loc_nplus1,epsilondev_loc,&
+ istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,deltat)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -50,7 +50,8 @@
kappavstore,muvstore,ibool,idoubling, &
c11store,c33store,c12store,c13store,c44store,R_memory,one_minus_sum_beta,deltat,veloc_inner_core,&
alphaval,betaval,gammaval,factor_common, &
- vx,vy,vz,vnspec,PARTIAL_PHYS_DISPERSION_ONLY)
+ vx,vy,vz,vnspec,PARTIAL_PHYS_DISPERSION_ONLY,&
+ istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL)
implicit none
@@ -210,6 +211,11 @@
integer, dimension(NSPEC2D_BOTTOM_INNER_CORE) :: ibelm_bottom_inner_core
real(kind=CUSTOM_REAL) templ
+!for LDDRK
+ integer :: istage
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_lddrk
+ real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -339,6 +345,11 @@
else
ispec_strain = ispec
endif
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
@@ -652,14 +663,33 @@
do i_SLS = 1,N_SLS
factor_common_use = factor_common(i_SLS,:,:,:,ispec)
- do i_memory = 1,5
- R_memory(i_memory,i_SLS,:,:,:,ispec) = &
- alphaval(i_SLS) * &
- R_memory(i_memory,i_SLS,:,:,:,ispec) + muvstore(:,:,:,ispec) * &
- factor_common_use * &
- (betaval(i_SLS) * &
- epsilondev_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
- enddo
+! do i_memory = 1,5
+! R_memory(i_memory,i_SLS,:,:,:,ispec) = &
+! alphaval(i_SLS) * &
+! R_memory(i_memory,i_SLS,:,:,:,ispec) + muvstore(:,:,:,ispec) * &
+! factor_common_use * &
+! (betaval(i_SLS) * &
+! epsilondev_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+! enddo
+ if(USE_LDDRK)then
+ do i_memory = 1,5
+ R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec) = ALPHA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec) + &
+ deltat * (muvstore(:,:,:,ispec)*factor_common_use*epsilondev_loc(i_memory,:,:,:) - &
+ R_memory(i_memory,i_SLS,:,:,:,ispec)*(1._CUSTOM_REAL/tau_sigma_CUSTOM_REAL(i_SLS)))
+ R_memory(i_memory,i_SLS,:,:,:,ispec) = R_memory(i_memory,i_SLS,:,:,:,ispec) + &
+ BETA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec)
+ enddo
+ else
+ do i_memory = 1,5
+ R_memory(i_memory,i_SLS,:,:,:,ispec) = &
+ alphaval(i_SLS) * &
+ R_memory(i_memory,i_SLS,:,:,:,ispec) + muvstore(:,:,:,ispec) * &
+ factor_common_use * &
+ (betaval(i_SLS) * &
+ epsilondev_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+ enddo
+ endif
+
enddo
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_Dev.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_Dev.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -44,7 +44,8 @@
hprime_xx,hprime_xxT, &
hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool,MOVIE_VOLUME)
+ ibool,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -166,11 +167,18 @@
integer :: computed_elements
+! for LDDRK
+ integer :: istage
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROTATION) :: &
+ A_array_rotation_lddrk,B_array_rotation_lddrk
+
! ****************************************************
! big loop over all spectral elements in the fluid
! ****************************************************
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. icall == 1) div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
+ if(istage == 1) then
+ if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. icall == 1) div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
+ endif
computed_elements = 0
@@ -394,16 +402,18 @@
+ dpotentialdzl * gzl)
endif
- ! divergence of displacement field with gravity on
- ! note: these calculations are only considered for SIMULATION_TYPE == 1 .and. SAVE_FORWARD
- ! and one has set MOVIE_VOLUME_TYPE == 4 when MOVIE_VOLUME is .true.;
- ! in case of SIMULATION_TYPE == 3, it gets overwritten by compute_kernels_outer_core()
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. MOVIE_VOLUME) then
- div_displfluid(i,j,k,ispec) = &
- minus_rho_g_over_kappa_fluid(int_radius) &
- * (dpotentialdx_with_rot * gxl &
- + dpotentialdy_with_rot * gyl &
- + dpotentialdzl * gzl)
+ if(istage == 1)then
+ ! divergence of displacement field with gravity on
+ ! note: these calculations are only considered for SIMULATION_TYPE == 1 .and. SAVE_FORWARD
+ ! and one has set MOVIE_VOLUME_TYPE == 4 when MOVIE_VOLUME is .true.;
+ ! in case of SIMULATION_TYPE == 3, it gets overwritten by compute_kernels_outer_core()
+ if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. MOVIE_VOLUME) then
+ div_displfluid(i,j,k,ispec) = &
+ minus_rho_g_over_kappa_fluid(int_radius) &
+ * (dpotentialdx_with_rot * gxl &
+ + dpotentialdy_with_rot * gyl &
+ + dpotentialdzl * gzl)
+ endif
endif
endif
@@ -474,9 +484,18 @@
! update rotation term with Euler scheme
if(ROTATION_VAL) then
- ! use the source saved above
- A_array_rotation(:,:,:,ispec) = A_array_rotation(:,:,:,ispec) + source_euler_A(:,:,:)
- B_array_rotation(:,:,:,ispec) = B_array_rotation(:,:,:,ispec) + source_euler_B(:,:,:)
+ if(USE_LDDRK)then
+ ! use the source saved above
+ A_array_rotation_lddrk(:,:,:,ispec) = ALPHA_LDDRK(istage) * A_array_rotation_lddrk(:,:,:,ispec) + source_euler_A(:,:,:)
+ A_array_rotation(:,:,:,ispec) = A_array_rotation(:,:,:,ispec) + BETA_LDDRK(istage) * A_array_rotation_lddrk(:,:,:,ispec)
+
+ B_array_rotation_lddrk(:,:,:,ispec) = ALPHA_LDDRK(istage) * B_array_rotation_lddrk(:,:,:,ispec) + source_euler_B(:,:,:)
+ B_array_rotation(:,:,:,ispec) = B_array_rotation(:,:,:,ispec) + BETA_LDDRK(istage) * B_array_rotation_lddrk(:,:,:,ispec)
+ else
+ ! use the source saved above
+ A_array_rotation(:,:,:,ispec) = A_array_rotation(:,:,:,ispec) + source_euler_A(:,:,:)
+ B_array_rotation(:,:,:,ispec) = B_array_rotation(:,:,:,ispec) + source_euler_B(:,:,:)
+ endif
endif
enddo ! spectral element loop
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_noDev.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_noDev.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -44,7 +44,8 @@
hprime_xx,hprime_yy,hprime_zz, &
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
- ibool,MOVIE_VOLUME)
+ ibool,MOVIE_VOLUME,&
+ istage,A_array_rotation_lddrk,B_array_rotation_lddrk)
implicit none
@@ -143,12 +144,18 @@
integer :: computed_elements
+! for LDDRK
+ integer :: istage
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROTATION) :: &
+ A_array_rotation_lddrk,B_array_rotation_lddrk
+
! ****************************************************
! big loop over all spectral elements in the fluid
! ****************************************************
+ if(istage == 1) then
+ if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. icall == 1) div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
+ endif
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. icall == 1) div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
-
computed_elements = 0
do ispec = 1,NSPEC_OUTER_CORE
@@ -339,14 +346,16 @@
dpotentialdy_with_rot * gyl + dpotentialdzl * gzl)
endif
- ! divergence of displacement field with gravity on
- ! note: these calculations are only considered for SIMULATION_TYPE == 1 .and. SAVE_FORWARD
- ! and one has set MOVIE_VOLUME_TYPE == 4 when MOVIE_VOLUME is .true.;
- ! in case of SIMULATION_TYPE == 3, it gets overwritten by compute_kernels_outer_core()
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. MOVIE_VOLUME ) then
- div_displfluid(i,j,k,ispec) = &
- minus_rho_g_over_kappa_fluid(int_radius) * (dpotentialdx_with_rot * gxl + &
- dpotentialdy_with_rot * gyl + dpotentialdzl * gzl)
+ if(istage == 1) then
+ ! divergence of displacement field with gravity on
+ ! note: these calculations are only considered for SIMULATION_TYPE == 1 .and. SAVE_FORWARD
+ ! and one has set MOVIE_VOLUME_TYPE == 4 when MOVIE_VOLUME is .true.;
+ ! in case of SIMULATION_TYPE == 3, it gets overwritten by compute_kernels_outer_core()
+ if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. MOVIE_VOLUME ) then
+ div_displfluid(i,j,k,ispec) = &
+ minus_rho_g_over_kappa_fluid(int_radius) * (dpotentialdx_with_rot * gxl + &
+ dpotentialdy_with_rot * gyl + dpotentialdzl * gzl)
+ endif
endif
endif
@@ -386,9 +395,18 @@
! update rotation term with Euler scheme
if(ROTATION_VAL) then
- ! use the source saved above
- A_array_rotation(:,:,:,ispec) = A_array_rotation(:,:,:,ispec) + source_euler_A(:,:,:)
- B_array_rotation(:,:,:,ispec) = B_array_rotation(:,:,:,ispec) + source_euler_B(:,:,:)
+ if(USE_LDDRK)then
+ ! use the source saved above
+ A_array_rotation_lddrk(:,:,:,ispec) = ALPHA_LDDRK(istage) * A_array_rotation_lddrk(:,:,:,ispec) + source_euler_A(:,:,:)
+ A_array_rotation(:,:,:,ispec) = A_array_rotation(:,:,:,ispec) + BETA_LDDRK(istage) * A_array_rotation_lddrk(:,:,:,ispec)
+
+ B_array_rotation_lddrk(:,:,:,ispec) = ALPHA_LDDRK(istage) * B_array_rotation_lddrk(:,:,:,ispec) + source_euler_B(:,:,:)
+ B_array_rotation(:,:,:,ispec) = B_array_rotation(:,:,:,ispec) + BETA_LDDRK(istage) * B_array_rotation_lddrk(:,:,:,ispec)
+ else
+ ! use the source saved above
+ A_array_rotation(:,:,:,ispec) = A_array_rotation(:,:,:,ispec) + source_euler_A(:,:,:)
+ B_array_rotation(:,:,:,ispec) = B_array_rotation(:,:,:,ispec) + source_euler_B(:,:,:)
+ endif
endif
enddo ! spectral element loop
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -108,6 +108,12 @@
one_minus_sum_beta = one_minus_sum_beta - beta(i)
enddo
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we have put the 1/2 factor there we need to remove it
+!ZN from the expression in which we use the strain here in the code.
+!ZN This is why here Brian Savage multiplies beta(:) * tauinv(:) by 2.0 to compensate for the 1/2 factor used before
factor_common(:) = 2.0d0 * beta(:) * tauinv(:)
end subroutine get_attenuation_property_values
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -548,7 +548,7 @@
c11store_inner_core,c12store_inner_core,c13store_inner_core, &
c33store_inner_core,c44store_inner_core, &
alphaval,betaval,gammaval,b_alphaval,b_betaval,b_gammaval, &
- deltat,b_deltat,LOCAL_PATH)
+ deltat,b_deltat,LOCAL_PATH,tau_sigma_dble)
! precomputes attenuation factors
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2013-06-28 21:32:46 UTC (rev 22446)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2013-06-28 21:43:25 UTC (rev 22447)
@@ -940,6 +940,21 @@
logical :: you_can_start_doing_IOs
#endif
+!this is for LDDRK
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
+ displ_crust_mantle_lddrk,veloc_crust_mantle_lddrk
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
+ displ_outer_core_lddrk,veloc_outer_core_lddrk
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
+ displ_inner_core_lddrk,veloc_inner_core_lddrk
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ A_array_rotation_lddrk,B_array_rotation_lddrk
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:,:,:), allocatable :: R_memory_crust_mantle_lddrk
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:,:,:), allocatable :: R_memory_inner_core_lddrk
+ integer :: NSTAGE_TIME_SCHEME,istage
+ double precision, dimension(N_SLS) :: tau_sigma_dble
+ real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
+
integer msg_status(MPI_STATUS_SIZE)
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: b_displ_crust_mantle_store_buffer
@@ -955,6 +970,7 @@
undo_att_sim_type_3 = .false.
+
! *************************************************
! ************** PROGRAM STARTS HERE **************
! *************************************************
@@ -1921,7 +1937,12 @@
c11store_inner_core,c12store_inner_core,c13store_inner_core, &
c33store_inner_core,c44store_inner_core, &
alphaval,betaval,gammaval,b_alphaval,b_betaval,b_gammaval, &
- deltat,b_deltat,LOCAL_PATH)
+ deltat,b_deltat,LOCAL_PATH,tau_sigma_dble)
+ if(CUSTOM_REAL == SIZE_REAL) then
+ tau_sigma_CUSTOM_REAL(:) = sngl(tau_sigma_dble(:))
+ else
+ tau_sigma_CUSTOM_REAL(:) = tau_sigma_dble(:)
+ endif
if(UNDO_ATTENUATION) then
b_alphaval = alphaval
@@ -1966,6 +1987,97 @@
veloc_inner_core(:,:) = 0._CUSTOM_REAL
accel_inner_core(:,:) = 0._CUSTOM_REAL
+ if(USE_LDDRK)then
+ if(SIMULATION_TYPE /= 1 .or. SAVE_FORWARD .or. NOISE_TOMOGRAPHY /= 0) &
+ stop 'error: LDDRK is not implemented for adjoint tomography'
+ NSTAGE_TIME_SCHEME = 6
+ allocate(displ_crust_mantle_lddrk(NDIM,NGLOB_CRUST_MANTLE),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array displ_crust_mantle_lddrk'
+ allocate(veloc_crust_mantle_lddrk(NDIM,NGLOB_CRUST_MANTLE),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array veloc_crust_mantle_lddrk'
+ allocate(displ_outer_core_lddrk(NGLOB_OUTER_CORE),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array displ_outer_core_lddrk'
+ allocate(veloc_outer_core_lddrk(NGLOB_OUTER_CORE),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array veloc_outer_core_lddrk'
+ allocate(displ_inner_core_lddrk(NDIM,NGLOB_INNER_CORE),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array displ_inner_core_lddrk'
+ allocate(veloc_inner_core_lddrk(NDIM,NGLOB_INNER_CORE),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array veloc_inner_core_lddrk'
+
+ displ_crust_mantle_lddrk(:,:) = 0._CUSTOM_REAL
+ veloc_crust_mantle_lddrk(:,:) = 0._CUSTOM_REAL
+ displ_outer_core_lddrk(:) = 0._CUSTOM_REAL
+ veloc_outer_core_lddrk(:) = 0._CUSTOM_REAL
+ displ_inner_core_lddrk(:,:) = 0._CUSTOM_REAL
+ veloc_inner_core_lddrk(:,:) = 0._CUSTOM_REAL
+
+ allocate(A_array_rotation_lddrk(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROTATION),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array A_array_rotation_lddrk'
+ allocate(B_array_rotation_lddrk(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROTATION),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array B_array_rotation_lddrk'
+
+ A_array_rotation_lddrk(:,:,:,:) = 0._CUSTOM_REAL
+ B_array_rotation_lddrk(:,:,:,:) = 0._CUSTOM_REAL
+
+ allocate(R_memory_crust_mantle_lddrk(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array R_memory_crust_mantle_lddrk'
+ allocate(R_memory_inner_core_lddrk(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array R_memory_inner_core_lddrk'
+
+ R_memory_crust_mantle_lddrk(:,:,:,:,:,:) = 0._CUSTOM_REAL
+ R_memory_inner_core_lddrk(:,:,:,:,:,:) = 0._CUSTOM_REAL
+
+ if(FIX_UNDERFLOW_PROBLEM) then
+ displ_crust_mantle_lddrk(:,:) = VERYSMALLVAL
+ veloc_crust_mantle_lddrk(:,:) = VERYSMALLVAL
+ displ_outer_core_lddrk(:) = VERYSMALLVAL
+ veloc_outer_core_lddrk(:) = VERYSMALLVAL
+ displ_inner_core_lddrk(:,:) = VERYSMALLVAL
+ veloc_inner_core_lddrk(:,:) = VERYSMALLVAL
+ A_array_rotation_lddrk(:,:,:,:) = VERYSMALLVAL
+ B_array_rotation_lddrk(:,:,:,:) = VERYSMALLVAL
+ R_memory_crust_mantle_lddrk(:,:,:,:,:,:) = VERYSMALLVAL
+ R_memory_inner_core_lddrk(:,:,:,:,:,:) = VERYSMALLVAL
+ endif
+ else
+ NSTAGE_TIME_SCHEME = 1
+ allocate(displ_crust_mantle_lddrk(NDIM,1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array displ_crust_mantle_lddrk'
+ allocate(veloc_crust_mantle_lddrk(NDIM,1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array veloc_crust_mantle_lddrk'
+ allocate(displ_outer_core_lddrk(1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array displ_outer_core_lddrk'
+ allocate(veloc_outer_core_lddrk(1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array veloc_outer_core_lddrk'
+ allocate(displ_inner_core_lddrk(NDIM,1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array displ_inner_core_lddrk'
+ allocate(veloc_inner_core_lddrk(NDIM,1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array veloc_inner_core_lddrk'
+
+ displ_crust_mantle_lddrk(:,:) = 0._CUSTOM_REAL
+ veloc_crust_mantle_lddrk(:,:) = 0._CUSTOM_REAL
+ displ_outer_core_lddrk(:) = 0._CUSTOM_REAL
+ veloc_outer_core_lddrk(:) = 0._CUSTOM_REAL
+ displ_inner_core_lddrk(:,:) = 0._CUSTOM_REAL
+ veloc_inner_core_lddrk(:,:) = 0._CUSTOM_REAL
+
+ allocate(A_array_rotation_lddrk(NGLLX,NGLLY,NGLLZ,1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array A_array_rotation_lddrk'
+ allocate(B_array_rotation_lddrk(NGLLX,NGLLY,NGLLZ,1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array B_array_rotation_lddrk'
+
+ A_array_rotation_lddrk(:,:,:,:) = 0._CUSTOM_REAL
+ B_array_rotation_lddrk(:,:,:,:) = 0._CUSTOM_REAL
+
+ allocate(R_memory_crust_mantle_lddrk(5,N_SLS,NGLLX,NGLLY,NGLLZ,1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array R_memory_crust_mantle_lddrk'
+ allocate(R_memory_inner_core_lddrk(5,N_SLS,NGLLX,NGLLY,NGLLZ,1),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array R_memory_inner_core_lddrk'
+
+ R_memory_crust_mantle_lddrk(:,:,:,:,:,:) = 0._CUSTOM_REAL
+ R_memory_inner_core_lddrk(:,:,:,:,:,:) = 0._CUSTOM_REAL
+ endif
+
! put negligible initial value to avoid very slow underflow trapping
if(FIX_UNDERFLOW_PROBLEM) then
displ_crust_mantle(:,:) = VERYSMALLVAL
@@ -2161,7 +2273,11 @@
!! DK DK this first part handles the cases SIMULATION_TYPE == 1 and SIMULATION_TYPE == 2
!! DK DK it also handles the cases NOISE_TOMOGRAPHY == 1 and NOISE_TOMOGRAPHY == 2
!! DK DK
- include "part1_classical.f90"
+ if(USE_LDDRK)then
+ include "part1_undo_att.f90"
+ else
+ include "part1_classical.f90"
+ endif
!! DK DK
!! DK DK this first part handles the case SIMULATION_TYPE == 3
@@ -2397,7 +2513,6 @@
displ_inner_core,veloc_inner_core,accel_inner_core, &
displ_outer_core,veloc_outer_core,accel_outer_core, &
R_memory_crust_mantle,R_memory_inner_core, &
-!ZN epsilondev_crust_mantle,epsilondev_inner_core, &
A_array_rotation,B_array_rotation,LOCAL_PATH)
! synchronize all processes
More information about the CIG-COMMITS
mailing list