      SUBROUTINE GRID
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
************************************************************************
*
*  GRID CALCULATES THE ENERGY-SURFACE RESULTING FROM VARIATION OF
*       TWO COORDINATES. THE STEP-SIZE IS STEP1 AND STEP2, AND A 11
*       BY 11 GRID OF POINTS IS GENERATED
*
*       This subroutine is extensively modified by Manyin Yi, Aug 1989.
*
*       New features:
*      1. The input geometry definition should always be the upper-left
*          corner(smallest coordinates) instead of the middle point;
*       2. The starting point for calculation can be one of the four
*          corners by setting "+/-" STEP1/2;
*       3. The grid size(max 23*23) is controlled by POINT1/2,
*          if POINT1/2 is omitted, then a size of 11*11 is assumed;
*          kwd MAX sets the max size;
*       4. The upper-left corner of the plotting grid always corresponds
*          to the smallest coordinates, the lower-right corner to the
*          largest, no matter where the calculation starts;
*       5. Restartable.
*       6. Write out UNIMAP irregular data UMP.DAT
*
************************************************************************
      common /GEOM  / GEO(3,NUMATM)
      common /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
      common /GRADNT/ GRAD(MAXPAR),GNORM
      common /GRAVEC/ COSINE
      common /MESH  / LATOM1, LPARA1, LATOM2, LPARA2
      common /GPARAM/ CURRT1,CURRT2
      common /IJLP  / IJLP, ILP, JLP, JLP1, IONE
      common /SURF  / SURF
      common /KEYWRD/ KEYWRD
      common /TITLES/ KOMENT, TITLE
      character KEYWRD*241, KOMENT*81, TITLE*81, GETNAM*80
      logical RESTRT
      logical FIRST /.true./                            !A.M. 07.1993
      dimension GD(MAXPAR),XLAST(MAXPAR),MDFP(20),XDFP(20)
      dimension SURFAC(23,23)
      dimension SURF(23*23)
      dimension UMPX(23),UMPY(23),UMPZ(23*23)
C
      SAVE FIRST                        !A.M. 07.1991

      STEP1 = READA(KEYWRD,INDEX(KEYWRD,'STEP1')+6)
      STEP2 = READA(KEYWRD,INDEX(KEYWRD,'STEP2')+6)
      NPTS1=11
      NPTS2=11
      if (INDEX(KEYWRD,' MAX').NE.0) then
         NPTS1=23
         NPTS2=23
         go to 10
      end if
      if (INDEX(KEYWRD,'POINT1').NE.0)
     1     NPTS1=ABS(READA(KEYWRD,INDEX(KEYWRD,'POINT1')+7))
      if (INDEX(KEYWRD,'POINT2').NE.0)
     1     NPTS2=ABS(READA(KEYWRD,INDEX(KEYWRD,'POINT2')+7))
   10 RESTRT=(INDEX(KEYWRD,'RESTART').NE.0)
C
C  THE TOP-LEFT VALUE OF THE FIRST AND SECOND dimensionS ARE
C      GEO(LPARA1,LATOM1) AND GEO(LPARA2,LATOM2)
C
      UMPY(1)=GEO(LPARA1,LATOM1)
      UMPX(1)=GEO(LPARA2,LATOM2)
      DEGREE=180.D0/3.14159265359D0
      if (LPARA1.NE.1) STEP1=STEP1/DEGREE
      if (LPARA2.NE.1) STEP2=STEP2/DEGREE
C
C  NOW SET THE STARTING POINT TO THE DESIRED CORNER
C
      if (STEP1.GT.0.0.AND.STEP2.GT.0.0) then
         START1=GEO(LPARA1,LATOM1)
         START2=GEO(LPARA2,LATOM2)
      end if
C BOTTOM-LEFT
      if (STEP1.LT.0.0.AND.STEP2.GT.0.0) then
         START1=GEO(LPARA1,LATOM1)+(NPTS1-1)*(ABS(STEP1))
         START2=GEO(LPARA2,LATOM2)
      end if
