SUBROUTINE propags (f1, f3, ig)
! ----------------------------------------------------------------------
!**** *PROPAGS* - COMPUTATION OF A PROPAGATION TIME STEP.
!     S.D. HASSELMANN.
!     OPTIMIZED BY: L. ZAMBRESKY AND H. GUENTHER
!     MODIFIED BY   H. GUNTHER   01/06/90    -   LAND POINTS ARE TAKEN
!                             OUT OF BLOCKS AND REFRACTION INTEGRATION
!                             CORRECTED FOR N-S AND S-N PROPAGATION.
!     K.P. HUBBERT                /07/89    -   DEPTH AND CURRENT
!     S. HASSELMANN   MPIFM       /04/90        REFRACTION SHALLOW
!     H. GUNTHER   GKSS/ECMWF   17/01/91    -   MODIFIED FOR CYCLE_4
!*    PURPOSE.
!     --------
!       COMPUTATION OF A PROPAGATION TIME STEP.
!**   INTERFACE.
!     ----------
!       *CALL* *PROPAGS(F1, F3, IG)*
!          *F1* - SPECTRUM AT TIME T.
!          *F3* - SPECTRUM AT TIME T+DELT.
!          *IG* - BLOCK NUMBER.
!     METHOD.
!     -------
!       FIRST ORDER FLUX SCHEME.
!     EXTERNALS.
!     ----------
!       *DOTDC*     - READ DOT TERMS FOR REFRACTION AND SCATTER TABLE.
!     REFERENCE.
!     ----------
!       NONE.
! ----------------------------------------------------------------------
USE params ; USE current ; USE fredir ; USE gridpar ; USE map
USE refdot ; USE shallow ; USE source ; USE stat    ; USE ubuf
IMPLICIT NONE

REAL   , INTENT(IN) :: f1(0:niblo,nang,nfre)
REAL   , INTENT(OUT):: f3(0:niblo,nang,nfre)
INTEGER, INTENT(IN) :: ig
INCLUDE'globals.h'
REAL   :: DCO(NIBLO), DP1(NIBLO), DP2(NIBLO)
REAL   :: DPN(NIBLO), DPS(NIBLO), DPH(0:NIBLO)
REAL   :: DLE(NIBLO), DLW(NIBLO), DLA(0:NIBLO)
REAL   :: DOP(NIBLC), DOM(NIBLC)
REAL   :: DTP(NIBLO), DTM(NIBLO), DRGP(NIBLO), DRGM(NIBLO)
REAL   :: DRDP(NIBLD), DRDM(NIBLD), DRCP(NIBLC), DRCM(NIBLC)
REAL   :: DTC(NIBLO), CGOND(0:NIBLO)
INTEGER:: K, M, IJ, JH, KM1, KP1, MM1, MP1, IJLA, IJPH
REAL   :: CD, SD, SM, SP
REAL   :: CGC, CGS, DEA, DFM, DFP, DNO, DTH, DTT, CDA, SDA
REAL   :: DLWE, DLEA, DPSO, DPNO, DTHM, DTHP
REAL   :: DELPRO, DELPH0, DELLA0, DELTH0, DELFR0, TANPH
! ----------------------------------------------------------------------
!*    0. SPECTRUM AT LAND TO ZERO.
!        -------------------------
DO M = 1, NFRE
   DO K = 1, NANG
      F3 (0, K, M) = 0.
   ENDDO
ENDDO
!*    0.1 READ REFRACTION DOT TERMS.
!         --------------------------
IF (IREFRA.NE.0) CALL DOTDC (IJS (IG), IJL (IG), IG, IGL, ISHALLO, IREFRA)
!*    0.2 SPHERICAL OR CARTESIAN GRID?
!         ----------------------------
IF (ICASE.EQ.1) THEN
!*    0.2.1 SPHERICAL GRID.
!           ---------------
!*    0.2.1.1 COSINE OF LATITUDE.
!             -------------------
   DO IJ = 1, IJLT (IG)
      JH = KXLT (IJ, IG)
      DCO (IJ) = 1. / COSPH (JH)
   ENDDO
!*    0.2.1.2 COMPUTE COS PHI FACTOR FOR ADJOINING GRID POINT.
!             ------------------------------------------------
   DO IJ = IJS (IG), IJL (IG)
      JH = KLAT (IJ, 1)
      IF (JH.LE.0) THEN
         DP1 (IJ) = 1.
      ELSE
         DP1 (IJ) = DCO (IJ) / DCO (JH)
      ENDIF
      JH = KLAT (IJ, 2)
      IF (JH.LE.0) THEN
         DP2 (IJ) = 1.
      ELSE
         DP2 (IJ) = DCO (IJ) / DCO (JH)
      ENDIF
   ENDDO
   IF (IREFRA.NE.2) THEN
