SUBROUTINE snonlin (f, fl, ijs, ijl, ig)
! ----------------------------------------------------------------------
!**** *SNONLIN* - COMPUTATION OF NONLINEAR TRANSFER RATE AND ITS
!****             FUNCTIONAL DERIVATIVE (DIAGONAL TERMS ONLY) AND
!****             ADDITION TO CORRESPONDING NET EXPRESSIONS.
!     S.D. HASSELMANN.  MPI
!     G. KOMEN, P. JANSSEN   KNMI             MODIFIED TO SHALLOW WATER
!     H. GUENTHER, L. ZAMBRESKY               OPTIMIZED
!     H. GUENTHER       GKSS/ECMWF  JUNE 1991 INTERACTIONS BETWEEN DIAG-
!                                             AND PROGNOSTIC PART.
!*    PURPOSE.
!     --------
!       SEE ABOVE.
!**   INTERFACE.
!     ----------
!       *CALL* *SNONLIN (F, FL, IJS, IJL, IG)*
!          *F*   - SPECTRUM.
!          *FL*  - DIAGONAL MATRIX OF FUNCTIONAL DERIVATIVE
!          *IJS* - INDEX OF FIRST GRIDPOINT
!          *IJL* - INDEX OF LAST GRIDPOINT
!          *IG*  - BLOCK NUMBER.
!     METHOD.
!     -------
!       NONE.
!     EXTERNALS.
!     ----------
!       NONE.
!     REFERENCE.
!     ----------
!       NONE.
! ----------------------------------------------------------------------
USE params; USE indnl ; USE source ; USE stat
!SHALLOW
USE meanpa ; USE shallow
!SHALLOW
IMPLICIT NONE

REAL   , INTENT(IN) :: f(0:niblo,nang,nfre)
REAL   , INTENT(OUT):: fl(0:niblo,nang,nfre)
INTEGER, INTENT(IN) :: ijs, ijl, ig
INTEGER :: k, k1, k2, k11, k21, kh, ic, ij, im, im1, ip, ip1
INTEGER :: mc, mm, mm1, mp, mp1
REAL :: ffacp, ffacp1, ftail, ffacm1
REAL :: fklamp, fklamp1, fklamp2, fklampa, fklapa2, fklampb, fklapb2
REAL :: fklap12, fklap22, fklamm, fklamm1, fklamma, fklammb, fklamm2
REAL :: fklama2, fklamb2, fklam12, fklam22
REAL :: gw1, gw2, gw3, gw4, gw5, gw6, gw7, gw8, fad1, fad2, fcen, fij, sam, sap
REAL,ALLOCATABLE :: ftemp(:), ad(:), delad(:), delap(:), delam(:), enh(:)


ALLOCATE(ftemp(ijs:ijl)) ; ALLOCATE(ad(ijs:ijl))    ; ALLOCATE(delad(ijs:ijl))
ALLOCATE(delap(ijs:ijl)) ; ALLOCATE(delam(ijs:ijl)) ; ALLOCATE(enh(ijs:ijl))

!SHALLOW
!*    1. SHALLOW WATER INITIALISATION.
!        -----------------------------
IF( ishallo .NE. 1 ) THEN
   DO ij = ijs, ijl
      enh(ij) = 0.75 * depth(ij,ig) * akmean(ij)
      enh(ij) = MAX (enh(ij), .5)
      enh(ij) = 1. + (5.5 / enh(ij)) * (1. - .833 * enh(ij)) * EXP( - 1.25 * enh(ij))
   ENDDO
ENDIF
!SHALLOW

