      SUBROUTINE LOWLAY                                                         
c These parameter statements are to support a test run of this code.
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT= 100)

      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,          
     +ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),           
     +FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),           
     +AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),        
     +GL(4000),POT(37,73),                                                      
     +GCNK(37,37),GSNK(37,37),TES(26011)                                        
     +  ,GOO,F(37,73)                                                           
     +  ,LMX,ZSIG(37,73)                                                        
     +   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO                                
      COMMON/PHIFT/ PHINF                                                       
      COMMON /CURT/ TCUR,SUMCUP,SUMCDN                                          
      COMMON/TRANS/ DSIGZ                                                       
      COMMON / ANGLE / ANGX , ANGY                                              
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      

      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      EXTERNAL COFZ                                                             
      NLAT=NLATS                                                           
      NLONG=NLONGS                                                          
      NVRT=NVERT                                                                
      NZLV=NZLVL                                                                
      C2=1.E-4*ALOG(10.0)                                                       
      PII=4.0*ATAN(1.0)                                                         
      DLAT=PII/(NLAT-1)                                                         
      DLONG=(PII+PII)/(NLONG-1)                                                 

      DO 19999 J=1,NLAT                                                         
       LATMSH(J)=(J-1)*DLAT                                                      
19999 CONTINUE                                                                  

      DO 19998 K=1,NLONG                                                        
        LNGMSH(K)=(K-1)*DLONG                                                     
19998 CONTINUE                                                                  

	open (unit=8, file= "nzlvl.data")
      READ (8,100) NZLV,(ZLEVEL(I),I=1,NZLV)                                        
  100 FORMAT(I5/(16F5.0))                                                       
	close (unit=8)

      CALL IOUT (NZLV, "NZLV")                                                   
      CALL VECOUT (ZLEVEL,"ZLEVEL",NZLV)                                        
      ZSGMIN=1.E+100                                                            
      ORGMAX=0.0                                                                

      DO 19997 K=1,NLONG                                                        
        DO 19996 J=1,NLAT                                                         
          ZSGMIN=AMIN1(ZSIG (J,K),ZSGMIN)                                           
          ORGMAX=AMAX1(OROG(J,K),ORGMAX)                                            
19996   CONTINUE                                                                  
19997 CONTINUE                                                                  

      IF (.NOT.(ZLEVEL(1).LT.ORGMAX.OR.ZLEVEL(NZLV).GT.ZSGMIN))                 
     +GO TO 19995                                                               

      WRITE (6,200)                                                             
  200 FORMAT('WARNING**, * PRESCRIBED LEVELS OUTSIDE GROUND-ZSIG BN*')        

      CALL VOUT(ZSGMIN, "**ZSGMIN")                                              
      CALL VOUT(ORGMAX, "**ORGMAX")                                              
      CALL VOUT(ZLEVEL(1), "ZLEVELMN")                                           
      CALL VOUT(ZLEVEL(NZLV), "ZLEVELMX")                                        
19995 CONTINUE                                                                  
      DO 19994 K=1,NLONG                                                        
      XLONG=LNGMSH(K)                                                           
      DO 19993 J=1,NLAT                                                         
      XLAT=LATMSH(J)                                                            
      ZBOT = OROG(J,K)                                                          
      ZTOP=ZSIG(J,K)                                                            
      C3=EXP(-C2*ZTOP)-1.0                                                      
      ZB = ZF(ZBOT)                                                             
      ZT = ZF(ZTOP)                                                             
      DLZ = (ZT-ZB)/(NVERT-1)                                                   
      FPRM(1)=ZB                                                                
      FPRM(2)=ZT                                                                
      FPRM(3)=0.0                                                               
      FPRM(4)=POT(J,K)+PHINF                                                    
      FPRM(5)=0.0                                                               
      FPRM(6)=1.0                                                               
      FPRM(7)=0.0                                                               
      FPRM(8)=1.0                                                               
      INTL = 0                                                                  
      IORDER=4                                                                  
      DO 19992 IZ=1,NVERT                                                       
      WRK(IZ) = 0.0                                                             
