#ifdef propagation1

SUBROUTINE propags1 (isub, ijss, ijll, ijsst, ijllt)
USE params ; USE cbound ; USE current ; USE fredir  ; USE gridpar
USE map    ; USE ndtf   ; USE refdot  ; USE shallow ; USE source
USE stat   ; USE testo  ; USE ubuf
USE spe,only:fl1,fl3
IMPLICIT NONE

INTEGER, INTENT(IN) :: isub, ijss, ijll, ijsst, ijllt
INCLUDE'globals.h'
REAL,PARAMETER    :: co = 1.1
REAL, ALLOCATABLE :: dco(:)   , dp1(:)   , dp2(:)    , dpn(:) , dps(:)
REAL, ALLOCATABLE :: dle(:)   , dlw(:)   , dla(:)    , dphd(:), dlad(:)
REAL, ALLOCATABLE :: dop(:)   , dom(:)   , dtp(:)    , dtm(:) , drgp(:)
REAL, ALLOCATABLE :: drdp(:)  , drdm(:)  , drcp(:)   , drcm(:)
REAL, ALLOCATABLE :: dtc(:)   , cgond(:) , dtcd(:)   , dph(:)
REAL, ALLOCATABLE :: dlsw(:)  , dlse(:)  , dlnw(:)   , dlne(:)
REAL, ALLOCATABLE :: sdtt(:)  , dnoo(:)  , deao(:)   , drgm(:)
REAL, ALLOCATABLE :: coloio(:), colain(:), colaout(:)
REAL, ALLOCATABLE :: dnoi(:)  , deai(:)  , scross(:)
REAL, ALLOCATABLE :: sdiago(:)

INTEGER :: ig, m, k, jh, ij, KP1, km1 , ijla, ijph, mm1, mp1
REAL    :: sd, cd, delfr0, delpro, delph0, dea, dno, dtt, dth, dpno, dpso, dlea, dlwe
REAL    :: cgc, cgs, dfp, sda, cda, sp, sm,  dfm, dthp, dthm, della0, delth0, tanph
! ----------------------------------------------------------------------
!*    0. SPECTRUM AT LAND TO ZERO.
!        -------------------------
ig = 1

DO m = 1, nfre
   DO k = 1, nang
      fl3(0, k, m) = 0.
   ENDDO
ENDDO

!*    0.1 READ REFRACTION DOT TERMS.
!         --------------------------
!
IF( (irefra .NE. 0) .AND. (isub .EQ. 1) )  &
                            CALL dotdc (ijss, ijll, ig, igl, ishallo, irefra)
!
!*    0.2 SPHERICAL OR CARTESIAN GRID?
!         ----------------------------
ALLOCATE(dco(ijsst:ijllt) ) ; ALLOCATE(dp1(ijss:ijll) ) ; ALLOCATE(dp2(ijss:ijll) )

IF( icase .EQ. 1 ) THEN
!*    0.2.1 SPHERICAL GRID.
!           ---------------
!*    0.2.1.1 COSINE OF LATITUDE.
!             -------------------
   DO IJ = IJSST, IJLLT
      JH = KXLT (IJ, IG)
      DCO (IJ) = 1. / COSPH (JH)
   ENDDO
!*    0.2.1.2 COMPUTE COS PHI FACTOR FOR ADJOINING GRID POINT.
!             ------------------------------------------------
   DO IJ = IJSS, IJLL
      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)
DELPRO = FLOAT (IDTMIN)
DELPH0 = DELPRO / DELPHI
DELLA0 = DELPRO / DELLAM
DELTH0 = 0.25 * DELPRO / DELTH

