Merge pull request #841 from danielpeter/devel

.github/workflows/CI.yml

@@ -11,7 +11,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Check for changes
       id: diff
@@ -86,7 +86,7 @@ jobs:
     needs: changesCheck
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Install dependencies
       run: |
@@ -150,7 +150,7 @@ jobs:
         os: [ubuntu-latest,ubuntu-20.04]
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Install packages
       run: ./.github/scripts/run_install.sh
@@ -169,7 +169,7 @@ jobs:
     needs: changesCheck
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Cache Intel oneapi packages
       id: cache-intel-oneapi
@@ -285,7 +285,7 @@ jobs:
     needs: [linuxCheck]
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Install packages
       run: ./.github/scripts/run_install.sh
@@ -303,7 +303,7 @@ jobs:
      needs: [linuxCheck]
 
      steps:
-     - uses: actions/checkout@v3
+     - uses: actions/checkout@v4
 
      - name: Install packages
        run: ./.github/scripts/run_install.sh
@@ -328,7 +328,7 @@ jobs:
     needs: [linuxCheck]
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Install packages
       run: ./.github/scripts/run_install.sh
@@ -352,7 +352,7 @@ jobs:
     needs: [linuxCheck]
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Install packages
       run: ./.github/scripts/run_install.sh
@@ -376,7 +376,7 @@ jobs:
     needs: [linuxCheck]
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Install packages
       run: ./.github/scripts/run_install.sh
@@ -401,7 +401,7 @@ jobs:
     needs: [linuxCheck]
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Install packages
       run: ./.github/scripts/run_install.sh
@@ -425,7 +425,7 @@ jobs:
     needs: [linuxCheck]
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
 
     - name: Install packages
       run: ./.github/scripts/run_install.sh

src/auxiliaries/create_movie_GMT_global.f90

@@ -304,6 +304,10 @@ program create_movie_GMT_global
     print *
   endif
 
+  ! initializes
+  mute_factor = 0.0_CUSTOM_REAL
+  t0 = 0.0_CUSTOM_REAL
+
   if (MUTE_SOURCE) then
     ! initializes
     LAT_SOURCE = -1000.0

src/auxiliaries/write_profile.f90

@@ -79,7 +79,7 @@ program xwrite_profile
     CRUSTAL,ONE_CRUST
 
   ! ellipticity
-  use meshfem_models_par, only: nspl,rspl,ellipicity_spline,ellipicity_spline2
+  use meshfem_models_par, only: nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2
 
   use shared_parameters, only: &
     doubling_index,MAX_NUMBER_OF_MESH_LAYERS,rmaxs,rmins, &
@@ -523,7 +523,7 @@ program xwrite_profile
             ! ellipticity
             if (ELLIPTICITY) then
               ! adds ellipticity factor to radius