!*       BRANCH TO 3. IF WITHOUT REFRACTION OR DEPTH.
!        --------------------------------------------
      GOTO 3000
   ELSE
!*       BRANCH TO 4. IF DEPTH AND CURRENT REFRACTION.
!        ---------------------------------------------
      GOTO 4000
   ENDIF
ELSE
!*    0.2.2 CARTESIAN GRID.
!           ---------------
!*    0.2.2.1 BRANCH TO 2. IF DEPTH AND CURRENT REFRACTION.
!             ---------------------------------------------
   IF (IREFRA.EQ.2) GOTO 2000
ENDIF
! ----------------------------------------------------------------------
!*    1. PROPAGATION FOR CARTESIAN GRID
!*       WITHOUT REFRACTION OR DEPTH REFRATION.
!        --------------------------------------
DELPRO = FLOAT (IDELPRO)
DELPH0 = DELPRO / DELPHI
DELLA0 = DELPRO / DELLAM
DELTH0 = 0.25 * DELPRO / DELTH
!*    1.1 LOOP OVER DIRECTIONS.
!         ---------------------
DO K = 1, NANG
   SD = SINTH (K) * DELLA0
   CD = COSTH (K) * DELPH0
!*    1.1.1 INDEX FOR ADJOINING POINTS.
!           ---------------------------
   IF (SD.LT.0) THEN
      IJLA = 2
   ELSE
      IJLA = 1
   ENDIF
   IF (CD.LT.0) THEN
      IJPH = 2
   ELSE
      IJPH = 1
   ENDIF
   IF (ISHALLO.EQ.1) THEN
!*    1.1.2 DEEP WATER.
!           -----------
      SD = ABS (SD)
      CD = ABS (CD)
      DTH = SD+CD
!*    1.1.2.1 LOOP OVER FREQUENCIES.
!             ----------------------
      DO M = 1, NFRE
!*    1.1.2.1.1 LOOP OVER GRIDPOINTS.
!               ---------------------
         DTT = 1. - DTH * GOM (M)
         DNO = CD * GOM (M)
         DEA = SD * GOM (M)
         DO IJ = IJS (IG), IJL (IG)
            F3(IJ,K,M) = DTT * F1(IJ,K,M) + DNO * F1(KLAT(IJ,IJPH),K,M) &
                        +DEA * F1(KLON(IJ,IJLA),K,M)
         ENDDO
!*    BRANCH BACK TO 1.1.2.1 FOR NEXT FREQUENCY.
      ENDDO
   ELSE
!SHALLOW
!*    1.1.3 SHALLOW WATER.
!           --------------
      SD = 0.5 * SD
      CD = 0.5 * CD
!*    1.1.3.1 DEPTH REFRACTION.
!             -----------------
      IF (IREFRA.EQ.1) THEN
         KP1 = K + 1
         IF (KP1.GT.NANG) KP1 = 1
         KM1 = K - 1
         IF (KM1.LT.1) KM1 = NANG
         DO IJ = IJS (IG), IJL (IG)
            DRDP (IJ) = (THDD (IJ, K) + THDD (IJ, KP1) ) * DELTH0
            DRDM (IJ) = (THDD (IJ, K) + THDD (IJ, KM1) ) * DELTH0
         ENDDO
      ENDIF
!*    1.1.3.2 LOOP OVER FREQUENCIES.
!             ----------------------
      DO M = 1, NFRE
!*    1.1.3.2.1 GROUP VELOCITIES.
!               -----------------
         CGOND (0) = TCGOND (NDEPTH, M)
         DO IJ = 1, IJLT (IG)
            CGOND (IJ) = TCGOND (INDEP (IJ), M)
         ENDDO
!*    1.1.3.2.2 WEIGHTS IN INTEGRATION SCHEME.
!               ------------------------------
         IF (SD.GE.0.) THEN
            DO IJ = IJS (IG), IJL (IG)
               DLA (IJ) = SD * (CGOND (KLON (IJ, 1) ) + CGOND (IJ) )
               DTC (IJ) = SD * (CGOND (KLON (IJ, 2) ) + CGOND (IJ) )
            ENDDO
         ELSE
            DO IJ = IJS (IG), IJL (IG)
               DLA (IJ) = - SD * (CGOND (KLON (IJ, 2) ) + CGOND (IJ) )
               DTC (IJ) = - SD * (CGOND (KLON (IJ, 1) ) + CGOND (IJ) )
            ENDDO
         ENDIF
         IF (CD.GE.0.) THEN
            DO IJ = IJS (IG), IJL (IG)
               DPH (IJ) = CD * (CGOND (KLAT (IJ, 1) ) + CGOND (IJ) )
               DTC (IJ) = DTC (IJ) + CD * (CGOND (KLAT (IJ, 2) ) + CGOND (IJ) )
            ENDDO
         ELSE
            DO IJ = IJS (IG), IJL (IG)
               DPH (IJ) = - CD * (CGOND (KLAT (IJ, 2) ) + CGOND (IJ) )
               DTC (IJ) = DTC (IJ) - CD * (CGOND (KLAT (IJ, 1) ) + CGOND (IJ) )
            ENDDO
         ENDIF
         IF (IREFRA.EQ.1) THEN
            DO IJ = IJS (IG), IJL (IG)
               DTHP = FCONST (IJ, M) * DRDP (IJ)
               DTHM = FCONST (IJ, M) * DRDM (IJ)
               DTC (IJ) = DTC (IJ) + DTHP + ABS (DTHP) - DTHM + ABS (DTHM)
               DTP (IJ) = - DTHP + ABS (DTHP)
               DTM (IJ) = DTHM + ABS (DTHM)
            ENDDO
         ENDIF
