Merge pull request #6 from IvanaEscobar/eigenray

Eigenray data storage

Author: IvanaEscobar

mitgcm_input/data.ihop

@@ -2,7 +2,7 @@
 # IHOP package parameters
 # ***********************
  &IHOP_PARM01 
- IHOP_fileroot='baro',
+ IHOP_fileroot='baroA',
  IHOP_title='baroclinic_gyre+ihop',
  &
 
@@ -19,12 +19,12 @@
 # rr + ihop_step < ihop_ranges(end)
  IHOP_rr=2542.0,
 
- IHOP_runopt='A',
+ IHOP_runopt='e',
  IHOP_nalpha=4,
  IHOP_alpha=-20.0, 20.0,
  &
 
  &IHOP_PARM03
  IHOP_interpfile = 'bc-gyre_ihop-grid.nc',
- ihop_iter = 1,2,3,4,5,6,7,8,9,
+ ihop_iter = 1,2,3,4,5,6,7,8,
  &

src/bellhop.F90

@@ -295,15 +295,17 @@ SUBROUTINE IHOP_INIT ( myTime, myIter, myThid )
     CALL PRINT_MESSAGE(msgBuf, PRTFile, SQUEEZE_RIGHT, myThid)
 
     ! close all files
-    SELECT CASE ( Beam%RunType( 1 : 1 ) )
+    SELECT CASE ( Beam%RunType( 1:1 ) )
     CASE ( 'C', 'S', 'I' )  ! TL calculation
        CLOSE( SHDFile )
     CASE ( 'A', 'a' )       ! arrivals calculation
        CLOSE( ARRFile )
-       CLOSE( RAYFile )
-       IF ( writeDelay ) CLOSE( DELFile )
     CASE ( 'R', 'E' )       ! ray and eigen ray trace
        CLOSE( RAYFile )
+    CASE ( 'e' )
+       CLOSE( RAYFile ) 
+       CLOSE( ARRFile )
+       IF ( writeDelay ) CLOSE( DELFile )
     END SELECT
 
     if (numberOfProcs.gt.1) then
@@ -389,7 +391,7 @@ SUBROUTINE BellhopCore( myThid )
     END SELECT
 
       IF (ALLOCATED(U)) DEALLOCATE(U)
-      SELECT CASE ( Beam%RunType( 1 : 1 ) )
+      SELECT CASE ( Beam%RunType( 1:1 ) )
       ! for a TL calculation, allocate space for the pressure matrix
       CASE ( 'C', 'S', 'I' )        ! TL calculation
           ALLOCATE ( U( NRz_per_range, Pos%NRr ), Stat = iAllocStat )
@@ -401,13 +403,13 @@ SUBROUTINE BellhopCore( myThid )
 #endif /* IHOP_WRITE_OUT */
               STOP 'ABNORMAL END: S/R BellhopCore'
           END IF
-      CASE ( 'A', 'a', 'R', 'E' )   ! Arrivals calculation
+      CASE ( 'A', 'a', 'R', 'E', 'e' )   ! Arrivals calculation
           ALLOCATE ( U( 1, 1 ), Stat = iAllocStat )   ! open a dummy variable
       END SELECT
 
       ! for an arrivals run, allocate space for arrivals matrices
-      SELECT CASE ( Beam%RunType( 1 : 1 ) )
-      CASE ( 'A', 'a' )
+      SELECT CASE ( Beam%RunType( 1:1 ) )
+      CASE ( 'A', 'a', 'e' )
           ! allow space for at least MinNArr arrivals
           MaxNArr = MAX( ArrivalsStorage / ( NRz_per_range * Pos%NRr ), MinNArr )
 #ifdef IHOP_WRITE_OUT
@@ -446,10 +448,10 @@ SUBROUTINE BellhopCore( myThid )
     SourceDepth: DO is = 1, Pos%NSz
        xs = [ zeroRL, Pos%Sz( is ) ]   ! source coordinate, assuming source @ r=0
 
-       SELECT CASE ( Beam%RunType( 1 : 1 ) )
+       SELECT CASE ( Beam%RunType( 1:1 ) )
        CASE ( 'C', 'S', 'I' ) ! TL calculation, zero out pressure matrix
           U = 0.0
-       CASE ( 'A', 'a' )      ! Arrivals calculation, zero out arrival matrix
+       CASE ( 'A', 'a', 'e' )   ! Arrivals calculation, zero out arrival matrix
           NArr = 0
        END SELECT
 