C TOP-RIGHT
      if (STEP1.GT.0.0.AND.STEP2.LT.0.0) then
         START1=GEO(LPARA1,LATOM1)
         START2=GEO(LPARA2,LATOM2)+ABS((NPTS2-1)*STEP2)
      end if
C BOTTOM-RIGHT
      if (STEP1.LT.0.0.AND.STEP2.LT.0.0) then
         START1=GEO(LPARA1,LATOM1)+ABS((NPTS1-1)*STEP1)
         START2=GEO(LPARA2,LATOM2)+ABS((NPTS2-1)*STEP2)
      end if
C
C  NOW TO SWEEP THROUGH THE GRID OF POINTS LEFT TO RIGHT then RIGHT
C  TO LEFT OR VISA VERSA. THIS SHOULD AVOID THE GEOMETRY OR SCF GETTING
C  MESSED UP.
C
      if (LPARA1.NE.1) then
         C1=DEGREE
      else
         C1=1.D0
      end if
      if (LPARA2.NE.1) then
         C2=DEGREE
      else
         C2=1.D0
      end if
C   THESE PARAMETERS NEED TO BE DUMPED IN '.RES'
      CURRT1=START1
      CURRT2=START2
      IONE=-1.D0
      CPUTOT=0.0
      IJLP=0
      ILP=1
      JLP=1
      JLP1=1
      SURF(1)=0.D0
C
      if (RESTRT) then
         MDFP(9)=0
         call DFPSAV(CPUTOT,XPARAM,GD,XLAST,ESCF,MDFP,XDFP)
      end if
C
      GEO(LPARA1,LATOM1)=CURRT1
      GEO(LPARA2,LATOM2)=CURRT2
      do 30 ILOOP=ILP,NPTS1
         IONE=-IONE
         do 20 JLOOP=JLP,NPTS2
            JLOOP1=0
            if (IONE.LT.0)JLOOP1=NPTS2+1
            if (RESTRT) then
               JLOOP1=JLP1
               IONE=-IONE
               RESTRT=.FALSE.
            else
               JLOOP1=JLOOP1+IONE
               JLP1=JLOOP1
            end if
            CPU1=SECOND()
            CURRT1=GEO(LPARA1,LATOM1)
            CURRT2=GEO(LPARA2,LATOM2)
            call FLEPO(XPARAM, NVAR, ESCF)
            CPU2=SECOND()
            CPU3=CPU2-CPU1
            CPUTOT=CPUTOT+CPU3
            JLP=JLP+1
            IJLP=IJLP+1
            SURF(IJLP)=ESCF
            if (FIRST) then               !A.M. 07.1993
               FIRST=.false.
               write(*,
     %'(/''       FIRST VARIABLE   SECOND VARIABLE '')')
               write(*,'(''  '',F16.5,F16.5)')
     %             GEO(LPARA1,LATOM1)*C1,GEO(LPARA2,LATOM2)*C2
            endif
            write(*,'(''+ '',F16.5,F16.5)')
     %             GEO(LPARA1,LATOM1)*C1,GEO(LPARA2,LATOM2)*C2

            write(MIOB,
     %'(/''       FIRST VARIABLE   SECOND VARIABLE         function'')')
            write(MIOB,'('' :'',F16.5,F16.5,F16.6)')
     %             GEO(LPARA1,LATOM1)*C1,GEO(LPARA2,LATOM2)*C2,ESCF
            call GEOUT(MIOB)
            GEO(LPARA2,LATOM2)=GEO(LPARA2,LATOM2)+STEP2*IONE
   20    continue
         GEO(LPARA1,LATOM1)=GEO(LPARA1,LATOM1)+STEP1
         GEO(LPARA2,LATOM2)=GEO(LPARA2,LATOM2)-STEP2*IONE
         ILP=ILP+1
         JLP=1
   30 continue
      write(MIOB,'(/10X,''HORIZONTAL: VARYING SECOND PARAMETER,'',
     %               /10X,''VERTICAL:   VARYING FIRST PARAMETER'')')
      write(MIOB,'(/10X,''WHOLE OF GRID, SUITABLE FOR PLOTTING'',//)')