!*    1.1.3.2.3 LOOP OVER GRIDPOINTS.
!               ---------------------
         DO IJ = IJS (IG), IJL (IG)
            F3 (IJ, K, M) = (1. - DTC (IJ) ) * F1 (IJ, K, M) &
             + DPH (IJ) * F1 (KLAT (IJ, IJPH), K, M) + DLA (IJ) &
             * F1 (KLON (IJ, IJLA), K, M)
         ENDDO
         IF (IREFRA.EQ.1) THEN
            DO IJ = IJS (IG), IJL (IG)
               F3(IJ,K,M) = F3(IJ,K,M) + DTP(IJ) * F1(IJ,KP1,M) + DTM(IJ) * F1(IJ,KM1,M)
            ENDDO
         ENDIF
!*    BRANCH BACK TO 1.1.3.2 FOR NEXT FREQUENCY.
      ENDDO
!SHALLOW
   ENDIF
!*    BRANCH BACK TO 1.1 FOR NEXT DIRECTION.
ENDDO
!*    1.2 END OF PROPAGATION FOR CARTESIAN GRID
!*        WITHOUT REFRACTION OR DEPTH REFRACTION, RETURN.
!         -----------------------------------------------
RETURN
! ----------------------------------------------------------------------
!*    2. PROPAGATION FOR CARTESIAN GRID
!*       WITH DEPTH AND CURRENT REFRACTION.
!        ----------------------------------
2000 CONTINUE
DELPRO = FLOAT (IDELPRO)
DELPH0 = 0.25 * DELPRO / DELPHI
DELTH0 = 0.25 * DELPRO / DELTH
DELLA0 = 0.25 * DELPRO / DELLAM
DELFR0 = 0.25 * DELPRO / (0.1 * ZPI)
!*    2.1 LOOP OVER DIRECTIONS.
!         ---------------------
DO K = 1, NANG
   KP1 = K + 1
   IF (KP1.GT.NANG) KP1 = 1
   KM1 = K - 1
   IF (KM1.LT.1) KM1 = NANG
   SD = SINTH (K) * DELLA0
   CD = COSTH (K) * DELPH0
!*    2.1.1 DEPTH REFRACTION IF SHALLOW WATER.
!           ----------------------------------
   IF (ISHALLO.NE.1) THEN
      DO IJ = IJS (IG), IJL (IG)
         DRDP (IJ) = (THDD (IJ, K) + THDD (IJ, KP1) ) * DELTH0
         DRDM (IJ) = (THDD (IJ, K) + THDD (IJ, KM1) ) * DELTH0
      ENDDO
   ENDIF
!*    2.1.2 CURRENT REFRACTION.
!           -------------------
   DO IJ = IJS (IG), IJL (IG)
      DRCP (IJ) = (THDC (IJ, K) + THDC (IJ, KP1) ) * DELTH0
      DRCM (IJ) = (THDC (IJ, K) + THDC (IJ, KM1) ) * DELTH0
   ENDDO
!*    2.1.3 LOOP OVER FREQUENCIES.
!           ----------------------
   DO M = 1, NFRE
      IF (ISHALLO.EQ.1) THEN
!*    2.1.3.1 DEEP WATER.
!             -----------
         MP1 = MIN (NFRE, M + 1)
         MM1 = MAX (1, M - 1)
         DFP = PI * 2.1 * DELFR0
!*    2.1.3.1.1 GROUP VELOCITIES.
!               -----------------
         CGS = GOM (M) * SD
         CGC = GOM (M) * CD
