diff --git a/src/meshfem3D/model_crust_berkeley.f90 b/src/meshfem3D/model_crust_berkeley.f90 index 8f545c50e..78154c0e6 100644 --- a/src/meshfem3D/model_crust_berkeley.f90 +++ b/src/meshfem3D/model_crust_berkeley.f90 @@ -148,7 +148,7 @@ subroutine model_berkeley_crust(x,theta,phi,vp,vs,rho,moho,found_crust,elem_in_c depth = (1.d0 - x) * EARTH_R_KM - call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth) + call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth,moho_only) ! using crustal values if (USE_OLD_VERSION_FORMAT) then @@ -226,7 +226,7 @@ subroutine model_berkeley_crust_aniso(x,theta,phi,vpv,vph,vsv,vsh,eta_aniso,rho, depth = (1.d0 - x) * EARTH_R_KM - call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth) + call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth,moho_only) ! using crustal values if (USE_OLD_VERSION_FORMAT) then @@ -266,7 +266,7 @@ end subroutine model_berkeley_crust_aniso !-------------------------------------------------------------------------------------------------- ! - subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth) + subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth,moho_only) use model_crust_berkeley_par @@ -274,6 +274,7 @@ subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth) double precision,intent(in) :: theta,phi,z double precision,intent(out) :: rho,vp,vsv,vsh,moho_depth + logical,intent(in) :: moho_only ! local parameters ! 4-th order GLL positions @@ -290,6 +291,12 @@ subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth) double precision,external :: moho_filtre double precision,external :: lagrange + ! initialize + rho = 0.d0 + vp = 0.d0 + vsv = 0.d0 + vsh = 0.d0 + ! get moho depth moho_depth = moho1D_depth - moho_filtre(theta,phi) @@ -297,6 +304,9 @@ subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth) !print *,"debug: [get_crust_val_csem] Moho depth:",moho_depth, & ! "moho1D_depth:",moho1D_depth,"moho_filtre:",moho_filtre(theta,phi) + ! check if anything further to do or return only moho + if (moho_only) return + ! ! horizontal interpolation for all registered depths ! diff --git a/src/shared/get_model_parameters.F90 b/src/shared/get_model_parameters.F90 index fc87ab5c3..d329a60dc 100644 --- a/src/shared/get_model_parameters.F90 +++ b/src/shared/get_model_parameters.F90 @@ -1190,7 +1190,8 @@ subroutine get_model_parameters_radii() RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS use shared_parameters, only: & - HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL,REFERENCE_1D_MODEL + HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL,REFERENCE_1D_MODEL, & + NCHUNKS,NEX_XI,NEX_ETA ! reference models use model_prem_par @@ -1439,17 +1440,17 @@ subroutine get_model_parameters_radii() ! (same as values in model_ccrem.f90) CCREM_RSURFACE = R_PLANET ROCEAN = CCREM_RSURFACE ! no ocean - RMIDDLE_CRUST = CCREM_RSURFACE - 20000.d0 ! depth = 20 km - RMOHO = CCREM_RSURFACE - 35000.d0 ! depth = 35 km - R80 = CCREM_RSURFACE - 80000.d00 ! depth = 80 km - R220 = CCREM_RSURFACE - 220000.d0 ! depth = 220 km - R400 = CCREM_RSURFACE - 410000.d0 ! depth = 410 km - CCREM depth 410km discontinuity - R600 = CCREM_RSURFACE - 600000.d0 ! depth = 600 km - R670 = CCREM_RSURFACE - 660000.d0 ! depth = 660 km - CCREM depth 660km discontinuity - R771 = CCREM_RSURFACE - 771000.d0 ! depth = 771 km (PREM) + RMIDDLE_CRUST = CCREM_RSURFACE - 20000.d0 ! depth = 20 km + RMOHO = CCREM_RSURFACE - 35000.d0 ! depth = 35 km + R80 = CCREM_RSURFACE - 80000.d00 ! depth = 80 km + R220 = CCREM_RSURFACE - 220000.d0 ! depth = 220 km + R400 = CCREM_RSURFACE - 410000.d0 ! depth = 410 km - CCREM depth 410km discontinuity + R600 = CCREM_RSURFACE - 600000.d0 ! depth = 600 km + R670 = CCREM_RSURFACE - 660000.d0 ! depth = 660 km - CCREM depth 660km discontinuity + R771 = CCREM_RSURFACE - 771000.d0 ! depth = 771 km (PREM) RTOPDDOUBLEPRIME = CCREM_RSURFACE - 2741000.d0 ! depth = 2741 km (PREM) - RCMB = CCREM_RSURFACE - 2891000.d0 ! depth = 2891 km (PREM) - RICB = CCREM_RSURFACE - 5153500.d0 ! depth = 5153.5 km + RCMB = CCREM_RSURFACE - 2891000.d0 ! depth = 2891 km (PREM) + RICB = CCREM_RSURFACE - 5153500.d0 ! depth = 5153.5 km RHO_TOP_OC = 9.9131d0 * 1000.d0 / RHOAV RHO_BOTTOM_OC = 12.1478d0 * 1000.d0 / RHOAV @@ -1457,15 +1458,15 @@ subroutine get_model_parameters_radii() case (REFERENCE_MODEL_SEMUCB) ! Berkeley SEMUCB Model - discontinuities ROCEAN = 6368000.d0 - RMIDDLE_CRUST = 6356000.d0 - RMOHO = 6341000.d0 ! moho depth = 30 km - R80 = 6291000.d0 - R120 = -1.d0 ! no d120 discontinuity, set to fictitious value - R220 = 6151000.d0 - R400 = 5961000.d0 - R600 = 5771000.d0 - R670 = 5721000.d0 - R771 = 5600000.d0 + RMIDDLE_CRUST = 6356000.d0 ! depth = 15 km + RMOHO = 6341000.d0 ! moho depth = 30 km + R80 = 6291000.d0 ! depth = 80 km + R120 = -1.d0 ! no d120 discontinuity, set to fictitious value + R220 = 6151000.d0 ! depth = 220 km + R400 = 5961000.d0 ! depth = 410 km + R600 = 5771000.d0 ! depth = 600 km + R670 = 5721000.d0 ! depth = 650 km + R771 = 5600000.d0 ! depth = 771 km RTOPDDOUBLEPRIME = 3630000.d0 RCMB = 3480000.d0 RICB = 1221500.d0 @@ -1605,6 +1606,14 @@ subroutine get_model_parameters_radii() ! moves fictitious moho closer to this 1d moho depth RMOHO_FICTITIOUS_IN_MESHER = R_PLANET - 29000.000 ! down at 29 km depth R80_FICTITIOUS_IN_MESHER = R_PLANET - 130000.000 ! down at 130 km depth + + ! coarse global mesh modification (to allow for small testing examples) + ! to avoid Jacobian errors around lat/lon ~ 80 deg / 100 deg due to stretching + if (NCHUNKS == 6 .and. min(NEX_XI,NEX_ETA) <= 48) then + RMOHO_FICTITIOUS_IN_MESHER = R_PLANET - 35000.000 ! down at 35 km depth + R80_FICTITIOUS_IN_MESHER = R_PLANET - 140000.000 ! down at 100 km depth + endif + endif endif diff --git a/src/shared/read_compute_parameters.f90 b/src/shared/read_compute_parameters.f90 index 12158836a..bd20db37d 100644 --- a/src/shared/read_compute_parameters.f90 +++ b/src/shared/read_compute_parameters.f90 @@ -524,6 +524,10 @@ subroutine rcp_check_parameters() ! for flat topography, NEX = 32 setting still okay nex_minimum = 32 endif + + ! Berkeley coarse global mesh modification (to allow for small testing examples) + if (REFERENCE_1D_MODEL == REFERENCE_MODEL_SEMUCB) nex_minimum = 32 + ! checks nex if (NEX_XI < nex_minimum) & stop 'NEX_XI must be greater to cut the sphere into slices with positive Jacobian'