ALLOCATE(DRDP(IJSS:IJLL) ) ; ALLOCATE(DRDM(IJSS:IJLL) ) ; ALLOCATE(CGOND(0:NIBLO) )
ALLOCATE(DTC(IJSS:IJLL) )  ; ALLOCATE(DTP(IJSS:IJLL) )  ; ALLOCATE(DTM(IJSS:IJLL) )
ALLOCATE(DPH(IJSS:IJLL) )  ; ALLOCATE(DLA (IJSS:IJLL) )
!*    1.1 LOOP OVER DIRECTIONS.
!         ---------------------
!
DO K = 1, NANG
   SD = SINTH (K) * DELLA0
   CD = COSTH (K) * DELPH0

!     1.1.1 INDEX FOR ADJOINING POINTS
!           --------------------------
!      COEFFICIENTS INITIALIZATION
   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
         IF (MOD (ISUB, NDTADF (M) ) .EQ.0) THEN
            IDELFD = NDTADF (M)
!*    1.1.2.1.1 LOOP OVER GRIDPOINTS.
!               ---------------------
            DTT = 1. - NDTADF (M) * DTH * GOM (M)
            DNO = CD * GOM (M) * NDTADF (M)
            DEA = SD * GOM (M) * NDTADF (M)
            DO IJ = IJSS, IJLL
!          ADVECTION MACHINE
               FL3 (IJ, K, M) = DTT * FL1 (IJ, K, M) + DNO * FL1 (KLAT (IJ, IJPH), K, M) &
                              + DEA * FL1 (KLON (IJ, IJLA), K, M)
!
            ENDDO
!*    BRANCH BACK TO 1.1.2.1 FOR NEXT FREQUENCY.
         ENDIF
      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 = IJSS, IJLL
            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
         IF (MOD (ISUB, NDTADF (M) ) .EQ.0) THEN
            IDELFD = NDTADF (M)
!*    1.1.3.2.1 GROUP VELOCITIES.
!               -----------------
            CGOND (0) = TCGOND (NDEPTH, M)
            DO IJ = IJSST, IJLLT
               CGOND (IJ) = TCGOND (INDEP (IJ), M)
            ENDDO
!
!*    1.1.3.2.2 WEIGHTS IN INTEGRATION SCHEME.
!               ------------------------------
            IF (SD.GE.0.) THEN
               DO IJ = IJSS, IJLL
                  DLA (IJ) = SD * (CGOND (KLON (IJ, 1) ) + CGOND (IJ) )
                  DTC (IJ) = SD * (CGOND (KLON (IJ, 2) ) + CGOND (IJ) )
                ENDDO
            ELSE
               DO IJ = IJSS, IJLL
                  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 = IJSS, IJLL
                  DPH (IJ) = CD * (CGOND (KLAT (IJ, 1) ) + CGOND (IJ) )
                  DTC (IJ) = DTC (IJ) + CD * (CGOND (KLAT (IJ, 2)) + CGOND (IJ) )
               ENDDO
            ELSE
               DO IJ = IJSS, IJLL
                  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 = IJSS, IJLL
                  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 = IJSS, IJLL
               FL3 (IJ, K, M) = (1. - DTC (IJ) * NDTADF (M) ) * FL1 (IJ, K, M)       &
                              + DPH (IJ) * NDTADF (M) * FL1 (KLAT (IJ, IJPH), K, M)  &
                              + DLA (IJ) * NDTADF (M) * FL1 (KLON (IJ, IJLA), K, M)
            ENDDO
            IF (IREFRA.EQ.1) THEN
               DO IJ = IJSS, IJLL
                  FL3(IJ, K, M) = FL3(IJ, K, M) + DTP(IJ) * NDTADF(M) * FL1(IJ, KP1, M) &
                                +DTM(IJ) * NDTADF(M) * FL1(IJ, KM1, M)
               ENDDO
            ENDIF
!*    BRANCH BACK TO 1.1.3.2 FOR NEXT FREQUENCY.
         ENDIF
      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.