!*    2.1.3.1.2 WEIGHTS IN INTEGRATION SCHEME.
!               ------------------------------
         DLA (0) = CGS
         DPH (0) = CGC
         DO IJ = 1, IJLT (IG)
            DLA (IJ) = U (IJ, IG) * DELLA0 + CGS
            DPH (IJ) = V (IJ, IG) * DELPH0 + CGC
         ENDDO
         DO IJ = IJS (IG), IJL (IG)
            DLWE = DLA (IJ) + DLA (KLON (IJ, 1) )
            DLEA = DLA (IJ) + DLA (KLON (IJ, 2) )
            DLE (IJ) = - DLEA + ABS (DLEA)
            DLW (IJ) = DLWE+ABS (DLWE)
            DTC (IJ) = DLEA + ABS (DLEA) - DLWE+ABS (DLWE)
            DPSO = DPH (IJ) + DPH (KLAT (IJ, 1) )
            DPNO = DPH (IJ) + DPH (KLAT (IJ, 2) )
            DPN (IJ) = - DPNO + ABS (DPNO)
            DPS (IJ) = DPSO + ABS (DPSO)
            DTC (IJ) = DTC (IJ) + DPNO + ABS (DPNO) - DPSO + ABS (DPSO)
            DTHP = DRCP (IJ)
            DTHM = DRCM (IJ)
            DTP (IJ) = - DTHP + ABS (DTHP)
            DTM (IJ) = DTHM + ABS (DTHM)
            DTC (IJ) = DTC (IJ) + DTHP + ABS (DTHP) - DTHM + ABS (DTHM)
            DTHP = SL (IJ, K, NFRE) * DFP
            DTC (IJ) = DTC (IJ) + 2. * ABS (DTHP)
            DOP (IJ) = ( - DTHP + ABS (DTHP) ) / 1.1
            DOM (IJ) = (DTHP + ABS (DTHP) ) * 1.1
         ENDDO
      ELSE
!SHALLOW
!*    2.1.3.2 SHALLOW WATER.
!             --------------
         MP1 = MIN (NFRE, M + 1)
         MM1 = MAX (1, M - 1)
         DFP = DELFR0 / FR (M)
         DFM = DELFR0 / FR (MM1)
!*    2.1.3.2.1 GROUP VELOCITIES.
!               -----------------
         CGOND (0) = TCGOND (NDEPTH, M)
         DO IJ = 1, IJLT (IG)
            CGOND (IJ) = TCGOND (INDEP (IJ), M)
         ENDDO
!*    2.1.3.2.2 WEIGHTS IN INTEGRATION SCHEME.
!               ------------------------------
         DLA (0) = SD * CGOND (0)
         DPH (0) = CD * CGOND (0)
         DO IJ = 1, IJLT (IG)
            DLA (IJ) = U (IJ, IG) * DELLA0 + SD * CGOND (IJ)
            DPH (IJ) = V (IJ, IG) * DELPH0 + CD * CGOND (IJ)
         ENDDO
         DO IJ = IJS (IG), IJL (IG)
            DLWE = DLA (IJ) + DLA (KLON (IJ, 1) )
            DLEA = DLA (IJ) + DLA (KLON (IJ, 2) )
            DLE (IJ) = - DLEA + ABS (DLEA)
            DLW (IJ) = DLWE+ABS (DLWE)
            DTC (IJ) = DLEA + ABS (DLEA) - DLWE+ABS (DLWE)
            DPSO = DPH (IJ) + DPH (KLAT (IJ, 1) )
            DPNO = DPH (IJ) + DPH (KLAT (IJ, 2) )
            DPN (IJ) = - DPNO + ABS (DPNO)
            DPS (IJ) = DPSO + ABS (DPSO)
            DTC (IJ) = DTC (IJ) + DPNO + ABS (DPNO) - DPSO + ABS (DPSO)
            DTHP = FCONST (IJ, M) * DRDP (IJ) + DRCP (IJ)
            DTHM = FCONST (IJ, M) * DRDM (IJ) + DRCM (IJ)
            DTC (IJ) = DTC (IJ) + DTHP + ABS (DTHP) - DTHM + ABS (DTHM)
            DTP (IJ) = - DTHP + ABS (DTHP)
            DTM (IJ) = DTHM + ABS (DTHM)
            DTHP = (SL (IJ, K, M) + SL (IJ, K, MP1) ) * DFP
            DTHM = (SL (IJ, K, M) + SL (IJ, K, MM1) ) * DFM
            DTC (IJ) = DTC (IJ) + DTHP + ABS (DTHP) - DTHM + ABS (DTHM)
            DOP (IJ) = ( - DTHP + ABS (DTHP) ) / 1.1
            DOM (IJ) = (DTHM + ABS (DTHM) ) * 1.1
         ENDDO
!SHALLOW
      ENDIF
!*    2.1.3.3 LOOP OVER GRIDPOINTS.
!             ---------------------
      DO IJ = IJS (IG), IJL (IG)
         F3 (IJ, K, M) = (1. - DTC (IJ) ) * F1 (IJ, K, M) + DPN ( &
          IJ) * F1 (KLAT (IJ, 2), K, M) + DPS (IJ) * F1 (KLAT (IJ, &
          1), K, M) + DLE (IJ) * F1 (KLON (IJ, 2), K, M) + DLW (IJ) &
          * F1 (KLON (IJ, 1), K, M) + DTP (IJ) * F1 (IJ, KP1, M) &
          + DTM (IJ) * F1 (IJ, KM1, M) + DOP (IJ) * F1 (IJ, K, MP1) &
          + DOM (IJ) * F1 (IJ, K, MM1)
      ENDDO
