C****************************************************************
C*****     PGPLOT making Contour and Vector map for RAMS43x   ***
C*****         THREE DIMENSIONAL VARIABLES VERSION            ***
C*****                 by M.KANDA                             ***
C****                  oroginal version 1.0 2001.7 June       ***
C****************************************************************
                   IMPLICIT REAL(A-H,O-Z)
C------------------------
      CALL VAR_INT
      CALL OPEN_FILE
      CALL RVDATA
      CALL DRIVER
      CALL CLOSE_FILE
C------------------------
      STOP
      END
C-------------------------------------------------------------
C
C***************************************************
      SUBROUTINE VAR_INT
C***************************************************
                    IMPLICIT REAL(A-H,O-Z)
      include 'ramspg.h'
C--------------------------------------------------
C     X and Y scale ( should be referred to RAMSIN )
C--------------------------------------------------
      XSCALE = 10.0
      YSCALE = 10.0
      WSCALE = 5.0
      IYEAR  = 98
      IMONTH = 8
      IDAY   = 8
C------------------------------
C  IFLAG1 =1 : OUTPUT PS FILES
C  IFLAG1 =0 : DISPLAY XWINDOW
C------------------------------
        IFLAG1 = 1
C
      RETURN
      END
C***************************************************
      SUBROUTINE OPEN_FILE
C***************************************************
                    IMPLICIT REAL(A-H,O-Z)
C----------------------------------------------------------
C-- Open variable files -------
      OPEN (79,FILE='RDATA',STATUS='OLD',FORM='FORMATTED')
C-- Open variable files -------
      OPEN (78,FILE='RTOPOG2',STATUS='OLD',FORM='FORMATTED')
C-----------------------------------------------------------
      RETURN
      END
C***************************************************
      SUBROUTINE CLOSE_FILE
C***************************************************
        IMPLICIT REAL(A-H,O-Z)
C----------------------------------------------------------
C-- Close files.
       CLOSE(79)
       CLOSE(78)
      RETURN
      END
C***************************************************
      SUBROUTINE RVDATA
C***************************************************
        IMPLICIT REAL(A-H,O-Z)
      include 'ramspg.h'
C---------------------------
C  3-dimensional variable
C---------------------------
            DO 100 M= 1,MMAX
       WRITE(*,*) 'M=',M
C--------
C READ U 
C--------
            DO 111 K=1,KMAX
               READ (79,*) 
               READ (79,777)
     +                ZLEVS(K),ITIME(M),ICCC,ICCC
               READ (79,888)
     +              ((U(M,I,J,K),I=1,IMAX),J=1,JMAX)
 111        CONTINUE
C--------
C READ V
C--------
            DO 112 K=1,KMAX
               READ(79,*)
               READ (79,777)
     +                ZLEVS(K),ICCC,ICCC,ICCC
               READ (79,888)
     +              ((V(M,I,J,K),I=1,IMAX),J=1,JMAX)
 112        CONTINUE
C--------
C READ W 
C--------
            DO 113 K=1,KMAX
               READ(79,*)
               READ (79,777)
     +                ZLEVS(K),ICCC,ICCC,ICCC
               READ (79,888)
     +              ((W(M,I,J,K),I=1,IMAX),J=1,JMAX)
 113         CONTINUE
C----------
C READ THT 
C----------
            DO 114 K=1,KMAX
               READ(79,*)
               READ (79,777)
     +                ZLEVS(K),ICCC,ICCC,ICCC
               READ (79,888)
     +              ((TP(M,I,J,K),I=1,IMAX),J=1,JMAX)
 114         CONTINUE
C
C-----------
C READ P1 
C-----------
            DO 115 K=1,KMAX
               READ(79,*)
               READ (79,777)
     +                ZLEVS(K),ICCC,ICCC,ICCC
               READ (79,888)
     +              ((P1(M,I,J,K),I=1,IMAX),J=1,JMAX)
 115        CONTINUE