19992 CONTINUE                                                                  
      CALL LTPBVP(INTL,IORDER,NVRT,FPRM,COFZ,WORK,ELEVRT,WRK(NVRT+1),           
     +PRTB,IERR)                                                                
      IF (.NOT.(IERR.NE.0)) GO TO 19991                                         
      CALL IOUT (IERR, "ERROR:")                                               
      CALL IOUT(K, "K")                                                          
      CALL IOUT(J, "J")                                                          
      CALL EXIT(0)
19991 CONTINUE                                                                  
      GRADGR(K,J) = (-3.*ELEVRT(1)+4.*ELEVRT(2)-ELEVRT(3))/(DLZ+DLZ)            
     +*DZDM(ZB)                                                                 
      JLAT=J                                                                    
      KLONG=K                                                                   
      CALL SETELE(JLAT,KLONG)                                                   
19993 CONTINUE                                                                  
19994 CONTINUE                                                                  
      CALL PROCESS                                                              
      RETURN                                                                    
      END                                                                       

      SUBROUTINE SETELE (JLAT,KLONG)                                            
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,          
     +ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),           
     +FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),           
     +AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),        
     +GL(4000),POT(37,73),                                                      
     +GCNK(37,37),GSNK(37,37),TES(26011)                                        
     +  ,GOO,F(37,73)                                                           
     +  ,LMX,ZSIG(37,73)                                                        
     +   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO                                
      COMMON/PHIFT/ PHINF                                                       
      COMMON /CURT/ TCUR,SUMCUP,SUMCDN                                          
      COMMON/TRANS/ DSIGZ                                                       
      COMMON / ANGLE / ANGX , ANGY                                              
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      J=JLAT                                                                    
      K=KLONG                                                                   
      DO 19999 L=1,NZLV                                                         
      ZLF = ZF(ZLEVEL(L))                                                       
      LL = INT((ZLF-ZB)/DLZ+0.5)+1                                              
      IF (LL.EQ.NVERT) LL = NVERT-1                                             
      ZLL = ZB+(LL-1)*DLZ                                                       
      LU = LL+1                                                                 
      IF (.NOT.(ZLF.GT.ZB.AND.ZLF.LT.ZT)) GO TO 19998                           
      ELEPOT(K,J,L) = ELEVRT(LL)+(ZLF-ZLL)/DLZ*(ELEVRT(LU)-ELEVRT(LL))          
      GO TO 19997                                                               
19998 CONTINUE                                                                  
      IF (.NOT.(ZLF.LE.ZB)) GO TO 19996                                         
      ELEPOT(K,J,L) =0.0                                                        
      GO TO 19997                                                               
19996 CONTINUE                                                                  
      IF (.NOT.(ZLF.GE.ZT)) GO TO 19995                                         
      ELEPOT(K,J,L) = POT(J,K)+PHINF                                            