!*    BRANCH BACK TO 2.1.3 FOR NEXT FREQUENCY.
   ENDDO
!*    BRANCH BACK TO 2.1 FOR NEXT DIRECTION.
ENDDO
!*    2.2 END OF PROPAGATION FOR CARTESIAN GRID
!*        WITH DEPTH AND CURRENT REFRACTION, RETURN.
!         ------------------------------------------
RETURN
! ----------------------------------------------------------------------
!*    3. PROPAGATION FOR SPHERICAL LATITUDE/LONGITUDE GRID
!*       WITHOUT OR DEPTH REFRACTION.
!        -------------------------------------------------
3000 CONTINUE
DELPRO = FLOAT (IDELPRO)
DELTH0 = 0.25 * DELPRO / DELTH
DELPH0 = 0.5 * DELPRO / DELPHI
IF (ISHALLO.EQ.1) THEN
   DELLA0 = DELPRO / DELLAM
ELSE
   DELLA0 = 0.5 * DELPRO / DELLAM
ENDIF
!*    3.1 LOOP OVER DIRECTIONS.
!         ---------------------
DO K = 1, NANG
   KP1 = K + 1
   IF (KP1.GT.NANG) KP1 = 1
   KM1 = K - 1
   IF (KM1.LT.1) KM1 = NANG
   SD = SINTH (K) * DELLA0
   CD = COSTH (K) * DELPH0
   SDA = ABS (SD)
   CDA = ABS (CD)
!*    3.1.1 COMPUTE GRID REFRACTION.
!           ------------------------
   SP = DELTH0 * (SINTH (K) + SINTH (KP1) ) / R
   SM = DELTH0 * (SINTH (K) + SINTH (KM1) ) / R
   DO IJ = IJS (IG), IJL (IG)
      JH = KXLT (IJ, IG)
      TANPH = SINPH (JH) * DCO (IJ)
      DRGP (IJ) = TANPH * SP
      DRGM (IJ) = TANPH * SM
   ENDDO
!*    3.1.2 INDEX FOR ADJOINING POINTS.
!           ---------------------------
   IF (SD.LT.0) THEN
      IJLA = 2
   ELSE
      IJLA = 1
   ENDIF
   IF (CD.LT.0) THEN
      IJPH = 2
   ELSE
      IJPH = 1
   ENDIF

   IF (ISHALLO.EQ.1) THEN
!*    3.1.3 DEEP WATER.
!           -----------
!*    3.1.3.1 LAT / LONG WEIGHTS IN INTEGRATION SCHEME.
!             -----------------------------------------
      DO IJ = IJS (IG), IJL (IG)
         DLE (IJ) = DCO (IJ) * SDA
      ENDDO
      IF (CD.GT.0.) THEN
         DO IJ = IJS (IG), IJL (IG)
            DTC (IJ) = DLE (IJ) + CDA * (DP2 (IJ) + 1.)
            DPN (IJ) = CDA * (DP1 (IJ) + 1.)
         ENDDO
      ELSE
         DO IJ = IJS (IG), IJL (IG)
            DTC (IJ) = DLE (IJ) + CDA * (DP1 (IJ) + 1.)
            DPN (IJ) = CDA * (DP2 (IJ) + 1.)
         ENDDO
      ENDIF
!*    3.1.3.2 REFRACTION WEIGHTS IN INTEGRATION SCHEME.
!             -----------------------------------------
      DO IJ = IJS (IG), IJL (IG)
         DTHP = DRGP (IJ)
         DTHM = DRGM (IJ)
         DTC (IJ) = DTC (IJ) + DTHP + ABS (DTHP) - DTHM + ABS (DTHM)
         DTP (IJ) = - DTHP + ABS (DTHP)
         DTM (IJ) = DTHM + ABS (DTHM)
      ENDDO
!*    3.1.3.3 LOOP OVER FREQUENCIES.
!             ----------------------
      DO M = 1, NFRE
!*    3.1.3.3.1 LOOP OVER GRIDPOINTS.
!               ---------------------
         DO IJ = IJS (IG), IJL (IG)
            DTT = 1. - DTC (IJ) * GOM (M)
            F3 (IJ, K, M) = DTT * F1 (IJ, K, M) + GOM (M) * &
             (DPN (IJ) * F1 (KLAT (IJ, IJPH), K, M) + DLE (IJ) &
             * F1 (KLON (IJ, IJLA), K, M) + DTP (IJ) * F1 (IJ, KP1, &
             M) + DTM (IJ) * F1 (IJ, KM1, M) )
         ENDDO
