C
      SUBROUTINE GRIDS(FUNC,XVLO,XVEC,XVHI,NVEC)
C
C Subroutine GRIDS: Grid-search optimization to find minimum of FUNC
C
C Update history:
C    V1.0    Original version                            8-Oct-1992
C    V2.0    Tidy up a bit                              13-Oct-1992
C
      IMPLICIT NONE
C
C Arguments:
C    FUNC = R*4 function = quantity to minimize
C    XVLO = R*4(1)       = array of lower limits to XVEC
C    XVEC = R*4(1)       = array of quantities to find at minimum
C    XVHI = R*4(1)       = array of upper limits to XVEC
C    NVEC = I*4          = number of parameters
C
      REAL*4    FUNC
      EXTERNAL  FUNC
C
      REAL*4    XVLO(1),XVEC(1),XVHI(1)
      INTEGER*4 NVEC
C
C Local parameter:
C    NVMAX = I*4 = maximum size of vectors that can be dealt with
C
      INTEGER*4  NVMAX
      PARAMETER (NVMAX=   128)
C
C Local variables:
C    FMIN    = R*4        = minimum of function
C    FVAL    = R*4        = current function value
C    I1      = I*4        = counting index
C    I2      = I*4        = counting index
C    I3      = I*4        = counting index
C    ILOOP   = I*4        = counter of loops
C    INDX    = I*4        = counter of varying parameters
C    ISUM    = I*4        = counter in parameter indexing
C    IVEC    = I*4        = index of position in vector
C    NINC    = I*4        = number of values for each parameter
C    NLA     = I*4        = actual number of loops to use
C    NLOOP   = I*4        = number of grid values to try
C    NVAR    = I*4        = number of parameters that vary
C    XINC(i) = R*4(NVMAX) = vector of increments in parameters
C    XINI(i) = R*4(NVMAX) = initial guess vector
C    XMIN(i) = R*4(NVMAX) = current minimum vector
C    VALUE   = R*4        = value of parameter
C
      REAL*4    FMIN,FVAL,VALUE
      REAL*4    XINC(NVMAX),XINI(NVMAX),XMIN(NVMAX)
      INTEGER*4 I1,I2,I3,ILOOP,INDX,ISUM,IVEC,NINC,NLA,NLOOP,NVAR
C
C Check that vector length isn't too large
C
      IF (NVEC .GT. NVMAX) THEN
        WRITE (6,699)
        RETURN
      ENDIF
C
C Get value of controlling parameter
C
      CALL GETOP('GRIDS.NLOOP',VALUE)
      NLOOP=IFIX(VALUE+0.5)
C
C Get number of times each variable will be incremented. Do by
C    counting only those variables for which XVLO =/= XVHI, and
C    set to a minimum of two values per variable.
C
      NVAR=0
      DO I1=1,NVEC
        IF (XVLO(I1) .NE. XVHI(I1)) THEN
          NVAR=NVAR+1
        ENDIF
      ENDDO
C
      NINC=IFIX(10.0**(ALOG10(FLOAT(NLOOP))/FLOAT(NVAR)))
      IF (NINC .LE. 1) THEN
        NINC=2
      ENDIF
C
C Calculate actual number of loops to use, and comment on the
C    need to truncate at the maximum permitted number (NLOOP)
C    if too large.
C
      NLA=NINC**NVAR
      IF (NLA .GT. NLOOP) THEN
        WRITE (6,698) NLA,NLOOP
        NLA=NLOOP
      ENDIF
C
C Get array of increments of each variable
C
      DO I1=1,NVEC
        XINC(I1)=(XVHI(I1)-XVLO(I1))/FLOAT(NINC-1)
      ENDDO
C
C Get starting value of FUNC, save vector. Save starting vector.
C
      FMIN=FUNC(XVEC,NVEC)
      DO I1=1,NVEC
        XINI(I1)=XVEC(I1)
        XMIN(I1)=XVEC(I1)
      ENDDO
C
C Debug: show initial value
C
      WRITE (6,999) FMIN
      WRITE (6,997) NVAR,NLA,NINC
C
C Now, loop over the following NLA times. In each loop the variables
C    are assigned within the allowed range.
C
      DO ILOOP=1,NLA
C
C For each variable, calculate the corresponding value that it
C    should have in a NINC**NVAR hypercubic grid at this value of 
C    ILOOP. This means that the routine only varies those
C    parameters with XVLO =/= XVHI.
C
        ISUM=ILOOP
        INDX=0
C
        DO IVEC=1,NVEC
C
          IF (XVLO(IVEC) .NE. XVHI(IVEC)) THEN
C
            INDX=INDX+1
            I2  =NINC**(NVAR-INDX)
            I3  =1+IFIX((ISUM-1)/I2)
            ISUM=ISUM-(I3-1)*(NINC**(NVAR-INDX))
C
            VALUE=XINI(IVEC)+FLOAT(I3-1)*XINC(IVEC)
            IF     (VALUE .GT. XVHI(IVEC)) THEN
              VALUE=VALUE-XVHI(IVEC)+XVLO(IVEC)
            ELSEIF (VALUE .LT. XVLO(IVEC)) THEN
              VALUE=VALUE+XVHI(IVEC)-XVLO(IVEC)
            ENDIF
C
            XVEC(IVEC)=VALUE
          ENDIF
C
        ENDDO
C
C Get new value of function, save it and the vector if smaller than
C    the previous value.
C
        FVAL=FUNC(XVEC,NVEC)
C 
        IF (FVAL .LT. FMIN) THEN
          FMIN=FVAL
          DO I2=1,NVEC
            XMIN(I2)=XVEC(I2)
          ENDDO
C
C Debug: show progress
C
        WRITE (6,998) FMIN,ILOOP
C
        ENDIF
C
      ENDDO
C
C Copy best vector back to XVEC
C
      DO I1=1,NVEC
        XVEC(I1)=XMIN(I1)
      ENDDO
C
C End
C
      RETURN
C
 698  FORMAT (' '/
     +        ' *** Warning: ',I6,' > ',I6,' steps needed:',
     +        ' truncating grid ***')
 699  FORMAT (' '/
     +        ' *** Error: parameter vector too large ***'/)
C
 997  FORMAT (' GRIDS:         no. of variables = ',I8/
     +        '                    no. of loops = ',I8/
     +        '          no. of values/variable = ',I8)
 998  FORMAT (' GRIDS:   improved minimum value = ',
     +        1PE15.5,' at iteration = ',I6)
 999  FORMAT (' '/
     +        ' GRIDS:   initial function value = ',
     +        1PE15.5)
C
      END