[cig-commits] r20731 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Wed Sep 19 05:40:14 PDT 2012
Author: surendra
Date: 2012-09-19 05:40:14 -0700 (Wed, 19 Sep 2012)
New Revision: 20731
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Processors without fault nodes now only call the assembling routines
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-09-19 12:39:55 UTC (rev 20730)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-09-19 12:40:14 UTC (rev 20731)
@@ -297,28 +297,28 @@
read(IIN_BIN) bc%coord(1,:)
read(IIN_BIN) bc%coord(2,:)
read(IIN_BIN) bc%coord(3,:)
- endif
- bc%dt = dt_tmp
+ bc%dt = dt_tmp
- allocate( bc%B(bc%nglob) )
- bc%B = 0e0_CUSTOM_REAL
- allocate( nx(bc%nglob),ny(bc%nglob),nz(bc%nglob) )
- nx = 0e0_CUSTOM_REAL
- ny = 0e0_CUSTOM_REAL
- nz = 0e0_CUSTOM_REAL
- do e=1,bc%nspec
- do ij = 1,NGLLSQUARE
- k = ibool1(ij,e)
- nx(k) = nx(k) + normal(1,ij,e)
- ny(k) = ny(k) + normal(2,ij,e)
- nz(k) = nz(k) + normal(3,ij,e)
- bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
+ allocate( bc%B(bc%nglob) )
+ bc%B = 0e0_CUSTOM_REAL
+ allocate( nx(bc%nglob),ny(bc%nglob),nz(bc%nglob) )
+ nx = 0e0_CUSTOM_REAL
+ ny = 0e0_CUSTOM_REAL
+ nz = 0e0_CUSTOM_REAL
+ do e=1,bc%nspec
+ do ij = 1,NGLLSQUARE
+ k = ibool1(ij,e)
+ nx(k) = nx(k) + normal(1,ij,e)
+ ny(k) = ny(k) + normal(2,ij,e)
+ nz(k) = nz(k) + normal(3,ij,e)
+ bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
+ enddo
enddo
- enddo
+ endif
if(PARALLEL_FAULT) then
accel=0._CUSTOM_REAL
- accel(1,bc%ibulk1) = bc%B(:)
+ if (bc%nspec>0) accel(1,bc%ibulk1) = bc%B(:)
! assembles with other MPI processes
call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,accel, &
@@ -326,12 +326,14 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
my_neighbours_ext_mesh)
- bc%B(:) = accel(1,bc%ibulk1)
+ if (bc%nspec>0) bc%B(:) = accel(1,bc%ibulk1)
accel=0._CUSTOM_REAL
- accel(1,bc%ibulk1) = nx(:)
- accel(2,bc%ibulk1) = ny(:)
- accel(3,bc%ibulk1) = nz(:)
+ if (bc%nspec>0) then
+ accel(1,bc%ibulk1) = nx(:)
+ accel(2,bc%ibulk1) = ny(:)
+ accel(3,bc%ibulk1) = nz(:)
+ endif
! assembles with other MPI processes
call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,accel, &
@@ -339,278 +341,282 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
my_neighbours_ext_mesh)
- nx(:) = accel(1,bc%ibulk1)
- ny(:) = accel(2,bc%ibulk1)
- nz(:) = accel(3,bc%ibulk1)
+ if (bc%nspec>0) then
+ nx(:) = accel(1,bc%ibulk1)
+ ny(:) = accel(2,bc%ibulk1)
+ nz(:) = accel(3,bc%ibulk1)
+ endif
endif
- do k=1,bc%nglob
- norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
- nx(k) = nx(k) / norm
- ny(k) = ny(k) / norm
- nz(k) = nz(k) / norm
- enddo
+ if (bc%nspec>0) then
+ do k=1,bc%nglob
+ norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
+ nx(k) = nx(k) / norm
+ ny(k) = ny(k) / norm
+ nz(k) = nz(k) / norm
+ enddo
- allocate( bc%R(3,3,bc%nglob) )
- call compute_R(bc%R,bc%nglob,nx,ny,nz)
+ allocate( bc%R(3,3,bc%nglob) )
+ call compute_R(bc%R,bc%nglob,nx,ny,nz)
- ! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
- allocate(bc%invM1(bc%nglob))
- allocate(bc%invM2(bc%nglob))
- bc%invM1 = Minv(bc%ibulk1)
- bc%invM2 = Minv(bc%ibulk2)
+ ! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
+ allocate(bc%invM1(bc%nglob))
+ allocate(bc%invM2(bc%nglob))
+ bc%invM1 = Minv(bc%ibulk1)
+ bc%invM2 = Minv(bc%ibulk2)
- ! Fault impedance, Z in : Trac=T_Stick-Z*dV
- ! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
- ! T_stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
- ! NOTE: same Bi on both sides, see note above
- allocate(bc%Z(bc%nglob))
- bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
- ! WARNING: In non-split nodes at fault edges M is assembled across the fault.
- ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
- ! In a symmetric mesh (M1=M2) Z will be twice its intended value
+ ! Fault impedance, Z in : Trac=T_Stick-Z*dV
+ ! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
+ ! T_stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
+ ! NOTE: same Bi on both sides, see note above
+ allocate(bc%Z(bc%nglob))
+ bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
+ ! WARNING: In non-split nodes at fault edges M is assembled across the fault.
+ ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
+ ! In a symmetric mesh (M1=M2) Z will be twice its intended value
- allocate(bc%T(3,bc%nglob))
- allocate(bc%D(3,bc%nglob))
- allocate(bc%V(3,bc%nglob))
- bc%T = 0e0_CUSTOM_REAL
- bc%D = 0e0_CUSTOM_REAL
- bc%V = 0e0_CUSTOM_REAL
+ allocate(bc%T(3,bc%nglob))
+ allocate(bc%D(3,bc%nglob))
+ allocate(bc%V(3,bc%nglob))
+ bc%T = 0e0_CUSTOM_REAL
+ bc%D = 0e0_CUSTOM_REAL
+ bc%V = 0e0_CUSTOM_REAL
- ! Set initial fault stresses
- allocate(bc%T0(3,bc%nglob))
- S1 = 0e0_CUSTOM_REAL
- S2 = 0e0_CUSTOM_REAL
- S3 = 0e0_CUSTOM_REAL
- n1=0
- n2=0
- n3=0
- read(IIN_PAR, nml=INIT_STRESS)
- bc%T0(1,:) = S1
- bc%T0(2,:) = S2
- bc%T0(3,:) = S3
+ ! Set initial fault stresses
+ allocate(bc%T0(3,bc%nglob))
+ S1 = 0e0_CUSTOM_REAL
+ S2 = 0e0_CUSTOM_REAL
+ S3 = 0e0_CUSTOM_REAL
+ n1=0
+ n2=0
+ n3=0
+ read(IIN_PAR, nml=INIT_STRESS)
+ bc%T0(1,:) = S1
+ bc%T0(2,:) = S2
+ bc%T0(3,:) = S3
- call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
- call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
- call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
+ call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
+ call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
+ call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
- !WARNING : Quick and dirty free surface condition at z=0
- ! do k=1,bc%nglob
- ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
- ! end do
+ !WARNING : Quick and dirty free surface condition at z=0
+ ! do k=1,bc%nglob
+ ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
+ ! end do
- if (TPV16) then
+ if (TPV16) then
- allocate(inp_nx(bc%nglob))
- allocate(inp_nz(bc%nglob))
- allocate(loc_str(bc%nglob))
- allocate(loc_dip(bc%nglob))
- allocate(sigma0(bc%nglob))
- allocate(tau0_str(bc%nglob))
- allocate(tau0_dip(bc%nglob))
- allocate(Rstress_str(bc%nglob))
- allocate(Rstress_dip(bc%nglob))
- allocate(static_fc(bc%nglob))
- allocate(dyn_fc(bc%nglob))
- allocate(swcd(bc%nglob))
- allocate(cohes(bc%nglob))
- allocate(tim_forcedRup(bc%nglob))
+ allocate(inp_nx(bc%nglob))
+ allocate(inp_nz(bc%nglob))
+ allocate(loc_str(bc%nglob))
+ allocate(loc_dip(bc%nglob))
+ allocate(sigma0(bc%nglob))
+ allocate(tau0_str(bc%nglob))
+ allocate(tau0_dip(bc%nglob))
+ allocate(Rstress_str(bc%nglob))
+ allocate(Rstress_dip(bc%nglob))
+ allocate(static_fc(bc%nglob))
+ allocate(dyn_fc(bc%nglob))
+ allocate(swcd(bc%nglob))
+ allocate(cohes(bc%nglob))
+ allocate(tim_forcedRup(bc%nglob))
- open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
- read(IIN_NUC,*) relz_num,sub_relz_num
- read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
- read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
- do ipar=1,bc%nglob
- read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
- Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
- enddo
- close(IIN_NUC)
+ open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
+ read(IIN_NUC,*) relz_num,sub_relz_num
+ read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
+ read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ do ipar=1,bc%nglob
+ read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
+ Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
+ enddo
+ close(IIN_NUC)
- allocate( bc%swf )
- allocate( bc%swf%mus(bc%nglob) )
- allocate( bc%swf%mud(bc%nglob) )
- allocate( bc%swf%Dc(bc%nglob) )
- allocate( bc%swf%theta(bc%nglob) )
- allocate( bc%swf%T(bc%nglob) )
- allocate( bc%swf%C(bc%nglob) )
+ allocate( bc%swf )
+ allocate( bc%swf%mus(bc%nglob) )
+ allocate( bc%swf%mud(bc%nglob) )
+ allocate( bc%swf%Dc(bc%nglob) )
+ allocate( bc%swf%theta(bc%nglob) )
+ allocate( bc%swf%T(bc%nglob) )
+ allocate( bc%swf%C(bc%nglob) )
- minX = minval(bc%coord(1,:))
+ minX = minval(bc%coord(1,:))
- do iLoc=1,bc%nglob
+ do iLoc=1,bc%nglob
- ipar = minloc( (minX+loc_str(:)-bc%coord(1,iLoc))**2 + (-loc_dip(:)-bc%coord(3,iLoc))**2 , 1) !iloc_dip is negative of Z-coord
- bc%T0(3,iLoc) = -sigma0(ipar)
- bc%T0(1,iLoc) = tau0_str(ipar)
- bc%T0(2,iLoc) = tau0_dip(ipar)
-
- bc%swf%mus(iLoc) = static_fc(ipar)
- bc%swf%mud(iLoc) = dyn_fc(ipar)
- bc%swf%Dc(iLoc) = swcd(ipar)
- bc%swf%C(iLoc) = cohes(ipar)
- bc%swf%T(iLoc) = tim_forcedRup(ipar)
- enddo
+ ipar = minloc( (minX+loc_str(:)-bc%coord(1,iLoc))**2 + (-loc_dip(:)-bc%coord(3,iLoc))**2 , 1) !iloc_dip is negative of Z-coord
+ bc%T0(3,iLoc) = -sigma0(ipar)
+ bc%T0(1,iLoc) = tau0_str(ipar)
+ bc%T0(2,iLoc) = tau0_dip(ipar)
+
+ bc%swf%mus(iLoc) = static_fc(ipar)
+ bc%swf%mud(iLoc) = dyn_fc(ipar)
+ bc%swf%Dc(iLoc) = swcd(ipar)
+ bc%swf%C(iLoc) = cohes(ipar)
+ bc%swf%T(iLoc) = tim_forcedRup(ipar)
+ enddo
- endif
+ endif
- ! Set friction parameters and initialize friction variables
- ! Slip weakening friction
- if(.not. Rate_AND_State) then
- allocate( bc%swf )
- allocate( bc%swf%mus(bc%nglob) )
- allocate( bc%swf%mud(bc%nglob) )
- allocate( bc%swf%Dc(bc%nglob) )
- allocate( bc%swf%theta(bc%nglob) )
- allocate( bc%swf%C(bc%nglob) )
- allocate( bc%swf%T(bc%nglob) )
- ! WARNING: if V_HEALING is negative we turn off healing
- bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+ ! Set friction parameters and initialize friction variables
+ ! Slip weakening friction
+ if(.not. Rate_AND_State) then
+ allocate( bc%swf )
+ allocate( bc%swf%mus(bc%nglob) )
+ allocate( bc%swf%mud(bc%nglob) )
+ allocate( bc%swf%Dc(bc%nglob) )
+ allocate( bc%swf%theta(bc%nglob) )
+ allocate( bc%swf%C(bc%nglob) )
+ allocate( bc%swf%T(bc%nglob) )
+ ! WARNING: if V_HEALING is negative we turn off healing
+ bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
- mus = 0.6e0_CUSTOM_REAL
- mud = 0.1e0_CUSTOM_REAL
- dc = 1e0_CUSTOM_REAL
- nmus = 0
- nmud = 0
- ndc = 0
+ mus = 0.6e0_CUSTOM_REAL
+ mud = 0.1e0_CUSTOM_REAL
+ dc = 1e0_CUSTOM_REAL
+ nmus = 0
+ nmud = 0
+ ndc = 0
- C = 0._CUSTOM_REAL
- T = HUGEVAL
- nC = 0
- nForcedRup = 0
+ C = 0._CUSTOM_REAL
+ T = HUGEVAL
+ nC = 0
+ nForcedRup = 0
- read(IIN_PAR, nml=SWF)
- bc%swf%mus = mus
- bc%swf%mud = mud
- bc%swf%Dc = dc
- bc%swf%C = C
- bc%swf%T = T
- call init_2d_distribution(bc%swf%mus,bc%coord,IIN_PAR,nmus)
- call init_2d_distribution(bc%swf%mud,bc%coord,IIN_PAR,nmud)
- call init_2d_distribution(bc%swf%Dc ,bc%coord,IIN_PAR,ndc)
- call init_2d_distribution(bc%swf%C ,bc%coord,IIN_PAR,nC)
- call init_2d_distribution(bc%swf%T ,bc%coord,IIN_PAR,nForcedRup)
+ read(IIN_PAR, nml=SWF)
+ bc%swf%mus = mus
+ bc%swf%mud = mud
+ bc%swf%Dc = dc
+ bc%swf%C = C
+ bc%swf%T = T
+ call init_2d_distribution(bc%swf%mus,bc%coord,IIN_PAR,nmus)
+ call init_2d_distribution(bc%swf%mud,bc%coord,IIN_PAR,nmud)
+ call init_2d_distribution(bc%swf%Dc ,bc%coord,IIN_PAR,ndc)
+ call init_2d_distribution(bc%swf%C ,bc%coord,IIN_PAR,nC)
+ call init_2d_distribution(bc%swf%T ,bc%coord,IIN_PAR,nForcedRup)
- bc%swf%theta = 0e0_CUSTOM_REAL
- allocate(bc%MU(bc%nglob))
- bc%MU = swf_mu(bc%swf)
+ bc%swf%theta = 0e0_CUSTOM_REAL
+ allocate(bc%MU(bc%nglob))
+ bc%MU = swf_mu(bc%swf)
- ! Rate and state friction
- else
- allocate( bc%rsf )
- allocate( bc%rsf%V0(bc%nglob) )
- allocate( bc%rsf%f0(bc%nglob) )
- allocate( bc%rsf%a(bc%nglob) )
- allocate( bc%rsf%b(bc%nglob) )
- allocate( bc%rsf%L(bc%nglob) )
- allocate( bc%rsf%V_init(bc%nglob) )
- allocate( bc%rsf%theta(bc%nglob) )
- allocate( bc%rsf%C(bc%nglob) )
- allocate( bc%rsf%T(bc%nglob) )
+ ! Rate and state friction
+ else
+ allocate( bc%rsf )
+ allocate( bc%rsf%V0(bc%nglob) )
+ allocate( bc%rsf%f0(bc%nglob) )
+ allocate( bc%rsf%a(bc%nglob) )
+ allocate( bc%rsf%b(bc%nglob) )
+ allocate( bc%rsf%L(bc%nglob) )
+ allocate( bc%rsf%V_init(bc%nglob) )
+ allocate( bc%rsf%theta(bc%nglob) )
+ allocate( bc%rsf%C(bc%nglob) )
+ allocate( bc%rsf%T(bc%nglob) )
- V0 =1.e-6_CUSTOM_REAL
- f0 =0.6_CUSTOM_REAL
- a =0.0080_CUSTOM_REAL !0.0080_CUSTOM_REAL
- b =0.0040_CUSTOM_REAL !0.0120_CUSTOM_REAL
- L =0.0135_CUSTOM_REAL
- V_init =1.e-12_CUSTOM_REAL
- theta_init =1.084207680000000e+09_CUSTOM_REAL
+ V0 =1.e-6_CUSTOM_REAL
+ f0 =0.6_CUSTOM_REAL
+ a =0.0080_CUSTOM_REAL !0.0080_CUSTOM_REAL
+ b =0.0040_CUSTOM_REAL !0.0120_CUSTOM_REAL
+ L =0.0135_CUSTOM_REAL
+ V_init =1.e-12_CUSTOM_REAL
+ theta_init =1.084207680000000e+09_CUSTOM_REAL
- nV0 =0
- nf0 =0
- na =0
- nb =0
- nL =0
- nV_init =0
- ntheta_init =0
+ nV0 =0
+ nf0 =0
+ na =0
+ nb =0
+ nL =0
+ nV_init =0
+ ntheta_init =0
- C = 0._CUSTOM_REAL
- T = HUGEVAL
- nC = 0
- nForcedRup = 0
+ C = 0._CUSTOM_REAL
+ T = HUGEVAL
+ nC = 0
+ nForcedRup = 0
- read(IIN_PAR, nml=RSF)
- bc%rsf%V0 = V0
- bc%rsf%f0 = f0
- bc%rsf%a = a
- bc%rsf%b = b
- bc%rsf%L = L
- bc%rsf%V_init = V_init
- bc%rsf%theta = theta_init
- bc%rsf%C = C
- bc%rsf%T = T
- call init_2d_distribution(bc%rsf%V0,bc%coord,IIN_PAR,nV0)
- call init_2d_distribution(bc%rsf%f0,bc%coord,IIN_PAR,nf0)
- call init_2d_distribution(bc%rsf%a,bc%coord,IIN_PAR,na)
- call init_2d_distribution(bc%rsf%b,bc%coord,IIN_PAR,nb)
- call init_2d_distribution(bc%rsf%L,bc%coord,IIN_PAR,nL)
- call init_2d_distribution(bc%rsf%V_init,bc%coord,IIN_PAR,nV_init)
- call init_2d_distribution(bc%rsf%theta,bc%coord,IIN_PAR,ntheta_init)
- call init_2d_distribution(bc%rsf%C ,bc%coord,IIN_PAR,nC)
- call init_2d_distribution(bc%rsf%T ,bc%coord,IIN_PAR,nForcedRup)
+ read(IIN_PAR, nml=RSF)
+ bc%rsf%V0 = V0
+ bc%rsf%f0 = f0
+ bc%rsf%a = a
+ bc%rsf%b = b
+ bc%rsf%L = L
+ bc%rsf%V_init = V_init
+ bc%rsf%theta = theta_init
+ bc%rsf%C = C
+ bc%rsf%T = T
+ call init_2d_distribution(bc%rsf%V0,bc%coord,IIN_PAR,nV0)
+ call init_2d_distribution(bc%rsf%f0,bc%coord,IIN_PAR,nf0)
+ call init_2d_distribution(bc%rsf%a,bc%coord,IIN_PAR,na)
+ call init_2d_distribution(bc%rsf%b,bc%coord,IIN_PAR,nb)
+ call init_2d_distribution(bc%rsf%L,bc%coord,IIN_PAR,nL)
+ call init_2d_distribution(bc%rsf%V_init,bc%coord,IIN_PAR,nV_init)
+ call init_2d_distribution(bc%rsf%theta,bc%coord,IIN_PAR,ntheta_init)
+ call init_2d_distribution(bc%rsf%C ,bc%coord,IIN_PAR,nC)
+ call init_2d_distribution(bc%rsf%T ,bc%coord,IIN_PAR,nForcedRup)
- W1=15000._CUSTOM_REAL
- W2=7500._CUSTOM_REAL
- w=3000._CUSTOM_REAL
- hypo_z = -7500._CUSTOM_REAL
- do iglob=1,bc%nglob
- x=bc%coord(1,iglob)
- z=bc%coord(3,iglob)
- c1=abs(x)<W1+w
- c2=abs(x)>W1
- c3=abs(z-hypo_z)<W2+w
- c4=abs(z-hypo_z)>W2
- if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
+ W1=15000._CUSTOM_REAL
+ W2=7500._CUSTOM_REAL
+ w=3000._CUSTOM_REAL
+ hypo_z = -7500._CUSTOM_REAL
+ do iglob=1,bc%nglob
+ x=bc%coord(1,iglob)
+ z=bc%coord(3,iglob)
+ c1=abs(x)<W1+w
+ c2=abs(x)>W1
+ c3=abs(z-hypo_z)<W2+w
+ c4=abs(z-hypo_z)>W2
+ if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
- if (c1 .and. c2) then
- b11 = w/(abs(x)-W1-w)
- b12 = w/(abs(x)-W1)
- B1 = 0.5 * (ONE + tanh(b11 + b12))
- elseif(abs(x)<=W1) then
- B1 = 1._CUSTOM_REAL
- else
- B1 = 0._CUSTOM_REAL
- endif
+ if (c1 .and. c2) then
+ b11 = w/(abs(x)-W1-w)
+ b12 = w/(abs(x)-W1)
+ B1 = 0.5 * (ONE + tanh(b11 + b12))
+ elseif(abs(x)<=W1) then
+ B1 = 1._CUSTOM_REAL
+ else
+ B1 = 0._CUSTOM_REAL
+ endif
- if (c3 .and. c4) then
- b21 = w/(abs(z-hypo_z)-W2-w)
- b22 = w/(abs(z-hypo_z)-W2)
- B2 = 0.5 * (ONE + tanh(b21 + b22))
- elseif(abs(z-hypo_z)<=W2) then
- B2 = 1._CUSTOM_REAL
+ if (c3 .and. c4) then
+ b21 = w/(abs(z-hypo_z)-W2-w)
+ b22 = w/(abs(z-hypo_z)-W2)
+ B2 = 0.5 * (ONE + tanh(b21 + b22))
+ elseif(abs(z-hypo_z)<=W2) then
+ B2 = 1._CUSTOM_REAL
+ else
+ B2 = 0._CUSTOM_REAL
+ endif
+
+
+ bc%rsf%a(iglob) = 0.008 + 0.008 * (ONE - B1*B2)
+ elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
+ bc%rsf%a(iglob) = 0.008
else
- B2 = 0._CUSTOM_REAL
+ bc%rsf%a(iglob) = 0.016
endif
-
- bc%rsf%a(iglob) = 0.008 + 0.008 * (ONE - B1*B2)
- elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
- bc%rsf%a(iglob) = 0.008
- else
- bc%rsf%a(iglob) = 0.016
- endif
-
- bc%rsf%theta(iglob) = L/V0 * exp( (bc%rsf%a(iglob) * log(2.0*sinh(-S1/S3/bc%rsf%a(iglob))) - f0 - bc%rsf%a(iglob)*log(V_init/V0) )/b )
+ bc%rsf%theta(iglob) = L/V0 * exp( (bc%rsf%a(iglob) * log(2.0*sinh(-S1/S3/bc%rsf%a(iglob))) - f0 - bc%rsf%a(iglob)*log(V_init/V0) )/b )
- enddo
+ enddo
- allocate(bc%MU(bc%nglob))
+ allocate(bc%MU(bc%nglob))
- bc%V(1,:) = bc%rsf%V_init
- allocate( bc%asp )
- allocate( bc%asp%Fload(bc%nglob) )
+ bc%V(1,:) = bc%rsf%V_init
+ allocate( bc%asp )
+ allocate( bc%asp%Fload(bc%nglob) )
- Fload =0.e0_CUSTOM_REAL
- nFload =0
+ Fload =0.e0_CUSTOM_REAL
+ nFload =0
- read(IIN_PAR, nml=ASP)
- bc%asp%Fload =Fload
- call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
+ read(IIN_PAR, nml=ASP)
+ bc%asp%Fload =Fload
+ call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
+ endif
endif
@@ -1428,10 +1434,12 @@
npoin = DataXZ%npoin
allocate(DataXZ%stg(nglob))
- if(.not. Rate_AND_State) then
- DataXZ%sta => bc%swf%theta
- else
- DataXZ%sta => bc%rsf%theta
+ if(bc%nglob > 0) then
+ if(.not. Rate_AND_State) then
+ DataXZ%sta => bc%swf%theta
+ else
+ DataXZ%sta => bc%rsf%theta
+ endif
endif
DataXZ%d1 => bc%d(1,:)
DataXZ%d2 => bc%d(2,:)
More information about the CIG-COMMITS
mailing list