C
C-----------
C READ P2 
C-----------
C            DO 116 K=1,KMAX
C               READ (79,777)
C     +                CVAR,ZLEVS(K),ICCC,ICCC,ICCC
C               READ (79,888)
C     +              ((P2(M,I,J,K),I=1,IMAX),J=1,JMAX)
C 116        CONTINUE
C
C-----------
C READ P3 
C-----------
C            DO 118 K=1,KMAX
C               READ (79,777) 
C     +                CVAR,ZLEVS(K),ICCC,ICCC,ICCC
C               READ (79,888) 
C     +              ((P3(M,I,J,K),I=1,IMAX),J=1,JMAX)
C 118        CONTINUE   
C
C-----------
C READ P4 
C-----------
C            DO 119 K=1,KMAX
C               READ (79,777) 
C     +                CVAR,ZLEVS(K),ICCC,ICCC,ICCC
C               READ (79,888) 
C     +              ((P4(M,I,J,K),I=1,IMAX),J=1,JMAX)
C 119        CONTINUE   
C
C-----------
C READ P5 
C-----------
C            DO 120 K=1,KMAX
C               READ (79,777) 
C     +                CVAR,ZLEVS(K),ICCC,ICCC,ICCC
C               READ (79,888) 
C     +              ((P5(M,I,J,K),I=1,IMAX),J=1,JMAX)
C 120        CONTINUE         
C
 100               CONTINUE          
C---------------------------------------------------
C READ topo  2D variable
C---------------------------------------------------
      READ (78,*)
      READ (78,777) 
     &                cccc,ICCC,ICCC,ICCC
      READ (78,888) 
     &              ((TOPT(I,J),I=1,IMAX),J=1,JMAX)
      ZLEVS(1) = 0.0
                   DO 333 I = 1,IMAX
                   DO 333 J = 1,JMAX
                   DO 333 K = 1,KMAX
      X(I,J,K) = XSCALE*(I-1)
      Y(I,J,K) = YSCALE*(J-1)
      Z(I,J,K) = (ZLEVS(K)-ZLEVS(1))/(ZLEVS(KMAX)-ZLEVS(1))*
     1             ((ZLEVS(KMAX)-ZLEVS(1))-TOPT(I,J)) + TOPT(I,J)
 333               CONTINUE
C---------------------------------------------------
C     DATA FORMAT
C---------------------------------------------------
  777    FORMAT(F10.2,3I3) 
  888    FORMAT(5E14.6) 
C---------------------------------------------------
      RETURN
      END
C***************************************************
      SUBROUTINE DRIVER
C***************************************************
        IMPLICIT REAL(A-H,O-Z)
      include 'ramspg.h'
      include 'conx.h'
C-------------------------------------- 
 210                         CONTINUE 
C--------------------------------------
C (ipg1,jpg1,kpg1) - (ipg2,jpg2,kpg2)  
C specification of the analyzed region 
C--------------------------------------
      ipg1 = 1
      ipg2 = IMAX 
      jpg1 = 1
      jpg2 = JMAX
      kpg1 = 1
      kpg2 = KMAX 
C
      WRITE(*,*) ' Do you want to specify plotted region ?'
      WRITE(*,*) ' PRESS 0 : Default, all region plotted '
      READ(*,*) ISC
      IF (ISC.NE.0) THEN 
      WRITE(*,*) ' if you chnage horizontal window, check source code'
      WRITE(*,*) ' only kpg2 can be changable here '
      READ(*,*) kpg2
                    ENDIF
C-------------------------------------
C     IPL : PLANE NUMBER
C     XY-PLANE 1
C     YZ-PLANE 2
C     XZ-PLANE 3
C-------------------------------------
      WRITE(*,*) 'SELECT PLANE TYPE'
      WRITE(*,*) 'XY-PLANE : TYPE 1'
      WRITE(*,*) 'YZ-PLANE : TYPE 2'
      WRITE(*,*) 'XZ-PLANE : TYPE 3'
      READ(*,*) IPL
C------------------------
C     XY-PLANE
C------------------------
      IF (IPL.EQ.1) THEN
      WRITE(*,15) kpg2
  15  FORMAT('SELECT SLICE NUMBER 1 to KMAX : KMAX=',I2)
      READ(*,*) ICUT 
      IF ((ICUT.LT.1).OR.(ICUT.GT.Kpg2)) THEN 
          WRITE(*,*) 'WRONG ICUT NUMBER ! PROGRAM STOP'
                      STOP
                                       ENDIF
      IPG = ipg2-ipg1+1
      JPG = jpg2-jpg1+1
      CALL D3_GRAPH
                                        ENDIF
