	PROGRAM SOUND
C
C   ***   This program accepts input from the sounding contained   ***
C   ***      in the file <<sounding>> and calculates various       ***
C   ***       properties of samples of air raised or lowered       ***
C   ***                   to different levels.                     ***
C   ***      Output is tabulated in the file <<calcsound.out>>     ***
C   ***                 and in NCAR graphics files                 ***
C
	PARAMETER (NA=200,NK=20,NS=71)
	REAL T(NA), P(NA), R(NA), TLR(NK,NA), TLP(NK,NA)
	REAL EV(NA), ES(NA), TVRDIF(NK,NA), TVPDIF(NK,NA)
	REAL TLVR(NK,NA), TLVP(NK,NA), LW(NK,NA), TVD(NA)
	REAL CAPER(NS,NK), CAPEP(NS,NK), DCAPE(NS,NK), PAR(NS,NK)
	REAL PAP(NS,NK), NAR(NS,NK), NAP(NS,NK), PL(NK)
	REAL ATVRD(NS,NK,NA),ATVPD(NS,NK,NA)
	REAL TRDBAR(NK,NA),TPDBAR(NK,NA),STDR(NK,NA),STDP(NK,NA)
	INTEGER JJIR(NK,NA),JJIP(NK,NA)
	CHARACTER*12 AA
	CHARACTER*12 PLLAB,PJLAB,TRLAB,TPLAB,SRLAB,SPLAB
	OPEN(UNIT=12, FILE='cape.out', STATUS='NEW')
C
C   ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
C
        CPD=1005.7
	CPV=1870.0
        CL=4190.0
        CPVMCL=2320.0
        RV=461.5
        RD=287.04
        EPS=RD/RV
        ALV0=2.501E6
	PLLAB='Origin Pressure (mb)'
	PJLAB='p (mb)'
	TRLAB='Tv_rev diff (K)'
	TPLAB='Tp_pse diff (K)'
	SRLAB='Std. Dev_rev (K)'
	SPLAB='Std. Dev_pse (K)'
	P1BAR=0.0
	NPC=0
	NMAX=0
	DO 8 I=1,NK
	 DO 7 J=1,NA
	  JJIR(I,J)=0
	  JJIP(I,J)=0
	  TRDBAR(I,J)=0.0
	  TPDBAR(I,J)=0.0
	  STDR(I,J)=0.0
	  STDP(I,J)=0.0
	  DO 6 K=1,NS
	   ATVRD(K,I,J)=0.0
	   ATVPD(K,I,J)=0.0
    6     CONTINUE
    7	 CONTINUE
    8   CONTINUE
	OPEN(UNIT=15,FILE='list',STATUS='OLD')
	DO 900 JJ=1,NS
	READ(15,10)AA
   10	FORMAT(A12)
C
C   ***   Read in the sounding from the file <<sounding>>          ***
c
	OPEN(UNIT=11,FILE=AA,STATUS='OLD')	
	DO 12 J=1,15
	 READ(11,14)
   14	 FORMAT(5X)
   12	CONTINUE
	KK=182
	I=0
	DO 20 KI=1,KK
	 READ(11,15,ERR=21,END=21)PK,TK,RK
   15	 FORMAT(7X,F6.1,1X,F5.1,1X,F5.1)
         IF(KI.GT.1.AND.PK.GT.1000.01)GOTO 20
         IF(PK.GT.98.0)THEN
          I=I+1
          P(I)=PK
          T(I)=TK
          R(I)=RK
	  IF(T(I).GT.100.0.OR.T(I).LT.-200.0)T(I)=-200.0
	  IF(R(I).GT.100.0)R(I)=-60.0
          EV(I)=6.112*EXP(17.67*R(I)/(243.5+R(I)))
          ES(I)=6.112*EXP(17.67*T(I)/(243.5+T(I)))
          IF(I.NE.1)THEN