@@ -458,7 +460,7 @@ SUBROUTINE BellhopCore( myThid )
 
        !!IESCO22: BEAM stuff !!
        RadMax = 5 * c / IHOP_freq  ! 5 wavelength max radius IEsco22: unused
-       IF ( Beam%RunType( 1 : 1 ) == 'C' ) THEN ! for Coherent TL Run
+       IF ( Beam%RunType( 1:1 ) == 'C' ) THEN ! for Coherent TL Run
        ! Are there enough rays?
           DalphaOpt = SQRT( c / ( 6.0 * IHOP_freq * Pos%Rr( Pos%NRr ) ) )
           NalphaOpt = 2 + INT( ( Angles%alpha( Angles%Nalpha ) &
@@ -496,7 +498,7 @@ SUBROUTINE BellhopCore( myThid )
              ! IEsco22: When a beam pattern isn't specified, Amp0 = 0
 
              ! Lloyd mirror pattern for semi-coherent option
-             IF ( Beam%RunType( 1 : 1 ) == 'S' ) &
+             IF ( Beam%RunType( 1:1 ) == 'S' ) &
                 Amp0 = Amp0 * SQRT( 2.0 ) * ABS( SIN( afreq / c * xs( 2 ) &
                        * SIN( Angles%alpha( ialpha ) ) ) )
              !!IESCO22: end BEAM stuff !!
@@ -517,11 +519,10 @@ SUBROUTINE BellhopCore( myThid )
              ! Write the ray trajectory to RAYFile
              IF ( Beam%RunType(1:1) == 'R') THEN   
                 CALL WriteRay2D( SrcDeclAngle, Beam%Nsteps )
-             ELSE ! Compute the contribution to the field
-                CALL WriteRay2D( SrcDeclAngle, Beam%Nsteps )
                 IF (writeDelay) CALL WriteDel2D( SrcDeclAngle, Beam%Nsteps )
+             ELSE ! Compute the contribution to the field
 
-                SELECT CASE ( Beam%Type( 1 : 1 ) )
+                SELECT CASE ( Beam%Type( 1:1 ) )
                 CASE ( 'g' )
                    CALL InfluenceGeoHatRayCen(    U, Angles%alpha( ialpha ), &
                                                   Angles%Dalpha, myThid )
@@ -541,7 +542,7 @@ SUBROUTINE BellhopCore( myThid )
 
        ! write results to disk
 
-       SELECT CASE ( Beam%RunType( 1 : 1 ) )
+       SELECT CASE ( Beam%RunType( 1:1 ) )
        CASE ( 'C', 'S', 'I' )   ! TL calculation
           CALL ScalePressure( Angles%Dalpha, ray2D( 1 )%c, Pos%Rr, U, &
                               NRz_per_range, Pos%NRr, Beam%RunType, IHOP_freq )
@@ -550,7 +551,7 @@ SUBROUTINE BellhopCore( myThid )
              IRec = IRec + 1
              WRITE( SHDFile, REC = IRec ) U( Irz1, 1 : Pos%NRr )
           END DO RcvrDepth
-       CASE ( 'A' )             ! arrivals calculation, ascii
+       CASE ( 'A', 'e' )             ! arrivals calculation, ascii
           CALL WriteArrivalsASCII(  Pos%Rr, NRz_per_range, Pos%NRr, &
                                     Beam%RunType( 4:4 ) )
        CASE ( 'a' )             ! arrivals calculation, binary
@@ -599,8 +600,8 @@ SUBROUTINE TraceRay2D( xs, alpha, Amp0, myThid )
     ray2D( 1 )%p         = [ 1.0, 0.0 ]   ! IESCO22: slowness vector
     ! second component of qv is not supported in geometric beam tracing
     ! set I.C. to 0 in hopes of saving run time
-    IF ( Beam%RunType( 2:2 ) == 'G' ) THEN ! IESCO22: geometric hat in Cartesian
-        ray2D( 1 )%q = [ 0.0, 0.0 ]        ! IESCO22: ray centered coords
+    IF ( Beam%RunType( 2:2 ) == 'G' .or. Beam%RunType( 2:2 ) == 'B')  THEN 
+        ray2D( 1 )%q = [ 0.0, 0.0 ]   ! IESCO22: geometric beam in Cartesian
     ELSE
         ray2D( 1 )%q = [ 0.0, 1.0 ]   ! IESCO22: ray centered coords 
     END IF

src/influence.F90

@@ -16,7 +16,7 @@ MODULE influence
 ! sspMod used to construct image beams in the Cerveny style beam routines
   USE ssp_mod,       only: Bdry 
   USE arr_mod,       only: WriteArrivalsASCII, WriteArrivalsBinary, AddArr
-  USE writeRay,     only: WriteRay2D
+  USE writeRay,     only: WriteRay2D, WriteDel2D
 
 ! ! USES
   implicit none
@@ -474,9 +474,16 @@ SUBROUTINE InfluenceGeoGaussianCart( U, alpha, Dalpha, myThid )
   SUBROUTINE ApplyContribution( U )
     COMPLEX, INTENT( INOUT ) :: U
 
-    SELECT CASE( Beam%RunType( 1 : 1 ) )
+    SELECT CASE( Beam%RunType( 1:1 ) )
     CASE ( 'E' )                ! eigenrays
        CALL WriteRay2D( SrcDeclAngle, iS )
+    CASE ( 'e' )                ! eigenrays AND arrivals
+       CALL WriteRay2D( SrcDeclAngle, iS )
+       IF (writeDelay) CALL WriteDel2D( SrcDeclAngle, iS )
+
+       CALL AddArr( afreq, iz, ir, Amp, phaseInt, delay, SrcDeclAngle, &
+                    RcvrDeclAngle, ray2D( iS )%NumTopBnc, &
+                    ray2D( iS )%NumBotBnc )
     CASE ( 'A', 'a' )           ! arrivals
        CALL AddArr( afreq, iz, ir, Amp, phaseInt, delay, SrcDeclAngle, &
                     RcvrDeclAngle, ray2D( iS )%NumTopBnc, &
@@ -556,10 +563,11 @@ SUBROUTINE InfluenceSGB( U, alpha, Dalpha )
              deltaz =  Pos%Rz( iz ) - x( 2 )   ! ray to rcvr distance
              ! Adeltaz    = ABS( deltaz )
              ! IF ( Adeltaz < RadiusMax ) THEN
-             SELECT CASE( Beam%RunType( 1 : 1 ) )
-             CASE ( 'E' )         ! eigenrays
+             SELECT CASE( Beam%RunType( 1:1 ) )
+             CASE ( 'E', 'e' )         ! eigenrays
                 SrcDeclAngle = rad2deg * alpha   ! take-off angle in degrees
                 CALL WriteRay2D( SrcDeclAngle, iS )
+                IF (writeDelay) CALL WriteDel2D( SrcDeclAngle, iS )
              CASE DEFAULT         ! coherent TL
                 CPA    = ABS( deltaz * ( rB - rA ) ) / SQRT( ( rB - rA )**2 &
                          + ( ray2D( iS )%x( 2 ) - ray2D( iS-1 )%x( 2 ) )**2  )

src/readenvihop.F90

@@ -614,7 +614,7 @@ SUBROUTINE ReadRunType( RunType, PlotType, myThid )
     CALL PRINT_MESSAGE( msgbuf, PRTFile, SQUEEZE_RIGHT, myThid )
 #endif /* IHOP_WRITE_OUT */
 
-    SELECT CASE ( RunType( 1 : 1 ) )
+    SELECT CASE ( RunType( 1:1 ) )
     CASE ( 'R' )
 #ifdef IHOP_WRITE_OUT
        WRITE(msgBuf,'(A)') 'Ray trace run'
@@ -649,6 +649,11 @@ SUBROUTINE ReadRunType( RunType, PlotType, myThid )
 #ifdef IHOP_WRITE_OUT
        WRITE(msgBuf,'(A)') 'Arrivals calculation, binary file output'
        CALL PRINT_MESSAGE( msgbuf, PRTFile, SQUEEZE_RIGHT, myThid )
+#endif /* IHOP_WRITE_OUT */
+    CASE ( 'e' )
+#ifdef IHOP_WRITE_OUT
+       WRITE(msgBuf,'(A)') 'Eigenrays + Arrivals run, ASCII file output'
+       CALL PRINT_MESSAGE( msgbuf, PRTFile, SQUEEZE_RIGHT, myThid )
 #endif /* IHOP_WRITE_OUT */
     CASE DEFAULT
 #ifdef IHOP_WRITE_OUT
@@ -960,7 +965,7 @@ SUBROUTINE OpenOutputFiles( fName, myTime, myIter, myThid )
         fullName = fName
     ENDIF
 
-    SELECT CASE ( Beam%RunType( 1 : 1 ) )
+    SELECT CASE ( Beam%RunType( 1:1 ) )
     CASE ( 'R', 'E' )   ! Ray trace or Eigenrays
 #ifdef IHOP_WRITE_OUT
        OPEN ( FILE = TRIM( fullName ) // '.ray', UNIT = RAYFile, &
@@ -977,9 +982,9 @@ SUBROUTINE OpenOutputFiles( fName, myTime, myIter, myThid )
 #else /* IHOP_THREED */
        WRITE( RAYFile, * ) '''rz'''
 #endif /* IHOP_THREED */
+       FLUSH( RAYFile )
 #endif /* IHOP_WRITE_OUT */
-
-    CASE ( 'A' )        ! arrival file in ascii format
+    CASE ( 'e' )        ! eigenrays + arrival file in ascii format
 #ifdef IHOP_WRITE_OUT
        OPEN ( FILE = TRIM( fullName ) // '.arr', UNIT = ARRFile, &
               FORM = 'FORMATTED' )
@@ -1008,7 +1013,7 @@ SUBROUTINE OpenOutputFiles( fName, myTime, myIter, myThid )
        WRITE( ARRFile, * ) Pos%Ntheta, Pos%theta( 1 : Pos%Ntheta )
 # endif /* IHOP_THREED */
 
-       ! IEsco22: add to arrivals output
+       ! IEsco22: add erays to arrivals output
        OPEN ( FILE = TRIM( fullName ) // '.ray', UNIT = RAYFile, &
               FORM = 'FORMATTED' )
        WRITE( RAYFile, * ) '''', Title( 1 : 50 ), ''''
@@ -1040,6 +1045,38 @@ SUBROUTINE OpenOutputFiles( fName, myTime, myIter, myThid )
         WRITE( DELFile, * ) '''rz'''
 # endif /* IHOP_THREED */
        ENDIF
+
+       FLUSH( RAYFile )
+       IF (writeDelay) FLUSH( DELFile )
+       FLUSH( ARRFile )
+#endif /* IHOP_WRITE_OUT */
+    CASE ( 'A' )        ! arrival file in ascii format
+#ifdef IHOP_WRITE_OUT
+       OPEN ( FILE = TRIM( fullName ) // '.arr', UNIT = ARRFile, &
+              FORM = 'FORMATTED' )
+
+# ifdef IHOP_THREED
+       WRITE( ARRFile, * ) '''3D'''
+# else /* IHOP_THREED */
+       WRITE( ARRFile, * ) '''2D'''
+# endif /* IHOP_THREED */
+
+       WRITE( ARRFile, * ) IHOP_freq
+
+       ! write source locations
+# ifdef IHOP_THREED
+       WRITE( ARRFile, * ) Pos%NSx,    Pos%Sx(    1 : Pos%NSx )
+       WRITE( ARRFile, * ) Pos%NSy,    Pos%Sy(    1 : Pos%NSy )
+# endif /* IHOP_THREED */
+       WRITE( ARRFile, * ) Pos%NSz,    Pos%Sz(    1 : Pos%NSz )
+
+       ! write receiver locations
+       WRITE( ARRFile, *    ) Pos%NRz,    Pos%Rz(    1 : Pos%NRz )
+       WRITE( ARRFile, *    ) Pos%NRr,    Pos%Rr(    1 : Pos%NRr )
+# ifdef IHOP_THREED
+       WRITE( ARRFile, * ) Pos%Ntheta, Pos%theta( 1 : Pos%Ntheta )
+# endif /* IHOP_THREED */
+       FLUSH( ARRFile )
 #endif /* IHOP_WRITE_OUT */
     CASE ( 'a' )        ! arrival file in binary format
 #ifdef IHOP_WRITE_OUT
@@ -1069,8 +1106,8 @@ SUBROUTINE OpenOutputFiles( fName, myTime, myIter, myThid )
 # ifdef IHOP_THREED
        WRITE( ARRFile    ) Pos%Ntheta, Pos%theta( 1 : Pos%Ntheta )
 # endif /* IHOP_THREED */
+       FLUSH( ARRFile )
 #endif /* IHOP_WRITE_OUT */
-
     CASE DEFAULT
        atten = 0.0
 
@@ -1263,7 +1300,7 @@ SUBROUTINE openPRTFile ( myTime, myIter, myThid )
     CALL PRINT_MESSAGE( msgbuf, PRTFile, SQUEEZE_RIGHT, myThid )
     WRITE(msgbuf,'(A)')
     CALL PRINT_MESSAGE( msgbuf, PRTFile, SQUEEZE_RIGHT, myThid )
-    WRITE( msgBuf, '(A,I,A,F20.2,A)') 'GCM iter ', myIter,' at time = ', &
+    WRITE( msgBuf, '(A,I10,A,F20.2,A)') 'GCM iter ', myIter,' at time = ', &
         myTime,' [sec]'
     CALL PRINT_MESSAGE( msgbuf, PRTFile, SQUEEZE_RIGHT, myThid )
 # ifdef ALLOW_CAL
@@ -1294,6 +1331,7 @@ SUBROUTINE resetMemory()
     USE bdry_mod,   only: Top,Bot
     USE angle_mod,  only: Angles
     USE arr_mod,    only: Narr, Arr
+    USE ihop_mod,   only: ray2D, MaxN, iStep
 
     ! From angle_mod
     IF (ALLOCATED(Angles%alpha))DEALLOCATE(Angles%alpha)
@@ -1322,6 +1360,17 @@ SUBROUTINE resetMemory()
     IF (ALLOCATED(SSP%cMat3))   DEALLOCATE(SSP%cMat3)
     IF (ALLOCATED(SSP%czMat3))  DEALLOCATE(SSP%czMat3)
     IF (ALLOCATED(SSP%Seg%r))   DEALLOCATE(SSP%Seg%r)
+    ! From ihop_mod
+    DO iStep = 1,MaxN
+        ray2D(iStep)%x = [zeroRL, zeroRL]
+        ray2D(iStep)%t = [zeroRL, zeroRL]
+        ray2D(iStep)%p = [zeroRL, zeroRL]
+        ray2D(iStep)%q = [zeroRL, zeroRL]
+        ray2D(iStep)%c = zeroRL
+        ray2D(iStep)%Amp = zeroRL
+        ray2D(iStep)%Phase = zeroRL
+        ray2D(iStep)%tau = (zeroRL, zeroRL)
+    END DO
 
   END ! SUBROUTINE resetMemory
 

src/writeray.F90

@@ -28,7 +28,7 @@ MODULE writeRay
 
 !=======================================================================
 
-  INTEGER, PRIVATE :: MaxNRayPoints = 500000   ! this is the maximum length of 
+  INTEGER, PRIVATE :: MaxNRayPoints = 50000   ! this is the maximum length of 
   ! the ray vector that is written out
   INTEGER, PRIVATE :: is, N2, iSkip
 
@@ -58,13 +58,14 @@ SUBROUTINE WriteRay2D( alpha0, Nsteps1 )
 
     ! write to ray file
 
+    WRITE(*, '(A,G,3I)') "Escobar:",alpha0, N2, ray2D(N2)%NumTopBnc, ray2D(N2)%NumBotBnc
 #ifdef IHOP_WRITE_OUT
-    WRITE( RAYFile, * ) alpha0
-    WRITE( RAYFile, * ) N2, ray2D( Nsteps1 )%NumTopBnc, &
+    WRITE( RAYFile, '(G)') alpha0
+    WRITE( RAYFile, '(3I)' ) N2, ray2D( Nsteps1 )%NumTopBnc, &
                         ray2D( Nsteps1 )%NumBotBnc
 
     DO is = 1, N2
-       WRITE( RAYFile, * ) ray2D( is )%x
+       WRITE( RAYFile, '(2G)' ) ray2D( is )%x
     END DO
 #endif /* IHOP_WRITE_OUT */
 
@@ -91,18 +92,19 @@ SUBROUTINE WriteDel2D( alpha0, Nsteps1 )
             MOD( is, iSkip ) == 0 .OR. is == Nsteps1 ) THEN
           N2 = N2 + 1
           ray2D( N2 )%x = ray2D( is )%x
+          ray2D( N2 )%tau = ray2D( is )%tau
        END IF
     END DO Stepping
 
     ! write to delay file
 
 #ifdef IHOP_WRITE_OUT
-    WRITE( DELFile, * ) alpha0
-    WRITE( DELFile, * ) N2, ray2D( Nsteps1 )%NumTopBnc, &
+    WRITE( DELFile, '(G)' ) alpha0
+    WRITE( DELFile, '(3I)' ) N2, ray2D( Nsteps1 )%NumTopBnc, &
                         ray2D( Nsteps1 )%NumBotBnc
 
     DO is = 1, N2
-       WRITE( DELFile, * ) REAL(ray2D( is )%tau), ray2D( is )%x(2)
+       WRITE( DELFile, '(2G)' ) REAL(ray2D( is )%tau), ray2D( is )%x(2)
     END DO
 #endif /* IHOP_WRITE_OUT */