      BLOCK DATA                                                        
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      COMMON /VERIFY/ VDAT(15),CDAT(15)                                 
      COMMON /PI/ PI                                                    
      COMMON /DEBUG/ NUMBER,TLAST                                       
C                                                                       
      DATA            TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT  
     1              /     0.,    0.,    0.,    0.,    0.,    0.,    0./ 
      DATA            TXSHUF,TUSHUF,TXPRNT,TXLOG2                       
     1              /     0.,    0.,    0.,    0./                      
      DATA            EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT                
     1              /     0.,    0.,    0.,    0.,    0./               
      DATA            NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT  
     1              /     0 ,    0 ,    0 ,    0 ,    0 ,    0 ,    0 / 
      DATA            NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
     1              /     0 ,    0 ,    0 ,    0 /                      
      DATA VDAT/                                                        
     * 0.10000E+03,-0.10000E+03, 0.10000E-04,-0.10000E+03,-0.15708E+01, 
     * 0.99983E+02,-0.18147E+01,-0.20427E+01,-0.22452E+01,-0.24181E+01, 
     * 0.95976E+02,-0.95977E+02, 0.10000E-04,-0.95978E+02,-0.25589E+01/ 
C                                                                       
C PI IS A USEFULL CONSTANT                                              
C                                                                       
      DATA  PI /3.141592653589793/                                      
      DATA NUMBER /1/, TLAST /0.0/                                      
      END                                                               