C	   IF(I.EQ.800)THEN
           RH=EV(I)/ES(I)
           IF(RH.LE.0.9)THEN
            RHN=RH+0.14286*RH-0.10204*RH*RH
           ELSE
            RHN=RH+0.4592*(1.-RH)
           END IF
           EV(I)=RHN*ES(I)
          END IF
          R(I)=0.622*EV(I)/(P(I)-EV(I))
          T(I)=T(I)+273.15
         END IF
   20	CONTINUE
   21	CONTINUE
	N=I
	CLOSE(UNIT=11,DISPOSE='SAVE')
	IF(P(1).LT.1000.0.OR.T(1).LT.293.0.OR.T(1).GT.306.0.OR.
     1    R(1).LT.0.01.OR.N.LT.160.OR.P(1).GT.1025.0)GOTO 900
        NMAX=MAX(N,NMAX)
C
C   ***  Begin outer loop, which cycles through parcel origin levels I ***
C  
	DO 500 I=1,NK
C
C   ***  Define various conserved parcel quantities: reversible   ***
C   ***        entropy, S, pseudo-adiabatic entropy, SP,          *** 
C   ***                   and enthalpy, AH                        ***
C
	PL(I)=P(I)
	RS=EPS*ES(I)/(P(I)-ES(I))
	ALV=ALV0-CPVMCL*(T(I)-273.15)
	EM=MAX(EV(I),1.0E-6) 
	S=(CPD+R(I)*CL)*LOG(T(I))-RD*LOG(P(I)-EV(I))+
     1    ALV*R(I)/T(I)-R(I)*RV*LOG(EM/ES(I))
	SP=CPD*LOG(T(I))-RD*LOG(P(I)-EV(I))+
     1    ALV*R(I)/T(I)-R(I)*RV*LOG(EM/ES(I))
	AH=(CPD+R(I)*CL)*T(I)+ALV*R(I)
C      
C   ***  Find the temperature and mixing ratio of the parcel at   ***
C   ***    level I saturated by a wet bulb process                ***
C
	SLOPE=CPD+ALV*ALV*RS/(RV*T(I)*T(I))
	TG=T(I)
	RG=RS  
	DO 100 J=1,20 
	 ALV1=ALV0-CPVMCL*(TG-273.15)
	 AHG=(CPD+CL*RG)*TG+ALV1*RG
	 TG=TG+(AH-AHG)/SLOPE
	 TC=TG-273.15
	 ENEW=6.112*EXP(17.67*TC/(243.5+TC))
	 RG=EPS*ENEW/(P(I)-ENEW)
  100	CONTINUE
C   
C   ***  Calculate conserved variable at top of downdraft   ***
C
	EG=RG*P(I)/(EPS+RG)
	SPD=CPD*LOG(TG)-RD*LOG(P(I)-EG)+
     1    ALV1*RG/TG
	TVD(I)=TG*(1.+RG/EPS)/(1.+RG)-T(I)*(1.+R(I)/EPS)/
     1    (1.+R(I))
	IF(P(I).LT.100.0)TVD(I)=0.0
	RGD0=RG
	TGD0=TG
C
C   ***   Find lifted condensation pressure     ***
C
	RH=R(I)/RS
	RH=MIN(RH,1.0)
	CHI=T(I)/(1669.0-122.0*RH-T(I))
	PLCL=1.0
	IF(RH.GT.0.0)THEN
	 PLCL=P(I)*(RH**CHI)
	END IF
C
C   ***  Begin updraft loop   ***
C
	SUM=0.0
	RG0=R(I)
	TG0=T(I)
	DO 200 J=I,N
C
C   ***  Calculate estimates of the rates of change of the entropies  ***
C   ***           with temperature at constant pressure               ***
C  
	 RS=EPS*ES(J)/(P(J)-ES(J))
	 ALV=ALV0-CPVMCL*(T(J)-273.15)
	 SL=(CPD+R(I)*CL+ALV*ALV*RS/(RV*T(J)*T(J)))/T(J)
	 SLP=(CPD+RS*CL+ALV*ALV*RS/(RV*T(J)*T(J)))/T(J)