!         -----------------------------------------------
DEALLOCATE(DCO, DP1, DP2, DRDP, DRDM) ; DEALLOCATE(DLA, DPH, DTC, DTP, DTM, CGOND)
RETURN
! ----------------------------------------------------------------------
!*    2. PROPAGATION FOR CARTESIAN GRID
!*       WITH DEPTH AND CURRENT REFRACTION.
!        ----------------------------------
2000 CONTINUE

DELPRO = FLOAT (IDTMIN)
DELPH0 = 0.25 * DELPRO / DELPHI
DELTH0 = 0.25 * DELPRO / DELTH
DELLA0 = 0.25 * DELPRO / DELLAM
!wluco begin
!     DELFR0 = 0.25*DELPRO/(0.1*ZPI)
DELFR0 = 0.25 * DELPRO / ( (CO - 1.) * ZPI)
!wluco end
!
!*    2.1 LOOP OVER DIRECTIONS.
!         ---------------------
ALLOCATE(DRDP(IJSS:IJLL)) ; ALLOCATE(DRDM(IJSS:IJLL)) ; ALLOCATE(DRCP(IJSS:IJLL))
ALLOCATE(DRCM(IJSS:IJLL)) ; ALLOCATE(CGOND(0:NIBLO))
ALLOCATE(DLA(0:NIBLO))    ; ALLOCATE(DPH(0:NIBLO))    ; ALLOCATE(DLE(IJSS:IJLL))
ALLOCATE(DLW(IJSS:IJLL))  ; ALLOCATE(DTC(IJSS:IJLL))  ; ALLOCATE(DTP(IJSS:IJLL))
ALLOCATE(DTM(IJSS:IJLL))  ; ALLOCATE(DPN(IJSS:IJLL))  ; ALLOCATE(DPS(IJSS:IJLL))
ALLOCATE(DOP(IJSS:IJLL))  ; ALLOCATE(DOM(IJSS:IJLL))

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 = IJSS, IJLL
         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 = IJSS, IJLL
      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 (MOD (ISUB, NDTADF (M) ) .EQ.0) THEN
         IDELFD = NDTADF (M)
         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 = IJSST, IJLLT
               DLA (IJ) = U (IJ, IG) * DELLA0 + CGS
               DPH (IJ) = V (IJ, IG) * DELPH0 + CGC
            ENDDO
            DO IJ = IJSS, IJLL
               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)
!wluco begin
               DOP (IJ) = ( - DTHP + ABS (DTHP) ) / CO
               DOM (IJ) = (DTHP + ABS (DTHP) ) * CO
!wluco end
            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 = IJSST, IJLLT
               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 = IJSST, IJLLT
               DLA (IJ) = U (IJ, IG) * DELLA0 + SD * CGOND (IJ)
               DPH (IJ) = V (IJ, IG) * DELPH0 + CD * CGOND (IJ)
            ENDDO
            DO IJ = IJSS, IJLL
               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)
!wluco begin
               DOP (IJ) = ( - DTHP + ABS (DTHP) ) / CO
               DOM (IJ) = (DTHM + ABS (DTHM) ) * CO
!wluco end
            ENDDO
!SHALLOW
         ENDIF
!
!*    2.1.3.3 LOOP OVER GRIDPOINTS.
!             ---------------------
!
         DO IJ = IJSS, IJLL
            FL3 (IJ, K, M) = (1. - DTC (IJ) * NDTADF (M) ) * FL1 (IJ, K, M)    &
                           + DPN (IJ) * NDTADF (M) * FL1 (KLAT (IJ, 2), K, M) &
                           + DPS (IJ) * NDTADF (M) * FL1 (KLAT (IJ, 1), K, M) &
                           + DLE (IJ) * NDTADF (M) * FL1 (KLON (IJ, 2), K, M) &
                           + DLW (IJ) * NDTADF (M) * FL1 (KLON (IJ, 1), K, M) &
                           + DTP (IJ) * NDTADF (M) * FL1 (IJ, KP1, M)         &
                           + DTM (IJ) * NDTADF (M) * FL1 (IJ, KM1, M)         &
                           + DOP (IJ) * NDTADF (M) * FL1 (IJ, K, MP1)         &
                           + DOM (IJ) * NDTADF (M) * FL1 (IJ, K, MM1)
         ENDDO
