C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LSM_FOUR_VAL     **
C**************************
C
C **** LSM SERVICE ROUTINE ONLY ****
C       
       REAL FUNCTION LSM_FOUR_VAL (P, Q, N, M,  X, Y)
C
C
       INCLUDE 'implicit_none.finc'
C
C
C Purpose: Evaluate Fourier function value
C
C Author:  Steffi Arzt, John Campbell, July 1997
C
C Arguments:
C
       DOUBLE PRECISION P(*)
       DOUBLE PRECISION Q(*)
       INTEGER N
       INTEGER M
       DOUBLE PRECISION X
       DOUBLE PRECISION Y

C Parameters:
C  P         (R)    The Fourier 'P' coefficients. (N+1)*(M+1)
C                   numbers stored in a 1 D array with M being the
C                   faster moving index and N the slower moving index
C                   P(0,0), P(0,1) ... P(0,M)
C                   P(1,0), P(1,1) ... P(1,M)
C                   .
C                   .
C                   P(N,0), P(N,1) ... P(N,M)
C  Q         (R)    The Fourier 'Q' coefficients (stored as for P)
C  N         (R)    Polynomial order in 'x' direction
C  N         (R)    Polynomial order in 'y' direction
C  X         (R)    The value of 'x' at which the function is to be evaluated
C                   (in degrees)
C  Y         (R)    The value of 'y' at which the function is to be evaluated
C                   (in degrees)
C 
C
C====== PARAMETERS:
C
      INCLUDE 'limits.finc'
C
C====== GLOBALS:
C

C
C====== FUNCTION TYPES:
C
      DOUBLE PRECISION COSD, SIND
      EXTERNAL COSD, SIND
C
C====== LOCALS:
C
      INTEGER I, K
      INTEGER J
      
 
C
C-------------------------------------------------------------------------------
C
      LSM_FOUR_VAL = 0.0
      K = 0
      DO 20 I = 0, N
         DO 10 J = 0, M
            K = K + 1
            LSM_FOUR_VAL = LSM_FOUR_VAL + P(K)*SIND(I*X+J*Y)+
     +                     Q(K)*COSD(I*X+J*Y)
 10      CONTINUE 
 20   CONTINUE
      RETURN
      END