!*    BRANCH BACK TO 3.1.3.3 FOR NEXT FREQUENCY.
      ENDDO
   ELSE
!SHALLOW
!*    3.1.4 SHALLOW WATER.
!           --------------
!*    3.1.4.1 COMPUTE DEPTH REFRACTION.
!             -------------------------
      IF (IREFRA.EQ.1) THEN
         DO IJ = IJS (IG), IJL (IG)
            DRDP (IJ) = (THDD (IJ, K) + THDD (IJ, KP1) ) * DELTH0
            DRDM (IJ) = (THDD (IJ, K) + THDD (IJ, KM1) ) * DELTH0
         ENDDO
      ENDIF
!*    3.1.4.2 LOOP OVER FREQUENCIES.
!             ----------------------
      DO M = 1, NFRE
!*    3.1.4.2.1 GROUP VELOCITIES.
!               -----------------
         CGOND (0) = TCGOND (NDEPTH, M)
         DO IJ = 1, IJLT (IG)
            CGOND (IJ) = TCGOND (INDEP (IJ), M)
         ENDDO
!*    3.1.4.3.2 LAT / LONG WEIGHTS IN INTEGRATION SCHEME.
!               -----------------------------------------
         IF (SD.GT.0.) THEN
            DO IJ = IJS (IG), IJL (IG)
               DTC (IJ) = 1. - DCO (IJ) * SDA * (CGOND (KLON (IJ,2) ) + CGOND (IJ) )
               DLE (IJ) = DCO (IJ) * SDA * (CGOND (KLON (IJ, 1) ) + CGOND (IJ) )
            ENDDO
         ELSE
            DO IJ = IJS (IG), IJL (IG)
               DTC (IJ) = 1. - DCO (IJ) * SDA * (CGOND (KLON (IJ, 1) ) + CGOND (IJ) )
               DLE (IJ) = DCO (IJ) * SDA * (CGOND (KLON (IJ, 2) ) + CGOND (IJ) )
            ENDDO
         ENDIF
         IF (CD.GT.0.) THEN
            DO IJ = IJS (IG), IJL (IG)
               DTC(IJ) = DTC(IJ) - CDA * (CGOND(KLAT(IJ, 2) ) * DP2(IJ) + CGOND(IJ) )
               DPN (IJ) = CDA * (CGOND (KLAT (IJ, 1) ) * DP1 (IJ) + CGOND (IJ) )
            ENDDO
         ELSE
            DO IJ = IJS (IG), IJL (IG)
               DTC(IJ) = DTC(IJ) - CDA * (CGOND(KLAT(IJ,1) ) * DP1(IJ) + CGOND(IJ) )
               DPN (IJ) = CDA * (CGOND (KLAT (IJ, 2) ) * DP2 (IJ) + CGOND(IJ) )
            ENDDO
         ENDIF
!*    3.1.4.2.3 REFRACTION WEIGHTS IN INTEGRATION SCHEME.
!               -----------------------------------------
         IF (IREFRA.EQ.0) THEN
            DO IJ = IJS (IG), IJL (IG)
               DTHP = DRGP (IJ) * CGOND (IJ)
               DTHM = DRGM (IJ) * CGOND (IJ)
               DTC (IJ) = DTC (IJ) - DTHP - ABS (DTHP) + DTHM - ABS (DTHM)
               DTP (IJ) = - DTHP + ABS (DTHP)
               DTM (IJ) = DTHM + ABS (DTHM)
            ENDDO
         ELSE
            DO IJ = IJS (IG), IJL (IG)
               DTHP = DRGP (IJ) * CGOND (IJ) + FCONST (IJ, M) * DRDP (IJ)
               DTHM = DRGM (IJ) * CGOND (IJ) + FCONST (IJ, M) * DRDM (IJ)
               DTC (IJ) = DTC (IJ) - DTHP - ABS (DTHP) + DTHM - ABS (DTHM)
               DTP (IJ) = - DTHP + ABS (DTHP)
               DTM (IJ) = DTHM + ABS (DTHM)
            ENDDO
         ENDIF
!*    3.1.4.2.4 LOOP OVER GRIDPOINTS.
!               ---------------------
         DO IJ = IJS (IG), IJL (IG)
            F3 (IJ, K, M) = DTC (IJ) * F1 (IJ, K, M) + DPN (IJ) &
             * F1 (KLAT (IJ, IJPH), K, M) + DLE (IJ) * F1 (KLON ( &
             IJ, IJLA), K, M) + DTP (IJ) * F1 (IJ, KP1, M) + DTM ( &
             IJ) * F1 (IJ, KM1, M)
         ENDDO
!*    BRANCH BACK TO 3.1.4.2 FOR NEXT FREQUENCY.
      ENDDO