!*    2. FREQUENCY LOOP.
!        ---------------
DO mc = 1, nfre+4
   mp     = ikp(mc)
   mp1    = ikp1(mc)
   mm     = ikm(mc)
   mm1    = ikm1(mc)
   ffacp  = 1.
   ffacp1 = 1.
   ffacm1 = 1.
   ftail  = 1.
   ic     = mc
   ip     = mp
   ip1    = mp1
   im     = mm
   im1    = mm1
   IF( ip1 .GT. nfre ) THEN
      ffacp1 = frh(ip1 - nfre+1)
      ip1    = nfre
      IF( ip .GT. nfre ) THEN
         ffacp = frh(ip - nfre+1)
         ip    = nfre
         IF( ic .GT. nfre ) THEN
            ftail = frh(ic - nfre+1)
            ic    = nfre
            IF( im1 .GT. nfre) THEN
               ffacm1 = frh(im1 - nfre+1)
               im1    = nfre
            ENDIF
         ENDIF
      ENDIF
   ENDIF
   fklamp  = fklap(mc)
   fklamp1 = fklap1(mc)
   gw2     = fklamp1 * ffacp * dal1
   gw1     = gw2 * cl11
   gw2     = gw2 * acl1
   gw4     = fklamp * ffacp1 * dal1
   gw3     = gw4 * cl11
   gw4     = gw4 * acl1
   fklampa = fklamp * cl11
   fklampb = fklamp * acl1
   fklamp2 = fklamp1 * acl1
   fklamp1 = fklamp1 * cl11
   fklapa2 = fklampa**2
   fklapb2 = fklampb**2
   fklap12 = fklamp1**2
   fklap22 = fklamp2**2
   fklamm  = fklam(mc)
   fklamm1 = fklam1(mc)
   gw6     = fklamm1 * dal2
   gw5     = gw6 * cl21
   gw6     = gw6 * acl2
   gw8     = fklamm * ffacm1 * dal2
   gw7     = gw8 * cl21
   gw8     = gw8 * acl2
   fklamma = fklamm * cl21
   fklammb = fklamm * acl2
   fklamm2 = fklamm1 * acl2
   fklamm1 = fklamm1 * cl21
   fklama2 = fklamma**2
   fklamb2 = fklammb**2
   fklam12 = fklamm1**2
   fklam22 = fklamm2**2
   IF( ishallo .EQ. 1 ) THEN
      DO ij = ijs, ijl
         ftemp(ij) = af11(mc)
      ENDDO
   ELSE
!SHALLOW
      DO ij = ijs, ijl
         ftemp(ij) = af11(mc) * enh(ij)
      ENDDO
!SHALLOW
   ENDIF
!*    2.1 LOOP FOR ANLULAR SYMMETRY.
!         -------------------------
   DO kh = 1, 2
!*    2.1.1   ANGULAR LOOP.
!             -------------
      DO k = 1, nang
         k1  = k1w(k, kh)
         k2  = k2w(k, kh)
         k11 = k11w(k, kh)
         k21 = k21w(k, kh)