19997 CONTINUE                                                                  
19995 CONTINUE                                                                  
19999 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       
      FUNCTION ZF(ZMET)                                                         
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      ZF = (EXP(-C2*ZMET)-1.0)/C3                                               
      RETURN                                                                    
      END                                                                       
      FUNCTION ZMETF(Z)                                                         
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      ZMETF = -ALOG(C3*Z+1.0)/C2                                                
      RETURN                                                                    
      END                                                                       
      FUNCTION DZDM(Z)                                                          
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      DZDM = -C2*(C3*Z+1.0)/C3                                                  
      RETURN                                                                    
      END                                                                       
      FUNCTION D2ZDM2(Z)                                                        
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      D2ZDM2 = C2**2*(C3*Z+1.0)/C3                                              
      RETURN                                                                    
      END                                                                       
      SUBROUTINE COFZ(Z,AF,BF,CF)                                               
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      DZDM2 = DZDM(Z)**2                                                        
      AF = DZDM2                                                                
      CALL CONDUCT(Z,SIGMA,SIGMAZ)                                              
      BF = D2ZDM2(Z) + DZDM2* SIGMAZ/SIGMA                                      
      CF = 0.0                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE CONDUCT(Z,SIGMA,SIGMAZ)                                        
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,          
     +ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),           
     +FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),           
     +AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),        
     +GL(4000),POT(37,73),                                                      
     +GCNK(37,37),GSNK(37,37),TES(26011)                                        
     +  ,GOO,F(37,73)                                                           
     +  ,LMX,ZSIG(37,73)                                                        
     +   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO                                
      COMMON/PHIFT/ PHINF                                                       
      COMMON /CURT/ TCUR,SUMCUP,SUMCDN                                          
      COMMON/TRANS/ DSIGZ                                                       
      COMMON / ANGLE / ANGX , ANGY                                              
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      ZMET=ZMETF(Z)                                                             
      XLATD=180.*XLAT/PII                                                       
      XLONGD=180.*XLONG/PII                                                     
      CALL COND(ZMET,XLATD,XLONGD)                                              
      SIGMA=SIG                                                                 
      SIGMAZ=DSIGZ/DZDM(Z)                                                      
      RETURN                                                                    
      END                                                                       


      SUBROUTINE PROCESS                                                        
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      LOGICAL LAT,LONG                                                          

      CALL NAME("GRADGR**",2)                                                   

      DO 19999 J= 1,NLAT                                                         
        XLATD= LATMSH(J)*180./PII                                                  
        CALL VOUT(XLATD, "LATTITUDE")                                               
        WRITE(6,200) (GRADGR(K,J),K=1,NLONG)                                      
  200   FORMAT(1H 12E11.3)                                                        
19999 CONTINUE                                                                  

      LONG = .TRUE.                                                             
      LAT=.TRUE.                                                                
      XLAT=65.*PII/180.                                                         
      XLONG = PII                                                               

c Plot gradient at ground...
	print *,"Calling GRAD ..."
      CALL GRAD(LAT,LONG)                                                       
	print *,"Returned from GRAD..."

	print *, "Calling GROUND..."
      CALL GROUND                                                               
	print *, "Returned from GROUND..."

	print *, "Calling CNTPOT..."
      CALL CNTPOT                                                               
	print *, "Returned from CNTPOT..."

      RETURN                                                                    
      END                                                                       


      SUBROUTINE GROUND                                                         
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)

      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,          
     +ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),           
     +FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),           
     +AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),        
     +GL(4000),POT(37,73),                                                      
     +GCNK(37,37),GSNK(37,37),TES(26011)                                        
     +  ,GOO,F(37,73)                                                           
     +  ,LMX,ZSIG(37,73)                                                        
     +   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO                                
      COMMON/PHIFT/ PHINF                                                       
      COMMON /CURT/ TCUR,SUMCUP,SUMCDN                                          
      COMMON/TRANS/ DSIGZ                                                       
      COMMON / ANGLE / ANGX , ANGY                                              
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      

      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      DIMENSION GRDGRT(NLATS,NLONGS)                                   

	character*50 title, xlabel, ylabel

        COMMON/POSNEG/CURPOS,CURNEG                                               
        COMMON/GRDDCC/GRDCUR                                                      
	parameter (rearth= 6.371e6)

c	write (title, "('Ground Gradient')")
c	call plchhq(400,900,TITLE,15,2,0)


      DO 19999 J=1,NLAT                                                         
        DO 19998 K=1,NLONG                                                        
          GRDGRT(J,K)= GRADGR(K,J)                                                   
19998   CONTINUE                                                                  
19999 CONTINUE                                                                  

	title= "Ground Gradient"
	xlabel= "Longitude"
	ylabel= "N. Latitude"
c	print *,"Calling scaleit..."
c	call scaleit(1.0, NLONG, NLAT, GRDGRT, scld, Zmin, Zmax)
c	print *, "Returned from scaleit..."
c	write (6, '("Zmin= ",f12.2,", Zmax= ",f12.2)') Zmin, Zmax
c	write (6, '("Calling mbsrfc...")')
c	call mbsrfc(scld, title, xlabel, ylabel, Zmin, Zmax, NLONG, NLAT)
      CALL EZSRFC(GRDGRT,NLAT,NLONG,ANGX,ANGY,WORK)                             

      DO 19997 J= 1, NLAT                                                         
        NLMJ= NLAT + 1 - J                                                             
        DO 19996 K= 1, NLONG                                                        