!SHALLOW
   ENDIF
!*    BRANCH BACK TO 3.1 FOR NEXT DIRECTION.
ENDDO
!*    3.2 END OF PROPAGATION FOR SPHERICAL GRID
!*        WITHOUT REFRACTION OR DEPTH REFRACTION, RETURN.
!         -----------------------------------------------
RETURN
! ----------------------------------------------------------------------
!*    4. PROPAGATION FOR SPHERICAL LATITUDE/LONGITUDE GRID
!*       WITH DEPTH AND CURRENT REFRACTION.
!        -------------------------------------------------
4000 CONTINUE
DELPRO = FLOAT (IDELPRO)
DELPH0 = 0.25 * DELPRO / DELPHI
DELTH0 = 0.25 * DELPRO / DELTH
DELLA0 = 0.25 * DELPRO / DELLAM
DELFR0 = 0.25 * DELPRO / (0.1 * ZPI)
!*    4.1 LOOP OVER DIRECTIONS.
!         ---------------------
DO K = 1, NANG
   KP1 = K + 1
   IF (KP1.GT.NANG) KP1 = 1
   KM1 = K - 1
   IF (KM1.LT.1) KM1 = NANG
   SD = SINTH (K) * DELLA0
   CD = COSTH (K) * DELPH0
!*    4.1.1 COMPUTE GRID REFRACTION.
!           ------------------------
   SP = DELTH0 * (SINTH (K) + SINTH (KP1) ) / R
   SM = DELTH0 * (SINTH (K) + SINTH (KM1) ) / R
   DO IJ = IJS (IG), IJL (IG)
      JH = KXLT (IJ, IG)
      TANPH = SINPH (JH) * DCO (IJ)
      DRGP (IJ) = TANPH * SP
      DRGM (IJ) = TANPH * SM
   ENDDO
!*    4.1.2 COMPUTE DEPTH REFRACTION.
!           -------------------------
   IF (ISHALLO.NE.1) THEN
      DO IJ = IJS (IG), IJL (IG)
         DRDP (IJ) = (THDD (IJ, K) + THDD (IJ, KP1) ) * DELTH0
         DRDM (IJ) = (THDD (IJ, K) + THDD (IJ, KM1) ) * DELTH0
      ENDDO
   ENDIF
!*    4.1.3 COMPUTE CURRENT REFRACTION.
!           ---------------------------
   DO IJ = IJS (IG), IJL (IG)
      DRCP (IJ) = (THDC (IJ, K) + THDC (IJ, KP1) ) * DELTH0
      DRCM (IJ) = (THDC (IJ, K) + THDC (IJ, KM1) ) * DELTH0
   ENDDO
!*    4.1.4 LOOP OVER FREQUENCIES.
!           ----------------------
   DO M = 1, NFRE
      MP1 = MIN (NFRE, M + 1)
      MM1 = MAX (1, M - 1)
      IF (ISHALLO.EQ.1) THEN
!*    4.1.4.1 DEEP WATER.
!             -----------
!*    4.1.4.1.1 GROUP VELOCITIES.
!               -----------------
         DFP = PI * 2.1 * DELFR0
         CGS = GOM (M) * SD
         CGC = GOM (M) * CD
!*    4.1.4.1.2 WEIGHTS IN INTEGRATION SCHEME.
!               ------------------------------
         DLA (0) = CGS
         DPH (0) = CGC
         DO IJ = 1, IJLT (IG)
            DLA (IJ) = (U (IJ, IG) * DELLA0 + CGS) * DCO (IJ)
            DPH (IJ) = V (IJ, IG) * DELPH0 + CGC
         ENDDO
         DO IJ = IJS (IG), IJL (IG)
            DLWE = DLA (IJ) + DLA (KLON (IJ, 1) )
            DLEA = DLA (IJ) + DLA (KLON (IJ, 2) )
            DLE (IJ) = - DLEA + ABS (DLEA)
            DLW (IJ) = DLWE+ABS (DLWE)
            DTC (IJ) = DLEA + ABS (DLEA) - DLWE+ABS (DLWE)
            DPSO = DPH (IJ) + DPH (KLAT (IJ, 1) ) * DP1 (IJ)
            DPNO = DPH (IJ) + DPH (KLAT (IJ, 2) ) * DP2 (IJ)
            DPN (IJ) = - DPNO + ABS (DPNO)
            DPS (IJ) = DPSO + ABS (DPSO)
            DTC (IJ) = DTC (IJ) + DPNO + ABS (DPNO) - DPSO + ABS (DPSO)
            DTHP = DRGP (IJ) * GOM (M) + DRCP (IJ)
            DTHM = DRGM (IJ) * GOM (M) + DRCM (IJ)
            DTC (IJ) = DTC (IJ) + DTHP + ABS (DTHP) - DTHM + ABS (DTHM)
            DTP (IJ) = - DTHP + ABS (DTHP)
            DTM (IJ) = DTHM + ABS (DTHM)
            DTHP = SL (IJ, K, NFRE) * DFP
            DTC (IJ) = DTC (IJ) + 2. * ABS (DTHP)
            DOP (IJ) = ( - DTHP + ABS (DTHP) ) / 1.1
            DOM (IJ) = (DTHP + ABS (DTHP) ) * 1.1
         ENDDO
      ELSE