C   
C   ***  Calculate lifted parcel temperature below its LCL   ***
C
	 IF(P(J).GE.PLCL)THEN
	  TLR(I,J)=T(I)*(P(J)/P(I))**(RD/CPD)
	  TLP(I,J)=TLR(I,J) 
	  LW(I,J)=0.0
	  TLVR(I,J)=TLR(I,J)*(1.+R(I)/EPS)/(1.+R(I))
	  TLVP(I,J)=TLVR(I,J)
	  TVRDIF(I,J)=TLVR(I,J)-T(J)*(1.+R(J)/EPS)/(1.+R(J))
	  TVPDIF(I,J)=TVRDIF(I,J)
	 ELSE
C
C   ***  Iteratively calculate lifted parcel temperature and mixing   ***
C   ***    ratios for both reversible and pseudo-adiabatic ascent     ***
C
	 TG=T(J)
	 RG=RS
	 DO 150 K=1,20
	  EM=RG*P(J)/(EPS+RG)
	  ALV=ALV0-CPVMCL*(TG-273.15)
	  SG=(CPD+R(I)*CL)*LOG(TG)-RD*LOG(P(J)-EM)+
     1      ALV*RG/TG
	  TG=TG+(S-SG)/SL  
	  TC=TG-273.15
	  ENEW=6.112*EXP(17.67*TC/(243.5+TC))
	  RG=EPS*ENEW/(P(J)-ENEW)           
  150	 CONTINUE
	 TLR(I,J)=TG
	 TLVR(I,J)=TG*(1.+RG/EPS)/(1.+R(I))
	 LW(I,J)=R(I)-RG
	 LW(I,J)=MAX(0.0,LW(I,J))
	 TVRDIF(I,J)=TLVR(I,J)-T(J)*(1.+R(J)/EPS)/(1.+R(J))
	 IF(TVRDIF(I,J).GE.7.0)GOTO 900
C
C   ***   Now do pseudo-adiabatic ascent   ***
C
	 TG=T(J)
	 RG=RS
	 DO 180 K=1,20 
	  CPW=0.0
	  IF(J.GT.1)THEN
	   CPW=SUM+CL*0.5*(RG0+RG)*(LOG(TG)-LOG(TG0))
	  END IF
	  EM=RG*P(J)/(EPS+RG)
	  ALV=ALV0-CPVMCL*(TG-273.15)
	  SPG=CPD*LOG(TG)-RD*LOG(P(J)-EM)+CPW+
     1      ALV*RG/TG
	  TG=TG+(SP-SPG)/SLP  
	  TC=TG-273.15
	  ENEW=6.112*EXP(17.67*TC/(243.5+TC))
	  RG=EPS*ENEW/(P(J)-ENEW)           
  180	 CONTINUE
	 TLP(I,J)=TG
	 TLVP(I,J)=TG*(1.+RG/EPS)/(1.+RG)
	 TVPDIF(I,J)=TLVP(I,J)-T(J)*(1.+R(J)/EPS)/(1.+R(J))
	 RG0=RG
	 TG0=TG
	 SUM=CPW
 	END IF
  200	CONTINUE
	IF(I.EQ.1)GOTO 500
C
C   ***  Begin downdraft loop   ***
C
	SUM2=0.0
	DO 300 J=I-1,1,-1
C
C   ***  Calculate estimate of the rate of change of entropy          ***
C   ***           with temperature at constant pressure               ***
C  
	 RS=EPS*ES(J)/(P(J)-ES(J))
	 ALV=ALV0-CPVMCL*(T(J)-273.15)
	 SLP=(CPD+RS*CL+ALV*ALV*RS/(RV*T(J)*T(J)))/T(J)
	 TG=T(J)
	 RG=RS
