[cig-commits] r21235 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Wed Jan 16 08:10:21 PST 2013
Author: xie.zhinan
Date: 2013-01-16 08:10:21 -0800 (Wed, 16 Jan 2013)
New Revision: 21235
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
Log:
clean the code a little
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90 2013-01-16 14:06:25 UTC (rev 21234)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90 2013-01-16 16:10:21 UTC (rev 21235)
@@ -222,6 +222,9 @@
! loop on all the standard linear solids
do i_sls = 1,N_SLS
+ phinu2 = phi_nu2(i,j,ispec,i_sls)
+ tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
+
! evolution e1 ! no need since we are just considering shear attenuation
! if(stage_time_scheme == 1) then
! Un = e11(i,j,ispec,i_sls)
@@ -239,23 +242,17 @@
! endif
if(stage_time_scheme == 1) then
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) &
- + deltat*e11_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e11_accel(i,j,ispec,i_sls)
+ e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) + deltat*e11_veloc(i,j,ispec,i_sls) &
+ + deltatsquareover2*e11_accel(i,j,ispec,i_sls)
e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- e11_accel(i,j,ispec,i_sls) = ZERO
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
- e11_accel(i,j,ispec,i_sls) = (dvx_dxl_n(i,j,ispec)-theta_n_v/TWO) * phinu2- &
- e11_veloc(i,j,ispec,i_sls)*tauinvnu2
- e11_accel(i,j,ispec,i_sls) = e11_accel(i,j,ispec,i_sls)/(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
+ e11_accel(i,j,ispec,i_sls) = ((dvx_dxl_n(i,j,ispec)-theta_n_v/TWO) * phinu2- &
+ e11_veloc(i,j,ispec,i_sls)*tauinvnu2) &
+ /(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
endif
if(stage_time_scheme == 6) then
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
+ deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) &
- deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
@@ -263,8 +260,6 @@
endif
if(stage_time_scheme == 4) then
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2- &
e11(i,j,ispec,i_sls) * tauinvnu2)
@@ -287,23 +282,17 @@
! evolution e13
if(stage_time_scheme == 1) then
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) &
- + deltat*e13_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
+ e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) + deltat*e13_veloc(i,j,ispec,i_sls) &
+ + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- e13_accel(i,j,ispec,i_sls) = ZERO
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
- e13_accel(i,j,ispec,i_sls) = (dvx_dzl_n(i,j,ispec) + dvz_dxl_n(i,j,ispec)) * phinu2- &
- e13_veloc(i,j,ispec,i_sls)*tauinvnu2
- e13_accel(i,j,ispec,i_sls) = e13_accel(i,j,ispec,i_sls)/(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
+ e13_accel(i,j,ispec,i_sls) = ((dvx_dzl_n(i,j,ispec) + dvz_dxl_n(i,j,ispec)) * phinu2- &
+ e13_veloc(i,j,ispec,i_sls)*tauinvnu2) &
+ /(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
endif
if(stage_time_scheme == 6) then
- phinu2=phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) &
+ deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) &
- deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
@@ -311,8 +300,6 @@
endif
if(stage_time_scheme == 4) then
- phinu2=phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2- &
e13(i,j,ispec,i_sls) * tauinvnu2)
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90 2013-01-16 14:06:25 UTC (rev 21234)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90 2013-01-16 16:10:21 UTC (rev 21235)
@@ -220,9 +220,13 @@
do i=1,NGLLX
theta_n_u = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
theta_n_v = dvx_dxl_n(i,j,ispec) + dvz_dzl_n(i,j,ispec)
+
! loop on all the standard linear solids
do i_sls = 1,N_SLS
+ phinu2 = phi_nu2(i,j,ispec,i_sls)
+ tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
+
! evolution e1 ! no need since we are just considering shear attenuation
! if(stage_time_scheme == 1) then
! Un = e1(i,j,ispec,i_sls)
@@ -241,22 +245,16 @@
! evolution e11
if(stage_time_scheme == 1) then
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) &
- + deltat*e11_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e11_accel(i,j,ispec,i_sls)
+ e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) + deltat*e11_veloc(i,j,ispec,i_sls) &
+ + deltatsquareover2*e11_accel(i,j,ispec,i_sls)
e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- e11_accel(i,j,ispec,i_sls) = ZERO
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
- e11_accel(i,j,ispec,i_sls) = (dvx_dxl_n(i,j,ispec)-theta_n_v/TWO) * phinu2- &
- e11_veloc(i,j,ispec,i_sls)*tauinvnu2
- e11_accel(i,j,ispec,i_sls) = e11_accel(i,j,ispec,i_sls)/(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
+ e11_accel(i,j,ispec,i_sls) = ((dvx_dxl_n(i,j,ispec)-theta_n_v/TWO) * phinu2- &
+ e11_veloc(i,j,ispec,i_sls)*tauinvnu2) &
+ /(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
endif
if(stage_time_scheme == 6) then
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
+ deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) &
- deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
@@ -264,8 +262,6 @@
endif
if(stage_time_scheme == 4) then
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2- &
e11(i,j,ispec,i_sls) * tauinvnu2)
@@ -288,23 +284,17 @@
! evolution e13
if(stage_time_scheme == 1) then
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) &
- + deltat*e13_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
+ e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) + deltat*e13_veloc(i,j,ispec,i_sls) &
+ + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- e13_accel(i,j,ispec,i_sls) = ZERO
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
- e13_accel(i,j,ispec,i_sls) = (dvx_dzl_n(i,j,ispec) + dvz_dxl_n(i,j,ispec)) * phinu2- &
- e13_veloc(i,j,ispec,i_sls)*tauinvnu2
- e13_accel(i,j,ispec,i_sls) = e13_accel(i,j,ispec,i_sls)/(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
+ e13_accel(i,j,ispec,i_sls) = ((dvx_dzl_n(i,j,ispec) + dvz_dxl_n(i,j,ispec)) * phinu2- &
+ e13_veloc(i,j,ispec,i_sls)*tauinvnu2) &
+ /(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
endif
if(stage_time_scheme == 6) then
- phinu2=phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) &
+ deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) &
- deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
@@ -312,8 +302,6 @@
endif
if(stage_time_scheme == 4) then
- phinu2=phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2- &
e13(i,j,ispec,i_sls) * tauinvnu2)
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-01-16 14:06:25 UTC (rev 21234)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-01-16 16:10:21 UTC (rev 21235)
@@ -276,31 +276,28 @@
! loop on all the standard linear solids
do i_sls = 1,N_SLS
+ phinu1 = phi_nu1(i,j,ispec,i_sls)
+ tauinvnu1 = inv_tau_sigma_nu1(i,j,ispec,i_sls)
+ phinu2 = phi_nu2(i,j,ispec,i_sls)
+ tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
+
! evolution e1
if(stage_time_scheme == 1) then
- e1(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls) &
- + deltat*e1_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e1_accel(i,j,ispec,i_sls)
+ e1(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls) + deltat*e1_veloc(i,j,ispec,i_sls) &
+ + deltatsquareover2*e1_accel(i,j,ispec,i_sls)
e1_veloc(i,j,ispec,i_sls) = e1_veloc(i,j,ispec,i_sls) + deltatover2*e1_accel(i,j,ispec,i_sls)
- phinu1 = phi_nu1(i,j,ispec,i_sls)
- tauinvnu1 = inv_tau_sigma_nu1(i,j,ispec,i_sls)
e1_accel(i,j,ispec,i_sls) = (theta_n_v * phinu1 - e1_veloc(i,j,ispec,i_sls) * tauinvnu1) / &
(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu1*deltat)
e1_veloc(i,j,ispec,i_sls) = e1_veloc(i,j,ispec,i_sls) + deltatover2*e1_accel(i,j,ispec,i_sls)
endif
if(stage_time_scheme == 6) then
- tauinvnu1 = inv_tau_sigma_nu1(i,j,ispec,i_sls)
- phinu1 = phi_nu1(i,j,ispec,i_sls)
e1_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls) + &
deltat * (theta_n_u * phinu1 - e1(i,j,ispec,i_sls) * tauinvnu1)
e1(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls) + beta_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls)
endif
if(stage_time_scheme == 4) then
-
- tauinvnu1 = inv_tau_sigma_nu1(i,j,ispec,i_sls)
- phinu1 = phi_nu1(i,j,ispec,i_sls)
e1_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (theta_n_u * phinu1 - e1(i,j,ispec,i_sls) * tauinvnu1)
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
@@ -326,12 +323,9 @@
! evolution e11
if(stage_time_scheme == 1) then
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) &
- + deltat*e11_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e11_accel(i,j,ispec,i_sls)
+ e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls) + deltat*e11_veloc(i,j,ispec,i_sls) &
+ + deltatsquareover2*e11_accel(i,j,ispec,i_sls)
e11_veloc(i,j,ispec,i_sls) = e11_veloc(i,j,ispec,i_sls) + deltatover2*e11_accel(i,j,ispec,i_sls)
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e11_accel(i,j,ispec,i_sls) = ((dvx_dxl_n(i,j,ispec)-theta_n_v/TWO) * phinu2- &
e11_veloc(i,j,ispec,i_sls)*tauinvnu2) / &
(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
@@ -339,8 +333,6 @@
endif
if(stage_time_scheme == 6) then
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
+ deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2) &
- deltat * (e11(i,j,ispec,i_sls) * tauinvnu2)
@@ -348,8 +340,6 @@
endif
if(stage_time_scheme == 4) then
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dxl_n(i,j,ispec)-theta_n_u/TWO) * phinu2- &
e11(i,j,ispec,i_sls) * tauinvnu2)
@@ -371,12 +361,9 @@
! evolution e13
if(stage_time_scheme == 1) then
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) &
- + deltat*e13_veloc(i,j,ispec,i_sls) &
- + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
+ e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls) + deltat*e13_veloc(i,j,ispec,i_sls) &
+ + deltatsquareover2*e13_accel(i,j,ispec,i_sls)
e13_veloc(i,j,ispec,i_sls) = e13_veloc(i,j,ispec,i_sls) + deltatover2*e13_accel(i,j,ispec,i_sls)
- phinu2 = phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e13_accel(i,j,ispec,i_sls) = ((dvx_dzl_n(i,j,ispec) + dvz_dxl_n(i,j,ispec)) * phinu2- &
e13_veloc(i,j,ispec,i_sls)*tauinvnu2) / &
(1._CUSTOM_REAL + 0.5_CUSTOM_REAL*tauinvnu2*deltat)
@@ -385,17 +372,13 @@
if(stage_time_scheme == 6) then
- phinu2=phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls) &
- + deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) &
- - deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
+ + deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2) &
+ - deltat * (e13(i,j,ispec,i_sls) * tauinvnu2)
e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)
endif
if(stage_time_scheme == 4) then
- phinu2=phi_nu2(i,j,ispec,i_sls)
- tauinvnu2 = inv_tau_sigma_nu2(i,j,ispec,i_sls)
e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * ((dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec))*phinu2- &
e13(i,j,ispec,i_sls) * tauinvnu2)
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
More information about the CIG-COMMITS
mailing list