!SHALLOW
!*    4.1.4.2 SHALLOW WATER.
!             --------------
!*    4.1.4.2.1 GROUP VELOCITIES.
!               -----------------
         DFP = DELFR0 / FR (M)
         DFM = DELFR0 / FR (MM1)
         CGOND (0) = TCGOND (NDEPTH, M)
         DO IJ = 1, IJLT (IG)
            CGOND (IJ) = TCGOND (INDEP (IJ), M)
         ENDDO
!*    4.1.4.2.2 LON/LAT/DIR WEIGHTS IN INTEGRATION SCHEME.
!               ------------------------------------------
         DLA (0) = SD * CGOND (0)
         DPH (0) = CD * CGOND (0)
         DO IJ = 1, IJLT (IG)
            DLA (IJ) = (U (IJ, IG) * DELLA0 + SD * CGOND (IJ) ) * DCO (IJ)
            DPH (IJ) = V (IJ, IG) * DELPH0 + CD * CGOND (IJ)
         ENDDO
         DO IJ = IJS (IG), IJL (IG)
            DLWE = DLA (IJ) + DLA (KLON (IJ, 1) )
            DLEA = DLA (IJ) + DLA (KLON (IJ, 2) )
            DLE (IJ) = - DLEA + ABS (DLEA)
            DLW (IJ) = DLWE+ABS (DLWE)
            DTC (IJ) = DLEA + ABS (DLEA) - DLWE+ABS (DLWE)
            DPSO = DPH (IJ) + DPH (KLAT (IJ, 1) ) * DP1 (IJ)
            DPNO = DPH (IJ) + DPH (KLAT (IJ, 2) ) * DP2 (IJ)
            DPN (IJ) = - DPNO + ABS (DPNO)
            DPS (IJ) = DPSO + ABS (DPSO)
            DTC (IJ) = DTC (IJ) + DPNO + ABS (DPNO) - DPSO + ABS (DPSO)
            DTHP = DRGP (IJ) * CGOND (IJ) + FCONST (IJ, M) * DRDP (IJ) + DRCP (IJ)
            DTHM = DRGM (IJ) * CGOND (IJ) + FCONST (IJ, M) * DRDM (IJ) + DRCM (IJ)
            DTC (IJ) = DTC (IJ) + DTHP + ABS (DTHP) - DTHM + ABS (DTHM)
            DTP (IJ) = - DTHP + ABS (DTHP)
            DTM (IJ) = DTHM + ABS (DTHM)
            DTHP = (SL (IJ, K, M) + SL (IJ, K, MP1) ) * DFP
            DTHM = (SL (IJ, K, M) + SL (IJ, K, MM1) ) * DFM
            DTC (IJ) = DTC (IJ) + DTHP + ABS (DTHP) - DTHM + ABS (DTHM)
            DOP (IJ) = ( - DTHP + ABS (DTHP) ) / 1.1
            DOM (IJ) = (DTHM + ABS (DTHM) ) * 1.1
         ENDDO
!SHALLOW
      ENDIF
!*    4.1.4.3 LOOP OVER GRIDPOINTS.
!             ---------------------
      DO IJ = IJS (IG), IJL (IG)
         F3 (IJ, K, M) = (1. - DTC (IJ) ) * F1 (IJ, K, M) + DPN ( &
          IJ) * F1 (KLAT (IJ, 2), K, M) + DPS (IJ) * F1 (KLAT (IJ, &
          1), K, M) + DLE (IJ) * F1 (KLON (IJ, 2), K, M) + DLW (IJ) &
          * F1 (KLON (IJ, 1), K, M) + DTP (IJ) * F1 (IJ, KP1, M) &
          + DTM (IJ) * F1 (IJ, KM1, M) + DOP (IJ) * F1 (IJ, K, MP1) &
          + DOM (IJ) * F1 (IJ, K, MM1)
      ENDDO
!*    BRANCH BACK TO 4.1.4 FOR NEXT FREQUENCY.
   ENDDO
!*    BRANCH BACK TO 4.2 FOR NEXT DIRECTION.
ENDDO
!*    4.4 END OF PROPAGATION FOR SPHERICAL GRID
!*        WITH DEPTH AND CURRENT REFRACTION, RETURN.
!         ------------------------------------------
RETURN
END SUBROUTINE propags