C------------------------
C     YZ-PLANE
C------------------------
      IF (IPL.EQ.2) THEN
      WRITE(*,16) ipg2
  16  FORMAT('SELECT SLICE NUMBER 1 to IMAX : IMAX=',I2)
      READ(*,*) ICUT 
      IF ((ICUT.LT.1).OR.(ICUT.GT.ipg2)) THEN 
          WRITE(*,*) 'WRONG ICUT NUMBER ! PROGRAM STOP'
                      STOP
                                       ENDIF
      IPG = jpg2-jpg1+1
      JPG = kpg2-kpg1+1
      CALL D3_GRAPH
                    ENDIF
C------------------------
C     XZ-PLANE
C------------------------
      IF (IPL.EQ.3) THEN
      WRITE(*,17) jpg2
  17  FORMAT('SELECT SLICE NUMBER 1 to JMAX : JMAX=',I2)
      READ(*,*) ICUT 
      IF ((ICUT.LT.1).OR.(ICUT.GT.jpg2)) THEN 
          WRITE(*,*) 'WRONG ICUT NUMBER ! PROGRAM STOP'
                      STOP
                                       ENDIF
      IPG = ipg2-ipg1+1
      JPG = kpg2-kpg1+1
      CALL D3_GRAPH
                    ENDIF
C-------------------------------------------------------
      WRITE(*,*) 'IF YOU MAKE ANOTHER GRAPH PRESS 9'
      READ(*,*)  ICON1
      IF (ICON1.EQ.9) GOTO 210 
C-------------------------------------------------------
      RETURN
      END
C********************************************************
      SUBROUTINE D3_GRAPH
C********************************************************
                   IMPLICIT REAL(A-H,O-Z)
      include 'ramspg.h'
      include 'conx.h'
      INTEGER  pgopen
      CHARACTER filename1*17,PP*2
      DIMENSION PGU(MMAX,IIMAX,IIMAX),PGV(MMAX,IIMAX,IIMAX)
     &         ,PGP(M1MAX,MMAX,IIMAX,IIMAX)
      DIMENSION CPP(1000),TR(6)
      EXTERNAL PLOT
C------------------------------------------------
C     PROJECTION 3D-DATA TO 2D-Graphic Plane
C------------------------------------------------
                   DO 390 I = ipg1,ipg2
                   DO 390 J = jpg1,jpg2
      IX = I-ipg1+1
      JX = J-jpg1+1
      PGZ(IX,JX) = TOPT(I,J)
 390               CONTINUE
      IF (IPL.EQ.1) THEN
                   DO 400 I = ipg1,ipg2
                   DO 400 J = jpg1,jpg2
      IX = I-ipg1+1
      JX = J-jpg1+1
      PGX(IX,JX) = X(I,J,ICUT)
      PGY(IX,JX) = Y(I,J,ICUT)
 400               CONTINUE
                   DO 401 M = 1,MMAX
                   DO 401 I = ipg1,ipg2
                   DO 401 J = jpg1,jpg2
      IX = I-ipg1+1
      JX = J-jpg1+1
      PGU(M,IX,JX)   = U(M,I,J,ICUT)
      PGV(M,IX,JX)   = V(M,I,J,ICUT)
      PGP(1,M,IX,JX) = U(M,I,J,ICUT)
      PGP(2,M,IX,JX) = V(M,I,J,ICUT)
      PGP(3,M,IX,JX) = W(M,I,J,ICUT)
      PGP(4,M,IX,JX) = TP(M,I,J,ICUT)
      PGP(5,M,IX,JX) = P1(M,I,J,ICUT)
 401               CONTINUE
                    ENDIF
C
      IF (IPL.EQ.2) THEN
                   DO 410 J = jpg1,jpg2
                   DO 410 K = kpg1,kpg2
      JX = J-jpg1+1
      KX = K-kpg1+1
      PGX(JX,KX) = Y(ICUT,J,K)
      PGY(JX,KX) = Z(ICUT,J,K)
 410              CONTINUE
                   DO 411 M = 1,MMAX
                   DO 411 J = jpg1,jpg2
                   DO 411 K = kpg1,kpg2
      JX = J-jpg1+1
      KX = K-kpg1+1
      PGU(M,JX,KX)   = V(M,ICUT,J,K)
      PGV(M,JX,KX)   = W(M,ICUT,J,K)*WSCALE
      PGP(1,M,JX,KX) = U(M,ICUT,J,K)
      PGP(2,M,JX,KX) = V(M,ICUT,J,K)
      PGP(3,M,JX,KX) = W(M,ICUT,J,K)
      PGP(4,M,JX,KX) = TP(M,ICUT,J,K)
      PGP(5,M,JX,KX) = P1(M,ICUT,J,K)
 411               CONTINUE
                    ENDIF
