C @(#)conts3.for	17.1.1.1 (ES0-DMD) 01/25/02 17:17:55
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      SUBROUTINE CONTS3(A,MX,NX,NY,H,D,K1G,M1)
C++++++++++++
C     ROUTINE FOR PLOTTING CONTOURS OF A TABULETED FUNCTION
C
C     A = INPUT MATRIX OF DIMENSION A(MX,NY)
C     NX= NUMBER OF POINTS USED IN X-DIRECTION, .LE. MX
C     H = LEVEL OF CONTOURS
C     D =MESH SIZE IN PLOTTER UNITS
C-----------------
      IMPLICIT  NONE
      REAL      A(1)
      INTEGER   MX
      INTEGER   NX
      INTEGER   NY
      REAL      H
      REAL      D
      INTEGER   K1G
      INTEGER   M1
C
      INCLUDE   'MID_INCLUDE:IDIDEV.INC'
      INCLUDE   'MID_INCLUDE:IDIMEM.INC'
C
      INTEGER   I, IN, INEXT, IST, IX, ITM, IKX
      INTEGER   NP, N1
      INTEGER   J, JFOR, JNB, JM, JX, JM1, JABS
      INTEGER   L, LDA, LDB
      INTEGER   K, KX
      INTEGER   JA(3,2)
      INTEGER   MXR
      INTEGER   LX(1000)
      INTEGER   IX2(1000), IY2(1000)
      INTEGER   ERRCOD,WHITE
      REAL      ADN, AH
      REAL      X, Y
      REAL      PX
      CHARACTER STRING*80
C
      DATA      PX/0.5/,WHITE/2/
  999 FORMAT('*** WARNING: Error in CONTS3 - PL',I5)
C
      IX  = 0
      MXR = MX*(NY-1)
      N1  = NX-1
C
C *** tabulating mesh points . . . begin
      DO 100 JX = 2,N1
         DO 200 JM = JX,MXR,MX
            J  = JM+MX
            AH = A(J)-H
            IF (AH) 60,60,61
 60         IF (A(JM)-H) 62,62,1
 61         IF (A(JM)-H)  1,1,63
  1         IX     =  IX+1
            LX(IX) = -J
            IF (AH) 62,62,63
 62         IF (A(J-1)-H) 200,200,4
 63         IF (A(J-1)-H) 4,4,200
  4          IX     =  IX+1
            LX(IX) = J
200      CONTINUE
100   CONTINUE
C
      DO 300 JM = NX,MXR,MX
         J  = JM+MX
         AH = A(J)-H
         IF (AH) 64,64,65
 64      IF (A(J-1)-H) 66,66,76
 65      IF (A(J-1)-H) 76,76,67
 76      IX     = IX+1
         LX(IX) = J
         IF (AH) 66,66,67
 66      IF (A(JM)-H) 300,300,78 
 67      IF (A(JM)-H) 78,78,300
 78      IX     =IX+1
         LX(IX) = -J
300   CONTINUE
C
      DO 400 JM=1,MXR,MX
         J = JM+MX
         IF (A(J)-H) 68,68,69
 68      IF (A(JM)-H) 400,400,71
 69      IF (A(JM)-H) 71,71,400
 71      IX     = IX+1
         LX(IX) = -J
400   CONTINUE
C
      DO 500 J=2,NX
         IF (A(J)-H) 58,58,59
 58      IF (A(J-1)-H) 500,500,72
 59      IF (A(J-1)-H) 72,72,500
 72      IX     = IX+1
         LX(IX) = J
500   CONTINUE
C
C *** tabulating mesh points . . . end
      IF (IX) 50,50,8
  8   IN=IX
C 
C *** start on a new line
      NP   = 0
      JFOR = 0
      LDA  = 1
      LDB  = 2
 30   JX   = LX(IN)
C
C *** interploation to neigbouring point
      IF (JX) 21,22,22
 21   JABS = -JX
      JNB  = JABS-MX
      GO TO 23
 22   JABS = JX
      JNB  = JABS-1
 23   ADN=A(JABS)-A(JNB)
      IF (ADN) 24,9,24
 24   PX   = (A(JABS)-H)/ADN
  9   KX   = (JABS-1)/MX
      IKX  = JABS-MX*KX-1
      IF (JX)25,26,26
C
C *** evaluation of coordinates
 25   X    = (IKX)*D
      Y    = D*(KX-PX)+K1G
      GO TO 27
 26   X    = D*(IKX-PX)
      Y    = D*KX+K1G
 27   CONTINUE
C
C *** plotting is done here
      NP      = NP+1
      IX2(NP) = X
      IY2(NP) = Y
C
C *** start looking for next point
      ITM     = 1
      JM1     = MX
      IF (JX) 10,10,11
  10  JM1     = -1
  11  JA(1,1) = JX+JM1
      JA(2,1) = -JA(1,1)
      JA(3,1) = -JX+1-MX
      JA(3,2) = -JX
      JA(1,2) = JX-JM1
      JA(2,2) = JA(3,1)+JM1
      IF (KX) 14,14,39
  39  IF (IKX) 14,14,36
  36  IF (IKX+1-NX) 35,37,37
  37  IF (JX)   38,38,35
  35  IF (JFOR) 28,29,28
  28  DO 600 I=1,3
         IF (JFOR-JA(I,2)) 600,14,600
 600  CONTINUE
C
  38  LDA = 2
      GO TO 15
  14  LDA =1
  15  LDB = LDA
  29  DO 700 K=1,3
         DO 800 L=LDA,LDB
            DO 900 I=1,IX
               IF (LX(I)-JA(K,L)) 900,17,900
 900        CONTINUE
            GO TO 800
  17        CONTINUE
            ITM=ITM+1
            INEXT=I
            IF (JFOR) 19,33,19
  33        CONTINUE
            GO TO (800,800,20), ITM
 800     CONTINUE
 700  CONTINUE
C
  19  CONTINUE
      LX(IN) = 0
      GO TO (6,20,20),ITM
  20  CONTINUE
      JFOR   = JX
      IN     = INEXT
      GO TO 30
C
C *** end of a line. look for another
   6  CONTINUE
      CALL IIGPLY(QDSPNO,QOVCH, IY2, IX2, NP, WHITE, 1, ERRCOD)
      IF (ERRCOD .NE. 0) THEN
         WRITE (STRING,999) ERRCOD
         CALL STTPUT(STRING,IST)
      END IF
C   
  31  IF (LX(IX)) 8,7,8
   7  IX = IX-1
      IF (IX) 51,51,31
  51  CONTINUE
C
C *** raise pen before exit from routine
  50  RETURN
      END