c         WORK(K,J)= GRADGR(K,NLMJ)
C The line below is a guess added by Bill Roberts, 6/26/90.
          WORK(K,J,1)= GRADGR(K,NLMJ)
19996   CONTINUE                                                                  
19997 CONTINUE                                                                  

	call mbcntr(WORK, title, xlabel, ylabel, NLONG, NLAT)
c*     call plchhq(400,600,TITLE,15,2,0)
c*      CALL EZCNTR(WORK,NLONG,NLAT)                                              

      GRDCUR= 0.0                                                              
      CURPOS= 0.0                                                                
      CURNEG= 0.0                                                                

      DO 19995 J=1,NLAT                                                         
        XLAT= LATMSH(J)                                                          
        SLT= SIN(XLAT)                                                           
        DO 19994 K=1,NLONG                                                        
          XLONG= LNGMSH(K)                                                         
      	  C3= EXP(-C2*ZSIG(J,K))-1.0                                                 
      	  Z= ZF(OROG(J,K))                                                           
      	  CALL CONDUCT(Z,SIGMA,SIGMAZ)                                              
          GRADGR(K,J)= SIGMA*GRADGR(K,J)                                             
          GRDCUR= GRDCUR + GRADGR(K,J)*SLT                                           
          GRDE= GRADGR(J,K)*SLT                                                      
          IF (.NOT.(GRDE.LT.0.)) GO TO 19993                                        
          CURNEG= CURNEG+GRDE                                                        
          GO TO 19992                                                               
19993     CONTINUE                                                                  
          CURPOS= CURPOS+GRDE                                                        
19992     CONTINUE                                                                  
19994   CONTINUE                                                                  
19995 CONTINUE                                                                  

      GRDCUR= DLAT*DLONG*GRDCUR*REARTH**2                                      
      CALL VOUT(CURPOS, "POS. CURR.")                                              

      CURNEG= DLAT*DLONG*CURNEG*REARTH**2                                        
      CALL VOUT(CURNEG, "NEG. CURR.")                                              

      CURPOS= DLAT*DLONG*CURPOS*REARTH**2                                        
      CALL VOUT(GRDCUR, "GRD. CURR.")                                                

      DO 19991 J= 1,NLAT                                                         
        DO 19990 K=1, NLONG                                                        
          GRDGRT(J,K)= GRADGR(K,J)                                                   
19990   CONTINUE                                                                  
19991 CONTINUE                                                                  

	title= "Ground Gradient"
	xlabel= "LONGITUDE"
	ylabel= "N. LATITUDE"
c	call scaleit(1.0, NLAT, NLONG, GRDGRT, scld, Zmin, Zmax)
c	call mbsrfc(scld, title, xlabel, ylabel, Zmin, Zmax, NLAT, NLONG)
      call plchhq(300,900,TITLE,23,2,0)
      CALL EZSRFC(GRDGRT,NLAT,NLONG,ANGX,ANGY,WORK)                             

      DO 19989 J= 1,NLAT                                                         
        NLMJ= NLAT + 1 - J                                                             
        DO 19988 K= 1, NLONG                                                        
      	  WORK(K,J,1)= GRADGR(K,NLMJ)
19988   CONTINUE                                                                  
19989 CONTINUE                                                                  

	title= "Ground Gradient"
	xlabel= "LONGITUDE"
	ylabel= "N. LATITUDE"
      call mbcntr(WORK, title, xlabel, ylabel, NLONG, NLAT)
