[cig-commits] r13307 - seismo/3D/SPECFEM3D_SESAME/trunk

[nlegoff at geodynamics.org]

Author: nlegoff
Date: 2008-11-14 05:31:42 -0800 (Fri, 14 Nov 2008)
New Revision: 13307

Modified:
   seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90
Log:
added regolith layer for asteroid (in comments).

Modified: seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90	2008-11-14 13:21:24 UTC (rev 13306)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90	2008-11-14 13:31:42 UTC (rev 13307)
@@ -485,6 +485,11 @@
   logical, dimension(:), allocatable :: iglob_is_inner_ext_mesh
   integer :: iinterface
 
+!!!! NL NL REGOLITH : regolith layer for asteroid
+!!$  double precision, external :: materials_ext_mesh
+!!$  logical, dimension(:), allocatable :: ispec_is_regolith  
+!!!! NL NL REGOLITH
+
 ! ************** PROGRAM STARTS HERE **************
 
 ! sizeprocs returns number of processes started
@@ -1234,6 +1239,80 @@
 
   endif !
 
+!!!! NL NL REGOLITH : runs at cines for asteroid simulations. Elements in contact with surface are part of the regolith layer.
+!!$  allocate(ispec_is_regolith(NSPEC_AB))
+!!$  ispec_is_regolith(:) = .false.
+!!$  do ispec = 1, NSPEC_AB
+!!$    do k = 1, NGLLZ
+!!$      do j = 1, NGLLY
+!!$        do i = 1, NGLLX
+!!$          iglob = ibool(i,j,k,ispec)
+!!$          if (iglob_is_surface_external_mesh(iglob)) then
+!!$            ispec_is_regolith(ispec) = .true.
+!!$          endif
+!!$        enddo
+!!$      enddo
+!!$    enddo
+!!$  enddo
+!!$
+!!$  do ispec = 1, NSPEC_AB
+!!$    if (ispec_is_regolith(ispec)) then
+!!$      do k = 1, NGLLZ
+!!$        do j = 1, NGLLY
+!!$          do i = 1, NGLLX
+!!$             kappastore(i,j,k,ispec) = materials_ext_mesh(1,2)* &
+!!$                  (materials_ext_mesh(2,2)*materials_ext_mesh(2,2) - &
+!!$                  4.d0*materials_ext_mesh(3,2)*materials_ext_mesh(3,2)/3.d0)
+!!$             mustore(i,j,k,ispec) = materials_ext_mesh(1,2)*materials_ext_mesh(3,2)*&
+!!$                  materials_ext_mesh(3,2)
+!!$             
+!!$          enddo
+!!$        enddo
+!!$      enddo
+!!$    endif
+!!$  enddo
+!!$
+!!$
+!!$  call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+!!$  call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+!!$  call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+!!$
+!!$  rmass(:) = 0._CUSTOM_REAL
+!!$
+!!$  do ispec=1,NSPEC_AB
+!!$  do k=1,NGLLZ
+!!$    do j=1,NGLLY
+!!$      do i=1,NGLLX
+!!$        weight=wxgll(i)*wygll(j)*wzgll(k)
+!!$        iglob=ibool(i,j,k,ispec)
+!!$
+!!$        jacobianl=jacobian(i,j,k,ispec)
+!!$
+!!$! distinguish between single and double precision for reals
+!!$        if (.not. ispec_is_regolith(ispec)) then
+!!$        if(CUSTOM_REAL == SIZE_REAL) then
+!!$          rmass(iglob) = rmass(iglob) + &
+!!$               sngl(dble(materials_ext_mesh(1,1)) * dble(jacobianl) * weight)
+!!$        else
+!!$          rmass(iglob) = rmass(iglob) + materials_ext_mesh(1,1) * jacobianl * weight
+!!$        endif
+!!$        else
+!!$        if(CUSTOM_REAL == SIZE_REAL) then
+!!$          rmass(iglob) = rmass(iglob) + &
+!!$               sngl(dble(materials_ext_mesh(1,2)) * dble(jacobianl) * weight)
+!!$        else
+!!$          rmass(iglob) = rmass(iglob) + materials_ext_mesh(1,2) * jacobianl * weight
+!!$        endif
+!!$        endif
+!!$
+!!$      enddo
+!!$    enddo
+!!$  enddo
+!!$  enddo
+
+  
+!!!! NL NL REGOLITH
+
   endif
 
 ! $$$$$$$$$$$$$$$$$$$$$$$$ SOURCES $$$$$$$$$$$$$$$$$
@@ -4041,3 +4120,40 @@
 
   end subroutine specfem3D
 
+
+!!!! NL NL REGOLITH
+!!$  double precision function materials_ext_mesh(i,j)
+!!$
+!!$    implicit none
+!!$
+!!$    integer :: i,j
+!!$
+!!$    select case (j)
+!!$      case (1)
+!!$        select case (i)
+!!$          case (1)
+!!$            materials_ext_mesh = 2700.d0
+!!$          case (2)
+!!$            materials_ext_mesh = 3000.d0
+!!$          case (3)
+!!$            materials_ext_mesh = 1732.051d0
+!!$          case default
+!!$            call stop_all()
+!!$          end select
+!!$      case (2)
+!!$        select case (i)
+!!$          case (1)
+!!$            materials_ext_mesh = 2000.d0
+!!$          case (2)
+!!$            materials_ext_mesh = 900.d0
+!!$          case (3)
+!!$            materials_ext_mesh = 500.d0
+!!$          case default
+!!$            call stop_all()
+!!$          end select
+!!$      case default
+!!$        call stop_all()
+!!$    end select
+!!$
+!!$  end function materials_ext_mesh
+!!!! NL NL REGOLITH