!*    BRANCH BACK TO 2.1.3 FOR NEXT FREQUENCY.
      ENDIF
   ENDDO
!*    BRANCH BACK TO 2.1 FOR NEXT DIRECTION.
ENDDO
!*    2.2 END OF PROPAGATION FOR CARTESIAN GRID
!*        WITH DEPTH AND CURRENT REFRACTION, RETURN.
!         ------------------------------------------
DEALLOCATE(DCO, DP1, DP2)
DEALLOCATE(DRDP, DRDM, DRCP, DRCM, DLA, DPH, DLE, DLW)
DEALLOCATE (DPN, DPS, DOP, DOM, DTC, DTP, DTM, CGOND)
RETURN
! ----------------------------------------------------------------------
!*    3. PROPAGATION FOR SPHERICAL LATITUDE/LONGITUDE GRID
!*       WITHOUT OR DEPTH REFRACTION.
!        -------------------------------------------------
 3000 CONTINUE

DELPRO = FLOAT (IDTMIN)
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.
!         ---------------------
ALLOCATE(DRGP(IJSS:IJLL)) ; ALLOCATE(DRGM(IJSS:IJLL)) ; ALLOCATE(DLE(IJSS:IJLL))
ALLOCATE(DTP(IJSS:IJLL))  ; ALLOCATE(DTC(IJSS:IJLL))  ; ALLOCATE(DTM(IJSS:IJLL))
ALLOCATE(DRDP(IJSS:IJLL)) ; ALLOCATE(DRDM(IJSS:IJLL)) ; ALLOCATE(CGOND(0:NIBLO))
ALLOCATE(DPN(IJSS:IJLL))

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 = IJSS, IJLL
      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 = IJSS, IJLL
         DLE (IJ) = DCO (IJ) * SDA
      ENDDO
      IF (CD.GT.0.) THEN
         DO IJ = IJSS, IJLL
            DTC (IJ) = DLE (IJ) + CDA * (DP2 (IJ) + 1.)
            DPN (IJ) = CDA * (DP1 (IJ) + 1.)
         ENDDO
      ELSE
         DO IJ = IJSS, IJLL
            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 = IJSS, IJLL
         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.1 LOOP OVER FREQUENCIES.
!             ----------------------
      DO M = 1, NFRE
         IF (MOD (ISUB, NDTADF (M) ) .EQ.0) THEN
            IDELFD = NDTADF (M)
!*    3.1.3.3.1 LOOP OVER GRIDPOINTS.
!               ---------------------
!               ADVECTION MACHINE
            DO IJ = IJSS, IJLL
               DTT = 1. - DTC (IJ) * GOM (M) * NDTADF (M)
               FL3 (IJ, K, M) = DTT * FL1(IJ, K, M)                                  &
                              + GOM(M) * NDTADF(M)                                 &
                               *( DPN(IJ) * FL1(KLAT(IJ,IJPH),K,M)                  &
                                 +DLE(IJ) * FL1(KLON(IJ,IJLA),K,M)                  &
                                 +DTP(IJ) * FL1(IJ,KP1,M) + DTM(IJ) * FL1(IJ,KM1,M))
            ENDDO
!*    BRANCH BACK TO 3.1.3.1 FOR NEXT FREQUENCY.
         ENDIF
      ENDDO
   ELSE
!SHALLOW
!*    3.1.4 SHALLOW WATER.
!           --------------
!*    3.1.4.1 COMPUTE DEPTH REFRACTION.
!             -------------------------
      IF (IREFRA.EQ.1) THEN
         DO IJ = IJSS, IJLL
            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
         IF (MOD (ISUB, NDTADF (M) ) .EQ.0) THEN
            IDELFD = NDTADF (M)
