updates mesher array initialization (and comments out unneeded parameters)

Author: danielpeter

src/generate_databases/model_vs30.f90

@@ -74,7 +74,10 @@ subroutine model_vs30_setup_check()
   ! model values
   real(kind=CUSTOM_REAL) :: vs,min_vs,max_vs
 
-  ! checks is USGS Vs30 model interface is provided in folder USGS_VS30/
+  ! initializes
+  USE_MODEL_LAYER_VS30 = .false.
+
+  ! checks if USGS Vs30 model interface is provided in folder USGS_VS30/
   !
   ! the script ./utils/scripts/run_get_USGS_Vs30.py can create an interface file for a given region
   ! and stores it in folder USGS_VS30/ as 'interface_vs30.dat'
@@ -83,7 +86,7 @@ subroutine model_vs30_setup_check()
 
   inquire(file=trim(filename),exist=file_exists)
 
-  ! checks is anything to do
+  ! checks if anything to do
   if (.not. file_exists) return
 
   ! user info

src/meshfem3D/compute_parameters.f90

@@ -30,7 +30,7 @@ subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
                                 NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
                                 NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
                                 NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
-                                NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB, &
+                                NGLOB_AB, &
                                 USE_REGULAR_MESH,NDOUBLINGS,ner_doublings)
 
   use constants
@@ -53,7 +53,6 @@ subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
   integer,intent(out) :: NSPEC2D_A_XI,NSPEC2D_B_XI,NSPEC2D_A_ETA,NSPEC2D_B_ETA
 
   integer,intent(out) :: NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP
-  integer,intent(out) :: NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX
 
   ! local parameters
   integer :: NGLOB_DOUBLING,NGLOB_DOUBLING_SINGLE
@@ -102,8 +101,9 @@ subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
   NSPEC2DMAX_YMIN_YMAX = 0
   NSPEC2D_BOTTOM = 0
   NSPEC2D_TOP = 0
-  NPOIN2DMAX_XMIN_XMAX = 0
-  NPOIN2DMAX_YMIN_YMAX = 0
+
+  !NPOIN2DMAX_XMIN_XMAX = 0 - not needed..
+  !NPOIN2DMAX_YMIN_YMAX = 0
 
   !
   !--- case of a regular mesh
@@ -152,8 +152,8 @@ subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
 
     ! 2-D addressing and buffers for summation between slices
     ! we add one to number of points because of the flag after the last point
-    NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*NGLLY_M*NGLLZ_M + 1
-    NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*NGLLX_M*NGLLZ_M + 1
+    !NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*NGLLY_M*NGLLZ_M + 1   - not needed..
+    !NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*NGLLX_M*NGLLZ_M + 1
 
     ! exact number of global points per slice
     NGLOB_AB = (NEX_PER_PROC_XI*(NGLLX_M-1)+1) * (NEX_PER_PROC_ETA*(NGLLY_M-1)+1) * (NER*(NGLLZ_M-1)+1)
@@ -402,8 +402,8 @@ subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
 
     ! 2-D addressing and buffers for summation between slices
     ! we add one to number of points because of the flag after the last point
-    NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX * NGLLY_M * NGLLZ_M + 1
-    NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX * NGLLX_M * NGLLZ_M + 1
+    !NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX * NGLLY_M * NGLLZ_M + 1  - not needed..
+    !NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX * NGLLX_M * NGLLZ_M + 1
 
     ! exact number of global points
     NGLOB_AB = NGLOB_DOUBLING + NGLOB_NO_DOUBLING
@@ -610,7 +610,6 @@ subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
 !      NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*NGLLY_M*NGLLZ_M + 1
 !      NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*NGLLX_M*NGLLZ_M + 1
 !
-!
 !      ! exact number of global points
 !
 !      ! case of the doubling regions

src/meshfem3D/create_meshfem_mesh.f90

@@ -45,8 +45,8 @@ module create_meshfem_par
 
   ! flag indicating whether point is in the sediments
   !  logical point_is_in_sediments
-  logical, dimension(:,:,:,:), allocatable :: flag_sediments
-  logical, dimension(:), allocatable :: not_fully_in_bedrock
+  !logical, dimension(:,:,:,:), allocatable :: flag_sediments  - not used yet..
+  !logical, dimension(:), allocatable :: not_fully_in_bedrock
 
   ! boundary parameters locator
   integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top
@@ -146,7 +146,7 @@ subroutine create_meshfem_mesh()
                            nspec,nglob, &
                            prname,nodes_coords,ibool,ispec_material_id)
 
-  !! setting up wavefield discontinuity interface
+  ! setting up wavefield discontinuity interface
   if (IS_WAVEFIELD_DISCONTINUITY) then
     call find_wavefield_discontinuity_elements()
   endif