!*    2.1.1.1 LOOP OVER GRIDPOINTS.. NONLINEAR TRANSFER AND
!*            DIAGONAL MATRIX OF FUNCTIONAL DERIVATIVE.
!             ----------------------------------------------
         IF( mc .GT. 4 ) THEN
            DO ij = ijs, ijl
               sap       =  gw1 * f(ij, k1, ip)  + gw2 * f(ij, k11, ip) &
                          + gw3 * f(ij, k1, ip1) + gw4 * f(ij, k11, ip1)
               sam       =  gw5 * f(ij, k2, im)  + gw6 * f(ij, k21, im) &
                          + gw7 * f(ij, k2, im1) + gw8 * f(ij, k21, im1)
               fij       = f(ij,k,ic) * ftail
               fad1      = fij * (sap + sam)
               fad2      = fad1 - 2. * sap * sam
               fad1      = fad1 + fad2
               fcen      = ftemp(ij) * fij
               ad(ij)    = fad2 * fcen
               delad(ij) = fad1 * ftemp (ij)
               delap(ij) = (fij - 2. * sam) * dal1 * fcen
               delam(ij) = (fij - 2. * sap) * dal2 * fcen
            ENDDO
            DO ij = ijs, ijl
               sl(ij, k2, mm)  = sl(ij, k2, mm)  + ad(ij) * fklamm1
               sl(ij, k21, mm) = sl(ij, k21, mm) + ad(ij) * fklamm2
               fl(ij, k2, mm)  = fl(ij, k2, mm)  + delam(ij) * fklam12
               fl(ij, k21, mm) = fl(ij, k21, mm) + delam(ij) * fklam22
            ENDDO
            IF( mm1 .LE. nfre ) THEN
               DO ij = ijs, ijl
                  sl(ij, k2, mm1)  = sl(ij, k2, mm1) + ad(ij) * fklamma
                  sl(ij, k21, mm1) = sl(ij, k21, mm1) + ad(ij) * fklammb
                  fl(ij, k2, mm1)  = fl(ij, k2, mm1)  + delam(ij) * fklama2
                  fl(ij, k21, mm1) = fl(ij, k21, mm1) + delam(ij) * fklamb2
               ENDDO
               IF( mc .LE. nfre ) THEN
                  DO ij = ijs, ijl
                     sl(ij, k, mc) = sl(ij, k, mc) - 2. * ad(ij)
                     fl(ij, k, mc) = fl(ij, k, mc) - 2. * delad(ij)
                  ENDDO
                  IF( mp .LE. nfre ) THEN
                     DO ij = ijs, ijl
                        sl(ij, k1, mp)  = sl(ij, k1, mp)  + ad(ij) * fklamp1
                        sl(ij, k11, mp) = sl(ij, k11, mp) + ad(ij) * fklamp2
                        fl(ij, k1, mp)  = fl(ij, k1, mp)  + delap(ij) * fklap12
                        fl(ij, k11, mp) = fl(ij, k11, mp) + delap(ij) * fklap22
                     ENDDO
                     IF( mp1 .LE. nfre ) THEN
                        DO ij = ijs, ijl
                           sl(ij, k1, mp1)  = sl(ij, k1, mp1)  + ad(ij) * fklampa
                           sl(ij, k11, mp1) = sl(ij, k11, mp1) + ad(ij) * fklampb
                           fl(ij, k1, mp1)  = fl(ij, k1, mp1)  + delap(ij) * fklapa2
                           fl(ij, k11, mp1) = fl(ij, k11, mp1) + delap(ij) * fklapb2
                        ENDDO
                     ENDIF
                  ENDIF
               ENDIF
            ENDIF
         ELSE
            DO ij = ijs, ijl
               sap       =  gw1 * f(ij, k1, ip)  + gw2 * f(ij, k11, ip) &
                          + gw3 * f(ij, k1, ip1) + gw4 * f(ij, k11, ip1)
!     changed in 13-feb-97
               fij       = f(ij, k, ic)
               fad2      = fij * sap
               fad1      = 2. * fad2
               fcen      = ftemp(ij) * fij
               ad(ij)    = fad2 * fcen
               delad(ij) = fad1 * ftemp(ij)
               delap(ij) = fij  * dal1 * fcen
            ENDDO
            DO ij = ijs, ijl
               sl(ij, k, mc)    = sl(ij, k, mc)    - 2. * ad(ij)
               sl(ij, k1, mp)   = sl(ij, k1, mp)   + ad(ij) * fklamp1
               sl(ij, k11, mp)  = sl(ij, k11, mp)  + ad(ij) * fklamp2
               sl(ij, k1, mp1)  = sl(ij, k1, mp1)  + ad(ij) * fklampa
               sl(ij, k11, mp1) = sl(ij, k11, mp1) + ad(ij) * fklampb

               fl(ij, k, mc)    = fl(ij, k, mc)    - 2. * delad(ij)
               fl(ij, k1, mp)   = fl(ij, k1, mp)   + delap(ij) * fklap12
               fl(ij, k11, mp)  = fl(ij, k11, mp)  + delap(ij) * fklap22
               fl(ij, k1, mp1)  = fl(ij, k1, mp1)  + delap(ij) * fklapa2
               fl(ij, k11, mp1) = fl(ij, k11, mp1) + delap(ij) * fklapb2
            ENDDO
         ENDIF
!*    BRANCH BACK TO 2.1.1 FOR NEXT DIRECTION.
      ENDDO
!*    BRANCH BACK TO 2.1 FOR MIRROR INTERACTIONS.
   ENDDO
!*    BRANCH BACK TO 2. FOR NEXT FREQUENCY.
ENDDO
DEALLOCATE(ENH, DELAM, DELAP, DELAD, AD, FTEMP)

RETURN
END SUBROUTINE snonlin