!*    3.1.4.2.1 GROUP VELOCITIES.
!               -----------------
            CGOND (0) = TCGOND (NDEPTH, M)
            DO IJ = IJSST, IJLLT
               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 = IJSS, IJLL
                  DTC (IJ) = 1. - DCO (IJ) * SDA * NDTADF (M)                      &
                                 *(CGOND (KLON (IJ, 2) ) + CGOND (IJ) )
                  DLE (IJ) = DCO (IJ) * SDA * ( CGOND(KLON(IJ,1)) + CGOND(IJ) )
               ENDDO
            ELSE
               DO IJ = IJSS, IJLL
                  DTC (IJ) = 1. - DCO (IJ) * SDA * NDTADF (M)                      &
                                 *(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 = IJSS, IJLL
                  DTC (IJ) = DTC (IJ) - CDA * NDTADF (M)                           &
                                  *(CGOND (KLAT (IJ, 2) ) * DP2 (IJ) + CGOND (IJ) )
                  DPN (IJ) = CDA * (CGOND (KLAT (IJ, 1) ) * DP1 (IJ) + CGOND (IJ) )
               ENDDO
            ELSE
               DO IJ = IJSS, IJLL
                  DTC (IJ) = DTC (IJ) - CDA * NDTADF (M)                           &
                                   *(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 = IJSS, IJLL
                  DTHP    = DRGP(IJ) * CGOND(IJ)
                  DTHM    = DRGM(IJ) * CGOND(IJ)
                  DTC(IJ) = DTC(IJ) - (DTHP + ABS(DTHP) - DTHM + ABS(DTHM)) * NDTADF(M)
                  DTP(IJ) = - DTHP + ABS (DTHP)
                  DTM(IJ) = DTHM + ABS (DTHM)
               ENDDO
            ELSE
               DO IJ = IJSS, IJLL
                  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) ) * NDTADF(M)
                  DTP(IJ) = - DTHP + ABS(DTHP)
                  DTM(IJ) = DTHM + ABS(DTHM)
               ENDDO
            ENDIF
!*    3.1.4.2.4 LOOP OVER GRIDPOINTS.
!               ---------------------
            DO IJ = IJSS, IJLL
               FL3(IJ,K,M) = DTC(IJ) * FL1 (IJ, K, M) &
                           +NDTADF(M) * (DPN(IJ) * FL1(KLAT(IJ,IJPH),K,M)    &
                                         +DLE(IJ) * FL1(KLON(IJ, IJLA),K,M)  &
                                         +DTP(IJ) * FL1(IJ,KP1,M)            &
                                         +DTM(IJ) * FL1(IJ,KM1,M) )
            ENDDO
!*    BRANCH BACK TO 3.1.4.2 FOR NEXT FREQUENCY.
         ENDIF
      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.
!         -----------------------------------------------
DEALLOCATE(DCO, DP1, DP2)
DEALLOCATE(DRGP, DRGM, DLE, DTC, DPN) ; DEALLOCATE(DTP, DTM, DRDP, DRDM, CGOND)
RETURN
! ----------------------------------------------------------------------
!*    4. PROPAGATION FOR SPHERICAL LATITUDE/LONGITUDE GRID
!*       WITH DEPTH AND CURRENT REFRACTION.
!        -------------------------------------------------
4000 CONTINUE

DELPRO = FLOAT (IDTMIN)
DELPH0 = 0.25 * DELPRO / DELPHI
DELTH0 = 0.25 * DELPRO / DELTH
DELLA0 = 0.25 * DELPRO / DELLAM
DELFR0 = 0.25 * DELPRO / ( (CO - 1.) * ZPI)