C
C   ***  Do iteration to find downdraft temperature   ***
C
	 DO 250 K=1,20
	  CPW=SUM2+CL*0.5*(RGD0+RG)*(LOG(TG)-LOG(TGD0))
	  EM=RG*P(J)/(EPS+RG)
	  ALV=ALV0-CPVMCL*(TG-273.15)
	  SPG=CPD*LOG(TG)-RD*LOG(P(J)-EM)+CPW+
     1      ALV*RG/TG
	  TG=TG+(SPD-SPG)/SLP  
	  TC=TG-273.15
	  ENEW=6.112*EXP(17.67*TC/(243.5+TC))
	  RG=EPS*ENEW/(P(J)-ENEW)           
  250	 CONTINUE
	 SUM2=CPW
	 TGD0=TG
	 RGD0=RG
	 TLP(I,J)=TG
	 TLVP(I,J)=TG*(1.+RG/EPS)/(1.+RG)
	 TVPDIF(I,J)=TLVP(I,J)-T(J)*(1.+R(J)/EPS)/(1.+R(J))
	 IF(P(I).LT.100.0)TVPDIF(I,J)=0.0
	 TVPDIF(I,J)=MIN(TVPDIF(I,J),0.0)
	 TLR(I,J)=T(J)
	 TLVR(I,J)=T(J)
	 TVRDIF(I,J)=0.0
	 LW(I,J)=0.0
  300	CONTINUE
  500   CONTINUE
C
C  ***  Begin loop to find CAPE, PA, and NA from reversible and ***
C  ***            pseudo-adiabatic ascent, and DCAPE            ***
C
	DO 800 I=1,NK
	 CAPER(JJ,I)=0.0
	 CAPEP(JJ,I)=0.0
	 DCAPE(JJ,I)=0.0
	 PAP(JJ,I)=0.0
	 PAR(JJ,I)=0.0
	 NAP(JJ,I)=0.0
	 NAR(JJ,I)=0.0
C
C   ***   Find lifted condensation pressure     ***
C
	RS=EPS*ES(I)/(P(I)-ES(I))
	RH=R(I)/RS
	RH=MIN(RH,1.0)
	CHI=T(I)/(1669.0-122.0*RH-T(I))
	PLCL=1.0
	IF(RH.GT.0.0)THEN
	 PLCL=P(I)*(RH**CHI)
	END IF
C
C   ***  Find lifted condensation level and maximum level   ***
C   ***               of positive buoyancy                  ***
C
	 ICB=N
	 INBR=1
	 INBP=1
	 DO 550 J=N,I,-1
	  IF(P(J).LT.PLCL)ICB=MIN(ICB,J)
	  IF(TVRDIF(I,J).GT.0.0)INBR=MAX(INBR,J)
	  IF(TVPDIF(I,J).GT.0.0)INBP=MAX(INBP,J)
  550	 CONTINUE
	  IMAX=MAX(INBR,I)
	   DO 555 J=IMAX,N
	    TVRDIF(I,J)=0.0
  555	   CONTINUE
	  IMAX=MAX(INBP,I)
	   DO 565 J=IMAX,N
	    TVPDIF(I,J)=0.0
  565	   CONTINUE
C
C   ***  Do updraft loops        ***
C
	 IF(INBR.GT.I)THEN
	  DO 600 J=I+1,INBR
	   TVM=0.5*(TVRDIF(I,J)+TVRDIF(I,J-1))
	   PM=0.5*(P(J)+P(J-1))
	   IF(TVM.LE.0.0)THEN
	    NAR(JJ,I)=NAR(JJ,I)-RD*TVM*(P(J-1)-P(J))/PM
	   ELSE
	    PAR(JJ,I)=PAR(JJ,I)+RD*TVM*(P(J-1)-P(J))/PM
	   END IF
  600	  CONTINUE
	  CAPER(JJ,I)=PAR(JJ,I)-NAR(JJ,I)
	 END IF
	 IF(INBP.GT.I)THEN
	  DO 650 J=I+1,INBP
	   TVM=0.5*(TVPDIF(I,J)+TVPDIF(I,J-1))
	   PM=0.5*(P(J)+P(J-1))
	   IF(TVM.LE.0.0)THEN
	    NAP(JJ,I)=NAP(JJ,I)-RD*TVM*(P(J-1)-P(J))/PM
	   ELSE
	    PAP(JJ,I)=PAP(JJ,I)+RD*TVM*(P(J-1)-P(J))/PM
	   END IF
  650	  CONTINUE
	  CAPEP(JJ,I)=PAP(JJ,I)-NAP(JJ,I)
	 END IF