C
C  ARCHIVE
      open(unit=MIOF,file=GETNAM('FOR012'),status='UNKNOWN')
c
c J/M     open(unit=20,file=GETNAM('FOR020'),status='NEW')  ! why #$@! NEW
c
      open(unit=MIOJ,file=GETNAM('FOR020'),status='UNKNOWN')
      write(MIOF,40)
      call WRTTXT(MIOF)
   40 format(' ARCHIVE FILE FOR GRID CALCULATION'/'GRID OF HEATS'/)
      write(MIOF,'(/'' TOTAL CPU TIME IN FLEPO : '',F10.3/)') CPUTOT
C
C  WRITE OUT THE GRIDS
      IONE=1.D0
      ILOOP=1
      JLOOP1=1
      do 50 IJ=1,NPTS1*NPTS2
         SURFAC(JLOOP1,ILOOP)=SURF(IJ)
         N=IJ-(IJ/NPTS2)*NPTS2
         if (N.EQ.0) then
            ILOOP=ILOOP+1
            JLOOP1=JLOOP1+IONE
            IONE=-IONE
         end if
         JLOOP1=JLOOP1+IONE
   50 continue
C
      do 60 I=2,NPTS1
   60 UMPY(I)=UMPY(1)+(I-1)*ABS(STEP1)
      do 70 I=2,NPTS2
   70 UMPX(I)=UMPX(1)+(I-1)*ABS(STEP2)
      N=0
      if (STEP1.GT.0.0.AND.STEP2.GT.0.0) then
         do 90 I=1,NPTS1
            do 80 J=1,NPTS2
               N=N+1
   80       UMPZ(N)=SURFAC(J,I)
            write(MIOB,'(11F7.2)') (SURFAC(J,I),J=1,NPTS2)
   90    write(MIOF,'(11F7.2)') (SURFAC(J,I),J=1,NPTS2)
      end if
      if (STEP1.LT.0.0.AND.STEP2.GT.0.0) then
         do 110 I=NPTS1,1,-1
            do 100 J=1,NPTS2
               N=N+1
  100       UMPZ(N)=SURFAC(J,I)
            write(MIOB,'(11F7.2)') (SURFAC(J,I),J=1,NPTS2)
  110    write(MIOF,'(11F7.2)') (SURFAC(J,I),J=1,NPTS2)
      end if
      if (STEP1.GT.0.0.AND.STEP2.LT.0.0) then
         do 130 I=1,NPTS1
            do 120 J=NPTS2,1,-1
               N=N+1
  120       UMPZ(N)=SURFAC(J,I)
            write(MIOB,'(11F7.2)') (SURFAC(J,I),J=NPTS2,1,-1)
  130    write(MIOF,'(11F7.2)') (SURFAC(J,I),J=NPTS2,1,-1)
      end if
      if (STEP1.LT.0.0.AND.STEP2.LT.0.0) then
         do 150 I=NPTS1,1,-1
            do 140 J=NPTS2,1,-1
               N=N+1
  140       UMPZ(N)=SURFAC(J,I)
            write(MIOB,'(11F7.2)') (SURFAC(J,I),J=NPTS2,1,-1)
  150    write(MIOF,'(11F7.2)') (SURFAC(J,I),J=NPTS2,1,-1)
      end if
      do 160 I=0,NPTS1-1
         do 160 J=1,NPTS2
            N=I*NPTS2+J
  160 write(MIOJ,'(3(1X,F8.3))') UMPX(J),UMPY(I+1),UMPZ(N)
      close (unit=MIOJ)
      end