C
      IF (IPL.EQ.3) THEN
                   DO 420 I = ipg1,ipg2
                   DO 420 K = kpg1,kpg2
      IX = I-ipg1+1
      KX = K-kpg1+1
      PGX(IX,KX) = X(I,ICUT,K)
      PGY(IX,KX) = Z(I,ICUT,K)
 420               CONTINUE
                   DO 421 M = 1,MMAX
                   DO 421 I = ipg1,ipg2
                   DO 421 K = kpg1,kpg2
      IX = I-ipg1+1
      KX = K-kpg1+1
      PGU(M,IX,KX)   = U(M,I,ICUT,K)
      PGV(M,IX,KX)   = W(M,I,ICUT,K)*WSCALE
      PGP(1,M,IX,KX) = U(M,I,ICUT,K)
      PGP(2,M,IX,KX) = V(M,I,ICUT,K)
      PGP(3,M,IX,KX) = W(M,I,ICUT,K)
      PGP(4,M,IX,KX) = TP(M,I,ICUT,K)
      PGP(5,M,IX,KX) = P1(M,I,ICUT,K)
 421               CONTINUE
                    ENDIF
C------------------------------------------------
C     Graphic type selector
C------------------------------------------------
      WRITE(*,*) 'SELECT Graph Type '
      WRITE(*,*) 'Vector  : 0 '
      READ(*,*)  INUM0
      WRITE(*,*) 'Contour : Variable Number'
      WRITE(*,*) '0:NONE 1:U  2:V  3:W  4:THT 5:VAP'
      READ(*,*)  INUM1 
      WRITE(*,*) ' DRAWING COSTAL LINE ? Y:1 '
      READ(*,*) ICOS
C------------------------------------------------
C  Searching World-Coordinate range
C------------------------------------------------
       XMIN = 0.0
       YMIN = 0.0
       XMAX = 0.0
       YMAX = 0.0
               DO 510 I = 1,IPG
               DO 510 J = 1,JPG
       XMAX = MAX(XMAX,PGX(I,J))
       YMAX = MAX(YMAX,PGY(I,J))
 510           CONTINUE
       PXX = XMAX/20.0
       PYY = YMAX/20.0
       XMIN = 0 - PXX
       YMIN = 0 - PYY
       XMAX = XMAX + PXX
       YMAX = YMAX + PYY

C--------------------------------------------------------
C--------------------------------------------------------
C ----------------------
C Serching Contour MAX
C ----------------------
       PMIN = 1000000.0
       PMAX = -1000000.0
               DO 1000 M = 1,MMAX
               DO 520  I = 1,IPG
               DO 520  J = 1,JPG
       PG1(I,J) = PGU(M,I,J)
       PG2(I,J) = PGV(M,I,J)
       PG3(I,J) = 0.0
       IF (INUM1.NE.0) THEN 
       PG3(I,J) = PGP(INUM1,M,I,J)
       PMIN = MIN(PMIN,PG3(I,J))
       PMAX = MAX(PMAX,PG3(I,J))
                       ENDIF
 520                     CONTINUE
 1000                    CONTINUE
C-----------------------
C REFERENCE FOR CONTOUR
C-----------------------
               IF (INUM1.NE.0) THEN
       WRITE(*,*) 'PMAX=',PMIN
       WRITE(*,*) 'PMAX=',PMAX
       WRITE(*,*) 'READ PMIN ='
       READ(*,*) PMIN
       WRITE(*,*) 'READ DP ='
       READ(*,*) DP
                               
       NCP = 10
                    DO 555 II = 1, NCP
       CPP(II) = PMIN + DP*(II-1)
  555               CONTINUE
                              ENDIF
C COASTAL LINE 
       CPC = 0.01
C
C-----------------------
C REFERENCE FOR VECTOR
C-----------------------
C
      NCV = 1
      CV  = Xscale/2.5