c*      call plchhq(300,900,TITLE,23,2,0)
c*      CALL EZCNTR(WORK,NLONG,NLAT)                                              

      RETURN                                                                    
      END                                                                       


      SUBROUTINE CNTPOT                                                         
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)

      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      COMMON RE,RO,ZM,SIGM,PI,RAD,SRPI,DPHI,JMX,DTHETA,IMX,ABC1,SR2PI,          
     +ABC2,THET(73),THETA(73),PHII(37),PHI(37),DTHET,DPH,OROG(37,73),           
     +FS(37,73),Z1,SIGSL,SIGO,ALP,SIG,DRL,SIG1,AC1(37,37),BC1(37,37),           
     +AS1(37,37),BS1(37,37),AC2(37,37),BC2(37,37),AS2(37,37),BS2(37,37),        
     +GL(4000),POT(37,73),                                                      
     +GCNK(37,37),GSNK(37,37),TES(26011)                                        
     +  ,GOO,F(37,73)                                                           
     +  ,LMX,ZSIG(37,73)                                                        
     +   ,UBC(37,73),USNK(37,37),UCNK(37,37),UOO                                
      COMMON/PHIFT/ PHINF                                                       
      COMMON /CURT/ TCUR,SUMCUP,SUMCDN                                          
      COMMON/TRANS/ DSIGZ                                                       
      COMMON / ANGLE / ANGX , ANGY                                              
      DIMENSION EL1POT(NLATS,NLONGS), scld(NLATS,NLONGS)
	character*50 title, xlabel, ylabel


      DO 19999 L=1,NZLVL                                                        
        ZKM = ZLEVEL(L)/1000.                                                     
	write (title, "('Electric Potential at ,',f5.1,' KM')") ZKM 

        DO 19998 J= 1, NLAT                                                         
      	  NLMJ= NLAT + 1 - J                                                             
      	  DO 19997 K= 1, NLONG                                                        
c         WORK(K,J)=ELEPOT(K,NLMJ,L)
c The line below was added by Bill Roberts, 6/26/90.
      	  WORK(K,J,1)= ELEPOT(K,NLMJ,L)
      	  EL1POT(J,K)= ELEPOT(K,J,L)                                                 
19997     CONTINUE                                                                  
19998   CONTINUE                                                                  


	xlabel= "LONGITUDE"
	ylabel= "N. LATITUDE"
	call mbcntr(WORK, title, xlabel, ylabel, NLONG, NLAT)
c*      call plchhq(200,900,TITLE,30,2,0)
c*      CALL EZCNTR(WORK,NLONG,NLAT)                                              

	call scaleit(1.e-2, NLAT, NLONG, EL1POT, scld, Zmin, Zmax)
	call mbsrfc(scld, title, xlabel, ylabel, Zmin, Zmax, NLAT, NLONG)

c      CALL EZSRFC(EL1POT,NLAT,NLONG,ANGX,ANGY,WORK)                             
c      CALL FRAME                                                                
19999 CONTINUE                                                                  

      RETURN                                                                    
      END                                                                       

      SUBROUTINE GRAD(LAT,LONG)                                                 
      parameter (NLONGS= 73)
      parameter (NLATS= 37)
      parameter (NZLVL= 5)
      parameter (NVERT = 100)
      COMMON / LOWERL / ELEPOT(NLONGS,NLATS,NZLVL),                    
     +ZLEVEL(NZLVL),GRADGR(NLONGS,NLATS),LATMSH(NLATS),           
     +LNGMSH(NLONGS),WORK(NLONGS,NLATS,3),ELEVRT(NVERT),           
     +FPRM(8),ZTOP,ZBOT,ZT,ZB,DLZ,C2,C3,PII,DLAT,DLONG,XLAT,XLONG,              
     +NLAT,NLONG,NVRT,NZLV                                                      
      REAL LATMSH,LNGMSH                                                        
      DIMENSION WRK(10*NVERT)                                                   
      EQUIVALENCE(WRK,WORK)                                                     
      DIMENSION POTG(NLONGS)                                                
c      DIMENSION LABG(4)
      character*40 labg
c      REAL LABG
      REAL LDASH                                                                
      LOGICAL LAT,LONG                                                          
c      DATA NDASH,LDASH / 0, 1777777777777777777777B /                         
	data ndash / 0 /
	data ldash / O'37777777777'/
      IF (.NOT.(LAT)) GO TO 19999                                               
      XLATD=180.*XLAT/PII                                                       

	write (labg, "('Fixed Lattitude = ',f5.1,' degrees$')") xlatd

      J= INT(XLAT/DLAT+0.5)+1                                                    
      DO 19998 K=1,NLONG                                                        
        POTG(K)=GRADGR(K,J)                                                       