C
C  ***       Find DCAPE     ***
C
	 IF(I.EQ.1)GOTO 800
	 DO 700 J=I-1,1,-1
	  TVDIFM=TVPDIF(I,J+1)
	  IF(I.EQ.(J+1))TVDIFM=TVD(I)
	  TVM=0.5*(TVPDIF(I,J)+TVDIFM)
	  PM=0.5*(P(J)+P(J+1))
	  IF(TVM.LT.0.0)THEN
	   DCAPE(JJ,I)=DCAPE(JJ,I)-RD*TVM*(P(J)-P(J+1))/PM
	  END IF
  700	 CONTINUE	  	
  800	CONTINUE
	DO 520 I=1,NK
	 INBR=1
	 INBP=1
	 DO 505 J=N,I,-1
	  IF(TVRDIF(I,J).GT.0.0)INBR=MAX(INBR,J)
	  IF(TVPDIF(I,J).GT.0.0)INBP=MAX(INBP,J)
  505	 CONTINUE
	  IMAX=MAX(INBR,I)
	IF(IMAX.GT.I)THEN
	DO 508 J=I,IMAX
	 JJIR(I,J)=JJIR(I,J)+1
	 ATVRD(JJIR(I,J),I,J)=TVRDIF(I,J)
  508	CONTINUE
	END IF
	IMAX=MAX(I,INBP)
	DO 510 J=I,IMAX
	 IF(TVPDIF(I,J).NE.0.0)THEN
	 JJIP(I,J)=JJIP(I,J)+1
	 ATVPD(JJIP(I,J),I,J)=TVPDIF(I,J)
	 END IF
  510	CONTINUE
  520	CONTINUE	
	NPC=NPC+1
	P1BAR=P1BAR+P(1)
  900	CONTINUE
	CLOSE(UNIT=15,DISPOSE='SAVE')
	P1BAR=P1BAR/FLOAT(NPC)
	P(1)=P1BAR
	PL(1)=P1BAR
	DO 920 I=1,NK     
	PARS=0.0
	PAPS=0.0
	ANAPS=0.0
	ANARS=0.0
	CAPERS=0.0
	CAPEPS=0.0
	DCAPES=0.0
	JPAR=0
	JPAP=0
	JNAR=0
	JNAP=0
	JCAPER=0
	JCAPEP=0
	JDCAPE=0
	DO 903 JJ=1,NS
	 IF(PAP(JJ,I).NE.0.0)THEN
	  PAPS=PAPS+PAP(JJ,I)
	  JPAP=JPAP+1
	 END IF
	 IF(PAR(JJ,I).NE.0.0)THEN
	  PARS=PARS+PAR(JJ,I)
	  JPAR=JPAR+1
	 END IF
	 IF(NAP(JJ,I).NE.0.0)THEN
	  ANAPS=ANAPS+NAP(JJ,I)
	  JNAP=JNAP+1
	 END IF
	 IF(NAR(JJ,I).NE.0.0)THEN
	  ANARS=ANARS+NAR(JJ,I)
	  JNAR=JNAR+1
	 END IF
	 IF(CAPER(JJ,I).NE.0.0)THEN
	  CAPERS=CAPERS+CAPER(JJ,I)
	  JCAPER=JCAPER+1
	 END IF
	 IF(CAPEP(JJ,I).NE.0.0)THEN
	  CAPEPS=CAPEPS+CAPEP(JJ,I)
	  JCAPEP=JCAPEP+1
	 END IF
	 IF(DCAPE(JJ,I).NE.0.0)THEN
	  DCAPES=DCAPES+DCAPE(JJ,I)
	  JDCAPE=JDCAPE+1
	 END IF
  903	CONTINUE 
	JNAP=MAX(JNAP,1)
	JNAR=MAX(JNAR,1)
	JPAP=MAX(JPAP,1)
	JPAR=MAX(JPAR,1)
	JCAPER=MAX(JCAPER,1)
	JCAPEP=MAX(JCAPEP,1)
	JDCAPE=MAX(JDCAPE,1)
	PAP(1,I)=PAPS/FLOAT(JPAP)
	PAR(1,I)=PARS/FLOAT(JPAR)
	NAP(1,I)=ANAPS/FLOAT(JNAP)
	NAR(1,I)=ANARS/FLOAT(JNAR)
	CAPEP(1,I)=CAPEPS/FLOAT(JCAPEP)
	CAPER(1,I)=CAPERS/FLOAT(JCAPER)
	DCAPE(1,I)=DCAPES/FLOAT(JDCAPE)
	DO 915 J=1,NMAX
	 IF(JJIR(I,J).GT.0)THEN
	 DO 905 JJ=1,JJIR(I,J)
	  TRDBAR(I,J)=TRDBAR(I,J)+ATVRD(JJ,I,J)
  905	 CONTINUE
	 TRDBAR(I,J)=TRDBAR(I,J)/FLOAT(JJIR(I,J))
	 DO 907 JJ=1,JJIR(I,J)
	  STDR(I,J)=STDR(I,J)+(ATVRD(JJ,I,J)-TRDBAR(I,J))**2
  907	 CONTINUE
	 STDR(I,J)=SQRT(STDR(I,J)/FLOAT(JJIR(I,J)))
	 END IF
	 IF(JJIP(I,J).GT.0)THEN
	 DO 908 JJ=1,JJIP(I,J)
	  TPDBAR(I,J)=TPDBAR(I,J)+ATVPD(JJ,I,J)
  908	 CONTINUE
	 TPDBAR(I,J)=TPDBAR(I,J)/FLOAT(JJIP(I,J))
	 DO 909 JJ=1,JJIP(I,J)
	  STDP(I,J)=STDP(I,J)+(ATVPD(JJ,I,J)-TPDBAR(I,J))**2
  909	 CONTINUE
	 STDP(I,J)=SQRT(STDP(I,J)/FLOAT(JJIP(I,J)))
	 END IF
  915	CONTINUE
  920	CONTINUE