ALLOCATE(DRGP(IJSS:IJLL)) ; ALLOCATE(DRGM(IJSS:IJLL)) ; ALLOCATE(DRDP(IJSS:IJLL))
ALLOCATE(DRDM(IJSS:IJLL)) ; ALLOCATE(DRCP(IJSS:IJLL)) ; ALLOCATE(DRCM(IJSS:IJLL))
ALLOCATE(DLA(0:NIBLO))    ; ALLOCATE(DPH(0:NIBLO))    ; ALLOCATE(DLE(IJSS:IJLL))
ALLOCATE(DLW(IJSS:IJLL))  ; ALLOCATE(DTC(IJSS:IJLL))  ; ALLOCATE(DPN(IJSS:IJLL))
ALLOCATE(DPS(IJSS:IJLL))  ; ALLOCATE(DTP(IJSS:IJLL))  ; ALLOCATE(DTM(IJSS:IJLL))
ALLOCATE(DOP(IJSS:IJLL))  ; ALLOCATE(DOM(IJSS:IJLL))  ; ALLOCATE(CGOND(0:NIBLO))

!*    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 = IJSS, IJLL
      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 = IJSS, IJLL
         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 = IJSS, IJLL
      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
      IF(MOD (ISUB, NDTADF (M) ) .EQ.0) THEN
         IDELFD = NDTADF (M)
         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 = IJSST, IJLLT
               DLA (IJ) = (U (IJ, IG) * DELLA0 + CGS) * DCO (IJ)
               DPH (IJ) = V (IJ, IG) * DELPH0 + CGC
            ENDDO
            DO IJ = IJSS, IJLL
               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) ) / CO
               DOM (IJ) = (DTHP + ABS (DTHP) ) * CO
            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 = IJSST, IJLLT
               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 = IJSST, IJLLT
               DLA (IJ) = (U (IJ, IG) * DELLA0 + SD * CGOND (IJ) ) * DCO (IJ)
               DPH (IJ) = V (IJ, IG) * DELPH0 + CD * CGOND (IJ)
            ENDDO
            DO IJ = IJSS, IJLL
               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) ) / CO
               DOM (IJ) = (DTHM + ABS (DTHM) ) * CO
!wluco end
            ENDDO
!SHALLOW
         ENDIF
!*    4.1.4.3 LOOP OVER GRIDPOINTS.
!             ---------------------
         DO IJ = IJSS, IJLL
            FL3 (IJ, K, M) = (1. - DTC (IJ) * NDTADF (M) ) * FL1 ( IJ, K, M)     &
                           +NDTADF (M) * (DPN (IJ) * FL1 (KLAT (IJ, 2), K, M)   &
                                         +DPS (IJ) * FL1 (KLAT (IJ, 1), K, M)   &
                                         +DLE (IJ) * FL1 (KLON (IJ, 2), K, M)   &
                                         +DLW (IJ) * FL1 (KLON (IJ, 1), K, M)   &
                                         +DTP (IJ) * FL1 (IJ, KP1, M)           &
                                         +DTM (IJ) * FL1 (IJ, KM1, M)           &
                                         +DOP (IJ) * FL1 (IJ, K, MP1)           &
                                         +DOM (IJ) * FL1 (IJ, K, MM1)          )
         ENDDO
!*    BRANCH BACK TO 4.1.4 FOR NEXT FREQUENCY.
      ENDIF
   ENDDO
!*    BRANCH BACK TO 4.2 FOR NEXT DIRECTION.
ENDDO
!*    4.4 END OF PROPAGATION FOR SPHERICAL GRID
!*        WITH DEPTH AND CURRENT REFRACTION.
!         -------------------------------------
!     5.0 INCLUDE ANGULAR DIFFUSION:
!         --------------------------
DEALLOCATE(DCO, DP1, DP2) ; DEALLOCATE(DRGP, DRGM, DRDP, DRDM, DRCP, DRCM, DLA, DPH)
DEALLOCATE(DLE, DLW, DTC, DPN, DPS, DTP) ; DEALLOCATE(DTM, DOP, DOM, CGOND)
RETURN

END SUBROUTINE PROPAGS1
#endif
