SUBROUTINE sinput (f, fl, ijs, ijl, ig)
!**** *SINPUT* - COMPUTATION OF INPUT SOURCE FUNCTION.
!     P.A.E.M. JANSSEN    KNMI      AUGUST    1990
!     OPTIMIZED BY : H. GUENTHER
!*    PURPOSE.
!     ---------
!       COMPUTE INPUT SOURCE FUNCTION AND STORE ADDITIVELY INTO NET
!       SOURCE FUNCTION ARRAY, ALSO COMPUTE FUNCTIONAL DERIVATIVE OF
!       INPUT SOURCE FUNCTION.
!**   INTERFACE.
!     ----------
!       *CALL* *SINPUT (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.
!     -------
!       SEE REFERENCE.
!     EXTERNALS.
!     ----------
!       NONE.
!     REFERENCE.
!     ----------
!       P. JANSSEN, J.P.O., 1989.
!       P. JANSSEN, J.P.O., 1991
! ----------------------------------------------------------------------
USE params ; USE coupl ; USE fredir ; USE meanpa ; USE source ; USE stat ; USE wind
!SHALLOW
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

REAL,ALLOCATABLE :: temp(:,:), temp1(:,:), uco(:), zco(:)
REAL,ALLOCATABLE :: ucn(:), zcn(:), ufac1(:), ufac2(:), cm(:)
INCLUDE'globals.h'
REAL,parameter :: roair = 1.225, rowater = 1000.
REAL,PARAMETER :: xeps = roair / rowater, xinveps = 1. / xeps
INTEGER        :: ij, k, m
REAL           :: x, fac, tkd, zarg, zlog, const, const1
! ----------------------------------------------------------------------
ALLOCATE(temp(ijs:ijl, nang)) ; ALLOCATE(temp1(ijs:ijl, nang))
ALLOCATE(uco(ijs:ijl))        ; ALLOCATE(zco(ijs:ijl))
ALLOCATE(ucn(ijs:ijl))        ; ALLOCATE(zcn(ijs:ijl))
ALLOCATE(ufac1(ijs:ijl))      ; ALLOCATE(ufac2(ijs:ijl)) ; ALLOCATE(cm(ijs:ijl))
!*    1. PRECALCULATED ANGULAR DEPENDENCE.
!        ---------------------------------

DO k=1,nang
  tkd=th(k)
  DO ij=ijs,ijl
    temp(ij,k)  = COS(tkd-thwold(ij,ig))
    temp1(ij,k) = COS(tkd-thwnew(ij))
  ENDDO
ENDDO
! ----------------------------------------------------------------------
!*    2. LOOP OVER FREQUENCIES.
!        ----------------------

const1  = xeps*betamax/xkappa**2
DO m=1,nfre
  fac  = zpi*fr(m)
  const=fac*const1

!*      INVERSE OF PHASE VELOCITIES.
!       ----------------------------
  IF( ishallo == 1 ) THEN
    DO ij=ijs,ijl
      cm(ij) = fac/g
    ENDDO
  ELSE
    DO ij=ijs,ijl
      cm(ij) = tfak(indep(ij),m)/fac
    ENDDO
  ENDIF

!*      PRECALCULATE FREQUENCY DEPENDENCE.
!       ----------------------------------
  DO ij=ijs,ijl
    uco(ij) = usold(ij,ig)*cm(ij) + zalp
    zco(ij) = ALOG(g*z0old(ij,ig)*cm(ij)**2)
    ucn(ij) = usnew(ij)*cm(ij) + zalp
    zcn(ij) = ALOG(g*z0new(ij)*cm(ij)**2)
  ENDDO

!*    2.1 LOOP OVER DIRECTIONS.
!         ---------------------
  DO k=1,nang
    DO ij=ijs,ijl
      ufac1(ij) = 0.
      ufac2(ij) = 0.
    ENDDO
    DO ij=ijs,ijl
      IF( temp(ij,k) > 0.01 ) THEN
        x    = temp(ij,k)*uco(ij)
        zarg = xkappa/x
        zlog = zco(ij) + zarg
        IF( zlog < 0. ) THEN
          ufac1(ij) = const*EXP(zlog)*zlog**4*x**2
        ENDIF
      ENDIF
    ENDDO

    DO ij=ijs,ijl
      IF( temp1(ij,k) > 0.01 ) THEN
        x    = temp1(ij,k)*ucn(ij)
        zarg = xkappa/x
        zlog = zcn(ij) + zarg
        IF( zlog < 0. ) THEN
          ufac2(ij) = const*EXP(zlog)*zlog**4*x**2
        ENDIF
      ENDIF
    ENDDO

!*    2.2 ADDING INPUT SOURCE TERM TO NET SOURCE FUNCTION.
!         ------------------------------------------------
    DO ij=ijs,ijl
      sl(ij,k,m) = 0.5*(ufac1(ij)+ufac2(ij))*f(ij,k,m)
      fl(ij,k,m) = ufac2(ij)
    ENDDO
  ENDDO
ENDDO
DEALLOCATE(cm, ufac2, ufac1, zcn, ucn, zco, uco, temp1, temp)

RETURN
END SUBROUTINE sinput