C
C  ***  Write values of PA, NA, CAPE and DCAPE   ***
C
	WRITE(12,921)NPC
  921	FORMAT(5X,'Number of soundings = ',I4)
	WRITE(12,801)
  801   FORMAT(20X,'ALL AREAS IN UNITS OF J/kg',//)
	WRITE(12,802)
  802	FORMAT(1X,'Origin',4X,'Rev.',4X,'P.A.',4X,'Rev.',4X,'P.A.',
     1   4X,'Rev.',4X,'P.A.')
	WRITE(12,804)
  804	FORMAT(1X,'p (mb)',4X,' PA ',4X,' PA ',4X,' NA ',4X,' NA ',
     1   4X,'CAPE',4X,'CAPE',4X,'DCAPE')
	WRITE(12,806)
  806	FORMAT(1X,'------',4X,'----',4X,'----',4X,'----',4X,'----',
     1   4X,'----',4X,'----',4X,'-----')
	DO 820 I=1,NK
	 WRITE(12,810)PL(I),PAR(1,I),PAP(1,I),NAR(1,I),
     1     NAP(1,I),CAPER(1,I),CAPEP(1,I),DCAPE(1,I)
  810	 FORMAT(1X,F6.1,7(F8.1))
  820	CONTINUE
C
C   ***  Use NCAR graphics to contour arrays of parcel buoyancy   ***
C   ***    (Delete these if NCAR graphics are not available)      ***
C
	XS=1.0
	YS=0.865
	NC=12
	NM=NMAX
	CALL OPNGKS
	CALL PLREG(PL,P,TRDBAR,NK,NA,NK,NM,NC,XS,YS,PLLAB,PJLAB,TRLAB)
	CALL PLREG(PL,P,TPDBAR,NK,NA,NK,NM,NC,XS,YS,PLLAB,PJLAB,TPLAB)
	NC=8
	CALL PLREG(PL,P,STDR,NK,NA,NK,NM,NC,XS,YS,PLLAB,PJLAB,SRLAB)
	CALL PLREG(PL,P,STDP,NK,NA,NK,NM,NC,XS,YS,PLLAB,PJLAB,SPLAB)
	CALL CLSGKS
	STOP
	END
C
	SUBROUTINE PLREG(X,Y,Z,NA,NB,NX,NY,NC,XSCA,YSCA,XLAB,YLAB,ZLAB)
c
c   Subroutine PLREG contours the array Z(NA,NB) as a function of X(NA) and
c   Y(NB), using NX points in the X direction and NY points in the Y 
c   direction. The aspect ratio of the plot can be set by changing the 
c   scaling of the X axis, XSCA, and/or the Y axis, YSCA. (0<XSCA,YSCA<1).
c   The X and Y axes are labelled by XLAB and YLAB (ideal dimensions
c   CHARACTER*12), and the contoured quantity is called ZLAB.
c         NC is the number of contours.
c
	REAL X(NA), Y(NB), Z(NA,NB)
	CHARACTER*(*) XLAB,YLAB,ZLAB
	CHARACTER*12 CHXMIN,CHXMAX,CHYMIN,CHYMAX
	CHARACTER*50 CHZ
	XSC=MIN(XSCA,1.0)
	XSC=MAX(XSCA,0.2)
	YSC=MIN(YSCA,1.0)
	YSC=MAX(YSCA,0.2)
	XSC=0.8*XSCA
	YSC=0.8*YSCA
	XMAX=X(NX)
	YMAX=Y(NY)
	XMIN=X(1)
	YMIN=Y(1)
	ZMAX=Z(1,1)
	ZMIN=Z(1,1)
	DO 100 I=1,NX
	 DO 100 J=1,NY
	  ZMAX=MAX(ZMAX,Z(I,J))
	  ZMIN=MIN(ZMIN,Z(I,J))
  100	CONTINUE
	DELZ=MAX((ZMAX-ZMIN),1.0E-16) 
	AINC=DELZ/FLOAT(NC)
	WRITE(CHZ,210)ZLAB,ZMIN,ZMAX
  210	FORMAT(1X,A10,'  From  ',F10.4,'  to  ',F10.4)
	WRITE(CHXMIN,220)XMIN
  220	FORMAT(1X,F8.2)
	WRITE(CHXMAX,220)XMAX
	WRITE(CHYMIN,220)YMIN
	WRITE(CHYMAX,220)YMAX
	SXMAX=0.5*(1.+XSC)+0.05
	SXMIN=0.5*(1.-XSC)+0.05
	SXMIN2=SXMIN+0.02
	SYMAX=0.5*(1.+YSC)
	SYMIN=0.5*(1.-YSC)
	SYMAX2=SYMAX+0.05
	SYMIN2=SYMIN-0.05
	CALL GSELNT(0)
        CALL SET(SXMIN,SXMAX,SYMIN,SYMAX,0.0,1.,0.,1.,1)
	CALL PERIM(NX,1,NY,1)
	CALL CONREC(Z,NA,NX,NY,0.,0.,AINC,-1,0,-750)
	CALL SET(0.,1.,0.,1.,0.,1.,0.,1.,1)
	CALL GSELNT(0)
	CALL PLCHHQ(0.55,SYMAX2,CHZ,0.012,0.,0.)
	CALL PLCHHQ(SXMIN2,SYMIN,CHYMIN,0.012,0.,1.)
	CALL PLCHHQ(SXMIN2,0.5,YLAB,0.012,0.,1.)
	CALL PLCHHQ(SXMIN2,SYMAX,CHYMAX,0.012,0.,1.)
	CALL PLCHHQ(SXMIN,SYMIN2,CHXMIN,0.012,0.,0.)
	CALL PLCHHQ(SXMAX,SYMIN2,CHXMAX,0.012,0.,0.)
	CALL PLCHHQ(0.55,SYMIN2,XLAB,0.012,0.,0.)
	CALL FRAME
	RETURN
	END