19998 CONTINUE                                                                  

      IBAC=1                                                                    
      ISET=1                                                                    

      CALL ANOTAT ("Longitude", "Gradient at Ground", IBAC,
     1  ISET,NDASH,LDASH)                                                    

      CALL VOUT(XLATD, "XLATD")                                                  
      CALL IOUT(NLONG, "NLONG")                                                  
      CALL VECOUT(LNGMSH, "LNGMSH",NLONG)                                        
      CALL VECOUT(POTG, "POTG",NLONG)                                            
      CALL EZXY(LNGMSH,POTG,NLONG,LABG)                                         
      CALL VOUT(XLATD, "XLATD")                                                  

19999 CONTINUE                                                                  
      IF (.NOT.(LONG)) GO TO 19997                                              
      XLONGD=180.*XLONG/PII                                                     

	write (labg, "('Fixed Longitude = ',f5.1,' degrees$')") xlongd

      K=INT(XLONG/DLONG+0.5)+1                                                  

      DO 19996 J=1,NLAT                                                         
        POTG(J)=GRADGR(K,J)                                                       
19996 CONTINUE                                                                  

      IBAC=1                                                                    
      ISET=1                                                                    

      CALL ANOTAT("Latitude", "Gradient at Ground", IBAC, 
     1   ISET, NDASH, LDASH)

      CALL VOUT(XLONGD,"XLONGD")                                                
      CALL IOUT(NLAT, "NLAT")                                                    
      CALL VECOUT(LATMSH, "LATMSH",NLAT)                                         
      CALL VECOUT(POTG, "POTG",NLAT)                                             
      CALL EZXY(LATMSH,POTG,NLAT,LABG)                                          

19997 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       


      SUBROUTINE VOUT(VAR,NAME)                                                 
      character*10 NAME                                                                 

      WRITE (6,100) NAME,VAR                                                    
  100 FORMAT(1x,a,"= ",e12.5)                                                 
      RETURN                                                                    
      END                                                                       


      SUBROUTINE IOUT(IVAR,NAME)                                                
      character*10 NAME                                                                 

      WRITE (6,100) NAME,IVAR                                                   
  100 FORMAT(1x,a,"= ",i8)                                                    
      RETURN                                                                    
      END                                                                       


      SUBROUTINE VECOUT(VEC,NAME,LENGTH)                                        
      character*10 NAME                                                                 
      DIMENSION VEC(LENGTH)                                                     

      WRITE (6,100) NAME                                                        
  100 FORMAT(1x,"*",A8)                                                             
      WRITE(6,200) (VEC(I),I=1,LENGTH)                                          
  200 FORMAT(1x,10e12.5)                                                        

      RETURN                                                                    
      END                                                                       


      SUBROUTINE IVECOUT(IVEC,NAME,LENGTH)                                      
      character*10 NAME                                                                 
      DIMENSION IVEC(LENGTH)                                                    

      WRITE(6,100) NAME                                                         
  100 FORMAT(1x,"*",A8)                                                             
      WRITE (6,200) (IVEC(I),I=1,LENGTH)                                        
  200 FORMAT(1x,10i10)                                                          
      RETURN                                                                    
      END                                                                       


      SUBROUTINE AOUT(ARRAY,NAME,IA,NROW,NCOL)                                  
      character*10 NAME                                                                 
      DIMENSION ARRAY(IA,NCOL)                                                  

      WRITE (6,100) NAME                                                        
  100 FORMAT(1x,"*",A8)                                                             
      DO 1 J=1,NCOL                                                             
      WRITE (6,200) J                                                           
  200 FORMAT(1x," COLUMN",I4)                                                    
      CALL VECOUT(ARRAY(1,J),NAME,NROW)                                          
    1 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       


      SUBROUTINE NAME(NAM,KIP)                                                  
      character*10 NAM                                                                  

      DO 1 I=1,KIP                                                              
        WRITE (6,50)                                                              
   50   FORMAT(1x," ")                                                               
    1 CONTINUE                                                                  
      WRITE(6,100) NAM                                                          
  100 FORMAT(1x,"*",A8)                                                             
      RETURN                                                                    
      END                                                                       