C                                                                       
C                                                                       
C******NOTE: USE OF MACHINE-DEPENDENT CONSTANT: ZERO IN SUBROUTINE RUNIN
C                                                                       
C            ZERO = 1.0E-200 (MAY HAVE TO REPLACE WITH YOUR MACHINES    
C                             RELATIVE PRECISION FOR REAL DATA TYPE)    
C                                                                       
C                    THIS PROGRAM VERSION HAS BEEN SET UP AT PRINCETON B
C                    FRED SEIDL (I.E., "FGPS") FOR USE AS A PERFECT-CLUB
C                    BENCHMARK.   IT MAY BE NOTED THAT TIMING MEASUREMEN
C                    MAKE USE OF (1) THE CYBER205 INTRINSIC SCALAR FUNCT
C                    SECOND() AND OF (2) THE CYBER205 SYSTEM-INTERFACE-L
C                    UAGE CALL Q5DCDMSC('MILLSEC=',T),  WHERE "T" IS THE
C                    WALL-CLOCK TIME IN MILLISECONDS.  SECOND() AND Q5DC
C                    ARE EMBEDDED IN SUBROUTINES CPUTIM AND ELAPSE RESPE
C                    TIVELY.  MOREOVER, IN ORDER TO DISPLAY THE OBJECT M
C                    ULE SIZE AND LOAD MODULE SIZE (BY WAY OF THE LOG FI
C                    THE CYBER205 CONTROL STATEMENT "FILES." IS INCLUDED
C                    THE JCL.                                           
C                                                                       
C     FGPS 12-MAY-88 WHEN CODE IS USED AS A BENCHMARK, SET "RLAMDA = 1.0
C                    IN SUBROUTINE DATAST                               
C     FGPS 11-MAY-88 REMOVED SOME INCONSISTENCIES IN NORMALIZING THE FFT
C                    AND RESET EPS = 100. IN SUBROUTINE DATAST          
C     FGPS  9-MAY-88 TIMING MEASUREMENTS ON DEALSE REPLACED BY MEASUREME
C                    ON FTRVMT (SEE "FGPS 10-MAR-88")                   
C     FGPS  4-MAY-88 MADE PROVISIONS TO MEASURE CPU EXECUTION TIMES DURI
C                    THE LAST TIMESTEP FOR SUBROUTINES RCS, CSR, SCSC, A
C     FGPS 29-APR-88 REPLACED FUNCTION LOG2 (NA) BY SUBROUTINE XLOG2 (NA
C     FGPS 15-APR-88 INCREASED THE INITIAL AMPLITUDE PRESCRIBED FOR THE 
C                    STREAMFUNCTION BY *50 (I.E., EPS = 50.0 IN DATAST) 
C     FGPS 15-APR-88 THIS VERSION OF THE SOURCE HAS,  IN A WAY CONSISTEN
C                    WITH THE PRESCRIBED PARAMETERS, ALL NONTRIVIAL UNRE
C                    ERENCED MODULES COMMENTED OUT USING THE FLAG "C2%2"
C                    WHERE C IS LOCATED IN COLUMN 1                     
C     FGPS 14-APR-88 WHILE IN USE AS A BENCHMARK-TEST CODE, CERTAIN STAT
C                    MENTS HAVE BEEN COMMENTED OUT.  IN PARTICULAR, CALL
C                    TO LUMULT2 AND SOLVE2 WHERE THESE ARE HANDLED AS DY
C                    NAMIC EXTERNALS ARE NOT TO BE MADE. TO FACILITATE R
C                    TORATION OF THE DETACHED PROGRAMMING,  THE CHARACTE
C                    STRING C2%2, WHERE C IS PLACED IN COLUMN 1, HAS BEE
C                    USED TO TRANSFORM STATEMENTS INTO COMMENT LINES    
C     FGPS  8-APR-88 INTRODUCED SUBROUTINE RPRNT2 (A MODIFIED VER-      
C                    SION OF RPRINT), WHICH AT INTERVALS PRINTS THE     
C                    TEMPERATURE AT THE CENTER OF THE FIRST CONVECTION  
C                    CELL SO THAT THIS QUANTITY CAN BE USED TO MONITOR  
C                    CODE EXECUTION                                     
C     FGPS  5-APR-88 FUNCTION LCMREQ (LARGE) CONTAINS                   
C                      PARAMETER (LCMSIZ=165117)                        
C                      COMMON /MEMORY/ MEMSIZ,RM(LCMSIZ)                
C                    ELSEWHERE RESET RM(25000)==>RM(165117)             
C     FGPS  4-APR-88 MOVED THE CALLS TO SUBROUTINE PSIPHY AND COR-      
C                    RECTED FOR ANY CHANGES MADE TO THE STREAMFUNC-     
C                    TION IN CALLING PSIPHY BY CALLING IN (CA,LPN,      
C                    NW) BEFORE AND AFTER CALLING PSIPHY                
C     FGPS  2-APR-88 REWROTE SUBROUTINE PSIPHY                          
C     FGPS 29-MAR-88 ADDED SUBROUTINES SHUF, UNSHUF, AND RPRINT TO      
C                    THE LIST SUBJECT TO MEASUREMENTS OF EXECUTION      
C                    TIME AND NUMBER-OF-CALLS-TO                        
C     FGPS 15-MAR-88 MODIFIED SUBROUTINE RPRINT(ARRAY,NAME,N1H,N2PD,IMOD
C                    AND REFERENCING SUBROUTINES SO AS TO ALLOW SELECTIO
C                    OF ONE OF THREE ALTERNATIVE PRINT FORMATS ACCORDING
C                    TO IMODE = -1,0,+1                                 
C     FGPS 10-MAR-88 MEASUREMENTS OF CPU EXECUTION TIME PER CALL AND    
C                    OF THE CUMULATIVE NUMBER OF CALLS WERE ADDED TO THE
C                    SUBROUTINES:  IN, OUT, RCS, CSR, SCSC, ACAC, DEALSE
C                    MEASURED VALUES ENTERED INTO THE MAIN PROGRAM VIA  
C                    COMMON /TIMING/, AND PRINTED OUT FROM THERE.  LATER
C                    (9-MAY-88), MEASUREMENTS ON DEALSE REPLACED BY FTRV
C     FGPS  8-MAR-88 REPROGRAMMED SUBROUTINE RPRINT SO AS TO PRINT OUT  
C                    QUANTITIES REPRESENTED IN PHYSICAL (OR REAL) SPACE 
C     FGPS 12-FEB-88 ONLY UNUSED MODULES IN THIS VERSION ARE DATE AND   
C                    TIME.  ALSO, SUBR. PLOTME COMMENTED "C" OUT        
C     FGPS 10-FEB-88 INSERTED "COMPLEX EX(1)" INTO SUBROUTINE TEMPHY    
C     FGPS 9-FEB-88  DELETED EXTRA LISTINGS OF SUBROUTINE EXSET AND     
C                    FUNCTION LOG2                                      
C     FGPS 4-DEC-87  TEMPORARILY COMMENTED OUT (I.E., "C") CALLS TO SUB-
C                    ROUTINE PLOTME                                     
C     FGPS PRIOR TO DEC-87 THE USE OF "*" AS A DELIMITER IN FORMAT      
C                    STATEMENTS WAS REPLACED BY THE SINGLE QUOTE        
C                                                                       
C***********************************************************************
C                        OCEAN FORTRAN PROGRAMS                         
C                        CREATION DATE: 3/22/81                         
C***********************************************************************
      PROGRAM OCEAN                                                     
C                                                                       
C  D(Z)/DT + (U,V).GRAD(Z) = RNU*DELSQ(Z) + ALPHA*D(TEMP)/DX            
C                                                                       
C  D(TEMP)/DT + (U,V).GRAD(TEMP) = RKAPA*DELSQ(TEMP) + BETA*V           
C                                                                       
C                                                                       
C                                                                       
C THIS PROGRAM SIMULATES TWO-DIMENSIONAL FLUID FLOW ON A   (N2+1) * N1  
C GRID. IT IS POSSIBLE TO SPECIFY NON-LINEAR DISTRIBUTION OF DATA POINTS
C THE Y-DIRECTION (THE FIRST INDEX CORRESPONDS TO Y DIRECTION, SECOND TO
C X DIRECTION). N1 AND N2 HAVE TO BE POWERS OF TWO GREATER OR EQUAL TO 8
C                                                                       
C THE PROGRAM USES LARGE CORE MEMORY. IT REQUIRES 8*(N2+1)*N1+(N2-1)*(N2
C WORDS OF IT. IF NOTEMP = 1 IT IS POSSIBLE TO SAVE 2*(N2+1)*N1 WORDS,  
C IF LMAPP = 1 YOU SAVE ADDITIONAL (N2-1)*(N2+1) WORDS.                 
C                                                                       
C IF  LMAPP = 2 , LCM REQUIREMENTS ARE AT LEAST  3 * (N2-1) * (N2-1)    
C MORE MEMORY COULD BE NEEDED IN THIS CASE BECAUSE OF EIGSET AND LCMRQR.
C                                                                       
C                                                                       
C FOR CDC7600 RUNS: DELETE ROUTINES LCMREQ, IN AND OUT FROM THE END OF  
C                   THIS PACKAGE. ON CDC7600 THE LCM IS USED THROUGH    
C                   THE ASSEMBLY ROUTINES.                              
C                                                                       
C                   LIBRARY.                                            
C FOR CRAY-1 RUNS : DELETE ROUTINES SECOND, DATE, TIME, EXSET AND       
C                   LOG2 FROM THE END OF THIS PACKAGE. THEY ARE         
C                   OBTAINED AT LOAD TIME FROM CRAY LIBRARIES AND       
C                   FFTLIB.                                             
C                   ALSO, ADJUST THE SIZE OF PARAMETER LCMSIZ IN        
C                   THE ROUTINE LCMREQ SO THAT THE COMMON IS LARGE      
C                   ENOUGH TO EMULATE CDC7600 LCM REQUIRED FOR THE RUN. 
C                                                                       
C THE FOLLOWING PARAMETERS DEFINE THE SIZE OF THE RUN:                  
C                                                                       
      PARAMETER (NX=128, NY=128)                                        
      PARAMETER (NA1=2*(NY+1),NA2=NX/2,NE=NY-1,NW1=NY+1)                
      PARAMETER (NW2=128,NEX=256)                                       
      PARAMETER (NI=2*NY-1,NEQ=3*NY-2,NCW2=2*NW1+1,NCW3=4*NW1+1)        
C                                                                       
C PLEASE REFER TO THE DIMENSIONING AND EQUIVALENCING INSTRUCTIONS IN    
C THE CASE YOUR COMPILER CANNOT HANDLE PARAMETER STATEMENTS...          
C COMMON BLOCKS:                                                        
C                                                                       
      COMMON /FT/ NPTS, NSKIP, MTRN, MSKIP, ISIGN, LOG, IEXSHF          
C                                                                       
C THE DIMENSION OF  GSOLV  HAS TO BE AT LEAST (N2-1)                    
C                                                                       
      COMMON /SOLV/ LVECT, LSKIP, NVECT, NVSKIP, GSOLV(127)             
C                                                                       
C THE DIMENSIONS OF  GMAPP  AND  G  HAVE TO BE AT LEAST (N2+1)          
C                                                                       
      COMMON /MAPING/ GMAPP(129), G(129)                                
      COMMON /MFLOW/ U(129)                                             
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1     KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,       
     2     KPRESS, NOALS                                                
      COMMON /PI/ PI                                                    
      COMMON /NMBRS/ N1H, N1Q, N1E, N2M, N2P, N2DP, N2DM, N2PD          
      COMMON /EXSIZE/ NMAX2                                             
      COMMON /LADDR/ LZN,LTN,LPN,LZO,LTO,LPY,LZX,LZY,LEIG,LE,LH,LEV     
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
C                                                                       
C                                                                       
C COMPLEX ARRAYS HAVE TO BE DIMENSIONED EXACTLY AS THE FORMULAS SAY:    
C            CA        (N2+1) X N1/2                                    
C            CWORK     N2+1                                             
C            CWORK2    N2+1                                             
C            CWORK3    N2+1                                             
C            EX        MAX(2*N2, N1)                                    
C                                                                       
C EX ARRAY MUST NOT BE THE FIRST ARRAY IN THE MEMORY.                   
C                                                                       
      COMPLEX CA(NW1,NA2), CWORK(NW1), CWORK2(NW1)                      
      COMPLEX CWORK3(NW1), EX(NEX), CF                                  
C                                                                       
C REAL ARRAYS HAVE TO BE DIMENSIONED EXACTLY AS THE FORMULAS SAY:       
C            A      2*(N2+1) X N1/2                                     
C            E      (N2-1)   X (N2-1)                                   
C            WORK   N2+1     X MAX(6, 2*MIN(64, N1/2))   (AT LEAST)     
C            EVAL   N2-1                                                
C            IR     N2-1                                                
C            IC     N2-1                                                
C                                                                       
      REAL A(NA1,NA2), E(NE,NE)                                         
      REAL WORK(NW1*NW2)                                                
      REAL EVAL(NE), IR(NE), IC(NE)                                     
      INTEGER VALID                                                     
C                                                                       
C EQUIVALENCE STATEMENTS ARE SAVING MEMORY SPACE. THE FOLLOWING PATTERN 
C HAS TO BE FOLLOWED EXACTLY:                                           
C                                                                       
C     A-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+           
C     CA+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+           
C     EVAL+-+IR-+-+-IC+-+-+E+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+             
C                                                                       
C                                                                       
C     WORK+-+-+-+-+-+-+-+-+-+-+-+-+-+-+                                 
C     CWORK-+-+-+CWORK2-+-+-CWORK3+-+-+                                 
C                                                                       
C      EQUIVALENCE (A,CA,EVAL), (A(NY),IR), (A(NI),IC), A(NEQ),E) - ORGI
C      EQUIVALENCE (A,CA,EVAL), (EVAL(NY),IR), (EVAL(NI),IC),           
C     *             (EVAL(NEQ),E) - MACFARLANE'S FIX. BAD.              
C POINTER'S FIX. ANSI REQUIRES THAT THE NUMBER OF SUBSCRIPTS MATCH      
C THE NUMBER OF DIMENSIONS OF THE ARRAY. THAT WAS THE PROBLEM WITH      
C THE ORGINAL STATEMENT. THE PROBLEM WITH MACFARLANE'S FIX IS THAT      
C NY,NI,NEQ>NE (THE DIMENSION OF EVAL). HERE 3*NE+1-NA1 IS THE 3*NE+1TH 
C ELEMENT OF A AND 3*NE+1 = NEQ.                                        
      EQUIVALENCE (A,CA,EVAL), (A(NY,1),IR), (A(NI,1),IC),              
     *            (A(3*NE+1-NA1,2),E)                                   
      EQUIVALENCE (WORK,CWORK), (WORK(NCW2),CWORK2), (WORK(NCW3),CWORK3)
      COMMON /VERIFY/ VDAT(15),CDAT(15)                                 
CX00001 CALL MAINIX (1)                                                 
CY00001 CALL MAINIY (1)                                                 
CX60000 CALL CNTINX ()                                                  
CY60000 CALL CNTINY ()                                                  
C                                                                       
C NOW WE CAN START WITH THE PROGRAM...                                  
C SUBROUTINE START TELLS OCEAN WHAT TO DO:                              
C                                                                       
      CALL CPUTIM (CPU1)                                                
      CALL ELAPSE (TCLCK1)                                              
      OPEN(1,FILE='OCV')                                                
      OPEN(6,FILE='OCO6')                                               
      WRITE(1,9991)                                                     
9991  FORMAT(1X,
     *  'OUTPUT FOR PERFECT CLUB BENCHMARK: OCEAN INSTRUMENTED VERSION'/  
     *       1X,'$Revision: 1.A.1 $ $Author: psharma $'/)                    
                                                                        
      CALL START (NX,NY)                                                
C                                                                       
C IF  LMAPP = 1 , H  (THE SIZE OF THE BOX IN Y DIRECTION)  SHOULD BE    
C SPECIFIED IN THE MAIN PROGRAM.                                        
C                                                                       
      H = PI                                                            
C                                                                       
C DEFINING CONSTANTS THAT WILL BE NEEDED:                               
C                                                                       
      N1H = N1/2                                                        
      N1Q = N1/4                                                        
      N1E = N1/8                                                        
      N2M = N2-1                                                        
      N2P = N2+1                                                        
      N2DP = 2*N2+1                                                     
      N2DM = 2*N2-1                                                     
      N2PD = 2*N2P                                                      
      NW = N2PD*N1H                                                     
      NWH = NW/2                                                        
      NWQ = NW/4                                                        
      NWEIG = N2M                                                       
      IF (LMAPP .NE. 1) NWEIG = (N2M+3)*N2M                             
      NMAX2 = MAX(2*N2,N1)                                              
      K1Q = N1E+1                                                       
      K2Q = K1Q+N1E                                                     
      K3Q = K2Q+N1E                                                     
      K1QM = K1Q-1                                                      
      K2QM = K2Q-1                                                      
      K3QM = K3Q-1                                                      
C                                                                       
C LCM ALLOCATION AND COMPUTATION OF ADDRESSES:                          
C                                                                       
      LARGE = 8*NW                                                      
      IF (NOTEMP .EQ. 1) LARGE = 6*NW                                   
      IF (LMAPP .NE. 2) GO TO 66                                        
C                                                                       
      IF (KPRESS .GT. 1) LARGE = MAX(LARGE, 3*N2M*N2M) +NWEIG+NWEIG     
      IF (KPRESS .LE. 1) LARGE = MAX(LARGE+NWEIG, 3*N2M*N2M)            
      GO TO 67                                                          
C                                                                       
   66 LARGE = LARGE + NWEIG                                             
C                                                                       
   67 LZN = LCMREQ (LARGE)                                              
      LTN = LZN+NW                                                      
      IF (NOTEMP .EQ. 1) LTN = LZN                                      
      LPN = LTN+NW                                                      
      LZO = LPN+NW                                                      
      LTO = LZO+NW                                                      
      IF (NOTEMP .EQ. 1) LTO = LZO                                      
      LPY = LTO+NW                                                      
      LZX = LPY+NW                                                      
      LZY = LZX+NW                                                      
      LPRESS = LZN + LARGE - NWEIG                                      
      LEIG = LPRESS                                                     
      IF ((KPRESS .GT. 1) .AND. (LMAPP .NE. 1)) LEIG = LPRESS - NWEIG   
C                                                                       
C THESE LCM ADDRESSES ARE USED BY EIGSET AND LCMRQR                     
C                                                                       
      LE = LZN                                                          
      LH = LE + N2M*N2M                                                 
      LEV = LH + N2M*N2M                                                
C                                                                       
      CALL EXSET (EX,NMAX2)                                             
CX10501 CALL DOINX( 501 )                                               
CY10501 CALL DOINY( 501 )                                               
C                                                                       
C SET ALL  U(I) TO 0. TO AVOID TROUBLE IF NO MEAN FLOW IS DESIRED       
C                                                                       
      DO 11 I = 1,N2P                                                   
11         U(I) = 0.                                                    
CX20502 CALL DOOUTX( 502 )                                              
CY20502 CALL DOOUTY( 502 )                                              
C                                                                       
C SET UP ARRAYS ASSOCIATED WITH THE MAPPING FUNCTION:                   
C                                                                       
      IF (LMAPP .NE. 1) GO TO 2                                         
C                                                                       
C LINEAR MAPPING - EIGENVALUES KNOWN, EIGENFUNCTIONS ARE SIMPLY SINES   
C AND COSINES, SO THAT NO MATRICES ARE NEEDED:                          
C                                                                       
           F = H/REAL(N2)                                               
           FZ = PI/H                                                    
           FZSQ = FZ*FZ                                                 
           GMAPP(1) = 0.                                                
           G(1) = FZ                                                    
CX10503 CALL DOINX( 503 )                                               
CY10503 CALL DOINY( 503 )                                               
           DO 1 I = 1,N2M                                               
                GMAPP(I+1) = F*REAL(I)                                  
                G(I+1) = FZ                                             
1               EVAL(I) = -(REAL(I)**2)*FZSQ                            
CX20504 CALL DOOUTX( 504 )                                              
CY20504 CALL DOOUTY( 504 )                                              
           GMAPP(N2P) = H                                               
           G(N2P) = FZ                                                  
           CALL OUT (EVAL,LEIG,NWEIG)                                   
           GO TO 4                                                      
2     IF (LMAPP .NE. 2) GO TO 3                                         
3     CONTINUE                                                          
4     CONTINUE                                                          
C                                                                       
C READ THE DATA:                                                        
C                                                                       
      IF ((KDATA .NE. 1) .AND. (KDATA .NE. 2)) GO TO 13                 
           CALL DATAST (A,N2P,N1,GMAPP)                                 
           GO TO 14                                                     
13    CONTINUE                                                          
14    CONTINUE                                                          
      FX = 2.0*PI/BOXX                                                  
      FX2 = 2.*FX                                                       
C                                                                       
C PREPARE THE INPUT FOR PROCESSING:                                     
C                                                                       
      IF (KDATA .NE. 1) GO TO 40                                        
C                                                                       
C THE INPUT WAS STREAMFUNCTION IN PHYSICAL SPACE - SHUFFLE THE DATA AND 
C GO TO FOURIER SPACE:                                                  
C                                                                       
           CALL IN (A,LPN,NW)                                           
           CALL SHUF (A,N2P,N1,WORK)                                    
           NPTS = N1                                                    
           NSKIP = N2P                                                  
           MTRN = N2P                                                   
           MSKIP = 1                                                    
           ISIGN = -1                                                   
           CALL XLOG2 (NPTS,LOG)                                        
           IEXSHF = NMAX2/NPTS                                          
           CALL RCS (A,EX)                                              
           NPTS = 2*N2                                                  
           NSKIP = 1                                                    
           MTRN = N1H                                                   
           MSKIP = N2P                                                  
           ISIGN = -1                                                   
           CALL XLOG2 (NPTS,LOG)                                        
           IEXSHF = NMAX2/NPTS                                          
           CALL ACAC (A,EX,WORK)                                        
C           F = 0.125/REAL(N1*N2)                                       
           F = 0.0625/REAL(N1*N2)                                       
CX10505 CALL DOINX( 505 )                                               
CY10505 CALL DOINY( 505 )                                               
           DO 30 J = 1,N1H                                              
                DO 30 I = 1,N2PD                                        
30                   A(I,J) = A(I,J)*F                                  
CX20506 CALL DOOUTX( 506 )                                              
CY20506 CALL DOOUTY( 506 )                                              
           CALL OUT (A,LPN,NW)                                          
           CALL OUT (A,LZN,NW)                                          
C                                                                       
CC*****C H E C K   F O U R I E R   C O E F F I C I E N T S**************
C                                                                       
C         WRITE (6,7010) A(1,1),A(1,2),A(1,3)                           
C         WRITE (6,7020) A(2,1),A(2,2),A(2,3)                           
C         WRITE (6,7030) A(3,1),A(3,2),A(3,3)                           
C 7010    FORMAT(/2X,26HCHECK FOURIER COEFFICIENTS/4X,8HA(1,1) =,1PE13.6
C     1   2X,8HA(1,2) =,1PE13.6,2X,8HA(1,3) =,1PE13.6)                  
C 7020    FORMAT(4X,8HA(2,1) =,1PE13.6,2X,8HA(2,2) =,1PE13.6,2X,        
C     1   8HA(2,3) =,1PE13.6)                                           
C 7030    FORMAT(4X,8HA(3,1) =,1PE13.6,2X,8HA(3,2) =,1PE13.6,2X,        
C     1   8HA(3,3) =,1PE13.6/)                                          
C                                                                       
C                                                                       
40    IF (KDATA .EQ. 3) GO TO 145                                       
C                                                     2                 
C STREAMFUNCTION KNOWN - FIND THE VORTICITY:  Z =  DEL  (PSI):          
C                                                                       
           CALL IN (EVAL,LEIG,NWEIG)                                    
           KZN = LZN                                                    
           LVECT = N2M                                                  
           LSKIP = 2                                                    
           NVECT = 2                                                    
           NVSKIP = 1                                                   
CX10507 CALL DOINX( 507 )                                               
CY10507 CALL DOINY( 507 )                                               
           DO 60 J = 1,N1H                                              
                CALL IN (CWORK,KZN,N2PD)                                
                F = REAL((J-1)*(J-1))*FX*FX                             
                DO 50 I = 2,N2                                          
50                   CWORK(I) = (EVAL(I-1)-F)*CWORK(I)                  
                CWORK(1) = -F*CWORK(1)                                  
                CWORK(N2P) = 0.0                                        
                CALL OUT (CWORK,KZN,N2PD)                               
60              KZN = KZN+N2PD                                          
CX20508 CALL DOOUTX( 508 )                                              
CY20508 CALL DOOUTY( 508 )                                              
           NPTS = 2*N2                                                  
           NSKIP = 1                                                    
           MTRN = N1H                                                   
           MSKIP = N2P                                                  
           ISIGN = 1                                                    
           CALL XLOG2 (NPTS,LOG)                                        
           IEXSHF = NMAX2/NPTS                                          
           CALL IN (A,LZN,NW)                                           
           CALL ACAC (A,EX,WORK)                                        
CX10509 CALL DOINX( 509 )                                               
CY10509 CALL DOINY( 509 )                                               
           DO 70 J = 1,N1H                                              
                DO 70 I = 1,N2PD                                        
70                   A(I,J) = 0.5*A(I,J)                                
CX20510 CALL DOOUTX( 510 )                                              
CY20510 CALL DOOUTY( 510 )                                              
           CALL OUT (A,LZO,NW)                                          
      IF (NOTEMP .EQ. 1) GO TO 140                                      
      IF (KDATA .NE. 1) GO TO 120                                       
C                                                                       
C TEMPERATURE IS SUPPLIED IN PHYSICAL SPACE - SHUFFLE THE DATA AND      
C GO TO FOURIER REPRESENTATION:                                         
C                                                                       
           CALL IN (A,LTN,NW)                                           
           CALL SHUF (A,N2P,N1,WORK)                                    
           NPTS = N1                                                    
           NSKIP = N2P                                                  
           MTRN = N2P                                                   
           MSKIP = 1                                                    
           ISIGN = -1                                                   
           CALL XLOG2 (NPTS,LOG)                                        
           IEXSHF = NMAX2/NPTS                                          
           CALL RCS (A,EX)                                              
           F = 0.5/REAL(N1)                                             
CX10511 CALL DOINX( 511 )                                               
CY10511 CALL DOINY( 511 )                                               
           DO 100 J = 1,N1H                                             
                DO 100 I = 1,N2PD                                       
100                  A(I,J) = F*A(I,J)                                  
CX20512 CALL DOOUTX( 512 )                                              
CY20512 CALL DOOUTY( 512 )                                              
           CALL OUT (A,LTO,NW)                                          
           NPTS = 2*N2                                                  
           NSKIP = 1                                                    
           MTRN = N1H                                                   
           MSKIP = N2P                                                  
           ISIGN = -1                                                   
           CALL XLOG2 (NPTS,LOG)                                        
           IEXSHF = NMAX2/NPTS                                          
           CALL ACAC (A,EX,WORK)                                        
C           F = 0.25/REAL(N2)                                           
           F = 0.125/REAL(N2)                                           
CX10513 CALL DOINX( 513 )                                               
CY10513 CALL DOINY( 513 )                                               
           DO 110 J = 1,N1H                                             
                DO 110 I = 1,N2PD                                       
110                  A(I,J) = F*A(I,J)                                  
CX20514 CALL DOOUTX( 514 )                                              
CY20514 CALL DOOUTY( 514 )                                              
           CALL OUT (A,LTN,NW)                                          
      GO TO 140                                                         
120        CALL IN (A,LTN,NW)                                           
C                                                                       
C NEW TEMPERATURE KNOWN - SET UP T(OLD):                                
C                                                                       
           NPTS = 2*N2                                                  
           NSKIP = 1                                                    
           MTRN = N1H                                                   
           MSKIP = N2P                                                  
           ISIGN = 1                                                    
           CALL XLOG2 (NPTS,LOG)                                        
           IEXSHF = NMAX2/NPTS                                          
           CALL ACAC (A,EX,WORK)                                        
CX10515 CALL DOINX( 515 )                                               
CY10515 CALL DOINY( 515 )                                               
           DO 130 J = 1,N1H                                             
                DO 130 I = 1,N2PD                                       
130                  A(I,J) = A(I,J)*0.5                                
CX20516 CALL DOOUTX( 516 )                                              
CY20516 CALL DOOUTY( 516 )                                              
           CALL OUT (A,LTO,NW)                                          
140   T = 0.                                                            
      DTT = DT                                                          
      NSTEPS = 1                                                        
145   CALL RUNINF                                                       
      CALL BUG ('AFT INITAL')                                           
C                                                                       
C THIS IS THE BEGINNING OF THE TIME STEP.                               
C                                                                       
150   IF ((MOD(NSTEPS, KTF) .EQ. 1) .AND. (NOTEMP .NE. 1))              
     1                             CALL TEMPFO (CA,N2P,N1H,CWORK)       
      IF (MOD(NSTEPS, KPF) .EQ. 1) CALL PSIFOU  (CA,N2P,N1H,CWORK)      
      IF (MOD(NSTEPS, KZF) .EQ. 1) CALL ZETAFO (CA,N2P,N1H,CWORK)       
      IF (MOD(NSTEPS, KVF) .EQ. 1) CALL VVELFO (CA,N2P,N1H,FX,CWORK)    
      IF (MOD (NSTEPS+1,KTA) .NE. 1) GO TO 211                          
C                                                                       
C TIME AVERAGE VORTICITY:                                               
C                                                                       
           DTT = DT                                                     
           T = T-0.5*DT                                                 
           CALL IN (A,LZO,NW)                                           
           NPTS = 2*N2                                                  
           NSKIP = 1                                                    
           MTRN = N1H                                                   
           MSKIP = N2P                                                  
           ISIGN = -1                                                   
           IEXSHF = NMAX2/NPTS                                          
           CALL XLOG2 (NPTS,LOG)                                        
           CALL ACAC (A,EX,WORK)                                        
           F = 0.125/REAL(N2)                                           
CX10517 CALL DOINX( 517 )                                               
CY10517 CALL DOINY( 517 )                                               
           DO 160 J = 1,N1H                                             
                DO 160 I = 1,N2PD                                       
160                  A(I,J) = A(I,J)*F                                  
CX20518 CALL DOOUTX( 518 )                                              
CY20518 CALL DOOUTY( 518 )                                              
           CALL OUT (A,LZO,NW)                                          
           ISIGN = 1                                                    
           MTRN = N1Q                                                   
           KZN = LZN                                                    
           KZO = LZO                                                    
CX10519 CALL DOINX( 519 )                                               
CY10519 CALL DOINY( 519 )                                               
           DO 180 K = 1,2                                               
                CALL IN (A,KZO,NWH)                                     
                CALL IN (A(1,K2Q),KZN,NWH)                              
                DO 170 J = 1,N1Q                                        
                     K2QJ = J+K2QM                                      
                     DO 170 I = 1,N2PD                                  
170                       A(I,J) = A(I,J)+0.5*A(I,K2QJ)                 
                CALL OUT (A,KZN,NWH)                                    
                CALL ACAC (A,EX,WORK)                                   
                DO 175 J = 1,N1Q                                        
                     DO 175 I = 1,N2PD                                  
175                       A(I,J) = 0.5*A(I,J)                           
                CALL OUT (A,KZO,NWH)                                    
                KZO = KZO+NWH                                           
180             KZN = KZN+NWH                                           
CX20520 CALL DOOUTX( 520 )                                              
CY20520 CALL DOOUTY( 520 )                                              
      IF (NOTEMP .EQ. 1) GO TO 211                                      
C                                                                       
C TIME AVERAGE TEMPERATURE:                                             
C                                                                       
           ISIGN = -1                                                   
           MTRN = N1H                                                   
           CALL IN (A,LTO,NW)                                           
           CALL ACAC (A,EX,WORK)                                        
           F = 0.125/REAL(N2)                                           
CX10521 CALL DOINX( 521 )                                               
CY10521 CALL DOINY( 521 )                                               
           DO 190 J = 1,N1H                                             
                DO 190 I = 1,N2PD                                       
190                  A(I,J) = F*A(I,J)                                  
CX20522 CALL DOOUTX( 522 )                                              
CY20522 CALL DOOUTY( 522 )                                              
           CALL OUT (A,LTO,NW)                                          
           ISIGN = 1                                                    
           MTRN = N1Q                                                   
           KTO = LTO                                                    
           KTN = LTN                                                    
CX10523 CALL DOINX( 523 )                                               
CY10523 CALL DOINY( 523 )                                               
           DO 210 K = 1,2                                               
                CALL IN (A,KTO,NWH)                                     
                CALL IN (A(1,K2Q),KTN,NWH)                              
                DO 200 J = 1,N1Q                                        
                     K2QJ = J+K2QM                                      
                     DO 200 I = 1,N2PD                                  
200                       A(I,J) = A(I,J)+0.5*A(I,K2QJ)                 
                CALL OUT (A,KTN,NWH)                                    
                CALL ACAC (A,EX,WORK)                                   
                DO 205 J = 1,N1Q                                        
                     DO 205 I = 1,N2PD                                  
205                       A(I,J) = 0.5*A(I,J)                           
                CALL OUT (A,KTO,NWH)                                    
                KTO = KTO+NWH                                           
210             KTN = KTN+NWH                                           
CX20524 CALL DOOUTX( 524 )                                              
CY20524 CALL DOOUTY( 524 )                                              
C                                                                       
C THE BEGINING OF A NORMAL TIME STEP                                    
C (U IS X VELOCITY AND V IS Y VELOCITY)                                 
C FIND (-U) = D(PSI)/DY                                                 
C                                                                       
211   CONTINUE                                                          
C                                                                       
      IF (NSTEPS.LT.3) THEN                                             
         CALL IN (CA,LPN,NW)                                            
         CALL PSIPHY (CA,N2P,N1H,EX,G,U,WORK,CWORK3,+1)                 
      ELSE                                                              
         CALL IN (CA,LPN,NW)                                            
         IF (MOD (NSTEPS , KPP) .EQ. 1) CALL PSIPHY (CA,N2P,N1H,EX,G,   
     1                                              U,WORK,CWORK3,-1)   
      END IF                                                            
C                                                                       
      NPTS = 2*N2                                                       
      NSKIP = 1                                                         
      MTRN = N1H                                                        
      MSKIP = N2P                                                       
      ISIGN = 1                                                         
      CALL XLOG2 (NPTS,LOG)                                             
      IEXSHF = NMAX2/NPTS                                               
      CALL IN (CA,LPN,NW)                                               
      F = 0.                                                            
CX10525 CALL DOINX( 525 )                                               
CY10525 CALL DOINY( 525 )                                               
      DO 230 I = 1,N2                                                   
           DO 220 J = 1,N1H                                             
220             CA(I,J) = CMPLX(-F*AIMAG(CA(I,J)), F*REAL(CA(I,J)))     
230        F = F+1.                                                     
CX20526 CALL DOOUTX( 526 )                                              
CY20526 CALL DOOUTY( 526 )                                              
CX10527 CALL DOINX( 527 )                                               
CY10527 CALL DOINY( 527 )                                               
      DO 235 J = 1,N1H                                                  
235        CA(N2P,J) = 0.                                               
CX20528 CALL DOOUTX( 528 )                                              
CY20528 CALL DOOUTY( 528 )                                              
      CALL SCSC (CA,EX,WORK)                                            
      CALL OUT (CA,LPY,NW)                                              
      IF (KSM .EQ. 1) CALL MEAN (CA,U,G)                                
      IF (MOD(NSTEPS, KUF) .EQ. 1) CALL UVELFO (CA,N2P,N1H,G,EX,CWORK)  
      CALL IN (CA,LPN,NW)                                               
C                                                                       
C FIND  V = D(PSI)/DX                                                   
C                                                                       
      F = 0.                                                            
CX10529 CALL DOINX( 529 )                                               
CY10529 CALL DOINY( 529 )                                               
      DO 250 J = 1,N1H                                                  
           DO 240 I = 1,N2P                                             
240             CA(I,J) = CMPLX(-F*AIMAG(CA(I,J)), F*REAL(CA(I,J)))     
250        F = F+FX                                                     
CX20530 CALL DOOUTX( 530 )                                              
CY20530 CALL DOOUTY( 530 )                                              
      CALL ACAC (CA,EX,WORK)                                            
      CALL OUT (CA,LPN,NW)                                              
      IF ((MOD(NSTEPS, KPP) .EQ. 1) .AND. (NOTEMP .NE. 1))              
     1                              CALL IN (CA,LPN,NW)                 
      NPTS = 2*N2                                                       
      NSKIP = 1                                                         
      MTRN = N1H                                                        
      MSKIP = N2P                                                       
      ISIGN = 1                                                         
      CALL XLOG2 (NPTS,LOG)                                             
      IEXSHF = NMAX2/NPTS                                               
      IF (NOTEMP .EQ. 1) GO TO 271                                      
C                                                                       
C T(TOTAL) = T-BETA*Y                                                   
C WHERE BETA < 0 MEANS   STABLE STRATIFICATION AND                      
C       BETA > 0 MEANS UNSTABLE STRATIFICATION                          
C                                                                       
           KTO = LTO                                                    
           DTB = 0.5*BETA*DTT                                           
CX10531 CALL DOINX( 531 )                                               
CY10531 CALL DOINY( 531 )                                               
           DO 270 J = 1,N1H                                             
                CALL IN (CWORK,KTO,N2PD)                                
                F = -REAL(J-1)*FX*DTT*0.25                              
                DO 260 I = 1,N2P                                        
260                  CWORK(I) = CWORK(I)*CMPLX(1., U(I)*F)+DTB*CA(I,J)  
                CALL OUT (CWORK,KTO,N2PD)                               
270             KTO = KTO+N2PD                                          
CX20532 CALL DOOUTX( 532 )                                              
CY20532 CALL DOOUTY( 532 )                                              
271   CALL IN (CA,LZN,NW)                                               
C                                                                       
C FIND DZ/DY:                                                           
C                                                                       
      F = 0.                                                            
CX10533 CALL DOINX( 533 )                                               
CY10533 CALL DOINY( 533 )                                               
      DO 290 I = 1,N2                                                   
           DO 280 J = 1,N1H                                             
280             CA(I,J) = CMPLX(-F*AIMAG(CA(I,J)), F*REAL(CA(I,J)))     
290        F = F+1.                                                     
CX20534 CALL DOOUTX( 534 )                                              
CY20534 CALL DOOUTY( 534 )                                              
CX10535 CALL DOINX( 535 )                                               
CY10535 CALL DOINY( 535 )                                               
      DO 295 J = 1,N1H                                                  
295        CA(N2P,J) = 0.                                               
CX20536 CALL DOOUTX( 536 )                                              
CY20536 CALL DOOUTY( 536 )                                              
      CALL SCSC (CA,EX,WORK)                                            
      CALL OUT  (CA,LZY,NW)                                             
      CALL IN   (CA,LZN,NW)                                             
      CALL ACAC (CA,EX,WORK)                                            
CX10537 CALL DOINX( 537 )                                               
CY10537 CALL DOINY( 537 )                                               
      DO 300 J = 1,N1H                                                  
           DO 300 I = 1,N2P                                             
300             CA(I,J) = CA(I,J)*0.5                                   
CX20538 CALL DOOUTX( 538 )                                              
CY20538 CALL DOOUTY( 538 )                                              
      CALL OUT (CA,LZN,NW)                                              
      NPTS = N1                                                         
      NSKIP = N2P                                                       
      MTRN = N2P                                                        
      MSKIP = 1                                                         
      ISIGN = 1                                                         
      CALL XLOG2 (NPTS,LOG)                                             
      IEXSHF = NMAX2/NPTS                                               
      IF (MOD(NSTEPS, KZP) .EQ. 1) CALL ZETAPH (CA,N2PD,N1H,EX,WORK,-1) 
      IF (MOD(NSTEPS, KZP) .EQ. 1) CALL IN (CA,LZN,NW)                  
C                                                                       
C FIND DZ/DX:                                                           
C                                                                       
      F = 0.                                                            
CX10539 CALL DOINX( 539 )                                               
CY10539 CALL DOINY( 539 )                                               
      DO 320 J = 1,N1H                                                  
           DO 310 I = 1,N2P                                             
310             CA(I,J) = CMPLX(-F*AIMAG(CA(I,J)), F*REAL(CA(I,J)))     
320        F = F+FX2                                                    
CX20540 CALL DOOUTX( 540 )                                              
CY20540 CALL DOOUTY( 540 )                                              
      CALL CSR (CA,EX)                                                  
      CALL OUT (CA,LZX,NW)                                              
      CALL IN  (CA,LZY,NW)                                              
      CALL CSR (CA,EX)                                                  
      CALL OUT (CA,LZY,NW)                                              
      CALL IN  (CA,LPN,NW)                                              
      CALL CSR (CA,EX)                                                  
      CALL OUT (CA,LPN,NW)                                              
      CALL IN  (CA,LPY,NW)                                              
      CALL CSR (CA,EX)                                                  
      CALL OUT (CA,LPY,NW)                                              
C                                                                       
C WE ARE IN PHYSICAL SPACE NOW - FIND  (U,V).GRAD (Z):                  
C                                                                       
      IF (MOD(NSTEPS, KUP) .EQ. 1) CALL UVELPH (A,N2PD,N1H,WORK,-1)     
      IF (MOD(NSTEPS, KVP) .EQ. 1) CALL VVELPH (A,N2PD,N1H,WORK,-1)     
      KPY = LPY                                                         
      KPN = LPN                                                         
      KZY = LZY                                                         
      KZX = LZX                                                         
      UKSM1F = REAL(2-KSM)                                              
CX10541 CALL DOINX( 541 )                                               
CY10541 CALL DOINY( 541 )                                               
      DO 340 K = 1,4                                                    
           CALL IN (A,KPY,NWQ)                                          
           CALL IN (A(1,K1Q),KPN,NWQ)                                   
           CALL IN (A(1,K2Q),KZY,NWQ)                                   
           CALL IN (A(1,K3Q),KZX,NWQ)                                   
           DO 330 I = 1,N2PD                                            
                IPH = (I+1)/2                                           
                GG = G(IPH)                                             
                UU = U(IPH)*UKSM1F                                      
                DO 330 J = 1,N1E                                        
                     A(I,J) = A(I,J)*GG + UU                            
                     A(I,K2QM+J) = A(I,J)*A(I,K3QM+J) -                 
     1                             A(I,K1QM+J)*A(I,K2QM+J)*GG           
330        CONTINUE                                                     
           CALL OUT (A,KPY,NWQ)                                         
           CALL OUT (A(1,K2Q),KZY,NWQ)                                  
           KPN = KPN+NWQ                                                
           KPY = KPY+NWQ                                                
           KZX = KZX+NWQ                                                
340        KZY = KZY+NWQ                                                
CX20542 CALL DOOUTX( 542 )                                              
CY20542 CALL DOOUTY( 542 )                                              
C                                                                       
C [Z(OLD) + (MEAN -U) * (-DZ(OLD)/DX)] + (U,V).GRAD(Z)                  
C                                                                       
      ISIGN = -1                                                        
      CALL IN  (A,LZY,NW)                                               
      CALL RCS (A,EX)                                                   
      CALL OUT (A,LZY,NW)                                               
      KZO = LZO                                                         
      KZN = LZN                                                         
      KZY = LZY                                                         
      DTFAC = 0.125*DTT/REAL(N1)                                        
      F1 = DTT*FX*0.25                                                  
      F = F1                                                            
CX10543 CALL DOINX( 543 )                                               
CY10543 CALL DOINY( 543 )                                               
      DO 360 K = 1,2                                                    
           CALL IN  (CA,KZO,NWH)                                        
           CALL IN  (CA(1,K2Q),KZN,NWH)                                 
           CALL OUT (CA(1,K2Q),KZO,NWH)                                 
           CALL IN  (CA(1,K2Q),KZY,NWH)                                 
           DO 350 J = 1,N1Q                                             
                K2QJ = J+K2QM                                           
                F = F-F1                                                
                DO 350 I = 1,N2P                                        
350                  CA(I,J) = CA(I,J)*CMPLX(1.,F*U(I))+DTFAC*CA(I,K2QJ)
           CALL OUT (CA,KZN,NWH)                                        
           KZN = KZN+NWH                                                
           KZO = KZO+NWH                                                
360        KZY = KZY+NWH                                                
CX20544 CALL DOOUTX( 544 )                                              
CY20544 CALL DOOUTY( 544 )                                              
      IF (NOTEMP .EQ. 1) GO TO 441                                      
C                                                                       
C FIND DT/DY:                                                           
C                                                                       
           NPTS = 2*N2                                                  
           NSKIP = 1                                                    
           MTRN = N1H                                                   
           MSKIP = N2P                                                  
           ISIGN = 1                                                    
           CALL XLOG2 (NPTS,LOG)                                        
           IEXSHF = NMAX2/NPTS                                          
           CALL IN (CA,LTN,NW)                                          
           F = 0.                                                       
CX10545 CALL DOINX( 545 )                                               
CY10545 CALL DOINY( 545 )                                               
           DO 380 I = 1,N2                                              
                DO 370 J = 1,N1H                                        
370                  CA(I,J) = CMPLX(-F*AIMAG(CA(I,J)), F*REAL(CA(I,J)))
380             F = F+1.                                                
CX20546 CALL DOOUTX( 546 )                                              
CY20546 CALL DOOUTY( 546 )                                              
CX10547 CALL DOINX( 547 )                                               
CY10547 CALL DOINY( 547 )                                               
           DO 385 J = 1,N1H                                             
385             CA(N2P,J) = 0.                                          
CX20548 CALL DOOUTX( 548 )                                              
CY20548 CALL DOOUTY( 548 )                                              
           CALL SCSC (CA,EX,WORK)                                       
           CALL OUT  (CA,LZY,NW)                                        
           CALL IN   (CA,LTN,NW)                                        
           CALL ACAC (CA,EX,WORK)                                       
           CALL OUT  (CA,LTN,NW)                                        
           NPTS = N1                                                    
           NSKIP = N2P                                                  
           MTRN = N2P                                                   
           MSKIP = 1                                                    
           CALL XLOG2 (NPTS,LOG)                                        
           IEXSHF = NMAX2/NPTS                                          
           IF (NSTEPS.EQ.1) THEN                                        
              CALL TEMPHY (A,N2PD,N1H,EX,WORK,0)                        
           ELSE                                                         
              IF ((MOD(NSTEPS, KTP).EQ.1).OR.(MOD(NSTEPS, KPP).EQ.1))   
     1        CALL TEMPHY (A,N2PD,N1H,EX,WORK,-1)                       
           END IF                                                       
C                                                                       
C TEMPERATURE INFLUENCE:                                                
C                                                                       
           F = -FX                                                      
           KZN = LZN                                                    
           KZX = LZX                                                    
           KTN = LTN                                                    
           DTA = 0.5*ALPHA*DTT                                          
CX10549 CALL DOINX( 549 )                                               
CY10549 CALL DOINY( 549 )                                               
           DO 400 K = 1,2                                               
                CALL IN (CA,KZN,NWH)                                    
                CALL IN (CA(1,K2Q),KTN,NWH)                             
                DO 390 J = 1,N1Q                                        
                     K2QJ = J+K2QM                                      
                     F = F+FX                                           
                     DO 390 I = 1,N2P                                   
                          CA(I,K2QJ) = CMPLX(-F*AIMAG(CA(I,K2QJ)),      
     1                                            F*REAL(CA(I,K2QJ)))   
390                       CA(I,J) = CA(I,J)+DTA*CA(I,K2QJ)              
                CALL OUT (CA,KZN,NWH)                                   
                CALL OUT (CA(1,K2Q),KZX,NWH)                            
                KZX = KZX+NWH                                           
                KTN = KTN+NWH                                           
400             KZN = KZN+NWH                                           
CX20550 CALL DOOUTX( 550 )                                              
CY20550 CALL DOOUTY( 550 )                                              
           CALL IN  (A,LZX,NW)                                          
           CALL CSR (A,EX)                                              
           CALL OUT (A,LZX,NW)                                          
           CALL IN  (A,LZY,NW)                                          
           CALL CSR (A,EX)                                              
           CALL OUT (A,LZY,NW)                                          
           KPN = LPN                                                    
           KPY = LPY                                                    
           KZX = LZX                                                    
           KZY = LZY                                                    
CX10551 CALL DOINX( 551 )                                               
CY10551 CALL DOINY( 551 )                                               
C                                                                       
C IN PHYSICAL SPACE - FIND  (U,V).GRAD(T):                              
C                                                                       
           DO 420 K = 1,4                                               
                CALL IN (A,KPY,NWQ)                                     
                CALL IN (A(1,K1Q),KPN,NWQ)                              
                CALL IN (A(1,K2Q),KZY,NWQ)                              
                CALL IN (A(1,K3Q),KZX,NWQ)                              
                DO 410 I = 1,N2PD                                       
                     IPH = (I+1)/2                                      
                     GG = G(IPH)                                        
                     DO 410 J = 1,N1E                                   
                          A(I,J) = A(I,J)*A(I,K3QM+J) -                 
     1                             A(I,K1QM+J)*A(I,K2QM+J)*GG           
410             CONTINUE                                                
                CALL OUT (A,KZY,NWQ)                                    
                KPN = KPN+NWQ                                           
                KPY = KPY+NWQ                                           
                KZX = KZX+NWQ                                           
420             KZY = KZY+NWQ                                           
CX20552 CALL DOOUTX( 552 )                                              
CY20552 CALL DOOUTY( 552 )                                              
           ISIGN = -1                                                   
           CALL IN  (A,LZY,NW)                                          
           CALL RCS (A,EX)                                              
           CALL OUT (A,LZY,NW)                                          
           KTN = LTN                                                    
           KTO = LTO                                                    
           KZY = LZY                                                    
CX10553 CALL DOINX( 553 )                                               
CY10553 CALL DOINY( 553 )                                               
           DO 440 K = 1,2                                               
                CALL IN (CA,KTO,NWH)                                    
                CALL IN (A(1,K2Q),KTN,NWH)                              
                DO 425 J = K2Q,N1H                                      
                     DO 425 I = 1,N2PD                                  
425                       A(I,J) = 0.5*A(I,J)                           
                CALL OUT (A(1,K2Q),KTO,NWH)                             
                CALL IN  (CA(1,K2Q),KZY,NWH)                            
                DO 430 J = 1,N1Q                                        
                     K2QJ = J+K2QM                                      
                     DO 430 I = 1,N2P                                   
430                       CA(I,J) = CA(I,J)+DTFAC*CA(I,K2QJ)            
                CALL OUT (CA,KTN,NWH)                                   
                KTN = KTN+NWH                                           
                KTO = KTO+NWH                                           
440             KZY = KZY+NWH                                           
CX20554 CALL DOOUTX( 554 )                                              
CY20554 CALL DOOUTY( 554 )                                              
441   DTF = 0.25*DTT*FX                                                 
      F = -DTF                                                          
      KZN = LZN                                                         
      KTN = LTN                                                         
CX10555 CALL DOINX( 555 )                                               
CY10555 CALL DOINY( 555 )                                               
      DO 460 K = 1,2                                                    
           IF (NOTEMP .NE. 1) CALL IN (CA,KTN,NWH)                      
           CALL IN (CA(1,K2Q),KZN,NWH)                                  
           DO 450 J = 1,N1Q                                             
                K2QJ = J+K2QM                                           
                F = F+DTF                                               
                IF (NOTEMP.EQ.1) GO TO 449                              
                DO 448 I = 1,N2P                                        
                     WORK(I) = F*U(I)                                   
                     CA(I,J) = CMPLX (1.0, -WORK(I)) * CA(I,J) /        
     1                               (1.0 + WORK(I)*WORK(I))            
448             CONTINUE                                                
  449 CONTINUE                                                          
      DO      450      I      =      1,N2P                              
                     WORK(I) = F*U(I)                                   
                     CA(I,K2QJ) = CMPLX (1.0, -WORK(I)) * CA(I,K2QJ) /  
     1                                  (1.0 + WORK(I)*WORK(I))         
450        CONTINUE                                                     
           IF (NOTEMP .NE. 1) CALL OUT (CA,KTN,NWH)                     
           CALL OUT (CA(1,K2Q),KZN,NWH)                                 
           KTN = KTN+NWH                                                
460        KZN = KZN+NWH                                                
CX20556 CALL DOOUTX( 556 )                                              
CY20556 CALL DOOUTY( 556 )                                              
      NPTS = 2*N2                                                       
      NSKIP = 1                                                         
      MTRN = N1H                                                        
      MSKIP = N2P                                                       
      ISIGN = -1                                                        
      IEXSHF = NMAX2/NPTS                                               
      CALL XLOG2 (NPTS,LOG)                                             
      CALL IN   (A,LZN,NW)                                              
      CALL ACAC (A,EX,WORK)                                             
      CALL OUT  (A,LZN,NW)                                              
      IF (NOTEMP .EQ. 1) GO TO 461                                      
           CALL IN   (A,LTN,NW)                                         
           CALL ACAC (A,EX,WORK)                                        
           CALL OUT  (A,LTN,NW)                                         
461   KZN = LZN                                                         
      KPN = LPN                                                         
      F2 = 1./(4.*REAL(N2))                                             
      F1 = RNU*DTT                                                      
C                                                                       
C FIND PRESSURE NOW OR NEVER (IF REQUIRED...)                           
C                                                                       
      CALL IN (EVAL,LEIG,NWEIG)                                         
      LVECT = N2M                                                       
      LSKIP = 2                                                         
      NVECT = 2                                                         
      NVSKIP = 1                                                        
      FXSQ = FX*FX                                                      
CX10557 CALL DOINX( 557 )                                               
CY10557 CALL DOINY( 557 )                                               
C                                                                       
C UPDATE VORTICITY AND FIND THE NEW STREAMFUNCTION:                     
C                                                                       
      DO 480 J = 1,N1H                                                  
           CALL IN (CWORK,KZN,N2PD)                                     
           F = REAL((J-1)*(J-1))*FXSQ                                   
           DO 470 I = 2,N2                                              
                CWORK(I) = CWORK(I)*F2 / (1.0-F1*(EVAL(I-1)-F))         
                CWORK2(I) = CWORK(I) / (EVAL(I-1)-F)                    
470        CONTINUE                                                     
           CWORK(1) = CWORK(1)*F2/(1.0+F1*F)                            
           CWORK2(1) = 0.0                                              
           IF (F.NE.0.0) CWORK2(1) = -CWORK(1)/F                        
           CWORK(N2P) = 0.0                                             
           CWORK2(N2P) = 0.0                                            
           CALL OUT (CWORK,KZN,N2PD)                                    
           CALL OUT (CWORK2,KPN,N2PD)                                   
           KPN = KPN+N2PD                                               
480        KZN = KZN+N2PD                                               
CX20558 CALL DOOUTX( 558 )                                              
CY20558 CALL DOOUTY( 558 )                                              
C                                                                       
C DE-ALIASE UPDATED FIELDS IF REQUESTED:                                
C                                                                       
      IF (NOALS.EQ.1) CALL DEALSE (CA,N2P,N1H,LZN)                      
      IF (NOALS.EQ.1) CALL DEALSE (CA,N2P,N1H,LPN)                      
      IF (NOTEMP .EQ. 1) GO TO 501                                      
C                                                                       
C UPDATE TEMPERATURE:                                                   
C                                                                       
           IF (NOALS.EQ.1) CALL IN (EVAL,LEIG,NWEIG)                    
           KTN = LTN                                                    
           F1 = RKAPA*DTT                                               
CX10559 CALL DOINX( 559 )                                               
CY10559 CALL DOINY( 559 )                                               
           DO 500 J = 1,N1H                                             
                CALL IN (CWORK,KTN,N2PD)                                
                F = REAL((J-1)*(J-1))*FXSQ                              
                DO 490 I = 2,N2                                         
                     CWORK(I) = CWORK(I)*F2 / (1.0-F1*(EVAL(I-1)-F))    
490             CONTINUE                                                
                CWORK(1) = CWORK(1)*F2 / (1.0+F1*F)                     
                CWORK(N2P) = 0.0                                        
                CALL OUT (CWORK,KTN,N2PD)                               
500             KTN = KTN+N2PD                                          
CX20560 CALL DOOUTX( 560 )                                              
CY20560 CALL DOOUTY( 560 )                                              
C                                                                       
C DE-ALIASE IF REQUESTED:                                               
C                                                                       
           IF (NOALS.EQ.1) CALL DEALSE (CA,N2P,N1H,LTN)                 
C                                                                       
C AT LAST, THE END OF THE TIME STEP:                                    
C                                                                       
501   DTT = DT+DT                                                       
      T = T+DT                                                          
      NSTEPS = NSTEPS+1                                                 
      IF (NUMBER .GE. NSTEPS) GO TO 150                                 
C                                                                       
      WRITE (6,1010) TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT,  
     1  TXSHUF,TUSHUF,TXPRNT,TXLOG2                                     
      WRITE (6,1012) EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT                 
      WRITE (6,1020) NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT,  
     1  NXSHUF,NUSHUF,NXPRNT,NXLOG2                                     
C                                                                       
      CALL BUG ('AFTER  ALL')                                           
C                                                                       
 1010 FORMAT(/2X,42HSUBROUTINE    IN, EX-TIME FOR FIRST CALL =,1PE13.6, 
     .        8H SECONDS/                                               
     1        2X,42HSUBROUTINE   OUT, EX-TIME FOR FIRST CALL =,1PE13.6/ 
     2        2X,42HSUBROUTINE   RCS, EX-TIME FOR FIRST CALL =,1PE13.6/ 
     3        2X,42HSUBROUTINE   CSR, EX-TIME FOR FIRST CALL =,1PE13.6/ 
     4        2X,42HSUBROUTINE  SCSC, EX-TIME FOR FIRST CALL =,1PE13.6/ 
     5        2X,42HSUBROUTINE  ACAC, EX-TIME FOR FIRST CALL =,1PE13.6/ 
     6        2X,42HSUBROUTINE FTRVMT EX-TIME FOR FIRST CALL =,1PE13.6/ 
     7        2X,42HSUBROUTINE  SHUF, EX-TIME FOR FIRST CALL =,1PE13.6/ 
     8        2X,42HSUBROUTINE UNSHUF EX-TIME FOR FIRST CALL =,1PE13.6/ 
     9        2X,42HSUBROUTINE RPRINT EX-TIME FOR FIRST CALL =,1PE13.6/ 
     1        2X,42HSUBROUTINE XLOG2, EX-TIME FOR FIRST CALL =,1PE13.6) 
 1012 FORMAT(/2X,42HSUBROUTINE   RCS, EX-TIME LAST TIMESTEP  =,1PE13.6, 
     .        8H SECONDS/                                               
     1        2X,42HSUBROUTINE   CSR, EX-TIME LAST TIMESTEP  =,1PE13.6/ 
     2        2X,42HSUBROUTINE  SCSC, EX-TIME LAST TIMESTEP  =,1PE13.6/ 
     3        2X,42HSUBROUTINE  ACAC, EX-TIME LAST TIMESTEP  =,1PE13.6/ 
     4        2X,42HSUBROUTINE FTRVMT EX-TIME LAST TIMESTEP  =,1PE13.6) 
 1020 FORMAT(/2X,42HSUBROUTINE    IN, CUMULATIVE CALLS TO    =,I12/     
     1        2X,42HSUBROUTINE   OUT, CUMULATIVE CALLS TO    =,I12/     
     2        2X,42HSUBROUTINE   RCS, CUMULATIVE CALLS TO    =,I12/     
     3        2X,42HSUBROUTINE   CSR, CUMULATIVE CALLS TO    =,I12/     
     4        2X,42HSUBROUTINE  SCSC, CUMULATIVE CALLS TO    =,I12/     
     5        2X,42HSUBROUTINE  ACAC, CUMULATIVE CALLS TO    =,I12/     
     6        2X,42HSUBROUTINE FTRVMT CUMULATIVE CALLS TO    =,I12/     
     7        2X,42HSUBROUTINE  SHUF, CUMULATIVE CALLS TO    =,I12/     
     8        2X,42HSUBROUTINE UNSHUF CUMULATIVE CALLS TO    =,I12/     
     9        2X,42HSUBROUTINE RPRINT CUMULATIVE CALLS TO    =,I12/     
     1        2X,42HSUBROUTINE XLOG2, CUMULATIVE CALLS TO    =,I12)     
C                                                                       
      CALL CPUTIM (CPU2)                                                
      CALL ELAPSE (TCLCK2)                                              
C                                                                       
      TWALL = MOD(TCLCK2-TCLCK1+1000000.0,1000000.0)                    
      WRITE (6,1040) TWALL                                              
 1040 FORMAT(/2X,50HMEASURE OF WALL-CLOCK TIME (ENVIRONMENT SENSITIVE)/ 
     1 2X,20HFOR EXECUTION ONLY =,1PE13.6,1X,7HSECONDS)                 
      WRITE(1,*)                                                        
      WRITE(1,*) 'ELAPSED CPU TIME IN SECONDS: ',CPU2-CPU1              
C MFLOPS RATE = MFLOP 1 CPU CRAY X-MP/CPU TIME                          
      WRITE(1,*) 'MFLOPS RATE: ', 1657.863/(CPU2-CPU1)                  
      WRITE(1,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ',TWALL            
      WRITE(1,9992)                                                     
      WRITE(1,9995) (CDAT(I),I=1,15)                                    
      VALID = 1                                                         
CX10561 CALL DOINX( 561 )                                               
CY10561 CALL DOINY( 561 )                                               
      DO  7  I=1,15                                                     
      IF (I.EQ.3 .OR. I.EQ.13) THEN                                     
         IF (ABS(CDAT(I)) .GT. ABS(VDAT(I))) VALID = 0                  
      ELSE                                                              
         IF (ABS(VDAT(I)-CDAT(I)) .GT. ABS(.01*VDAT(I))) VALID = 0      
      ENDIF                                                             
    7 CONTINUE                                                          
CX20562 CALL DOOUTX( 562 )                                              
CY20562 CALL DOOUTY( 562 )                                              
C                                                                       
      IF (VALID .EQ. 1) THEN                                            
         WRITE(1,9993)                                                  
      ELSE                                                              
         WRITE(1,9994)                                                  
      ENDIF                                                             
C                                                                       
9992  FORMAT(/1X,'VALIDATION PARAMETERS:'/)                             
9993  FORMAT(//1X,'RESULTS FOR THIS RUN ARE:   VALID')                  
9994  FORMAT(//1X,'RESULTS FOR THIS RUN ARE: INVALID')                  
9995  FORMAT(5E15.5)                                                    
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
      STOP                                                              
      END                                                               
C***********************************************************************
      SUBROUTINE DEALSE (CA,N2P,N1H,LOC)                                
      COMPLEX CA(N2P,1)                                                 
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
C                                                                       
C THIS FUNCTION GIVES THE INDEX OF THE HIGHEST MODE TO BE RETAINED      
C IN TERMS OF THE NUMBER OF MODES WHICH CAN BE RESOLVED:                
C                                                                       
      LSTMOD (MODEHI) = INT(0.666666666666666*REAL(MODEHI-1)+0.999)     
CX30003 CALL SBINX (3)                                                  
CY30003 CALL SBINY (3)                                                  
C                                                                       
      NW = 2*N2P*N1H                                                    
      CALL IN (CA,LOC,NW)                                               
      N1CUT = LSTMOD (N1H) +1                                           
      N2CUT = LSTMOD (N2P) +1                                           
CX10563 CALL DOINX( 563 )                                               
CY10563 CALL DOINY( 563 )                                               
C                                                                       
C THIS DEALIASING PROCEDURE HAS TO BE IMPROVED IF NONLINEAR MAPPING     
C IN Y-DIRECTION IS USED:                                               
C                                                                       
      DO 10 J=N1CUT,N1H                                                 
            DO 10 I=1,N2P                                               
   10             CA(I,J) = 0.0                                         
CX20564 CALL DOOUTX( 564 )                                              
CY20564 CALL DOOUTY( 564 )                                              
CX10565 CALL DOINX( 565 )                                               
CY10565 CALL DOINY( 565 )                                               
      DO 20 J=1,N1H                                                     
            DO 20 I=N2CUT,N2P                                           
   20              CA(I,J) = 0.0                                        
CX20566 CALL DOOUTX( 566 )                                              
CY20566 CALL DOOUTY( 566 )                                              
      CALL OUT (CA,LOC,NW)                                              
CX40004 CALL SBOUTX (4)                                                 
CY40004 CALL SBOUTY (4)                                                 
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C     FGPS  2-MAY-99 SET KPP = KZP = KUP = KVP = 249 AND NUMBER = 250   
C     FGPS 18-APR-88 SET KPP = KZP = KUP = KVP = 499 AND NUMBER = 500   
C                                                                       
      SUBROUTINE START (NX,NY)                                          
C                                                                       
C THIS SUBROUTINE IS USED TO SPECIFY OPTIONS FOR PROGRAM OCEAN.         
C                                                                       
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1       KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,     
     2       KPRESS, NOALS                                              
CX30005 CALL SBINX (5)                                                  
CY30005 CALL SBINY (5)                                                  
C                                                                       
C N1 AND N2 ARE THE DIMENSIONS OF THE MATRIX. THEY HAVE TO BE POWERS    
C OF  2  AND GREATER THAN OR EQUAL TO  8                                
C                                                                       
      N1 = NX                                                           
      N2 = NY                                                           
C                                                                       
C KDATA SPECIFIES THE WAY DATA IS SUPPLIED TO THE PROGRAM.              
C                                                                       
C            KDATA = 1     MEANS THAT IT WILL BE IN PHYSICAL SPACE      
C                  = 2     MEANS THAT IT WILL BE IN FOURIER SPACE       
C                  = 3     MEANS THAT THIS IS A RESTART RUN (->READ TAPE
C                                                                       
      KDATA = 1                                                         
C                                                                       
C LMAPP SPECIFIES THE MAPPING IN THE Y DIRECTION.                       
C                                                                       
C            LMAPP = 1     MEANS LINEAR MAPPING (-> EIGSET NOT CALLED)  
C                  = 2     MEANS NON-LINEAR MAPPING                     
C                  = 3     MEANS RESTART RUN, NON-LINEAR MAPPING        
C                                                                       
      LMAPP = 1                                                         
C                                                                       
C NOTEMP EXCLUDES TEMPERATURE FIELD FROM CALCULATIONS WHEN SET TO 1.    
C                                                                       
      NOTEMP = 0                                                        
C                                                                       
C NOALS TURNS ON DEALIASING AT EACH TIME STEP WHEN SET TO 1:            
C (TO BE USED ONLY WITH LINEAR MAPPING)                                 
C                                                                       
      NOALS = 1                                                         
C                                                                       
C TIME AVERAGING IS PERFORMED EACH KTA TIME STEPS:                      
C                                                                       
      KTA = 50                                                          
C                                                                       
C THE FOLLOWING PARAMETERS SPECIFY WHEN TO CALL CORRESPONDING OUTPUT    
C SUBROUTINES. IF YOU WANT, FOR EXAMPLE, PSIPHY TO BE CALLED EVERY      
C 50 TIME STEPS, YOU SHOULD SET KPP TO 50.                              
C                                                                       
C WARNING: SETTING A PARAMETER TO 1 TURNS OFF THE SUBROUTINE. IT DOES   
C          NOT CAUSE IT TO BE CALLED EVERY TIME STEP.                   
C          YOU SHOULD NEVER SET ANY OF THEM TO 0, SINCE THEY ARE USED   
C          WITH THE FUNCTION MOD(NSTEPS,KXYZ).                          
C                                                                       
      KPP = 249                                                         
      KZP = 249                                                         
      KTP = 50                                                          
      KUP = 249                                                         
      KVP = 249                                                         
      KPF = 1                                                           
      KZF = 1                                                           
      KTF = 1                                                           
      KUF = 1                                                           
      KVF = 1                                                           
      KPRESS = 1                                                        
C                                                                       
C KTAPE SPECIFIES HOW WILL HISTORY TAPE BE WRITTEN:                     
C                                                                       
C          KTAPE = 1  WRITE ONLY EIG-STUFF ON TAPE                      
C                = 2  WRITE ONLY DATA ON TAPE                           
C                = 3  WRITE BOTH EIG-STUFF AND DATA ON TAPE             
C                                                                       
C ANY OTHER VALUE WILL PREVENT THE HISTORY TAPE FROM BEING WRITTEN.     
C                                                                       
      KTAPE = 0                                                         
C                                                                       
C NUMBER SPECIFIES THE NUMBER OF TIME STEPS IN A SINGLE RUN:            
C                                                                       
      NUMBER = 250                                                      
C                                                                       
C KSM GOVERNS THE USE OF THE MEAN FLOW OPTION.                          
C                                                                       
C                   KSM = 0       MEANS THAT U(MEAN) IS ZERO            
C                         1       MEANS THAT IT SHOULD BE COMPUTED      
C                         2       MEANS THAT IT IS USER-SPECIFIED       
C                                                                       
      KSM = 1                                                           
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
C                                                                       
C THAT WAS ALL. RETURN TO MOTHER OCEAN.                                 
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE GSET (GMAPP)                                           
      DIMENSION GMAPP(1)                                                
      COMMON /PI/ PI                                                    
      COMMON /NMBRS/ NSKIP(4), N2P, NDUMY1(3)                           
CX30007 CALL SBINX (7)                                                  
CY30007 CALL SBINY (7)                                                  
      F = PI/REAL(N2P-1)                                                
CX10567 CALL DOINX( 567 )                                               
CY10567 CALL DOINY( 567 )                                               
C                                                                       
C GSET DEFINES THE MAPPING FUNCTION.                                    
C                                                                       
      DO 20 I = 1,N2P                                                   
20         GMAPP(I) = F*REAL(I-1)                                       
CX20568 CALL DOOUTX( 568 )                                              
CY20568 CALL DOOUTY( 568 )                                              
CX40008 CALL SBOUTX (8)                                                 
CY40008 CALL SBOUTY (8)                                                 
      RETURN                                                            
      END                                                               
C***********************************************************************
C     FGPS 12-MAY-88 RESET "RLAMDA = 1.00"                              
C     FGPS 18-APR-88 RESET "DT" TO DT = 0.0001, OTHERWISE THE INCREASED 
C                    FUNCTION RESULTS IN NUMERICAL INSTABILITY          
C     FGPS 15-APR-88 INCREASED THE INITIAL AMPLITUDE PRESCRIBED FOR THE 
C                    FUNCTION BY *50                                    
C     FGPS 16-FEB-88 RESET "DT" TO DT = 0.001                           
C                                                                       
      SUBROUTINE DATAST (A,N2P,N1,GMAPP)                                
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /PI/ PI                                                    
      COMMON /LADDR/ LZN,LTN,LPN,LZO,LTO,LPY,LZX,LZY,LEIG,LDUMY1(3)     
      DIMENSION A(N2P,1), GMAPP(1)                                      
C                                                                       
C PROGRAM TO INITIALIZE THE RUN (EITHER IN PHYSICAL OR FOURIER SPACE)   
C PHYSICAL SPACE: KDATA=1    --> A (I,J) ... POINT (I-1, J-1)           
C FOURIER SPACE : KDATA=2    --> CA(I,J) ... MODE  (I-1, J-1)           
C                                                                       
C ALSO, IF YOU WANT TO SPECIFY THE MEAN FLOW (KSM=2), DO IT HERE.       
C                                                                       
      COMMON /MFLOW/ U(129)                                             
CX30009 CALL SBINX (9)                                                  
CY30009 CALL SBINY (9)                                                  
C                                                                       
      ALPHA = 10.0                                                      
C                                                                       
C UNSTABLE CASE:                                                        
      BETA = 1.0                                                        
C      RLAMDA = 1.35                                                    
C     WHEN CODE IS USED AS A BENCHMARK, SET RLAMDA = 1.00               
      RLAMDA = 1.00                                                     
      BOXX = 2.0*PI*RLAMDA                                              
      FX = 2.0*PI/BOXX                                                  
C                                                                       
C ASSUMING THAT GSET IS CALLED. IF LMAPP = 1, H SHOULD BE SET TO        
C A DEFINED POSITIVE NUMBER.                                            
C                                                                       
      H = GMAPP(N2P)                                                    
C                                                                       
C COMPUTE VISCOSITY AND HEAT CONDUCTIVITY FROM PRANDTL NUMBER           
C AND THE RATIO OF RAYLEIGH NUMBER TO ITS CRITICAL VALUE:               
C                                                                       
      PRANDT = 10.0                                                     
      RRATIO = 20.0                                                     
      W = ALPHA*BETA*PRANDT*4.0/(27.0*RRATIO)                           
      RNU = SQRT (ABS(W))                                               
      RKAPA = RNU/PRANDT                                                
C                                                                       
C DETERMINE DIFFUSION LIMIT FOR THE TIME STEP AND COMPARE IT TO         
C THE CHOSEN VALUE:                                                     
C                                                                       
      DT = 0.0001                                                       
      X = FX*REAL(N1/2-1)                                               
      Y = REAL(N2P-1)                                                   
      DTMAX = 1.0/((X*X+Y*Y) * MAX(RNU,RKAPA))                          
      WRITE(6,170) DTMAX,DT                                             
  170 FORMAT (' DIFFUSION CONSTRAINT ON DT: DTMAX = ',F10.6,/           
     1        ' ACTUAL DT CHOSEN FOR THIS RUN: DT = ',F10.6)            
      NW = N2P*N1                                                       
C                                                                       
C                                                                       
      EPS = 100.                                                        
CX10569 CALL DOINX( 569 )                                               
CY10569 CALL DOINY( 569 )                                               
      DO 20 J=1,N1                                                      
         X = REAL(J-1)*2.0*PI/REAL(N1)                                  
         DO 10 I=1,N2P                                                  
            Y = REAL(I-1)*PI/REAL(N2P-1)                                
            A(I,J) = EPS*SIN(X)*SIN(Y)                                  
   10    CONTINUE                                                       
   20 CONTINUE                                                          
CX20570 CALL DOOUTX( 570 )                                              
CY20570 CALL DOOUTY( 570 )                                              
C                                                                       
      WRITE(6,100)                                                      
  100 FORMAT (///' INITIAL STREAMFUNCTION IS   A SIN(X) SIN(Y),'/       
     1           ' WHERE  0 <= X <= 2*PI  AND  0 <= Y <= PI.'/          
     2           ' THIS IMPLIES TWO CONVECTION CELLS INITIALLY.'/       
     3           ' NOTE THAT THE PRESCRIBED INITIAL STREAM-'/           
     4           ' FUNCTION VS X, FOR Y=YMAX/2, IS PRINTED OUT'/        
     5           ' BELOW FOR NSTEPS=1 (TIME ZERO) AND ALSO FOR'/        
     6           ' THE NEXT STEP, NSTEPS=2'//)                          
      CALL OUT (A,LPN,NW)                                               
CX10571 CALL DOINX( 571 )                                               
CY10571 CALL DOINY( 571 )                                               
      DO 40 J=1,N1                                                      
         DO 30 I=1,N2P                                                  
            A(I,J) = 0.0                                                
   30    CONTINUE                                                       
   40 CONTINUE                                                          
CX20572 CALL DOOUTX( 572 )                                              
CY20572 CALL DOOUTY( 572 )                                              
      WRITE(6,200)                                                      
  200 FORMAT (' INITIAL TEMPERATURE PERTURBATION IS ZERO.'/             
     1        ' NOTE THAT THE INITIAL TEMPERATURE DROP VS'/             
     2        ' Y IS PRINTED OUT BELOW, WHERE INITIALLY'/               
     3        ' THIS IS INDEPENDENT OF X,  I.E., 1 <= IX'/              
     4        ' <= 64, AND WHERE THE Y REPRESENTATION IN'/              
     5        ' PHYSICAL SPACE SPANS 1 <= JY <= 129.   A'/              
     6        ' KEY VARIABLE, THE EVOLUTION OF WHICH CAN'/              
     7        ' MONITOR CODE PERFORMANCE, IS  "TEMP PHY('/              
     8        ' JY=65,IX=33)"; I.E., A MEASURE OF TEMPER-'/             
     9        ' ATURE IN PHYSICAL SPACE AT THE CENTERPOINT'/            
     1        ' OF THE BOUNDARY SEPARATING THE TWO CONVEC-'/            
     2        ' TION CELLS.')                                           
      CALL OUT (A,LTN,NW)                                               
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE RUNINF                                                 
      COMMON /OPTONS/ N1,N2,NUMBER,KDATA,LMAPP,NOTEMP,KTA,KTAPE,KSM,    
     1       KPP,KZP,KTP,KUP,KVP,KPF,KZF,KTF,KUF,KVF,KPRESS,NOALS       
      COMMON /DATA/ ALPHA,BETA,BOXX,H,RNU,RKAPA,DT,DTT,T,NSTEPS         
CX30011 CALL SBINX (11)                                                 
CY30011 CALL SBINY (11)                                                 
C                                                                       
C THIS SUBROUTINE TELLS US WHAT  OCEAN  WILL BE DOING:                  
C                                                                       
C*****MACHINE-DEPENDENT CONSTANT (REPLACE FOR YOUR MACHINE)             
C                                                                       
C*********************                                                  
      ZERO = 1.0E-200                                                   
C*********************                                                  
C                                                                       
      WRITE(6,100)                                                      
C                                                                       
      CALL CPUTIM (CPUTM)                                               
C      CALL DATE   (TODAY)                                              
C      CALL TIME   (NOW)                                                
      WRITE(6,200) CPUTM                                                
C                                                                       
      MOREST = NUMBER-NSTEPS+1                                          
      WRITE(6,300) N1,N2,MOREST,KTA                                     
      IF (KSM.EQ.1) WRITE(6,301)                                        
      IF (KSM.EQ.2) WRITE(6,302)                                        
C                                                                       
      GO TO (1,2,3), KDATA                                              
      WRITE(6,400)                                                      
      GO TO 10                                                          
1     WRITE(6,401)                                                      
      GO TO 10                                                          
2     WRITE(6,402)                                                      
      GO TO 10                                                          
3     WRITE(6,403)                                                      
C                                                                       
10    GO TO (11,123,123), LMAPP                                         
      WRITE(6,500)                                                      
      GO TO 20                                                          
11    WRITE(6,501)                                                      
      GO TO 20                                                          
123   WRITE(6,523)                                                      
C                                                                       
20    GO TO (21), NOTEMP                                                
      WRITE(6,600)                                                      
      GO TO 30                                                          
21    WRITE(6,601)                                                      
C                                                                       
30    GO TO (31,32,33), KTAPE                                           
      WRITE(6,700)                                                      
      GO TO 40                                                          
31    WRITE(6,701)                                                      
      GO TO 40                                                          
32    WRITE(6,702)                                                      
      GO TO 40                                                          
33    WRITE(6,703)                                                      
C                                                                       
40    IF (KPP.NE.1) WRITE(6,801) KPP                                    
      IF (KZP.NE.1) WRITE(6,802) KZP                                    
      IF (KTP.NE.1) WRITE(6,803) KTP                                    
      IF (KUP.NE.1) WRITE(6,804) KUP                                    
      IF (KVP.NE.1) WRITE(6,805) KVP                                    
      IF (KPF.NE.1) WRITE(6,806) KPF                                    
      IF (KZF.NE.1) WRITE(6,807) KZF                                    
      IF (KTF.NE.1) WRITE(6,808) KTF                                    
      IF (KUF.NE.1) WRITE(6,809) KUF                                    
      IF (KVF.NE.1) WRITE(6,810) KVF                                    
      IF (KPRESS.NE.1) WRITE(6,811) KPRESS                              
C                                                                       
      IF (NOALS.EQ.1) WRITE(6,901)                                      
C                                                                       
      REY = 1.0/(RNU+ZERO)                                              
      RH =  H*H*H*H*ABS(ALPHA*BETA)/(RNU*RKAPA+ZERO)                    
C                                                                       
C CRITICAL RAYLEIGH NUMBER IS   (PI**4) * 27/4:                         
C                                                                       
      RHCRIT = 657.5113648                                              
      RRATIO = RH/RHCRIT                                                
      PR  = RNU/(RKAPA+ZERO)                                            
      WRITE(6,1000)  BOXX,H,RNU,RKAPA,ALPHA,BETA,REY,RH,RRATIO,PR,DT,T  
CX40012 CALL SBOUTX (12)                                                
CY40012 CALL SBOUTY (12)                                                
C                                                                       
      RETURN                                                            
C                                                                       
C FORMAT STATEMENTS:                                                    
C                                                                       
100   FORMAT(/12X,' OOOOOO     CCCCCC    EEEEEEE    AAAAAA    N     N'/ 
     2        12X,'O      O   C      C   E         A      A   NN    N'/ 
     3        12X,'O      O   C          E         A      A   N N   N'/ 
     4        12X,'O      O   C          EEEE      AAAAAAAA   N  N  N'/ 
     5        12X,'O      O   C          E         A      A   N   N N'/ 
     6        12X,'O      O   C      C   E         A      A   N    NN'/ 
     7        12X,' OOOOOO     CCCCCC    EEEEEEE   A      A   N     N') 
C                                                                       
C200   FORMAT (////' THIS RUN BEGINS ON DATE ',A8,' AT TIME ',A8,/,     
C     1            ' CPU TIME USED SO FAR IS',F9.6,' SECONDS.')         
200    FORMAT (////' CPU TIME USED SO FAR IS',F9.6,' SECONDS.'//        
     1      'W A R N I N G ******** POSSIBLE UNDERFLOW ********'/       
     2      'THE STATEMENT "ZERO = 1.0E-200" APPEARS IN RUNINF.')       
C                                                                       
  300 FORMAT (//' THIS RUN USES ',I3,' BY ',I3,' GRID, AND DOES',I5,    
     1    ' TIME STEPS',/,' TIME AVERAGING IS DONE EVERY',I3,' STEPS')  
  301 FORMAT (' MEAN FLOW IS INCLUDED IN CALCULATIONS')                 
  302 FORMAT (' MEAN FLOW IS SPECIFIED BY USER')                        
C                                                                       
  400 FORMAT (' KDATA HAS ILLEGAL VALUE')                               
  401 FORMAT (' INITIAL DATA IN PHYSICAL SPACE')                        
  402 FORMAT (' INITIAL DATA IN FOURIER SPACE')                         
  403 FORMAT (/' THIS IS A RESTART RUN'/)                               
C                                                                       
  500 FORMAT (' LMAPP HAS AN ILLEGAL VALUE')                            
  501 FORMAT (' LINEAR MAPPING IN BOTH DIRECTIONS')                     
  523 FORMAT (' NON-LINEAR MAPPING IN Y DIRECTION')                     
C                                                                       
  600 FORMAT (' TEMPERATURE FIELD INCLUDED IN CALCULATIONS')            
  601 FORMAT (' CALCULATIONS DONE WITHOUT TEMPERATURE FIELD')           
C                                                                       
  700 FORMAT (/' NO HISTORY TAPE WILL BE WRITTEN'/)                     
  701 FORMAT (/' HISTORY TAPE WILL INCLUDE OPERATOR MATRICES ONLY'/)    
  702 FORMAT (/' HISTORY TAPE WILL NOT INCLUDE OPERATOR MATRICES'/)     
  703 FORMAT (/' HISTORY TAPE WILL INCLUDE ALL IMPORTANT ARRAYS'/)      
C                                                                       
  801 FORMAT (' OUTPUT PHYSICAL STREAMFUNCTION EVERY',I4,' TIME STEPS') 
  802 FORMAT (' OUTPUT PHYSICAL VORTICITY      EVERY',I4,' TIME STEPS') 
  803 FORMAT (' OUTPUT PHYSICAL TEMPERATURE    EVERY',I4,' TIME STEPS') 
  804 FORMAT (' OUTPUT PHYSICAL X VELOCITY     EVERY',I4,' TIME STEPS') 
  805 FORMAT (' OUTPUT PHYSICAL Y VELOCITY     EVERY',I4,' TIME STEPS') 
  806 FORMAT (' OUTPUT FOURIER  STREAMFUNCTION EVERY',I4,' TIME STEPS') 
  807 FORMAT (' OUTPUT FOURIER  VORTICITY      EVERY',I4,' TIME STEPS') 
  808 FORMAT (' OUTPUT FOURIER  TEMPERATURE    EVERY',I4,' TIME STEPS') 
  809 FORMAT (' OUTPUT FOURIER  X VELOCITY     EVERY',I4,' TIME STEPS') 
  810 FORMAT (' OUTPUT FOURIER  V VELOCITY     EVERY',I4,' TIME STEPS') 
  811 FORMAT (' OUTPUT PHYSICAL PRESSURE       EVERY',I4,' TIME STEPS') 
C                                                                       
  901 FORMAT (' DE-ALIASING TURNED ON.')                                
C                                                                       
 1000 FORMAT (////' PHYSICAL PARAMETERS FOR THIS RUN ARE:',/,           
     1     ' BOX DIMENSIONS          : ',F10.6,' BY ',F10.6,/,          
     2     ' VISCOSITY (NU)          : ',E24.14,/,                      
     3     ' HEAT CONDUCTIVITY (KAPA): ',E24.14,/,                      
     4     ' ALPHA AND BETA          : ',F10.6,' AND',F10.6,/,          
     5     ' REYNOLDS NUMBER         : ',E24.14,/,                      
     6     ' RAYLEIGH NUMBER         : ',E24.14,/,                      
     7     ' RATIO RAY/RAY(CRITICAL) : ',E24.14,/,                      
     8     ' PRANDTL  NUMBER         : ',E24.14,//,' THE TIME STEP IS', 
     9     F10.6,/,' WE BEGIN AT  T =',F10.6,///,1H ,60(2H+-),//)       
CX40012 CALL SBOUTX (12)                                                
CY40012 CALL SBOUTY (12)                                                
      END                                                               
C***********************************************************************
      SUBROUTINE SHUF (A,N2P,N1,WORK)                                   
      DIMENSION A(N2P,1),WORK(1)                                        
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      LOGICAL   LFIRST                                                  
      DATA      LFIRST /.TRUE./                                         
      SAVE   LFIRST                                                     
CX30013 CALL SBINX (13)                                                 
CY30013 CALL SBINY (13)                                                 
C                                                                       
      IF (LFIRST) CALL CPUTIM (T1)                                      
CX10573 CALL DOINX( 573 )                                               
CY10573 CALL DOINY( 573 )                                               
C                                                                       
      DO 10 J = 1,N1,2                                                  
            DO 20 I = 1,N2P                                             
                  II = I+I                                              
                  WORK(II-1) = A(I,J)                                   
   20             WORK(II) = A(I,J+1)                                   
            DO 10 I = 1,N2P                                             
                  A(I,J) = WORK(I)                                      
   10             A(I,J+1) = WORK(I+N2P)                                
CX20574 CALL DOOUTX( 574 )                                              
CY20574 CALL DOOUTY( 574 )                                              
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TXSHUF = T2-T1                                                  
      END IF                                                            
C                                                                       
      NXSHUF = NXSHUF+1                                                 
CX40014 CALL SBOUTX (14)                                                
CY40014 CALL SBOUTY (14)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE UNSHUF (A,N2P,N1,WORK)                                 
      DIMENSION A(N2P,1),WORK(1)                                        
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      LOGICAL   LFIRST                                                  
      DATA      LFIRST /.TRUE./                                         
      SAVE   LFIRST                                                     
CX30015 CALL SBINX (15)                                                 
CY30015 CALL SBINY (15)                                                 
C                                                                       
      IF (LFIRST) CALL CPUTIM (T1)                                      
CX10575 CALL DOINX( 575 )                                               
CY10575 CALL DOINY( 575 )                                               
C                                                                       
      DO 10 J = 1,N1,2                                                  
            DO 20 I = 1,N2P                                             
                  WORK(I) = A(I,J)                                      
   20             WORK(I+N2P) = A(I,J+1)                                
            DO 10 I = 1,N2P                                             
                  II = I+I                                              
                  A(I,J) = WORK(II-1)                                   
   10             A(I,J+1) = WORK(II)                                   
CX20576 CALL DOOUTX( 576 )                                              
CY20576 CALL DOOUTY( 576 )                                              
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TUSHUF = T2-T1                                                  
      END IF                                                            
C                                                                       
      NUSHUF = NUSHUF+1                                                 
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE MEAN (CA,U,G)                                          
      DIMENSION U(1),G(1)                                               
      COMPLEX CA(1)                                                     
      COMMON /NMBRS/ NSKIP(4), N2P, NDUMY1(3)                           
CX30017 CALL SBINX (17)                                                 
CY30017 CALL SBINY (17)                                                 
CX10577 CALL DOINX( 577 )                                               
CY10577 CALL DOINY( 577 )                                               
      DO 10 I = 1,N2P                                                   
10         U(I) = -G(I)*REAL(CA(I))                                     
CX20578 CALL DOOUTX( 578 )                                              
CY20578 CALL DOOUTY( 578 )                                              
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
C                                                                       
C     FGPS  8-APR-88 SET UP RPRNT2 AS A MODIFICATION OF RPRINT.  THE IN 
C                    TENT IS TO PRINT THE VALUE OF A QUANTITY AT THE CEN
C                    OF THE FIRST VORTEX CELL IN ORDER TO MONITOR EXECUT
C                                                                       
      SUBROUTINE RPRNT2 (R,NAME,N1H,N2PD,IMODE)                         
C                                                                       
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      COMMON /VERIFY/ VDAT(15),CDAT(15)                                 
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1     KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,       
     2     KPRESS, NOALS                                                
      REAL            R(N2PD,1)                                         
      CHARACTER*8     NAME                                              
C                                                                       
      LOGICAL         LFIRST                                            
      DATA            LFIRST /.TRUE./                                   
      SAVE   LFIRST                                                     
CX30019 CALL SBINX (19)                                                 
CY30019 CALL SBINY (19)                                                 
C                                                                       
      IF (LFIRST)     CALL CPUTIM (T1)                                  
C                                                                       
      WRITE (6,1002) NSTEPS                                             
C                                                                       
C                                                                       
         IX1 = 1+N1H/4                                                  
         IX2 = 1+N1H/2                                                  
         JY1 = 1+N2PD/4                                                 
         WRITE (6,1010) NAME,JY1,IX1,R(JY1,IX1)                         
         WRITE (6,1020) NAME,JY1,IX2,R(JY1,IX2)                         
         IF (NSTEPS .EQ. 51) THEN                                       
            CDAT(7) = R(JY1,IX2)                                        
         ELSE IF (NSTEPS .EQ. 101) THEN                                 
            CDAT(8) = R(JY1,IX2)                                        
         ELSE IF (NSTEPS .EQ. 151) THEN                                 
            CDAT(9) = R(JY1,IX2)                                        
         ELSE IF (NSTEPS .EQ. 201) THEN                                 
            CDAT(10) = R(JY1,IX2)                                       
         ELSE IF (NSTEPS .EQ. NUMBER) THEN                              
            CDAT(15) = R(JY1,IX2)                                       
         ENDIF                                                          
C                                                                       
      IF (LFIRST) THEN                                                  
         LFIRST = .FALSE.                                               
         CALL CPUTIM (T2)                                               
         TXPRNT = T2-T1                                                 
      END IF                                                            
C                                                                       
      NXPRNT = NXPRNT+1                                                 
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
C                                                                       
      RETURN                                                            
C                                                                       
 1002 FORMAT(/2X,8HNSTEPS =,I7)                                         
 1010 FORMAT(2X,17HTEMPERATURE DROP:/                                   
     1 2X,41HAT THE CENTER OF FIRST CONVECTION CELL = ,A8,              
     2 4H(JY=,I3,4H,IX=,I2,3H) =,1PE13.6)                               
 1020 FORMAT(                                                           
     1 2X,41HAT CENTER OF CONVECTION-CELL BOUNDARY  = ,A8,              
     2 4H(JY=,I3,4H,IX=,I2,3H) =,1PE13.6/)                              
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
      END                                                               
C                                                                       
C***********************************************************************
C                                                                       
C     FGPS 25-MAR-88 CHANGED THE DIMENSIONING OF R(1,1) TO R(N2PD,1)    
C     FGPS 7-MAR-88 SUBROUTINE RPRINT REPROGRAMMED                      
C                                                                       
C     N O T E :  THE INDEXING OF ARRAYS HANDLING VARIABLES REPRESENTED I
C     PHYSICAL SPACE IS USUALLY OF THE FORM                             
C                                                                       
C                        A (1<=JY<=N2PD,1<=IX<=N1H)                     
C                                                                       
C     WHERE "JY" IS THE COUNT ALONG THE Y-AXIS AND "IX" RUNS ALONG THE X
C     AXIS, AND WHERE                                                   
C                                                                       
C                        N2PD = 2 * (N2+1)                 FOR THE Y-AXI
C     AND                                                               
C                         N1H = N1/2                       FOR THE X-AXI
C                                                                       
C     AND WHERE THE BASIC GRID IS DIMENSIONED 2*(N2+1) BY N1.           
C                                                                       
C                                                                       
      SUBROUTINE RPRINT (R,NAME,N1H,N2PD,IMODE)                         
C                                                                       
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      COMMON /VERIFY/ VDAT(15),CDAT(15)                                 
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1     KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,       
     2     KPRESS, NOALS                                                
      REAL            R(N2PD,1)                                         
      CHARACTER*8     NAME                                              
CX30021 CALL SBINX (21)                                                 
CY30021 CALL SBINY (21)                                                 
C                                                                       
C                                                                       
      WRITE (6,1002) NSTEPS                                             
C                                                                       
C     SELECT THE PRINT-OUT MODE                                         
C                                                                       
      IF (IMODE) 1000,2000,3000                                         
C                                                                       
 1000 CONTINUE                                                          
C                                                                       
         IX1 = 1+N1H/4                                                  
         IX2 = 1+N1H/2                                                  
         JY1 = 1+N2PD/4                                                 
         WRITE (6,1010) NAME,JY1,IX1,R(JY1,IX1)                         
C                                                                       
         WRITE (6,1020) NAME,JY1,IX2,R(JY1,IX2)                         
         IF (NSTEPS .EQ. 1) THEN                                        
            IF (NAME .EQ. 'ZETA PHY') THEN                              
               CDAT(2) = R(JY1,IX1)                                     
            ELSE IF (NAME .EQ. 'UVEL PHY') THEN                         
               CDAT(3) = R(JY1,IX2)                                     
            ELSE IF (NAME .EQ. 'VVEL PHY') THEN                         
               CDAT(4) = R(JY1,IX2)                                     
            ENDIF                                                       
         ELSE IF (NSTEPS .EQ. NUMBER) THEN                              
            IF (NAME .EQ. 'ZETA PHY') THEN                              
               CDAT(12) = R(JY1,IX1)                                    
            ELSE IF (NAME .EQ. 'UVEL PHY') THEN                         
               CDAT(13) = R(JY1,IX2)                                    
            ELSE IF (NAME .EQ. 'VVEL PHY') THEN                         
               CDAT(14) = R(JY1,IX2)                                    
            ENDIF                                                       
         ENDIF                                                          
C                                                                       
C                                                                       
      NXPRNT = NXPRNT+1                                                 
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
C                                                                       
      RETURN                                                            
C                                                                       
 1002 FORMAT(/2X,8HNSTEPS =,I7)                                         
 1010 FORMAT(2X,A8,4H(JY=,I3,4H,IX=,I3,3H) =,1PE13.6)                   
 1020 FORMAT(2X,A8,4H(JY=,I3,4H,IX=,I3,3H) =,1PE13.6/)                  
C                                                                       
 2000 CONTINUE                                                          
C                                                                       
C     VALUES OF THE VARIABLE ASSOCIATED WITH NAME ARE PRINTED OUT AT INT
C     VALS AS A FUNCTION OF Y (I.E., 1<=JY<=N2P WITH STRIDE = 4) AT CONS
C     X (I.E., IX = 1+N1H/2)                                            
C                                                                       
         IX1 = 1+N1H/2                                                  
         N2P = N2PD/2                                                   
         WRITE (6,2020) (JY,NAME,IX1,R(JY,IX1),JY=1,N2P,4)              
         WRITE (6,2030)                                                 
         IF (NSTEPS .EQ. 1 .AND. NAME .EQ. 'TEMP PHY') THEN             
            CDAT(5) = R(65,IX1)                                         
         ENDIF                                                          
C                                                                       
      NXPRNT = NXPRNT+1                                                 
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
C                                                                       
      RETURN                                                            
C                                                                       
 2020 FORMAT(34(2X,4HJY =,I5,2X,A8,7H(JY,IX=,I3,3H) =,1PE13.6/))        
 2030 FORMAT(///)                                                       
C                                                                       
 3000 CONTINUE                                                          
C                                                                       
C     VALUES OF THE VARIABLE ASSOCIATED WITH NAME ARE PRINTED OUT AT INT
C     VALS AS A FUNCTION OF X (I.E., 1<=IX<=N1H WITH STRIDE = 2) AT CONS
C     Y (I.E., JY = N2PD/4)                                             
C                                                                       
         JY1 = 3*N2PD/4                                                 
         WRITE (6,3020) (IX,NAME,JY1,R(JY1,IX),IX=1,N1H,2)              
C                                                                       
      NXPRNT = NXPRNT+1                                                 
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
C                                                                       
      RETURN                                                            
C                                                                       
 3020 FORMAT(/,32(2X,4HIX =,I5,2X,A8,4H(JY=,I3,6H,IX) =,1PE13.6,6X,     
     1   4HIX =,I5,2X,A8,4H(JY=,I3,6H,IX) =,1PE13.6/))                  
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
C                                                                       
      END                                                               
C***********************************************************************
C                                                                       
C     FGPS 29-MAR-88 MODIFIED RPRINT =====> CPRINT SO AS TO PRINT OUT   
C                    COMPLEX ARRAYS                                     
C                                                                       
C     N O T E :  THE INDEXING OF COMPLEX ARRAYS REPRESENTING PRINCIPAL V
C     IABLES IS OF THE FORM                                             
C                                                                       
C     COMPLEX  ======>  CA (1<=JY<=N2P,1<=IX<=N1H)                      
C                                                                       
C     WHERE "JY" IS THE COUNT ALONG THE Y-AXIS AND "IX" RUNS ALONG THE X
C     AXIS, AND WHERE                                                   
C                                                                       
C                          N2P = (N2+1)                    FOR THE Y-AXI
C     AND                                                               
C                           N1H = N1/2                     FOR THE X-AXI
C                                                                       
C     AND WHERE THE GRID IS NOMINALLY DIMENSIONED (N2+1) BY N1.  MOREOVE
C     THE INDEX PAIR  JY = (N2+1)/2, IX= N1H/2  IS ASSUMED TO IDENTIFY T
C     CENTER OF THE VORTEX CELL.  HOWEVER, IN ACTUAL USE HERE, THE PROGR
C     MING INVOLVES THE REAL EQUIVALENT TO CA; I.E.,                    
C                                                                       
C     REAL     ======>  RA (1<=JY<=N2P,1<=IX<=N1)                       
C                                                                       
C                                                                       
      SUBROUTINE CPRINT (RA,NAME,N1H,N2P,IMODE)                         
C                                                                       
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /VERIFY/ VDAT(15),CDAT(15)                                 
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1     KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,       
     2     KPRESS, NOALS                                                
      REAL            RA(N2P,1)                                         
      CHARACTER*8     NAME                                              
CX30023 CALL SBINX (23)                                                 
CY30023 CALL SBINY (23)                                                 
C                                                                       
      WRITE (6,1002) NSTEPS                                             
C                                                                       
C     SELECT THE PRINT-OUT MODE                                         
C                                                                       
      IF (IMODE) 1000,2000,3000                                         
C                                                                       
 1000 CONTINUE                                                          
C                                                                       
         IX1 = 1+N1H/4                                                  
         IX2 = 1+N1H/2                                                  
         JY1 = 1+N2P/2                                                  
         WRITE (6,1010) NAME,JY1,IX1,RA(JY1,IX1)                        
         WRITE (6,1020) NAME,JY1,IX2,RA(JY1,IX2)                        
         IF (NAME .EQ. 'PSI PHYS' .AND. NSTEPS .EQ. NUMBER) THEN        
            CDAT(11) = RA(JY1,IX2)                                      
         ENDIF                                                          
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
C                                                                       
      RETURN                                                            
C                                                                       
 1002 FORMAT(/2X,8HNSTEPS =,I7)                                         
 1010 FORMAT(2X,A8,4H(JY=,I3,4H,IX=,I3,3H) =,1PE13.6)                   
 1020 FORMAT(2X,A8,4H(JY=,I3,4H,IX=,I3,3H) =,1PE13.6/)                  
C                                                                       
 2000 CONTINUE                                                          
C                                                                       
C     VALUES OF THE VARIABLE ASSOCIATED WITH NAME ARE PRINTED OUT AT INT
C     VALS AS A FUNCTION OF Y (I.E., 1<=JY<=N2P WITH STRIDE = 4) AT CONS
C     X (E.G., IX = N1H/2)                                              
C                                                                       
         IX1 = N1H/2                                                    
         WRITE (6,2020) (JY,NAME,IX1,RA(JY,IX1),JY=1,N2P,4)             
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
C                                                                       
      RETURN                                                            
C                                                                       
 2020 FORMAT(/,34(2X,4HJY =,I5,2X,A8,7H(JY,IX=,I3,3H) =,1PE13.6,6X,     
     1   4HJY =,I5,2X,A8,7H(JY,IX=,I3,3H) =,1PE13.6/))                  
C                                                                       
 3000 CONTINUE                                                          
C                                                                       
C     VALUES OF THE VARIABLE ASSOCIATED WITH NAME ARE PRINTED OUT AT INT
C     VALS AS A FUNCTION OF X (I.E., 1<=IX<=N1 WITH STRIDE = 4) AT CONST
C     Y (E.G., JY = N2P/2 + 1)                                          
C                                                                       
         N1 = N1H + N1H                                                 
         JY1 = N2P/2 + 1                                                
         WRITE (6,3020) (IX,NAME,JY1,RA(JY1,IX),IX=1,N1,4)              
         IF (NAME .EQ. 'PSI PHYS') THEN                                 
            IF (NSTEPS .EQ. 1) THEN                                     
               CDAT(1) = RA(JY1,33)                                     
            ELSE IF (NSTEPS .EQ. 2) THEN                                
               CDAT(6) = RA(JY1,33)                                     
            ENDIF                                                       
         ENDIF                                                          
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
C                                                                       
      RETURN                                                            
C                                                                       
 3020 FORMAT(/,34(2X,4HIX =,I5,2X,A8,4H(JY=,I3,6H,IX) =,1PE13.6/))      
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
C                                                                       
      END                                                               
C***********************************************************************
C                                                                       
C     FGPS  2-APR-88 REWROTE SUBROUTINE SO AS TO RECEIVE THE COMPLEX ARR
C                    CA IN FOURIER SPACE, TRANSFORM THIS ONTO PHYSICAL S
C                    AND CALL CPRINT TO PRINT OUT THE RESULT            
C                                                                       
      SUBROUTINE PSIPHY (CA,N2P,N1H,EX,G,U,WORK,CWORK3,IMODE)           
      COMMON /EXSIZE/ NMAX2                                             
      COMMON /DATA/ ALPHA, BETA, FX, H, RNU, RKAPA, DT, DTT, T, NSTEPS  
      COMMON /FT/ NPTS, NSKIP, MTRN, MSKIP, ISIGN, LOG, IEXSHF          
      COMPLEX CA(N2P,1), CWORK3(1), EX(1)                               
      REAL    U(1), G(1), WORK(1)                                       
CX30025 CALL SBINX (25)                                                 
CY30025 CALL SBINY (25)                                                 
C                                                                       
      N2PD = N2P+N2P                                                    
      NPTS = N2PD - 2                                                   
      NSKIP = 1                                                         
      MTRN = N1H                                                        
      MSKIP = N2P                                                       
      ISIGN = 1                                                         
      CALL XLOG2 (NPTS,LOG)                                             
      IEXSHF = NMAX2/NPTS                                               
      CALL ACAC (CA,EX,WORK)                                            
C                                                                       
      N1 = N1H+N1H                                                      
      NPTS = N1                                                         
      NSKIP = N2P                                                       
      MTRN = N2P                                                        
      MSKIP = 1                                                         
      ISIGN = 1                                                         
      CALL XLOG2 (NPTS,LOG)                                             
      IEXSHF = NMAX2/NPTS                                               
      CALL CSR (CA,EX)                                                  
C                                                                       
      CALL TELTIM                                                       
      CALL UNSHUF (CA,N2P,N1,WORK)                                      
      CALL CPRINT (CA,'PSI PHYS',N1H,N2P,IMODE)                         
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
C      CALL PLOTME (CA,WORK,NCALL)                                      
      RETURN                                                            
C                                                                       
      END                                                               
C***********************************************************************
C                                                                       
C     FGPS 18-FEB-88 RELOCATED "CALL RPRINT ( , ,N1H,N2PD,IMODE)"       
C                                                                       
      SUBROUTINE ZETAPH (A,N2PD,N1H,EX,WORK,IMODE)                      
      DIMENSION A(N2PD,1), EX(1), WORK(1)                               
      SAVE NCALL                                                        
      DATA NCALL /0/                                                    
CX30027 CALL SBINX (27)                                                 
CY30027 CALL SBINY (27)                                                 
      CALL CSR (A,EX)                                                   
      N2P = N2PD/2                                                      
      N1 = N1H+N1H                                                      
      CALL UNSHUF (A,N2P,N1,WORK)                                       
      CALL TELTIM                                                       
CX10579 CALL DOINX( 579 )                                               
CY10579 CALL DOINY( 579 )                                               
      DO 3 J=1,N1H                                                      
    3    WORK(J) = 0.0                                                  
CX20580 CALL DOOUTX( 580 )                                              
CY20580 CALL DOOUTY( 580 )                                              
CX10581 CALL DOINX( 581 )                                               
CY10581 CALL DOINY( 581 )                                               
      DO 4 I=1,N2PD                                                     
         DO 4 J=1,N1H                                                   
    4       WORK(J) = WORK(J)+A(I,J)*A(I,J)                             
CX20582 CALL DOOUTX( 582 )                                              
CY20582 CALL DOOUTY( 582 )                                              
      SUMSQ=0.0                                                         
CX10583 CALL DOINX( 583 )                                               
CY10583 CALL DOINY( 583 )                                               
      DO 5 J=1,N1H                                                      
    5    SUMSQ = SUMSQ+WORK(J)                                          
CX20584 CALL DOOUTX( 584 )                                              
CY20584 CALL DOOUTY( 584 )                                              
      ZRMS = SQRT (SUMSQ/(REAL(N1)*REAL(N2P)))                          
      WRITE(6,6)  ZRMS                                                  
    6 FORMAT (/' SQRT (AVERAGE ZETASQUARED) =',F30.15)                  
      NCALL = NCALL+1                                                   
      CALL RPRINT (A,'ZETA PHY',N1H,N2PD,IMODE)                         
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
C      CALL PLOTME (A,WORK,NCALL)                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C                                                                       
C                                                                       
C     FGPS  9-APR-88 MADE TO CALL RPRNT2 (NOT RPRINT) WHENEVER IMODE =  
C     FGPS 18-FEB-88 RELOCATED "CALL RPRINT ( , ,N1H,N2PD,IMODE)"       
C     FGPS 10-FEB-88 ADDED "COMPLEX  EX(1)"                             
C                                                                       
      SUBROUTINE TEMPHY (A,N2PD,N1H,EX,WORK,IMODE)                      
      COMMON /LADDR/ LDUMMY, LTN, LDUMY1(10)                            
      DIMENSION A(N2PD,1), WORK(1)                                      
      COMPLEX   EX(1)                                                   
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /MAPING/ GMAPP(129), GDUMY1(129)                           
      SAVE NCALL                                                        
      DATA NCALL /0/                                                    
CX30029 CALL SBINX (29)                                                 
CY30029 CALL SBINY (29)                                                 
      NW = N2PD*N1H                                                     
      N2P = N2PD/2                                                      
      N1 = N1H+N1H                                                      
      CALL IN (A,LTN,NW)                                                
      CALL CSR (A,EX)                                                   
CX10585 CALL DOINX( 585 )                                               
CY10585 CALL DOINY( 585 )                                               
      DO 3 J=1,N1H                                                      
    3    WORK(J) = 0.0                                                  
CX20586 CALL DOOUTX( 586 )                                              
CY20586 CALL DOOUTY( 586 )                                              
CX10587 CALL DOINX( 587 )                                               
CY10587 CALL DOINY( 587 )                                               
      DO 4 I=1,N2PD                                                     
         IPH = (I+1)/2                                                  
         ADJUST = BETA*GMAPP(IPH)                                       
C WHEN INTERESTED IN TEMPERATURE DISTURBANCE ONLY, USE THIS INSTEAD:    
C        ADJUST = 0.0                                                   
         DO 4 J=1,N1H                                                   
            WORK(J) = WORK(J)+A(I,J)*A(I,J)                             
    4       A(I,J) = 0.5*A(I,J)-ADJUST                                  
CX20588 CALL DOOUTX( 588 )                                              
CY20588 CALL DOOUTY( 588 )                                              
      SUMSQ=0.0                                                         
CX10589 CALL DOINX( 589 )                                               
CY10589 CALL DOINY( 589 )                                               
      DO 5 J=1,N1H                                                      
    5    SUMSQ = SUMSQ+WORK(J)                                          
CX20590 CALL DOOUTX( 590 )                                              
CY20590 CALL DOOUTY( 590 )                                              
      ZRMS = 0.5*SQRT (SUMSQ/(REAL(N1)*REAL(N2P)))                      
      WRITE(6,6)  ZRMS                                                  
    6 FORMAT (/' SQRT (AVERAGE TEMPERATURESQUARED) =',F30.15)           
      NCALL = NCALL+1                                                   
      CALL UNSHUF (A,N2P,N1,WORK)                                       
      CALL TELTIM                                                       
      IF (IMODE.LT.0) THEN                                              
         CALL RPRNT2 (A,'TEMP PHY',N1H,N2PD,IMODE)                      
      ELSE                                                              
         CALL RPRINT  (A,'TEMP PHY',N1H,N2PD,IMODE)                     
      END IF                                                            
CX40030 CALL SBOUTX (30)                                                
CY40030 CALL SBOUTY (30)                                                
C      CALL PLOTME (A,WORK,NCALL)                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C                                                                       
C     FGPS 18-FEB-88 RELOCATED "CALL RPRINT ( , ,N1H,N2PD,IMODE)"       
C                                                                       
      SUBROUTINE UVELPH (A,N2PD,N1H,WORK,IMODE)                         
      COMMON /LADDR/ LDUMMY(5), LPY, LDUMY1(6)                          
      DIMENSION A(N2PD,1), WORK(1)                                      
      COMMON /DATA/ DUMY1(2), BOXX, DUMY2(6), IDUMY1                    
CX30031 CALL SBINX (31)                                                 
CY30031 CALL SBINY (31)                                                 
CX10591 CALL DOINX( 591 )                                               
CY10591 CALL DOINY( 591 )                                               
      DO 5 I=1,N2PD                                                     
    5      WORK(I) = 0.0                                                
CX20592 CALL DOOUTX( 592 )                                              
CY20592 CALL DOOUTY( 592 )                                              
      NW = N2PD*N1H                                                     
      N2P = N2PD/2                                                      
      N1 = N1H+N1H                                                      
      CALL TELTIM                                                       
      CALL IN (A,LPY,NW)                                                
CX10593 CALL DOINX( 593 )                                               
CY10593 CALL DOINY( 593 )                                               
      DO 10 J = 1,N1H                                                   
           DO 10 I = 1,N2PD                                             
   10           WORK(I) = MAX(WORK(I), ABS(A(I,J)) )                    
CX20594 CALL DOOUTX( 594 )                                              
CY20594 CALL DOOUTY( 594 )                                              
      WRKMAX = 0.0                                                      
CX10595 CALL DOINX( 595 )                                               
CY10595 CALL DOINY( 595 )                                               
      DO 20 I=1,N2PD                                                    
   20      WRKMAX = MAX(WRKMAX, WORK(I))                                
CX20596 CALL DOOUTX( 596 )                                              
CY20596 CALL DOOUTY( 596 )                                              
      CALL UNSHUF (A,N2P,N1,WORK)                                       
      CALL RPRINT (A,'UVEL PHY',N1H,N2PD,IMODE)                         
C                                                                       
C AT THIS POINT A(N2P,N1) CONTAINS 2.0*U                                
C                                                                       
      UMAX = 0.5*WRKMAX+1.0E-32                                         
      DTMAX = 0.2*(BOXX/REAL(N1))/UMAX                                  
      WRITE(6,100)  DTMAX                                               
  100 FORMAT (//,' DTMAX =',F16.10)                                     
CX40032 CALL SBOUTX (32)                                                
CY40032 CALL SBOUTY (32)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
C                                                                       
C     FGPS 18-FEB-88 RELOCATED "CALL RPRINT ( , ,N1H,N2PD,IMODE)"       
C                                                                       
      SUBROUTINE VVELPH (A,N2PD,N1H,WORK,IMODE)                         
      COMMON /LADDR/ LDUMMY(2), LPN, LDUMY1(9)                          
      DIMENSION A(N2PD,1), WORK(1)                                      
      COMMON /MAPING/ GMAPP(129),G(129)                                 
CX30033 CALL SBINX (33)                                                 
CY30033 CALL SBINY (33)                                                 
      NW = N2PD*N1H                                                     
      N2P = N2PD/2                                                      
      N1 = N1H+N1H                                                      
      CALL TELTIM                                                       
      CALL IN (A,LPN,NW)                                                
      CALL UNSHUF (A,N2P,N1,WORK)                                       
      CALL RPRINT (A,'VVEL PHY',N1H,N2PD,IMODE)                         
C                                                                       
C AT THIS POINT A(N2P,N1) CONTAINS 2.0*V                                
C                                                                       
      CALL TESTDT (A,N2P,N1,WORK,G)                                     
CX40034 CALL SBOUTX (34)                                                
CY40034 CALL SBOUTY (34)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE TESTDT (A,N2P,N1,WORK,G)                               
      COMMON /PI/ PI                                                    
      DIMENSION A(N2P,N1),WORK(N2P),G(N2P)                              
CX30035 CALL SBINX (35)                                                 
CY30035 CALL SBINY (35)                                                 
CX10597 CALL DOINX( 597 )                                               
CY10597 CALL DOINY( 597 )                                               
C                                                                       
C FIND MAXIMUM DT:                                                      
C                                                                       
      DO 10 I=1,N2P                                                     
   10    WORK(I) = 0.0                                                  
CX20598 CALL DOOUTX( 598 )                                              
CY20598 CALL DOOUTY( 598 )                                              
CX10599 CALL DOINX( 599 )                                               
CY10599 CALL DOINY( 599 )                                               
      DO 20 J=1,N1                                                      
         DO 20 I=1,N2P                                                  
   20       WORK(I) = MAX(WORK(I), ABS(G(I)*A(I,J)) )                   
CX20600 CALL DOOUTX( 600 )                                              
CY20600 CALL DOOUTY( 600 )                                              
CX10601 CALL DOINX( 601 )                                               
CY10601 CALL DOINY( 601 )                                               
C                                                                       
C G WAS 1.0/(DY/DX), AND X GOES FROM 0.0 TO PI, HENCE IF WE WANT        
C ((DY/DI)/V)=((DY/DX)(DX/DI)/V)=PI/(N1*G(I)*V(I))                      
C REMEMBER, A CONTAINS 2.0*V:                                           
      DO 30 I=2,N2P                                                     
   30    WORK(1) = MAX(WORK(1), WORK(I))                                
CX20602 CALL DOOUTX( 602 )                                              
CY20602 CALL DOOUTY( 602 )                                              
C                                                                       
C SO, WORK(1) = MAX (G(I)*2*V(I,J)):                                    
C DT(MAX) = 0.2* MIN (DY/V) = 0.2*2.0*PI/(N1*WORK(1)):                  
CSRD  DTMAX = 0.4*PI/(REAL(N1)*WORK(1)+1.0E-1000)                       
      DTMAX = 0.4*PI/(REAL(N1)*WORK(1)+1.0E-300)                        
      WRITE(6,200)  DTMAX                                               
  200 FORMAT (//,' DTMAX =',F16.10)                                     
CX40036 CALL SBOUTX (36)                                                
CY40036 CALL SBOUTY (36)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE PSIFOU (CA,N2P,N1H,CWORK)                              
      COMMON /LADDR/ LDUMMY(2), LPN, LDUMY1(9)                          
      COMPLEX CA(N2P,1),CWORK(1)                                        
CX30037 CALL SBINX (37)                                                 
CY30037 CALL SBINY (37)                                                 
      NW = N2P*N1H*2                                                    
      CALL IN (CA,LPN,NW)                                               
      CALL CPRINT (CA,'PSI FOU',N1H,N2P,-1)                             
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE ZETAFO (CA,N2P,N1H,CWORK)                              
      COMMON /LADDR/ LZN, LDUMY1(11)                                    
      COMPLEX CA(N2P,1),CWORK(1)                                        
CX30039 CALL SBINX (39)                                                 
CY30039 CALL SBINY (39)                                                 
      NW = N2P*N1H*2                                                    
      CALL IN (CA,LZN,NW)                                               
      CALL CPRINT (CA,'ZETA FOU',N1H,N2P,-1)                            
CX40040 CALL SBOUTX (40)                                                
CY40040 CALL SBOUTY (40)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE TEMPFO (CA,N2P,N1H,CWORK)                              
      COMMON /LADDR/ LDUMMY, LTN, LDUMY1(10)                            
      COMPLEX CA(N2P,1), CWORK(1)                                       
CX30041 CALL SBINX (41)                                                 
CY30041 CALL SBINY (41)                                                 
      NW = N2P*N1H*2                                                    
      CALL IN (CA,LTN,NW)                                               
      CALL CPRINT (CA,'TEMP FOU',N1H,N2P,-1)                            
CX40042 CALL SBOUTX (42)                                                
CY40042 CALL SBOUTY (42)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE UVELFO (CA,N2P,N1H,G,EX,CWORK)                         
      COMMON /FT/ NDUMMY(4), ISIGN, NDUMY2(2)                           
      COMPLEX CA(N2P,1),CWORK(1),EX(1)                                  
      DIMENSION G(1)                                                    
CX30043 CALL SBINX (43)                                                 
CY30043 CALL SBINY (43)                                                 
CX10603 CALL DOINX( 603 )                                               
CY10603 CALL DOINY( 603 )                                               
      DO 10 I = 1,N2P                                                   
           DO 10 J = 1,N1H                                              
10              CA(I,J) = G(I)*CA(I,J)                                  
CX20604 CALL DOOUTX( 604 )                                              
CY20604 CALL DOOUTY( 604 )                                              
      ISIGN = -1                                                        
      CALL SCSC (CA,EX,CWORK)                                           
      ISIGN = 1                                                         
      F = 0.125/REAL(N2P-1)                                             
CX10605 CALL DOINX( 605 )                                               
CY10605 CALL DOINY( 605 )                                               
      DO 20 J = 1,N1H                                                   
           DO 20 I = 1,N2P                                              
20              CA(I,J) = F*CA(I,J)                                     
CX20606 CALL DOOUTX( 606 )                                              
CY20606 CALL DOOUTY( 606 )                                              
      CALL CPRINT (CA,'UVEL FOU',N1H,N2P,-1)                            
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE VVELFO (CA,N2P,N1H,FX,CWORK)                           
      COMMON /LADDR/ LDUMMY(2), LPN, LDUMY1(9)                          
      COMPLEX CA(N2P,1),CWORK(1)                                        
CX30045 CALL SBINX (45)                                                 
CY30045 CALL SBINY (45)                                                 
      NW = N2P*N1H*2                                                    
      CALL IN (CA,LPN,NW)                                               
      F = 0.                                                            
CX10607 CALL DOINX( 607 )                                               
CY10607 CALL DOINY( 607 )                                               
      DO 20 J = 1,N1H                                                   
           DO 10 I = 1,N2P                                              
10             CA(I,J) = CMPLX(-F*AIMAG(CA(I,J)), F*REAL(CA(I,J)))      
20         F = F+FX                                                     
CX20608 CALL DOOUTX( 608 )                                              
CY20608 CALL DOOUTY( 608 )                                              
      CALL CPRINT (CA,'VVEL FOU',N1H,N2P,-1)                            
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
C***********************************************************************
      SUBROUTINE BUG (POINT)                                            
      CHARACTER*10 POINT                                                
      COMMON /DEBUG/ NUMBER,TLAST                                       
CX30047 CALL SBINX (47)                                                 
CY30047 CALL SBINY (47)                                                 
      CALL CPUTIM (TNOW)                                                
      TSINCE = TNOW-TLAST                                               
      WRITE(6,100) NUMBER,TNOW,TSINCE,POINT                             
100   FORMAT (' NUMBER =',I3,'; TNOW =',F9.3,' S; TSINCE =',            
     1     F9.3,' S; LOCATION IS: ',A10)                                
      NUMBER = NUMBER+1                                                 
      TLAST = TNOW                                                      
CX40048 CALL SBOUTX (48)                                                
CY40048 CALL SBOUTY (48)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
      SUBROUTINE TELTIM                                                 
      COMMON /DATA/ DUMMY(8),TIME, IDUMY1                               
CX30049 CALL SBINX (49)                                                 
CY30049 CALL SBINY (49)                                                 
      WRITE(6,100) TIME                                                 
100   FORMAT (' TIME = ',F20.6)                                         
CX40050 CALL SBOUTX (50)                                                
CY40050 CALL SBOUTY (50)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
C                                                                       
C     FGPS 12-FEB-1988 SUBR. PLOTME TEMPORARILY COMMENTED "C" OUT       
C                                                                       
C      SUBROUTINE PLOTME (Z,WORK,NCALL)                                 
C      COMMON /DATA/ ALPHA,BETA,BOXX,H, RDUMY1(5), IDUMY1               
C      COMMON /OPTONS/ NX,NY                                            
C      DIMENSION Z(1), WORK(1)                                          
C      COMMON /PLTDEF/ IDPAT, DXA, DXB, DYA, DYB                        
CC                                                                      
CC THIS SUBROUTINE CAN BE USED FOR PLOTTING OF FIELDS IN PHYSICAL       
CC SPACE. CONTOUR PLOTS ARE USUALLY BEST.                               
CC                                                                      
C      ASPECT = BOXX/H                                                  
C      A = (DXB-DXA) / (DYB-DYA)                                        
C      XMIN = DXA                                                       
C      XMAX = DXB                                                       
C      YMIN = DYA                                                       
C      YMAX = DYB                                                       
C      IF (A.LT.ASPECT) GO TO 10                                        
C      SPANH = 0.5*(DYB-DYA) * ASPECT                                   
C      RMEAN = 0.5*(DXA+DXB)                                            
C      XMIN = RMEAN-SPANH                                               
C      XMAX = RMEAN+SPANH                                               
C      GO TO 30                                                         
C   10 SPANH = 0.5*(DXB-DXA) / ASPECT                                   
C      RMEAN = 0.5*(DYA+DYB)                                            
C      YMIN = RMEAN-SPANH                                               
C      YMAX = RMEAN+SPANH                                               
C   30 CALL SET (XMIN,XMAX,YMIN,YMAX,1.0,REAL(NX+1),1.0,REAL(NY+1),1)   
C      CALL EDGES (-1)                                                  
C      CALL DRAWBOX (1.0,1.0,REAL(NX),REAL(NY+1))                       
C      CALL CONREC (Z,NY+1,NY+1,NX,0.0,0.0,0.0,1,0,0)                   
C      CALL FRAME                                                       
C      RETURN                                                           
C      END                                                              
C***********************************************************************
C      SUBROUTINE DATE (TODAY)                                          
C      TODAY = DATEF(1.0)                                               
C      RETURN                                                           
C      END                                                              
C***********************************************************************
C      SUBROUTINE TIME (NOW)                                            
C      NOW = 8HNOW                                                      
C      RETURN                                                           
C      END                                                              
C***********************************************************************
C                           E - O - D                                   
C***********************************************************************
C                                                                       
C FFTFORTRANPACK3                                                       
C FFTFTN.5                                                              
C 10/16/80                                                              
C                                                                       
      SUBROUTINE EXSET( EX, NPTS )                                      
      COMPLEX EX( 1 )                                                   
      COMPLEX A1,A2                                                     
CX30051 CALL SBINX (51)                                                 
CY30051 CALL SBINY (51)                                                 
      W = 6.283185307179586477/REAL( NPTS )                             
      NH = NPTS/2                                                       
CVD$ NOVECTOR                                                           
CVD$ NOCONCUR                                                           
CX10609 CALL DOINX( 609 )                                               
CY10609 CALL DOINY( 609 )                                               
      DO 10 J = 1, NH                                                   
C           EX( J ) = EXP( CMPLX( 0., REAL(J - 1)*W ) )                 
C10         EX( J + NH ) = EXP( CMPLX( 0., REAL(1 - J)*W ) )            
           A1 = CMPLX(0., REAL(J-1)*W)                                  
           EX(J) = EXP(A1)                                              
10         EX(J+NH) = EXP(-1.*A1)                                       
CX20610 CALL DOOUTX( 610 )                                              
CY20610 CALL DOOUTY( 610 )                                              
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
      END                                                               
C***********************************************************************
C                                                                       
C     FGPS 29-APR-88 INTRODUCED "SUBROUTINE XLOG2 (NA,NB)" IN ORDER TO I
C                    PROVE EFFICIENCY.  THIS REPLACES "FUNCTION LOG2 (N)
C                                                                       
C      FUNCTION LOG2 (N)                                                
C      LOG2 = LOG(REAL(N))/LOG(2.)+0.01                                 
C      RETURN                                                           
C      END                                                              
C                                                                       
      SUBROUTINE XLOG2 (NA,NB)                                          
C                                                                       
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
C                                                                       
      LOGICAL    INITL1,INITL2                                          
      DATA       INITL1,INITL2 /.TRUE.,.TRUE./                          
       REAL T1, T2, F1, B                                               
       SAVE INITL1, INITL2, F1                                          
CX30053 CALL SBINX (53)                                                 
CY30053 CALL SBINY (53)                                                 
C                                                                       
      IF        (INITL1) CALL CPUTIM (T1)                               
C                                                                       
      IF        (INITL2) THEN                                           
         INITL2 = .FALSE.                                               
         F1     = 1.0/LOG(2.0)                                          
      END IF                                                            
C                                                                       
      B  = F1*LOG( REAL(NA) ) + 0.01                                    
      NB = INT(B)                                                       
C                                                                       
      IF        (INITL1) THEN                                           
         INITL1 = .FALSE.                                               
         CALL CPUTIM (T2)                                               
         TXLOG2 = T2-T1                                                 
      END IF                                                            
C                                                                       
      NXLOG2 = NXLOG2+1                                                 
CX40054 CALL SBOUTX (54)                                                
CY40054 CALL SBOUTY (54)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
      FUNCTION LCMREQ (LARGE)                                           
      PARAMETER (LCMSIZ = 165117)                                       
      COMMON /MEMORY/ MEMSIZ,RM(LCMSIZ)                                 
CX30055 CALL SBINX (55)                                                 
CY30055 CALL SBINY (55)                                                 
      IF (LARGE.GT.LCMSIZ) GO TO 10                                     
      MEMSIZ = LARGE                                                    
      LCMREQ = 1                                                        
      WRITE(6,100) LARGE                                                
100   FORMAT (//' LCM REQUESTED:',I7,' WORDS',//)                       
CX40056 CALL SBOUTX (56)                                                
CY40056 CALL SBOUTY (56)                                                
      RETURN                                                            
C                                                                       
   10 WRITE(6,11) MEMSIZ,LARGE                                          
   11 FORMAT (///' LCMREQ: REQUEST EXCEEDS MEMORY SIZE.'/               
     1          ,' MEMORY SIZE:',I7,' ;   WORDS REQUESTED:',I22)        
CX40056 CALL SBOUTX (56)                                                
CY40056 CALL SBOUTY (56)                                                
      STOP                                                              
      END                                                               
C***********************************************************************
C                                                                       
      SUBROUTINE IN (ARRAY,LOC,NW)                                      
      COMMON /MEMORY/ MEMSIZ,RM(165117)                                 
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      DIMENSION ARRAY(1)                                                
      LOGICAL   LFIRST                                                  
      DATA      LFIRST /.TRUE./                                         
      SAVE   LFIRST                                                     
CX30057 CALL SBINX (57)                                                 
CY30057 CALL SBINY (57)                                                 
C                                                                       
      IF (LFIRST) CALL CPUTIM (T1)                                      
C                                                                       
      IF ((LOC.LT.1).OR.((LOC+NW-1).GT.MEMSIZ).OR.(NW.LT.1)) GO TO 100  
      LOCM=LOC-1                                                        
CX10611 CALL DOINX( 611 )                                               
CY10611 CALL DOINY( 611 )                                               
      DO 10 I=1,NW                                                      
           ARRAY(I) = RM(LOCM+I)                                        
   10 CONTINUE                                                          
CX20612 CALL DOOUTX( 612 )                                              
CY20612 CALL DOOUTY( 612 )                                              
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TXXXIN = T2-T1                                                  
      END IF                                                            
C                                                                       
      NXXXIN = NXXXIN+1                                                 
CX40058 CALL SBOUTX (58)                                                
CY40058 CALL SBOUTY (58)                                                
C                                                                       
      RETURN                                                            
C                                                                       
  100 WRITE(6,101) LOC,NW                                               
  101 FORMAT (///' IN: ILLEGAL REQUEST. LOCATION:',I22,                 
     1           ' ;   NUMBER OF WORDS:',I22,' ;  MEMSIZ:',I7)          
CX40058 CALL SBOUTX (58)                                                
CY40058 CALL SBOUTY (58)                                                
      STOP                                                              
      END                                                               
C***********************************************************************
C                                                                       
      SUBROUTINE OUT (ARRAY,LOC,NW)                                     
      COMMON /MEMORY/ MEMSIZ,RM(165117)                                 
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      DIMENSION ARRAY(1)                                                
      LOGICAL   LFIRST                                                  
      DATA      LFIRST /.TRUE./                                         
      SAVE   LFIRST                                                     
CX30059 CALL SBINX (59)                                                 
CY30059 CALL SBINY (59)                                                 
C                                                                       
      IF (LFIRST) CALL CPUTIM (T1)                                      
C                                                                       
      IF ((LOC.LT.1).OR.((LOC+NW-1).GT.MEMSIZ).OR.(NW.LT.1)) GO TO 100  
      LOCM = LOC-1                                                      
CX10613 CALL DOINX( 613 )                                               
CY10613 CALL DOINY( 613 )                                               
      DO 10 I=1,NW                                                      
           RM(LOCM+I) = ARRAY(I)                                        
   10 CONTINUE                                                          
CX20614 CALL DOOUTX( 614 )                                              
CY20614 CALL DOOUTY( 614 )                                              
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TXXOUT = T2-T1                                                  
      END IF                                                            
C                                                                       
      NXXOUT = NXXOUT+1                                                 
CX40060 CALL SBOUTX (60)                                                
CY40060 CALL SBOUTY (60)                                                
C                                                                       
      RETURN                                                            
C                                                                       
  100 WRITE(6,101) LOC,NW,MEMSIZ                                        
  101 FORMAT (///' OUT: ILLEGAL REQUEST. LOCATION:',I22,                
     1           ' ;   NUMBER OF WORDS:',I22,' ;  MEMSIZ:',I7)          
CX40060 CALL SBOUTX (60)                                                
CY40060 CALL SBOUTY (60)                                                
      STOP                                                              
      END                                                               
C                                                                       
      SUBROUTINE RCS( A, EX )                                           
      COMMON /FT/ NPTS, NSKIP, MTRN, MSKIP, ISIGN, LOG, IEXSHF          
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1     KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,       
     2     KPRESS, NOALS                                                
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      COMPLEX A( 1 ), EX( 1 ), BJ, BJSTAR                               
C                                                                       
C   PERFORM A REAL TO CONJUGATE SYMMETRIC TRANFORM, SAVING ONLY HALF OF 
C   THE RESULTING CONJUGATE SYMMETRIC ARRAY, AND PACKING A(0) AND A(N/2)
C   INTO A(0) SINCE THEY ARE BOTH REAL                                  
C                                                                       
      LOGICAL   LFIRST,LXSTEP,LXTIME                                    
      DATA      LFIRST,LXSTEP,LXTIME /.TRUE.,.FALSE.,.TRUE./            
      SAVE   LFIRST,LXSTEP,LXTIME                                       
CX30061 CALL SBINX (61)                                                 
CY30061 CALL SBINY (61)                                                 
C                                                                       
      IF (NSTEPS.EQ.NUMBER) LXSTEP = .TRUE.                             
      IF (LFIRST) CALL CPUTIM (T1)                                      
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) CALL CPUTIM (T1)                                    
      END IF                                                            
C                                                                       
      N = NPTS                                                          
      NPTS = NPTS/2                                                     
      M = NPTS                                                          
      ISHIFT = 0                                                        
      IF (ISIGN .LT. 0)  ISHIFT=M*IEXSHF                                
      MH = M/2                                                          
      MHP = MH + 1                                                      
      LOG = LOG - 1                                                     
      IEXSHF = 2*IEXSHF                                                 
      CALL FTRVMT( A( 1 ), EX( 1 ) )                                    
      LOG = LOG + 1                                                     
      IEXSHF = IEXSHF/2                                                 
      NPTS = N                                                          
CX10615 CALL DOINX( 615 )                                               
CY10615 CALL DOINY( 615 )                                               
      DO 20 LOOP = 1, MTRN                                              
      IST = (LOOP - 1)*MSKIP                                            
      ISTP = IST + 1                                                    
      A( ISTP ) = CMPLX( 2.*( REAL( A(ISTP) ) + AIMAG( A( ISTP ) ) ),   
     1                 2.*( REAL( A(ISTP) ) - AIMAG( A( ISTP ) ) ) )    
      DO 10 I = 2, MH                                                   
      J = (I - 1)*NSKIP + ISTP                                          
      JSTAR = (M - I + 1)*NSKIP + ISTP                                  
      BJ = A( J )                                                       
      BJSTAR = A( JSTAR )                                               
      IEXL = IEXSHF*(I - 1) + 1 + ISHIFT                                
      IEXH = IEXSHF*(M - I + 1) + 1 + ISHIFT                            
      A( J ) = BJ + CONJG( BJSTAR ) +                                   
     1                    (0., -1.)*EX( IEXL )*(BJ - CONJG( BJSTAR ) )  
      A( JSTAR ) = BJSTAR + CONJG( BJ ) +                               
     1             (0., -1.)*EX( IEXH )*(BJSTAR - CONJG( BJ ) )         
10    CONTINUE                                                          
      J = MH*NSKIP + ISTP                                               
      A( J ) = 2.*A( J )                                                
      IF (ISIGN .LT. 0)  A(J) = CONJG( A(J) )                           
20    CONTINUE                                                          
CX20616 CALL DOOUTX( 616 )                                              
CY20616 CALL DOOUTY( 616 )                                              
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TXXRCS = T2-T1                                                  
      END IF                                                            
C                                                                       
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) THEN                                                
          LXTIME = .FALSE.                                              
          CALL CPUTIM (T2)                                              
          EXXRCS = T2-T1                                                
        END IF                                                          
      END IF                                                            
C                                                                       
      NXXRCS = NXXRCS+1                                                 
CX40062 CALL SBOUTX (62)                                                
CY40062 CALL SBOUTY (62)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE CSR( A, EX )                                           
      COMMON /FT/ NPTS, NSKIP, MTRN, MSKIP, ISIGN, LOG, IEXSHF          
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1     KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,       
     2     KPRESS, NOALS                                                
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      COMPLEX A( 1 ), EX( 1 ), BJ, BJSTAR                               
C                                                                       
C   PERFORM A CONJUGATE SYMMETRIC TO REAL TRANSFORM ON A PACKED         
C   CONJUGATE SYMMETRIC ARRAY.  ONLY HALF THE ARRAY IS STORED, AND A(0) 
C   AND A(N/2) ARE BOTH PACKED TOGETHER IN A(0) SINCE THEY ARE BOTH REAL
C                                                                       
      LOGICAL   LFIRST,LXSTEP,LXTIME                                    
      DATA      LFIRST,LXSTEP,LXTIME /.TRUE.,.FALSE.,.TRUE./            
      SAVE   LFIRST,LXSTEP,LXTIME                                       
CX30063 CALL SBINX (63)                                                 
CY30063 CALL SBINY (63)                                                 
C                                                                       
      IF (NSTEPS.EQ.NUMBER) LXSTEP = .TRUE.                             
      IF (LFIRST) CALL CPUTIM (T1)                                      
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) CALL CPUTIM (T1)                                    
      END IF                                                            
C                                                                       
      N = NPTS                                                          
      NPTS = NPTS/2                                                     
      M = NPTS                                                          
      ISHIFT = 0                                                        
      IF (ISIGN .LT. 0)  ISHIFT=M*IEXSHF                                
      MH = M/2                                                          
      MHP = MH + 1                                                      
CX10617 CALL DOINX( 617 )                                               
CY10617 CALL DOINY( 617 )                                               
      DO 20 LOOP = 1, MTRN                                              
      IST = (LOOP - 1)*MSKIP                                            
      ISTP = IST + 1                                                    
      A( ISTP ) = CMPLX( REAL( A(ISTP) ) + AIMAG( A(ISTP) ),            
     1                    REAL( A(ISTP) ) - AIMAG( A(ISTP) ) )          
      DO 10 I = 2, MH                                                   
      J = (I - 1)*NSKIP + ISTP                                          
      JSTAR = (M - I + 1)*NSKIP + ISTP                                  
      BJ = A( J )                                                       
      BJSTAR = A( JSTAR )                                               
      IEXL = IEXSHF*(I - 1) + 1 + ISHIFT                                
      IEXH = IEXSHF*(M - I + 1) + 1 + ISHIFT                            
      A( J ) = BJ + CONJG( BJSTAR ) +                                   
     1                (0., 1.)*EX( IEXL )*( BJ - CONJG( BJSTAR ) )      
      A( JSTAR ) = BJSTAR + CONJG( BJ ) +                               
     1                (0., 1.)*EX( IEXH )*( BJSTAR - CONJG( BJ ) )      
10    CONTINUE                                                          
      J = MH*NSKIP + ISTP                                               
      A(J) = 2. * A(J)                                                  
      IF (ISIGN .GT. 0)  A(J) = CONJG( A(J) )                           
20    CONTINUE                                                          
CX20618 CALL DOOUTX( 618 )                                              
CY20618 CALL DOOUTY( 618 )                                              
      LOG = LOG - 1                                                     
      IEXSHF = 2*IEXSHF                                                 
      CALL FTRVMT( A( 1 ), EX( 1 ) )                                    
      LOG = LOG + 1                                                     
      IEXSHF = IEXSHF/2                                                 
      NPTS=N                                                            
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TXXCSR = T2-T1                                                  
      END IF                                                            
C                                                                       
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) THEN                                                
          LXTIME = .FALSE.                                              
          CALL CPUTIM (T2)                                              
          EXXCSR = T2-T1                                                
        END IF                                                          
      END IF                                                            
C                                                                       
      NXXCSR = NXXCSR+1                                                 
CX40064 CALL SBOUTX (64)                                                
CY40064 CALL SBOUTY (64)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE SCSC(AC,EX,WORK)                                       
C                                                                       
C COMPLEX SYMMETRIC TO COMPLEX SYMMETRIC                                
C                                                                       
      COMMON/FT/NPTS,NSKIP,MTRN,MSKIP,ISIGN,LOG,IEXSHF                  
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1     KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,       
     2     KPRESS, NOALS                                                
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      COMPLEX AC(1),WORK(1),TEMP1,TEMP2,X,Y,BJ,BJSYM                    
      DIMENSION EX(1)                                                   
C                                                                       
      LOGICAL   LFIRST,LXSTEP,LXTIME                                    
      DATA      LFIRST,LXSTEP,LXTIME /.TRUE.,.FALSE.,.TRUE./            
      SAVE   LFIRST,LXSTEP,LXTIME                                       
CX30065 CALL SBINX (65)                                                 
CY30065 CALL SBINY (65)                                                 
C                                                                       
      IF (NSTEPS.EQ.NUMBER) LXSTEP = .TRUE.                             
      IF (LFIRST) CALL CPUTIM (T1)                                      
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) CALL CPUTIM (T1)                                    
      END IF                                                            
C                                                                       
      M=NPTS/2                                                          
      NPTS=M                                                            
      MH=M/2                                                            
      MHP=MH+1                                                          
      ISTOP=1+(MTRN-1)*MSKIP                                            
      MNSKIP=M*NSKIP                                                    
      LOG=LOG-1                                                         
      NSKIP2=NSKIP*2                                                    
      MP2=M+2                                                           
      IEXSHF=IEXSHF*2                                                   
      FACT=REAL(ISIGN)/2.                                               
      MNSH=MNSKIP/2                                                     
      IPROD=MNSH-NSKIP                                                  
      IEXSYM=4+M*IEXSHF                                                 
CX10619 CALL DOINX( 619 )                                               
CY10619 CALL DOINY( 619 )                                               
      DO 30 LOOP=1,ISTOP,MSKIP                                          
         ISUM=LOOP+MNSKIP                                               
         WORK(1)=AC(LOOP)                                               
         WORK(MHP)=AC(ISUM)                                             
         IAC=LOOP+NSKIP                                                 
         Y=AC(IAC)                                                      
         DO 10 I=2,MH                                                   
            IAC=IAC+NSKIP2                                              
            Y=Y+AC(IAC)                                                 
            TEMP1=CMPLX(0.,1.)*(AC(IAC)-AC(IAC-NSKIP2))                 
            TEMP2=AC(IAC-NSKIP)                                         
            WORK(I)=TEMP2+TEMP1                                         
            WORK(MP2-I)=TEMP2-TEMP1                                     
   10    CONTINUE                                                       
         IAC=LOOP                                                       
         DO 20 I=1,M                                                    
            AC(IAC)=WORK(I)                                             
            IAC=IAC+NSKIP                                               
   20    CONTINUE                                                       
         AC(ISUM)=4.*Y                                                  
   30 CONTINUE                                                          
CX20620 CALL DOOUTX( 620 )                                              
CY20620 CALL DOOUTY( 620 )                                              
      CALL FTRVMT(AC,EX)                                                
CX10621 CALL DOINX( 621 )                                               
CY10621 CALL DOINY( 621 )                                               
      DO 40 LOOP=1,ISTOP,MSKIP                                          
         ISUM=LOOP+MNSKIP                                               
         JSUM=ISUM+LOOP                                                 
         X=2.*AC(LOOP)                                                  
         Y=AC(ISUM)                                                     
         AC(LOOP)=X+Y                                                   
         AC(ISUM)=X-Y                                                   
         AC(LOOP+MNSH)=2.*AC(LOOP+MNSH)                                 
         ISTART=LOOP+NSKIP                                              
         IFIN=LOOP+IPROD                                                
         I=2                                                            
         DO 35 J=ISTART,IFIN,NSKIP                                      
            I=I+IEXSHF                                                  
            JSYM=JSUM-J                                                 
            BJ=AC(J)                                                    
            BJSYM=AC(JSYM)                                              
            TEMP1=BJ+BJSYM                                              
            TEMP2=CMPLX(FACT/EX(I),0.)*(BJ-BJSYM)                       
            AC(J)=TEMP1+TEMP2                                           
            AC(JSYM)=TEMP1-TEMP2                                        
   35    CONTINUE                                                       
   40 CONTINUE                                                          
CX20622 CALL DOOUTX( 622 )                                              
CY20622 CALL DOOUTY( 622 )                                              
      NPTS=M*2                                                          
      LOG=LOG+1                                                         
      IEXSHF=IEXSHF/2                                                   
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TXSCSC = T2-T1                                                  
      END IF                                                            
C                                                                       
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) THEN                                                
          LXTIME = .FALSE.                                              
          CALL CPUTIM (T2)                                              
          EXSCSC = T2-T1                                                
        END IF                                                          
      END IF                                                            
C                                                                       
      NXSCSC = NXSCSC+1                                                 
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE ACAC(AC,EX,WORK)                                       
C                                                                       
C COMPLEX ANTISYMMETRIC TO COMPLEX ANTISYMMETRIC                        
C                                                                       
      COMMON/FT/NPTS,NSKIP,MTRN,MSKIP,ISIGN,LOG,IEXSHF                  
      COMMON /DATA/ ALPHA, BETA, BOXX, H, RNU, RKAPA, DT, DTT, T, NSTEPS
      COMMON /OPTONS/ N1, N2, NUMBER, KDATA, LMAPP, NOTEMP, KTA, KTAPE, 
     1     KSM, KPP, KZP, KTP, KUP, KVP, KPF, KZF, KTF, KUF, KVF,       
     2     KPRESS, NOALS                                                
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
      COMPLEX AC(1),WORK(1),TEMP1,TEMP2,X,Y,BJ,BJSYM                    
      DIMENSION EX(1)                                                   
C                                                                       
      LOGICAL   LFIRST,LXSTEP,LXTIME                                    
      DATA      LFIRST,LXSTEP,LXTIME /.TRUE.,.FALSE.,.TRUE./            
CX30067 CALL SBINX (67)                                                 
CY30067 CALL SBINY (67)                                                 
C                                                                       
      IF (NSTEPS.EQ.NUMBER) LXSTEP = .TRUE.                             
      IF (LFIRST) CALL CPUTIM (T1)                                      
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) CALL CPUTIM (T1)                                    
      END IF                                                            
C                                                                       
      M=NPTS/2                                                          
      NPTS=M                                                            
      MH=M/2                                                            
      MHP=MH+1                                                          
      ISTOP=1+(MTRN-1)*MSKIP                                            
      MNSKIP=M*NSKIP                                                    
      LOG=LOG-1                                                         
      NSKIP2=NSKIP*2                                                    
      MP2=M+2                                                           
      IEXSHF=IEXSHF*2                                                   
      FACT=REAL(ISIGN)/2.                                               
      MNSH=MNSKIP/2                                                     
      IPROD=MNSH                                                        
CX10623 CALL DOINX( 623 )                                               
CY10623 CALL DOINY( 623 )                                               
      DO 30 LOOP=1,ISTOP,MSKIP                                          
         ISUM=LOOP+MNSKIP                                               
         IAC=LOOP+NSKIP                                                 
         WORK(1)=CMPLX(0.,2.)*AC(IAC)                                   
         DO 10 I=2,MH                                                   
            IAC=IAC+NSKIP2                                              
            TEMP1=CMPLX(0.,1.)*(AC(IAC)-AC(IAC-NSKIP2))                 
            TEMP2=AC(IAC-NSKIP)                                         
            WORK(I)=TEMP2+TEMP1                                         
            WORK(MP2-I)=TEMP1-TEMP2                                     
   10    CONTINUE                                                       
         WORK(MHP)=CMPLX(0.,-2.)*AC(IAC)                                
         IAC=LOOP                                                       
         DO 20 I=1,M                                                    
            AC(IAC)=WORK(I)                                             
            IAC=IAC+NSKIP                                               
   20    CONTINUE                                                       
   30 CONTINUE                                                          
CX20624 CALL DOOUTX( 624 )                                              
CY20624 CALL DOOUTY( 624 )                                              
      CALL FTRVMT(AC,EX)                                                
CX10625 CALL DOINX( 625 )                                               
CY10625 CALL DOINY( 625 )                                               
      DO 40 LOOP=1,ISTOP,MSKIP                                          
         ISUM=LOOP+MNSKIP                                               
         JSUM=ISUM+LOOP                                                 
         AC(LOOP)=0.                                                    
         ISTART=LOOP+NSKIP                                              
         IFIN=LOOP+IPROD                                                
         I=2                                                            
         DO 35 J=ISTART,IFIN,NSKIP                                      
            I=I+IEXSHF                                                  
            JSYM=JSUM-J                                                 
            BJ=AC(J)                                                    
            BJSYM=AC(JSYM)                                              
            TEMP1=BJ-BJSYM                                              
            TEMP2=CMPLX(FACT/EX(I),0.)*(BJ+BJSYM)                       
            AC(J)=TEMP1+TEMP2                                           
            AC(JSYM)=TEMP2-TEMP1                                        
   35    CONTINUE                                                       
   40 CONTINUE                                                          
CX20626 CALL DOOUTX( 626 )                                              
CY20626 CALL DOOUTY( 626 )                                              
      NPTS=M*2                                                          
      LOG=LOG+1                                                         
      IEXSHF=IEXSHF/2                                                   
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TXACAC = T2-T1                                                  
      END IF                                                            
C                                                                       
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) THEN                                                
          LXTIME = .FALSE.                                              
          CALL CPUTIM (T2)                                              
          EXACAC = T2-T1                                                
        END IF                                                          
      END IF                                                            
C                                                                       
      NXACAC = NXACAC+1                                                 
CX40068 CALL SBOUTX (68)                                                
CY40068 CALL SBOUTY (68)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE FTRVMT (DATA,EX)                                       
C                                                                       
C THIS IS A FORTRAN VERSION OF THE ASSEMBLY LANGUAGE SUBROUTINE         
C FFT2, TO BE USED WITH THE NSSL SUBROUTINE PACKAGE FFTPOW2.            
C                                                                       
      COMMON/FT/NPTS,NSKIP,MTRN,MSKIP,ISIGN,LOG,IEXSHF                  
C                                                                       
C DATA IS A COMPLEX ARRAY OF LENGTH NPTS WHICH SERVES FOR BOTH          
C INPUT AND OUTPUT.                                                     
C                                                                       
      COMMON /DATA/ ALPHA, BETA, BOXX, H1, RNU, RKAPA, DT, DTT, T, NSTEP
      COMMON /TIMING/ TXXXIN,TXXOUT,TXXRCS,TXXCSR,TXSCSC,TXACAC,TFTVMT, 
     1                TXSHUF,TUSHUF,TXPRNT,TXLOG2,                      
     2                EXXRCS,EXXCSR,EXSCSC,EXACAC,EFTVMT,               
     3                NXXXIN,NXXOUT,NXXRCS,NXXCSR,NXSCSC,NXACAC,NFTVMT, 
     4                NXSHUF,NUSHUF,NXPRNT,NXLOG2                       
C                                                                       
      COMPLEX DATA(99), EXP1, EXK, EXJ, H, EX(1)                        
      REAL HH(2)                                                        
      EQUIVALENCE (HH,H)                                                
      LOGICAL   LFIRST,LXSTEP,LXTIME                                    
      DATA      LFIRST,LXSTEP,LXTIME /.TRUE.,.FALSE.,.TRUE./            
      SAVE   LFIRST,LXSTEP,LXTIME                                       
CX30069 CALL SBINX (69)                                                 
CY30069 CALL SBINY (69)                                                 
C                                                                       
      IF (NSTEPS.EQ.NUMBER) LXSTEP = .TRUE.                             
      IF (LFIRST) CALL CPUTIM (T1)                                      
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) CALL CPUTIM (T1)                                    
      END IF                                                            
C                                                                       
      I2K=NPTS                                                          
      IF(I2K.EQ.1) GO TO 7777                                           
      SGN1=ISIGN                                                        
      EXK=EX(1+IEXSHF)                                                  
      IF (SGN1.LT.0.) EXK=CONJG(EXK)                                    
  101 CONTINUE                                                          
      I2KP=I2K                                                          
      I2K=I2K/2                                                         
      I2KS=I2K*NSKIP                                                    
      IF(2*I2K.NE.I2KP) GO TO 117                                       
      JLI=I2K/2+1                                                       
CX10627 CALL DOINX( 627 )                                               
CY10627 CALL DOINY( 627 )                                               
      DO 109 JL=1,I2K                                                   
         IF(JL-1) 102,102,104                                           
  102    EXJ=(1.,0.)                                                    
         DO 103 JJ=JL,NPTS,I2KP                                         
            DO 103 MM=1,MTRN                                            
               JS=(JJ-1)*NSKIP+(MM-1)*MSKIP+1                           
               H=DATA(JS)-DATA(JS+I2KS)                                 
               DATA(JS)=DATA(JS)+DATA(JS+I2KS)                          
               DATA(JS+I2KS)=H                                          
  103    CONTINUE                                                       
         GO TO 109                                                      
  104    IF(JL-JLI) 105,107,105                                         
C                                                                       
C INCREMENT JL-DEPENDENT EXPONENTIAL FACTOR                             
C                                                                       
  105    EXJ=EXJ*EXK                                                    
         DO 106 JJ=JL,NPTS,I2KP                                         
            DO 106 MM=1,MTRN                                            
               JS=(JJ-1)*NSKIP+(MM-1)*MSKIP+1                           
               H=DATA(JS)-DATA(JS+I2KS)                                 
               DATA(JS)=DATA(JS)+DATA(JS+I2KS)                          
               DATA(JS+I2KS)=H*EXJ                                      
  106    CONTINUE                                                       
         GO TO 109                                                      
  107    EXJ=CMPLX(0.,SGN1)                                             
         DO 108 JJ=JL,NPTS,I2KP                                         
            DO 108 MM=1,MTRN                                            
               JS=(JJ-1)*NSKIP+(MM-1)*MSKIP+1                           
               H=DATA(JS)-DATA(JS+I2KS)                                 
               DATA(JS)=DATA(JS)+DATA(JS+I2KS)                          
               DATA(JS+I2KS)=CMPLX(-SGN1*HH(2),SGN1*HH(1))              
  108    CONTINUE                                                       
  109 CONTINUE                                                          
CX20628 CALL DOOUTX( 628 )                                              
CY20628 CALL DOOUTY( 628 )                                              
C                                                                       
C INCREMENT K-DEPENDENT EXPONENTIAL FACTOR                              
C                                                                       
      EXK=EXK*EXK                                                       
      IF(I2K-1) 110,110,101                                             
  110 IF(NPTS.LE.2) GO TO 7777                                          
      NPTSD2=NPTS/2                                                     
      JMIN=0                                                            
      JMAX=NPTS-4                                                       
      JREV=0                                                            
CX10629 CALL DOINX( 629 )                                               
CY10629 CALL DOINY( 629 )                                               
      DO 116 J=JMIN,JMAX,2                                              
         I2K=NPTSD2                                                     
         JREV2=JREV+I2K                                                 
         IF(JREV2-(J+1)) 113,113,111                                    
  111 CONTINUE                                                          
      DO      1110      MM=1,MTRN                                       
            JREV2S=(MM-1)*MSKIP+JREV2*NSKIP                             
            JS=(MM-1)*MSKIP+(J+1)*NSKIP                                 
            H=DATA(JREV2S+1)                                            
            DATA(JREV2S+1)=DATA(JS+1)                                   
            DATA(JS+1)=H                                                
 1110    CONTINUE                                                       
         IF(JREV-J) 113,113,112                                         
  112 CONTINUE                                                          
      DO      1120      MM=1,MTRN                                       
            JREVS=JREV*NSKIP+(MM-1)*MSKIP                               
            JS=J*NSKIP+(MM-1)*MSKIP                                     
            H=DATA(JREVS+1)                                             
            DATA(JREVS+1)=DATA(JS+1)                                    
            DATA(JS+1)=H                                                
 1120    CONTINUE                                                       
  113    CONTINUE                                                       
  114    I2K=I2K/2                                                      
         JREV=JREV-I2K                                                  
         IF(JREV) 115,114,114                                           
  115    JREV=JREV+2*I2K                                                
  116 CONTINUE                                                          
CX20630 CALL DOOUTX( 630 )                                              
CY20630 CALL DOOUTY( 630 )                                              
      GO TO 7777                                                        
  117 IERR=101                                                          
C                                                                       
 7777 CONTINUE                                                          
C                                                                       
      IF (LFIRST) THEN                                                  
        LFIRST = .FALSE.                                                
        CALL CPUTIM (T2)                                                
        TFTVMT = T2-T1                                                  
      END IF                                                            
C                                                                       
      IF (LXSTEP) THEN                                                  
        IF (LXTIME) THEN                                                
          LXTIME = .FALSE.                                              
          CALL CPUTIM (T2)                                              
          EFTVMT = T2-T1                                                
        END IF                                                          
      END IF                                                            
C                                                                       
      NFTVMT = NFTVMT+1                                                 
CX40070 CALL SBOUTX (70)                                                
CY40070 CALL SBOUTY (70)                                                
C                                                                       
      RETURN                                                            
      END                                                               
C***********************************************************************
C                                                                       
C     "CPUTIM" RETURNS ELAPSED USER (CPU) TIME, T, IN SECONDS, WHERE "T"
C     A REAL VARIABLE                                                   
C                                                                       
C     ON THE CYBER205/ETA-10, THE INTRINSIC FUNCTION "SECOND()" IS REFER
C                                                                       
C      SUBROUTINE CPUTIM (T)                                            
C                                                                       
C         REAL T                                                        
C         T = SECOND()                                                  
C                                                                       
C      RETURN                                                           
C      END                                                              
C***********************************************************************
C                                                                       
C     "ELAPSE" RETURNS ELAPSED WALL-CLOCK TIME, T, IN SECONDS           
C                                                                       
C     ON THE CYBER205 SYSTEM USE IS MADE OF THE SIL SUBROUTINE "Q5DCDMSC
C                                                                       
C     FGPS 19-APR-88                                                    
C                                                                       
C      SUBROUTINE ELAPSE (T)                                            
C                                                                       
C         REAL T                                                        
C         CALL Q5DCDMSC ('MILLSEC=',TELAPSED)                           
C         T = TELAPSED*1.0E-03                                          
C                                                                       
C      RETURN                                                           
C      END                                                              