C 
C      TR(1) = XMIN
C      TR(2) = XMAX/(1.0*IPG)
C      TR(3) = 0.0
C      TR(4) = YMIN
C      TR(5) = 0.0
C      TR(6) = YMAX/(1.0*JPG)
      TR(1) = -XSCALE
      TR(2) = XSCALE
      TR(3) = 0.0
      TR(4) = -YSCALE
      TR(5) = 0.0
      TR(6) = YSCALE
      BLANK = 100.0
C
C-----------------------------------------------------------
C-----------------------------------------------------------
C-----------------------------------------------------------
                         DO 100 M=1,MMAX
C-----------------------------------------------------------
C-----------------------------------------------------------
C-----------------------------------------------------------
C-----------------------------------------------------------
C Temporal Extraporation
C------------------------
                         DO 600 I = 1,IPG
                         DO 600 J = 1,JPG
       PG1(I,J) = PGU(M,I,J)
       PG2(I,J) = PGV(M,I,J)
       PG3(I,J) = 0.0
       IF (INUM1.NE.0) PG3(I,J) = PGP(INUM1,M,I,J)
 600           CONTINUE
C-----------------------------------------------------------
C Difinition of filename
C-----------------------
      IF (IPL.EQ.1) PP='XY'
      IF (IPL.EQ.2) PP='YZ'
      IF (IPL.EQ.3) PP='XZ'
      ICDAY = INT(1.0*(ITIME(M)+9)/24.0)   
      IITIME = (ITIME(M)+9)-24*ICDAY
      IIDAY  = IDAY + ICDAY
      WRITE(*,*) ITIME(M)
      WRITE(filename1,22) IYEAR,IMONTH,IIDAY,IITIME,PP,ICUT
   22 FORMAT(4(i2.2),A2,'-',(i3.3),'.ps')
C-----------------------
C Open device and check
C-----------------------
      IF (IFLAG1.EQ.0) CALL pgopen('/XWINDOW')
      IF (IFLAG1.EQ.1) CALL pgopen(filename1//'/PS')
C-----------------------
C Set Arrow size   
C-----------------------
      CALL PGSCH(1.5)
C ----------------------
C Open viewport
C ----------------------
      CALL pgenv(XMIN, XMAX, YMIN, YMAX, 0, 0)
C-----------------------
C Drawing lines
C-----------------------
      IF (INUM1.NE.0) 
C    & CALL PGCONT(PG3,IPG,JPG,1,IPG,1,JPG,CPP,NCP,TR)
     & CALL PGCONX(PG3,IPG,JPG,1,IPG,1,JPG,CPP,NCP,PLOT)
      IF (ICOS.EQ.1) 
     & CALL PGCONX(PGZ,IPG,JPG,1,IPG,1,JPG,CPC,-1,PLOT)
      IF (INUM0.EQ.0) THEN
       CALL PGSCH(0.2)
       CALL PGSAH(2,20.0,0.3)
       CALL PGVECT(PG1,PG2,IPG,JPG,1,IPG,1,JPG,CV,NCV,TR,BLANK
     &            ,2,PGX,PGY)
                      ENDIF
       CALL PGSCH(1.5)
       CALL PGTEXT(XMIN+10,YMIN+10,filename1)
C ----------------------
C Close device
C ----------------------
      CALL pgclos
C------------------------
C--------------------------------------------
 100                CONTINUE
C--------------------------------------------
      RETURN
      END
C----------------------------------------------------
      SUBROUTINE PLOT (IVISBLE,XX,YY,ZZ)
                   IMPLICIT REAL(A-H,O-Z)
      include 'conx.h'
      INTEGER IVISBLE
C
      II  = INT(XX)
      JJ  = INT(YY)
      XXX = XX -1.0*II
      YYY = YY -1.0*JJ
C----------------------------
      IF (II.EQ.IPG) THEN
      XWORLD = PGX(II,JJ)
                     ELSE
      XWORLD = PGX(II,JJ)+(PGX(II+1,JJ)-PGX(II,JJ))*XXX
                     ENDIF
      IF (JJ.EQ.JPG) THEN
      YWORLD = PGY(II,JJ)
                     ELSE
      YWORLD = PGY(II,JJ)+(PGY(II,JJ+1)-PGY(II,JJ))*YYY
                     ENDIF
C----------------------------
      IF (IVISBLE.EQ.0) THEN
          CALL PGMOVE (XWORLD, YWORLD)
      ELSE
          CALL PGDRAW (XWORLD, YWORLD)
      END IF
      RETURN
      END
C------------------------------------------------------