-              call add_ellipticity_rtheta(r_prem,theta,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+              call add_ellipticity_rtheta(r_prem,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
             endif
 
             ! scale values read from routines back to true values
@@ -1151,7 +1151,7 @@ subroutine write_profile_ocean(r_prem,theta,elevation,iline,iline_ocean)
         ! ellipticity
         if (ELLIPTICITY) then
           ! adds ellipticity factor to radius
-          call add_ellipticity_rtheta(r_ocean,theta,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+          call add_ellipticity_rtheta(r_ocean,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
         endif
 
         ! ocean properties (salt water parameters from PREM)
@@ -1171,7 +1171,7 @@ subroutine write_profile_ocean(r_prem,theta,elevation,iline,iline_ocean)
         ! ellipticity
         if (ELLIPTICITY) then
           ! adds ellipticity factor to radius
-          call add_ellipticity_rtheta(r_ocean,theta,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+          call add_ellipticity_rtheta(r_ocean,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
         endif
         ! last line
         write(57,'(F10.0,7F12.2,F12.5)') sngl(r_ocean*R_PLANET),1020.0,1450.,0.0,57822.5,0.0,1450.,0.0,1.0

src/meshfem3D/compute_element_properties.f90

@@ -50,7 +50,7 @@ subroutine compute_element_properties(ispec,iregion_code,idoubling,ipass, &
     REGIONAL_MOHO_MESH
 
   ! ellipticity
-  use meshfem_models_par, only: nspl,rspl,ellipicity_spline,ellipicity_spline2
+  use meshfem_models_par, only: nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2
 
   use regions_mesh_par2, only: &
     xixstore,xiystore,xizstore, &
@@ -236,10 +236,10 @@ subroutine compute_element_properties(ispec,iregion_code,idoubling,ipass, &
       ! note: after adding ellipticity, the mesh becomes elliptical and geocentric and geodetic/geographic colatitudes differ.
       if (USE_GLL) then
         ! make the Earth's ellipticity, use GLL points
-        call get_ellipticity_gll(xstore,ystore,zstore,ispec,nspec,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+        call get_ellipticity_gll(xstore,ystore,zstore,ispec,nspec,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
       else
         ! make the Earth's ellipticity, use element anchor points
-        call get_ellipticity(xelm,yelm,zelm,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+        call get_ellipticity(xelm,yelm,zelm,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
 
         ! re-interpolates GLL point locations
         ! needed for get_model(..) routine to consider stretched locations in xstore,.. arrays
@@ -373,10 +373,10 @@ subroutine compute_element_properties(ispec,iregion_code,idoubling,ipass, &
     if (.not. is_model_with_ellipticity) then
       if (USE_GLL) then
         ! make the Earth's ellipticity, use GLL points
-        call get_ellipticity_gll(xstore,ystore,zstore,ispec,nspec,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+        call get_ellipticity_gll(xstore,ystore,zstore,ispec,nspec,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
       else
         ! make the Earth's ellipticity, use element anchor points
-        call get_ellipticity(xelm,yelm,zelm,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+        call get_ellipticity(xelm,yelm,zelm,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
       endif
 
       !debug

src/meshfem3D/gravity_integrals.F90

@@ -109,7 +109,7 @@ subroutine gravity_observation_surface()
   use meshfem_par, only: x_observation,y_observation,z_observation, &
                            lon_observation,lat_observation,ELLIPTICITY,TOPOGRAPHY,OUTPUT_FILES
 
-  use meshfem_models_par, only: nspl,rspl,ellipicity_spline,ellipicity_spline2,ibathy_topo
+  use meshfem_models_par, only: nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2,ibathy_topo
 
   implicit none
 
@@ -248,7 +248,7 @@ subroutine gravity_observation_surface()
       end select
 
       ! add ellipticity
-      if (ELLIPTICITY) call add_ellipticity(x_top,y_top,z_top,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+      if (ELLIPTICITY) call add_ellipticity(x_top,y_top,z_top,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
 
       ! converts geocentric coordinates x/y/z to geographic radius/latitude/longitude (in degrees)
       call xyz_2_rlatlon_dble(x_top,y_top,z_top,r,lat,lon,ELLIPTICITY)

src/meshfem3D/meshfem3D_models.F90

@@ -38,7 +38,7 @@ subroutine meshfem3D_models_broadcast()
   integer :: ier
 
   ! sets up spline coefficients for ellipticity
-  if (ELLIPTICITY) call make_ellipticity(nspl,rspl,ellipicity_spline,ellipicity_spline2)
+  if (ELLIPTICITY) call make_ellipticity(nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
 
   ! read topography and bathymetry file
   if (TOPOGRAPHY) then
@@ -1745,7 +1745,7 @@ subroutine meshfem3D_plot_VTK_crust_moho()
   use shared_parameters, only: LOCAL_PATH,R_PLANET_KM,R_PLANET, &
     TOPOGRAPHY,ELLIPTICITY
   use meshfem_models_par, only: ibathy_topo
-  use meshfem_models_par, only: nspl,rspl,ellipicity_spline,ellipicity_spline2
+  use meshfem_models_par, only: nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2
 
   implicit none
 
@@ -1826,7 +1826,7 @@ subroutine meshfem3D_plot_VTK_crust_moho()
       ! ellipticity
       if (ELLIPTICITY) then
         ! adds ellipticity factor to radius
-        call add_ellipticity_rtheta(r,theta,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+        call add_ellipticity_rtheta(r,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
       endif
 
       ! gets moho
@@ -1917,7 +1917,7 @@ subroutine meshfem3D_plot_VTK_topo_bathy()
   use shared_parameters, only: LOCAL_PATH,RESOLUTION_TOPO_FILE,R_PLANET,R_PLANET_KM, &
     TOPOGRAPHY,ELLIPTICITY
   use meshfem_models_par, only: ibathy_topo
-  use meshfem_models_par, only: nspl,rspl,ellipicity_spline,ellipicity_spline2
+  use meshfem_models_par, only: nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2
 
   implicit none
 
@@ -2019,7 +2019,7 @@ subroutine meshfem3D_plot_VTK_topo_bathy()
       ! ellipticity
       if (ELLIPTICITY) then
         ! adds ellipticity factor to radius
-        call add_ellipticity_rtheta(r,theta,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+        call add_ellipticity_rtheta(r,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
       endif
 
       ! gets point's position

src/meshfem3D/meshfem3D_par.f90

@@ -100,8 +100,8 @@ module meshfem_models_par
   integer, dimension(:,:),allocatable :: ibathy_topo
 
   ! for ellipticity
-  double precision,dimension(NR_DENSITY) :: rspl,ellipicity_spline,ellipicity_spline2
-  integer :: nspl
+  double precision,dimension(NR_DENSITY) :: rspl_ellip,ellipicity_spline,ellipicity_spline2
+  integer :: nspl_ellip
 
   end module meshfem_models_par
 

src/meshfem3D/model_EMC.f90

@@ -3138,7 +3138,7 @@ subroutine model_EMC_crustmantle(iregion_code,r,theta,phi,vpv,vph,vsv,vsh,eta_an
   use constants
   use shared_parameters, only: R_PLANET,R_PLANET_KM,RHOAV,TOPOGRAPHY,ELLIPTICITY
   use meshfem_models_par, only: ibathy_topo
-  use meshfem_models_par, only: rspl,ellipicity_spline,ellipicity_spline2,nspl
+  use meshfem_models_par, only: rspl_ellip,ellipicity_spline,ellipicity_spline2,nspl_ellip
 
   use model_emc_par
 
@@ -3216,7 +3216,7 @@ subroutine model_EMC_crustmantle(iregion_code,r,theta,phi,vpv,vph,vsv,vsh,eta_an
   ! ellipticity
   if (ELLIPTICITY) then
     ! adds ellipticity factor to radius (needs geocentric colatitude)
-    call add_ellipticity_rtheta(r0,theta,nspl,rspl,ellipicity_spline,ellipicity_spline2)
+    call add_ellipticity_rtheta(r0,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2)
   endif
 
   ! depth with respect to surface (in km)

src/meshfem3D/model_bkmns.f90

@@ -920,6 +920,17 @@ subroutine read_mantle_bkmns_model()
   if (ier /= 0) call exit_MPI(myrank,'Error allocating nFactors zone4 array')
   nFactors_zone4(:) = 0.d0
 
+  ! initializes
+  spline_degree_zone2 = 0
+  rmin_zone2 = 0.d0
+  rmax_zone2 = 0.d0
+  spline_degree_zone3 = 0
+  rmin_zone3 = 0.d0
+  rmax_zone3 = 0.d0
+  spline_degree_zone4 = 0
+  rmin_zone4 = 0.d0
+  rmax_zone4 = 0.d0
+
   ! zone range within [2,4]
   do izone = 2,4
     ! zone string

src/meshfem3D/model_ccrem.f90

@@ -305,7 +305,7 @@ end subroutine model_ccrem
 
   subroutine define_model_ccrem(CRUSTAL)
 
-  use constants, only: myrank,IIN,SUPPRESS_CRUSTAL_MESH
+  use constants, only: myrank,IIN,SUPPRESS_CRUSTAL_MESH,ATTENUATION_COMP_MAXIMUM
   use shared_parameters, only: RICB,RCMB,ROCEAN,ONE_CRUST
 
   use model_ccrem_par
@@ -474,6 +474,9 @@ subroutine define_model_ccrem(CRUSTAL)
     r = r_prem
 
     ! determines attenuation Q factors from corresponding PREM shell values
+    Qmu = ATTENUATION_COMP_MAXIMUM
+    Qkappa = ATTENUATION_COMP_MAXIMUM
+
     ! inner core
     if (r >= 0.d0 .and. r <= CCREM_RICB) then
       Qmu = 84.6d0

src/meshfem3D/moho_stretching.f90

@@ -81,6 +81,10 @@ subroutine moho_stretching_honor_crust(xelm,yelm,zelm, &
 
   ! note: at this point, the mesh is still perfectly spherical.
 
+  ! initializes
+  MOHO_MINIMUM_DEFAULT = 0.d0
+  MOHO_MAXIMUM_DEFAULT = 0.d0
+
   ! min/max defaults
   select case(PLANET_TYPE)
   case (IPLANET_MARS)

src/meshfem3D/write_AVS_DX_output.f90

@@ -39,7 +39,7 @@ subroutine write_AVS_DX_output(npointot,iregion_code)
 
   use meshfem_models_par, only: &
     ELLIPTICITY,MODEL_3D_MANTLE_PERTUBATIONS, &
-    nspl,rspl,ellipicity_spline,ellipicity_spline2
+    nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2
 
   use regions_mesh_par2
 
@@ -77,23 +77,23 @@ subroutine write_AVS_DX_output(npointot,iregion_code)
     call write_AVS_DX_global_faces_data(prname,nspec,iMPIcut_xi, &
                                         iMPIcut_eta,ibool, idoubling,xstore,ystore,zstore,num_ibool_AVS_DX, &
                                         mask_ibool,npointot, rhostore,kappavstore,muvstore, &
-                                        nspl,rspl,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
+                                        nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
                                         MODEL_3D_MANTLE_PERTUBATIONS, RICB,RCMB, &
                                         RTOPDDOUBLEPRIME,R670,R220,R771,R400,R120,R80,RMOHO, &
                                         RMIDDLE_CRUST,iregion_code)
 
     call write_AVS_DX_global_chunks_data(prname,nspec,iboun,ibool, &
                                          idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool, &
                                          npointot,rhostore,kappavstore,muvstore, &
-                                         nspl,rspl,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
+                                         nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
                                          MODEL_3D_MANTLE_PERTUBATIONS, &
                                          RICB,RCMB,RTOPDDOUBLEPRIME,R670,R220,R771,R400,R120,R80,RMOHO, &
                                          RMIDDLE_CRUST,iregion_code)
 
     call write_AVS_DX_surface_data(prname,nspec,iboun,ibool, &
                                    idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool,npointot, &
                                    rhostore,kappavstore,muvstore, &
-                                   nspl,rspl,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
+                                   nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
                                    MODEL_3D_MANTLE_PERTUBATIONS, &
                                    RICB,RCMB,RTOPDDOUBLEPRIME,R670,R220,R771,R400,R120,R80,RMOHO, &
                                    RMIDDLE_CRUST,iregion_code)

src/meshfem3D/write_AVS_DX_output_adios.f90

@@ -39,7 +39,7 @@ subroutine write_AVS_DX_output_adios(npointot,iregion_code, &
 
   use meshfem_models_par, only: &
     ELLIPTICITY,MODEL_3D_MANTLE_PERTUBATIONS, &
-    nspl,rspl,ellipicity_spline,ellipicity_spline2
+    nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2
 
   use regions_mesh_par2
 
@@ -125,7 +125,7 @@ subroutine write_AVS_DX_output_adios(npointot,iregion_code, &
                                               iMPIcut_xi,iMPIcut_eta, &
                                               ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool, &
                                               npointot,rhostore,kappavstore,muvstore, &
-                                              nspl,rspl,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
+                                              nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
                                               MODEL_3D_MANTLE_PERTUBATIONS, &
                                               RICB,RCMB,RTOPDDOUBLEPRIME,R670,R220,R771,R400,R120,R80,RMOHO, &
                                               RMIDDLE_CRUST,iregion_code, &
@@ -136,7 +136,7 @@ subroutine write_AVS_DX_output_adios(npointot,iregion_code, &
   call prepare_AVS_DX_global_chunks_data_adios(prname,nspec, &
                                                iboun,ibool, idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool, &
                                                npointot,rhostore,kappavstore,muvstore, &
-                                               nspl,rspl,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
+                                               nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
                                                MODEL_3D_MANTLE_PERTUBATIONS, &
                                                RICB,RCMB,RTOPDDOUBLEPRIME,R670,R220,R771,R400,R120,R80,RMOHO, &
                                                RMIDDLE_CRUST,iregion_code, &
@@ -147,7 +147,7 @@ subroutine write_AVS_DX_output_adios(npointot,iregion_code, &
   call prepare_AVS_DX_surfaces_data_adios(nspec,iboun, &
                                           ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool,npointot, &
                                           rhostore,kappavstore,muvstore, &
-                                          nspl,rspl,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
+                                          nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2,ELLIPTICITY, &
                                           MODEL_3D_MANTLE_PERTUBATIONS, &
                                           RICB,RCMB,RTOPDDOUBLEPRIME,R670,R220,R771,R400,R120,R80,RMOHO, &
                                           RMIDDLE_CRUST,iregion_code, &