Merge pull request #849 from danielpeter/devel

adds initial SEMUCB Berkeley model; fixes attenuation for 1d models with 1d storage arrays

EXAMPLES/regional_Berkeley/DATA/CMTSOLUTION

@@ -0,0 +1,13 @@
+PDEW 2008  1  6  5 14 23.70  37.2200   22.6900  75.0 6.1 6.2 SOUTHERN GREECE
+event name:     200801060514A
+time shift:      0.000
+half duration:   3.0000
+latitude:       36.9800
+longitude:      22.8700
+depth:          92.3900
+Mrr:       7.740000e+24
+Mtt:      -1.830000e+25
+Mpp:       1.060000e+25
+Mrt:       1.290000e+25
+Mrp:      -8.450000e+24
+Mtp:      -4.430000e+24

EXAMPLES/regional_Berkeley/DATA/STATIONS

@@ -0,0 +1,6 @@
+GRFO       IU       49.6909     11.2203     384.0   116.0
+BEKI       YL       41.31500    34.26300   1415.0     0.0
+BGIO       SR       31.72200    35.08780    752.0   100.0
+LIT        HT       40.10080    22.49000    480.0     0.0
+ZKR        GE       35.11470    26.21700    270.0     0.0
+S001       XS       37.28300    21.71800    156.0     0.0

EXAMPLES/regional_Berkeley/README

@@ -0,0 +1,41 @@
+----------------------------------------------------------
+README
+----------------------------------------------------------
+
+This example creates a regional mesh using transversely isotropic PREM
+and runs a spectral-element simulation.
+
+The simulation is setup to use 4 CPUs only, thus can be run locally on a desktop machine
+(needs to have MPI installed).
+
+
+1. forward simulation:
+----------------------
+
+type:
+./run_this_example.sh
+
+to setup this directory for the simulation and to run a corresponding
+script "run_mesher_solver.bash".
+
+The simulation is accurate down to a shortest period of ~ 15 s.
+(minimum period resolved ~ (256/NEX_XI) * (ANGULAR_WIDTH_XI_IN_DEGREES/90) * 17)
+
+Reference seismograms for comparisons are provided in the directory: REF_SEIS/
+
+
+2. kernel simulation:
+---------------------
+
+type:
+./run_this_example.kernel.sh
+
+to setup this directory for the simulation and to submit a corresponding
+script "run_mesher_solver.kernel.bash".
+
+The simulation will use the adjoint source given in directory: SEM/
+to compute a travel-time sensitivity kernel.
+
+
+Please modify these scripts according to your specifics.
+

EXAMPLES/regional_Berkeley/run_mesher_solver.bash

