adds flag for chunk center position

danielpeter · Mar 5, 2024 · 0d89563 · 0d89563
1 parent eab11e9
commit 0d89563
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 8 deletions.
diff --git a/src/shared/euler_angles.f90 b/src/shared/euler_angles.f90
@@ -41,13 +41,25 @@ subroutine euler_angles(rotation_matrix,CENTER_LONGITUDE_IN_DEGREES,CENTER_LATIT
   double precision :: alpha,beta,gamma
   double precision :: sina,cosa,sinb,cosb,sing,cosg
 
+  ! flag to move chunk center position to geographic lon/lat instead of geocentric lon/lat when ellipticity is on.
+  ! by default, the chunk is centered around geocentric position.
+  logical, parameter :: USE_GEOGRAPHIC_CENTER_POSITION = .false.
+
   ! compute colatitude and longitude and convert to radians
-  ! longitude
-  alpha = CENTER_LONGITUDE_IN_DEGREES * DEGREES_TO_RADIANS
-  ! converts geographic latitude (degrees) to geocentric colatitude theta (radians)
-  call lat_2_geocentric_colat_dble(CENTER_LATITUDE_IN_DEGREES,beta,ELLIPTICITY)
-  ! gamma rotation
-  gamma = GAMMA_ROTATION_AZIMUTH * DEGREES_TO_RADIANS
+  if (USE_GEOGRAPHIC_CENTER_POSITION) then
+    ! longitude
+    alpha = CENTER_LONGITUDE_IN_DEGREES * DEGREES_TO_RADIANS
+    ! converts geographic latitude (degrees) to geocentric colatitude theta (radians)
+    call lat_2_geocentric_colat_dble(CENTER_LATITUDE_IN_DEGREES,beta,ELLIPTICITY)
+    ! gamma rotation
+    gamma = GAMMA_ROTATION_AZIMUTH * DEGREES_TO_RADIANS
+  else
+    ! uses center lon/lat without ellipticity correction,
+    ! assuming a perfect spherical Earth where geocentric and geographic positions are the same
+    alpha = CENTER_LONGITUDE_IN_DEGREES * DEGREES_TO_RADIANS
+    beta = (90.0d0 - CENTER_LATITUDE_IN_DEGREES) * DEGREES_TO_RADIANS
+    gamma = GAMMA_ROTATION_AZIMUTH * DEGREES_TO_RADIANS
+  endif
 
   sina = dsin(alpha)
   cosa = dcos(alpha)

diff --git a/src/specfem3D/locate_sources.f90 b/src/specfem3D/locate_sources.f90
@@ -479,7 +479,7 @@ subroutine locate_sources()
         write(IMAIN,*)
         ! different output for force point sources
         if (USE_FORCE_POINT_SOURCE) then
-          write(IMAIN,*) '  using force point source: '
+          write(IMAIN,*) '  using force point source:'
           write(IMAIN,*) '    xi coordinate of source in that element: ',xi_source(isource)
           write(IMAIN,*) '    eta coordinate of source in that element: ',eta_source(isource)
           write(IMAIN,*) '    gamma coordinate of source in that element: ',gamma_source(isource)
@@ -503,7 +503,7 @@ subroutine locate_sources()
             xyz_found_subset(2,isource_in_this_subset),xyz_found_subset(3,isource_in_this_subset)
         else
           ! moment tensor
-          write(IMAIN,*) '  using moment tensor source: '
+          write(IMAIN,*) '  using moment tensor source:'
           write(IMAIN,*) '    xi coordinate of source in that element: ',xi_source(isource)
           write(IMAIN,*) '    eta coordinate of source in that element: ',eta_source(isource)
           write(IMAIN,*) '    gamma coordinate of source in that element: ',gamma_source(isource)