@@ -248,23 +248,27 @@ subroutine cmm_allocate_arrays()
   ! use dynamic allocation to allocate memory for arrays
   allocate(ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1257')
+  ibool(:,:,:,:) = 0
 
   allocate(xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1258')
   allocate(ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1259')
   allocate(zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1260')
+  xstore(:,:,:,:) = 0.d0; ystore(:,:,:,:) = 0.d0; zstore(:,:,:,:) = 0.d0
 
   ! flag indicating whether point is in the sediments
-  allocate(flag_sediments(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier)
-  if (ier /= 0) call exit_MPI(myrank,'error allocating array 1261')
-  allocate(not_fully_in_bedrock(nspec),stat=ier)
-  if (ier /= 0) call exit_MPI(myrank,'error allocating array 1262')
+  !allocate(flag_sediments(NGLLX_M,NGLLY_M,NGLLZ_M,nspec),stat=ier)
+  !if (ier /= 0) call exit_MPI(myrank,'error allocating array 1261')
+  !allocate(not_fully_in_bedrock(nspec),stat=ier)
+  !if (ier /= 0) call exit_MPI(myrank,'error allocating array 1262')
+  !flag_sediments(:,:,:,:) = .false.; not_fully_in_bedrock(:) = .false.
 
   ! boundary locator
   allocate(iboun(6,nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1263')
+  iboun(:,:) = .false.
 
   ! boundary parameters locator
   allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),stat=ier)
@@ -279,15 +283,19 @@ subroutine cmm_allocate_arrays()
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1268')
   allocate(ibelm_top(NSPEC2D_TOP),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1269')
+  ibelm_xmin(:) = 0; ibelm_xmax(:) = 0; ibelm_ymin(:) = 0; ibelm_ymax(:) = 0
+  ibelm_bottom(:) = 0; ibelm_top(:) = 0
 
   ! MPI cut-planes parameters along xi and along eta
   allocate(iMPIcut_xi(2,nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1270')
   allocate(iMPIcut_eta(2,nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1271')
+  iMPIcut_xi(:,:) = .false.; iMPIcut_eta(:,:) = .false.
 
   allocate(ispec_material_id(nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'error allocating array 1279')
+  ispec_material_id(:) = 0
 
   ! synchronize
   call synchronize_all()

src/meshfem3D/meshfem3D.F90

@@ -193,7 +193,7 @@ program xmeshfem3D
                           NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
                           NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
                           NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
-                          NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB, &
+                          NGLOB_AB, &
                           USE_REGULAR_MESH,NDOUBLINGS,ner_doublings)
 
   ! check that the code is running with the requested nb of processes

src/meshfem3D/meshfem3D_par.f90

@@ -109,8 +109,9 @@ module meshfem_par
   integer :: NSPEC2D_A_XI,NSPEC2D_B_XI, &
              NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
              NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
-             NSPEC2D_BOTTOM,NSPEC2D_TOP, &
-             NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX
+             NSPEC2D_BOTTOM,NSPEC2D_TOP
+
+  !integer :: NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX - not needed...
 
   ! interfaces parameters
   integer :: number_of_interfaces,number_of_layers

src/meshfem3D/store_boundaries.f90

@@ -59,49 +59,62 @@ subroutine store_boundaries(iboun,nspec, &
   ispecb5 = 0
   ispecb6 = 0
 
+  ! boundary element arrays
+  ibelm_xmin(:) = 0
+  ibelm_xmax(:) = 0
+  ibelm_ymin(:) = 0
+  ibelm_ymax(:) = 0
+  ibelm_bottom(:) = 0
+  ibelm_top(:) = 0
+
+  nspec2D_xmin = 0
+  nspec2D_xmax = 0
+  nspec2D_ymin = 0
+  nspec2D_ymax = 0
+
   ! determine if the element falls on a boundary
   do ispec = 1,nspec
 
     ! on boundary: xmin
     if (iboun(1,ispec)) then
-      ispecb1 = ispecb1+1
+      ispecb1 = ispecb1 + 1
       if (ispecb1 > NSPEC2DMAX_XMIN_XMAX) stop 'error NSPEC2DMAX_XMIN_XMAX too small'
-      ibelm_xmin(ispecb1)=ispec
+      ibelm_xmin(ispecb1) = ispec
     endif
 
     ! on boundary: xmax
     if (iboun(2,ispec)) then
-      ispecb2 = ispecb2+1
+      ispecb2 = ispecb2 + 1
       if (ispecb2 > NSPEC2DMAX_XMIN_XMAX) stop 'error NSPEC2DMAX_XMIN_XMAX too small'
-      ibelm_xmax(ispecb2)=ispec
+      ibelm_xmax(ispecb2) = ispec
     endif
 
     ! on boundary: ymin
     if (iboun(3,ispec)) then
-      ispecb3 = ispecb3+1
+      ispecb3 = ispecb3 + 1
       if (ispecb3 > NSPEC2DMAX_YMIN_YMAX) stop 'error NSPEC2DMAX_YMIN_YMAX too small'
-      ibelm_ymin(ispecb3)=ispec
+      ibelm_ymin(ispecb3) = ispec
     endif
 
     ! on boundary: ymax
     if (iboun(4,ispec)) then
-      ispecb4 = ispecb4+1
+      ispecb4 = ispecb4 + 1
       if (ispecb4 > NSPEC2DMAX_YMIN_YMAX) stop 'error NSPEC2DMAX_YMIN_YMAX too small'
-      ibelm_ymax(ispecb4)=ispec
+      ibelm_ymax(ispecb4) = ispec
     endif
 
     ! on boundary: bottom
     if (iboun(5,ispec)) then
-      ispecb5 = ispecb5+1
+      ispecb5 = ispecb5 + 1
       if (ispecb5 > NSPEC2D_BOTTOM) stop 'error NSPEC2D_BOTTOM too small'
-      ibelm_bottom(ispecb5)=ispec
+      ibelm_bottom(ispecb5) = ispec
     endif
 
     ! on boundary: top
     if (iboun(6,ispec)) then
-      ispecb6 = ispecb6+1
+      ispecb6 = ispecb6 + 1
       if (ispecb6 > NSPEC2D_TOP) stop 'error NSPEC2D_TOP too small'
-      ibelm_top(ispecb6)=ispec
+      ibelm_top(ispecb6) = ispec
     endif
 
   enddo