@@ -0,0 +1,69 @@
+#!/bin/bash
+
+# gets settings from Par_file
+
+BASEMPIDIR=`grep ^LOCAL_PATH DATA/Par_file | cut -d = -f 2 `
+
+# script to run the mesher and the solver
+# read DATA/Par_file to get information about the run
+# compute total number of nodes needed
+NPROC_XI=`grep ^NPROC_XI DATA/Par_file | cut -d = -f 2 `
+NPROC_ETA=`grep ^NPROC_ETA DATA/Par_file | cut -d = -f 2`
+NCHUNKS=`grep ^NCHUNKS DATA/Par_file | cut -d = -f 2 `
+
+# total number of nodes is the product of the values read
+numnodes=$(( $NCHUNKS * $NPROC_XI * $NPROC_ETA ))
+
+mkdir -p OUTPUT_FILES
+
+# backup files used for this simulation
+cp DATA/Par_file OUTPUT_FILES/
+cp DATA/STATIONS OUTPUT_FILES/
+cp DATA/CMTSOLUTION OUTPUT_FILES/
+
+
+##
+## mesh generation
+##
+sleep 2
+
+echo
+echo `date`
+echo "starting MPI mesher on $numnodes processors"
+echo
+
+mpirun -np $numnodes $PWD/bin/xmeshfem3D
+
+# checks exit code
+if [[ $? -ne 0 ]]; then exit 1; fi
+
+echo "  mesher done: `date`"
+echo
+
+# backup important files addressing.txt and list*.txt
+cp OUTPUT_FILES/*.txt $BASEMPIDIR/
+
+
+##
+## forward simulation
+##
+
+# set up addressing
+#cp $BASEMPIDIR/addr*.txt OUTPUT_FILES/
+#cp $BASEMPIDIR/list*.txt OUTPUT_FILES/
+
+sleep 2
+
+echo
+echo `date`
+echo starting run in current directory $PWD
+echo
+
+mpirun -np $numnodes $PWD/bin/xspecfem3D
+
+# checks exit code
+if [[ $? -ne 0 ]]; then exit 1; fi
+
+echo "finished successfully"
+echo `date`
+

EXAMPLES/regional_Berkeley/run_this_example.sh

@@ -0,0 +1,81 @@
+#!/bin/bash
+#
+# regional simulation example
+#
+# script runs mesher and solver using mpirun
+# on 4 CPUs
+#
+# synthetics have an approximate shortest period ~ 20 s
+#
+# modify accordingly for your own system specifics
+##################################################
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo "directory: $currentdir"
+echo "(will take about 15 minutes)"
+echo
+
+# sets up directory structure in current example directoy
+echo
+echo "   setting up example..."
+echo
+
+mkdir -p DATABASES_MPI
+mkdir -p OUTPUT_FILES
+
+rm -rf DATABASES_MPI/*
+rm -rf OUTPUT_FILES/*
+
+# checks if executables were compiled and available
+if [ ! -e ../../bin/xspecfem3D ]; then
+  echo "Compiling first all binaries in the root directory..."
+  echo
+
+  # compiles executables in root directory
+  # using default configuration
+  cd ../../
+
+  # only in case static compilation would have been set to yes in Makefile:
+  cp $currentdir/DATA/Par_file DATA/Par_file
+
+  # compiles code
+  make clean
+  make -j4 all
+
+  # checks exit code
+  if [[ $? -ne 0 ]]; then exit 1; fi
+
+  # backup of constants setup
+  cp setup/* $currentdir/OUTPUT_FILES/
+  if [ -e OUTPUT_FILES/values_from_mesher ]; then
+    cp OUTPUT_FILES/values_from_mesher.h $currentdir/OUTPUT_FILES/values_from_mesher.h.compilation
+  fi
+  cp DATA/Par_file $currentdir/OUTPUT_FILES/
+
+  cd $currentdir
+fi
+
+# copy executables
+mkdir -p bin
+rm -rf bin/*
+cp ../../bin/x* ./bin/
+
+# links data directories needed to run example in this current directory with s362ani
+cd DATA/
+ln -s ../../../DATA/SEMUCB_A3d
+ln -s ../../../DATA/topo_bathy
+cd ../
+
+# run mesher & solver
+echo
+echo "  running script..."
+echo
+./run_mesher_solver.bash
+
+# checks exit code
+if [[ $? -ne 0 ]]; then exit 1; fi
+
+echo `date`
+

setup/constants.h.in

@@ -48,13 +48,48 @@
 
 !*********************************************************************************************************
 
+!!-----------------------------------------------------------
+!!
+!! Gauss-Lobatto-Legendre resolution
+!!
+!!-----------------------------------------------------------
+! number of GLL points in each direction of an element (degree plus one)
+  integer, parameter :: NGLLX = 5
+  integer, parameter :: NGLLY = NGLLX
+  integer, parameter :: NGLLZ = NGLLX
+
+
+!!-----------------------------------------------------------
+!!
+!! Energy
+!!
+!!-----------------------------------------------------------
 !! added this temporarily here to make SPECFEM3D and SPECFEM3D_GLOBE much more similar
 !! in terms of the structure of their main time iteration loop; these are future features
 !! that are missing in this code but implemented in the other and that could thus be cut and pasted one day
   logical, parameter :: OUTPUT_ENERGY = .false.
   integer, parameter :: IOUT_ENERGY = 937  ! file number for the energy file
   logical, parameter :: GRAVITY_SIMULATION = .false.
 
+
+!!-----------------------------------------------------------
+!!
+!! new version compatibility
+!!
+!!-----------------------------------------------------------
+! for comparison against older versions
+! in version 7.0: tiso was constrained to below moho, new version allows also in crust
+! in version 8.0: conversions between geodetic & geocentric latitudes were always applied,
+!                 and additional differences in 1-chunk centering & source positioning
+  logical, parameter :: USE_OLD_VERSION_FORMAT = .false.
+
+
+!!-----------------------------------------------------------
+!!
+!! I/O
+!!
+!!-----------------------------------------------------------
+
 ! if files on a local path on each node are also seen as global with same path
 ! set to .true. typically on a machine with a common (shared) file system, e.g. LUSTRE, GPFS or NFS-mounted /home
 ! set to .false. typically on a cluster of nodes with local disks used to store the mesh (not very common these days)
@@ -204,6 +239,39 @@
   logical,parameter :: PLOT_PNM_IMAGE_TOPO_BATHY = .false.
 
 
+!!-----------------------------------------------------------
+!!
+!! for Berkeley model
+!!
+!!-----------------------------------------------------------
+
+! source-time function is a UCB-style filtered heaviside
+  logical, parameter :: STF_IS_UCB_HEAVISIDE = .false.
+
+! Save mirror files
+  logical, parameter :: SAVE_MIRRORS = .false.
+
+! Path to A3d model
+! Make sure this line ends on "/"
+  character (len=*), parameter :: A3d_folder = "DATA/SEMUCB_A3d/"
+
+! topography model
+!!--- ETOPO1 5 - Berkeley filtered
+! size of topography and bathymetry file
+  integer, parameter :: NX_BATHY_BERKELEY = 360,NY_BATHY_BERKELEY = 180
+! resolution of topography file in minutes
+  integer, parameter :: RESOLUTION_TOPO_FILE_BERKELEY = 2
+! pathname of the topography file (un-smoothed)
+  character (len=*), parameter :: PATHNAME_TOPO_FILE_BERKELEY = trim(A3d_folder)//"ETOPO5_1x1_filtre.dat"
+
+!! uncomment this to use Berkeley topography as default - and comment out corresponding parameters above !!
+!! Topography defaults to Berkeley
+!  integer, parameter :: EARTH_NX_BATHY = NX_BATHY_BERKELEY
+!  integer, parameter :: EARTH_NY_BATHY = NY_BATHY_BERKELEY
+!  double precision, parameter :: EARTH_RESOLUTION_TOPO_FILE = RESOLUTION_TOPO_FILE_BERKELEY
+!  character (len=*), parameter :: EARTH_PATHNAME_TOPO_FILE = PATHNAME_TOPO_FILE_BERKELEY
+
+
 !!-----------------------------------------------------------
 !!
 !! for crustal model
@@ -227,18 +295,6 @@
   double precision, parameter :: MINIMUM_SEDIMENT_THICKNESS = 2.d0 ! minimim thickness in km
 
 
-!!-----------------------------------------------------------
-!!
-!! Gauss-Lobatto-Legendre resolution
-!!
-!!-----------------------------------------------------------
-
-! number of GLL points in each direction of an element (degree plus one)
-  integer, parameter :: NGLLX = 5
-  integer, parameter :: NGLLY = NGLLX
-  integer, parameter :: NGLLZ = NGLLX
-
-
 !!-----------------------------------------------------------
 !!
 !! source/receiver setup
@@ -564,17 +620,6 @@
 ! be careful when changing this flag when computing classical kernel
   integer, parameter :: NTSTEP_BETWEEN_COMPUTE_KERNELS = 1
 
-!!-----------------------------------------------------------
-!!
-!! new version compatibility
-!!
-!!-----------------------------------------------------------
-! for comparison against older versions
-! in version 7.0: tiso was constrained to below moho, new version allows also in crust
-! in version 8.0: conversions between geodetic & geocentric latitudes were always applied,
-!                 and additional differences in 1-chunk centering & source positioning
-  logical, parameter :: USE_OLD_VERSION_FORMAT = .false.
-
 
 !!-----------------------------------------------------------
 !!
@@ -892,6 +937,7 @@
   integer, parameter :: REFERENCE_MODEL_JP1D       = 7
   integer, parameter :: REFERENCE_MODEL_SEA1D      = 8
   integer, parameter :: REFERENCE_MODEL_CCREM      = 9
+  integer, parameter :: REFERENCE_MODEL_SEMUCB     = 10   ! Berkeley reference model
   ! Mars
   integer, parameter :: REFERENCE_MODEL_SOHL       = 101
   integer, parameter :: REFERENCE_MODEL_SOHL_B     = 102
@@ -913,6 +959,7 @@
   integer, parameter :: ICRUST_BKMNS_GLAD  = 8    ! Block Mantle Spherical-Harmonics model
   integer, parameter :: ICRUST_SPIRAL      = 9    ! SPiRaL
   integer, parameter :: ICRUST_SH_MARS     = 10   ! SH mars models (define crust & mantle values)
+  integer, parameter :: ICRUST_BERKELEY    = 11   ! Berkeley crustal model
 
 ! define flag for 3D Earth model
   integer, parameter :: THREE_D_MODEL_S20RTS        = 101
@@ -935,6 +982,7 @@
   integer, parameter :: THREE_D_MODEL_SPIRAL        = 118
   integer, parameter :: THREE_D_MODEL_HETEROGEN_PREM = 119
   integer, parameter :: THREE_D_MODEL_SH_MARS       = 120
+  integer, parameter :: THREE_D_MODEL_BERKELEY      = 121
   ! inner core model
   integer, parameter :: THREE_D_MODEL_INNER_CORE_ISHII = 201
 

src/auxiliaries/rules.mk

@@ -578,6 +578,7 @@ xwrite_profile_OBJECTS += \
 	$O/lgndr.check.o \
 	$O/meshfem3D_models.check.o \
 	$O/meshfem3D_par.check_module.o \
+	$O/model_1dberkeley.check.o \
 	$O/model_1dref.check.o \
 	$O/model_1066a.check.o \
 	$O/model_ak135.check.o \
@@ -587,11 +588,13 @@ xwrite_profile_OBJECTS += \
 	$O/model_atten3D_QRFSI12.check.o \
 	$O/model_attenuation_gll.check.o \
 	$O/model_attenuation.check.o \
+	$O/model_berkeley.check.o \
 	$O/model_bkmns.check.o \
 	$O/model_case65TAY.check.o \
 	$O/model_ccrem.check.o \
 	$O/model_crust_1_0.check.o \
 	$O/model_crust_2_0.check.o \
+	$O/model_crust_berkeley.check.o \
 	$O/model_crustmaps.check.o \
 	$O/model_eucrust.check.o \
 	$O/model_epcrust.check.o \
@@ -651,6 +654,7 @@ xwrite_profile_SHARED_OBJECTS = \
 	$O/rthetaphi_xyz.shared.o \
 	$O/smooth_weights_vec.shared.o \
 	$O/sort_array_coordinates.shared.o \
+	$O/spl_A3d.cc.o \
 	$O/spline_routines.shared.o \
 	$O/write_VTK_file.shared.o \
 	$O/ylm.shared.o \