      SUBROUTINE AMHMTM(DEL,RPI,SV,P1,P2,AM,HM,TM)                      
      DIMENSION DEL( 12 ),RPI( 11 ),SV( 12 ),                           
     1 P1( 12 ),P2( 12 ),LLL( 12 ),MMM( 12 ),                           
     2 AM( 12 , 12 ),HM( 12 , 12 ),TM( 12 , 12 )                        
CX30001 CALL SBINX (1)                                                  
CY30001 CALL SBINY (1)                                                  
CX10501 CALL DOINX( 501 )                                               
CY10501 CALL DOINY( 501 )                                               
      DO 5 J=1, 12                                                      
      DO 5 I=1, 12                                                      
      HM(I,J) = 0.                                                      
5     TM(I,J) = 0.                                                      
CX20502 CALL DOOUTX( 502 )                                              
CY20502 CALL DOOUTY( 502 )                                              
CX10503 CALL DOINX( 503 )                                               
CY10503 CALL DOINY( 503 )                                               
      DO 6 LE=1, 11                                                     
      HM(LE,LE) = 1.                                                    
      TM(LE,LE) = 0.5 *  1005.  * (RPI(LE)-1.)                          
C                                                                       
      HM(LE,LE+1) = -1.                                                 
6     TM(LE,LE+1) = 0.5 *  1005.  * (1.-1./RPI(LE))                     
CX20504 CALL DOOUTX( 504 )                                              
CY20504 CALL DOOUTY( 504 )                                              
CX10505 CALL DOINX( 505 )                                               
CY10505 CALL DOINY( 505 )                                               
      DO 7 LE=1, 12                                                     
      HM( 12 ,LE) = DEL(LE)                                             
7     TM( 12 ,LE) =  287.05  * DEL(LE)                                  
CX20506 CALL DOOUTX( 506 )                                              
CY20506 CALL DOOUTY( 506 )                                              
      CALL IMINV (HM,  12 , DET, LLL, MMM)                              
CX10507 CALL DOINX( 507 )                                               
CY10507 CALL DOINY( 507 )                                               
      DO 8 I=1, 12                                                      
      DO 8 J=1, 12                                                      
      AM(I,J) = 0.                                                      
      DO 8 K=1, 12                                                      
8     AM(I,J) = AM(I,J) + HM(I,K)*TM(K,J)                               
CX20508 CALL DOOUTX( 508 )                                              
CY20508 CALL DOOUTY( 508 )                                              
CX10509 CALL DOINX( 509 )                                               
CY10509 CALL DOINY( 509 )                                               
C     HERE IS A GOOD PLACE TO DIVIDE BY A**2 FOR LAPLACIAN              
C     STORE AM IN TM AND DIVIDE AM                                      
      DO 10 I=1, 12                                                     
      DO 10 J=1, 12                                                     
      TM(I,J) = AM(I,J)                                                 
      HM(I,J) = AM(I,J)                                                 
10    AM(I,J) = AM(I,J) / ( 6370000.  *  6370000. )                     
CX20510 CALL DOOUTX( 510 )                                              
CY20510 CALL DOOUTY( 510 )                                              
      CALL IMINV(TM, 12 ,DET,LLL,MMM)                                   
CX10511 CALL DOINX( 511 )                                               
CY10511 CALL DOINY( 511 )                                               
      DO 9 LE=1, 12                                                     
9     SV(LE) = -DEL(LE)                                                 
CX20512 CALL DOOUTX( 512 )                                              
CY20512 CALL DOOUTY( 512 )                                              
CX10513 CALL DOINX( 513 )                                               
CY10513 CALL DOINY( 513 )                                               
      DO 11 LE=1, 11                                                    
      P1(LE) = 1. / RPI(LE)                                             
11    P2(LE+1) = RPI(LE)                                                
CX20514 CALL DOOUTX( 514 )                                              
CY20514 CALL DOOUTY( 514 )                                              
      P1( 12 ) = 0.                                                     
      P2 ( 1 ) = 0.                                                     
C     PRINT 333                                                         
      WRITE(6,333)                                                      
333   FORMAT(1H0,'  END SR AMHMTM')                                     
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX40002 CALL SBOUTX (2)                                                 
CY40002 CALL SBOUTY (2)                                                 
      RETURN                                                            
      END                                                               
      BLOCK DATA                                                        
      COMMON /FFTS/ TRIGS(144)                                          
       DATA TRIGS(  1) /  1.000000000000000000000000000000 /            
       DATA TRIGS(  2) /  0.0                              /            
       DATA TRIGS(  3) /  0.991444861373810411144557526929 /            
       DATA TRIGS(  4) /  0.130526192220051591548406227895 /            
       DATA TRIGS(  5) /  0.965925826289068286749743199729 /            
       DATA TRIGS(  6) /  0.258819045102520762348898837624 /            
       DATA TRIGS(  7) /  0.923879532511286756128183189397 /            
       DATA TRIGS(  8) /  0.382683432365089771728459984030 /            
       DATA TRIGS(  9) /  0.866025403784438646763723170753 /            
       DATA TRIGS( 10) /  0.500000000000000000000000000000 /            
       DATA TRIGS( 11) /  0.793353340291235164579776961501 /            
       DATA TRIGS( 12) /  0.608761429008720639416097542898 /            
       DATA TRIGS( 13) /  0.707106781186547524400844362105 /            
       DATA TRIGS( 14) /  0.707106781186547524400844362105 /            
       DATA TRIGS( 15) /  0.608761429008720639416097542898 /            
       DATA TRIGS( 16) /  0.793353340291235164579776961501 /            
       DATA TRIGS( 17) /  0.500000000000000000000000000000 /            
       DATA TRIGS( 18) /  0.866025403784438646763723170753 /            
       DATA TRIGS( 19) /  0.382683432365089771728459984030 /            
       DATA TRIGS( 20) /  0.923879532511286756128183189397 /            
       DATA TRIGS( 21) /  0.258819045102520762348898837624 /            
       DATA TRIGS( 22) /  0.965925826289068286749743199729 /            
       DATA TRIGS( 23) /  0.130526192220051591548406227895 /            
       DATA TRIGS( 24) /  0.991444861373810411144557526929 /            
       DATA TRIGS( 25) /  0.000000000000000000000000000000 /            
       DATA TRIGS( 26) /  1.000000000000000000000000000000 /            
       DATA TRIGS( 27) / -0.130526192220051591548406227895 /            
       DATA TRIGS( 28) /  0.991444861373810411144557526929 /            
       DATA TRIGS( 29) / -0.258819045102520762348898837624 /            
       DATA TRIGS( 30) /  0.965925826289068286749743199729 /            
       DATA TRIGS( 31) / -0.382683432365089771728459984030 /            
       DATA TRIGS( 32) /  0.923879532511286756128183189397 /            
       DATA TRIGS( 33) / -0.500000000000000000000000000000 /            
       DATA TRIGS( 34) /  0.866025403784438646763723170753 /            
       DATA TRIGS( 35) / -0.608761429008720639416097542898 /            
       DATA TRIGS( 36) /  0.793353340291235164579776961501 /            
       DATA TRIGS( 37) / -0.707106781186547524400844362105 /            
       DATA TRIGS( 38) /  0.707106781186547524400844362105 /            
       DATA TRIGS( 39) / -0.793353340291235164579776961501 /            
       DATA TRIGS( 40) /  0.608761429008720639416097542898 /            
       DATA TRIGS( 41) / -0.866025403784438646763723170753 /            
       DATA TRIGS( 42) /  0.500000000000000000000000000000 /            
       DATA TRIGS( 43) / -0.923879532511286756128183189397 /            
       DATA TRIGS( 44) /  0.382683432365089771728459984030 /            
       DATA TRIGS( 45) / -0.965925826289068286749743199729 /            
       DATA TRIGS( 46) /  0.258819045102520762348898837624 /            
       DATA TRIGS( 47) / -0.991444861373810411144557526929 /            
       DATA TRIGS( 48) /  0.130526192220051591548406227896 /            
       DATA TRIGS( 49) / -1.000000000000000000000000000000 /            
       DATA TRIGS( 50) /  0.000000000000000000000000000000 /            
       DATA TRIGS( 51) / -0.991444861373810411144557526929 /            
       DATA TRIGS( 52) / -0.130526192220051591548406227895 /            
       DATA TRIGS( 53) / -0.965925826289068286749743199729 /            
       DATA TRIGS( 54) / -0.258819045102520762348898837624 /            
       DATA TRIGS( 55) / -0.923879532511286756128183189397 /            
       DATA TRIGS( 56) / -0.382683432365089771728459984030 /            
       DATA TRIGS( 57) / -0.866025403784438646763723170753 /            
       DATA TRIGS( 58) / -0.500000000000000000000000000000 /            
       DATA TRIGS( 59) / -0.793353340291235164579776961501 /            
       DATA TRIGS( 60) / -0.608761429008720639416097542898 /            
       DATA TRIGS( 61) / -0.707106781186547524400844362105 /            
       DATA TRIGS( 62) / -0.707106781186547524400844362105 /            
       DATA TRIGS( 63) / -0.608761429008720639416097542898 /            
       DATA TRIGS( 64) / -0.793353340291235164579776961501 /            
       DATA TRIGS( 65) / -0.500000000000000000000000000000 /            
       DATA TRIGS( 66) / -0.866025403784438646763723170753 /            
       DATA TRIGS( 67) / -0.382683432365089771728459984030 /            
       DATA TRIGS( 68) / -0.923879532511286756128183189397 /            
       DATA TRIGS( 69) / -0.258819045102520762348898837624 /            
       DATA TRIGS( 70) / -0.965925826289068286749743199729 /            
       DATA TRIGS( 71) / -0.130526192220051591548406227896 /            
       DATA TRIGS( 72) / -0.991444861373810411144557526929 /            
       DATA TRIGS( 73) / -0.000000000000000000000000000000 /            
       DATA TRIGS( 74) / -1.000000000000000000000000000000 /            
       DATA TRIGS( 75) /  0.130526192220051591548406227895 /            
       DATA TRIGS( 76) / -0.991444861373810411144557526929 /            
       DATA TRIGS( 77) /  0.258819045102520762348898837624 /            
       DATA TRIGS( 78) / -0.965925826289068286749743199729 /            
       DATA TRIGS( 79) /  0.382683432365089771728459984030 /            
       DATA TRIGS( 80) / -0.923879532511286756128183189397 /            
       DATA TRIGS( 81) /  0.500000000000000000000000000000 /            
       DATA TRIGS( 82) / -0.866025403784438646763723170753 /            
       DATA TRIGS( 83) /  0.608761429008720639416097542898 /            
       DATA TRIGS( 84) / -0.793353340291235164579776961501 /            
       DATA TRIGS( 85) /  0.707106781186547524400844362105 /            
       DATA TRIGS( 86) / -0.707106781186547524400844362105 /            
       DATA TRIGS( 87) /  0.793353340291235164579776961501 /            
       DATA TRIGS( 88) / -0.608761429008720639416097542898 /            
       DATA TRIGS( 89) /  0.866025403784438646763723170753 /            
       DATA TRIGS( 90) / -0.500000000000000000000000000000 /            
       DATA TRIGS( 91) /  0.923879532511286756128183189397 /            
       DATA TRIGS( 92) / -0.382683432365089771728459984031 /            
       DATA TRIGS( 93) /  0.965925826289068286749743199729 /            
       DATA TRIGS( 94) / -0.258819045102520762348898837624 /            
       DATA TRIGS( 95) /  0.991444861373810411144557526929 /            
       DATA TRIGS( 96) / -0.130526192220051591548406227896 /            
       DATA TRIGS( 97) /  1.000000000000000000000000000000 /            
       DATA TRIGS( 98) /  0.0                              /            
       DATA TRIGS( 99) /  0.997858923238603506738069791273 /            
       DATA TRIGS(100) /  0.065403129230143066815315558775 /            
       DATA TRIGS(101) /  0.991444861373810411144557526929 /            
       DATA TRIGS(102) /  0.130526192220051591548406227895 /            
       DATA TRIGS(103) /  0.980785280403230449126182236134 /            
       DATA TRIGS(104) /  0.195090322016128267848284868477 /            
       DATA TRIGS(105) /  0.965925826289068286749743199729 /            
       DATA TRIGS(106) /  0.258819045102520762348898837624 /            
       DATA TRIGS(107) /  0.946930129495105664255804274854 /            
       DATA TRIGS(108) /  0.321439465303161580701057624079 /            
       DATA TRIGS(109) /  0.923879532511286756128183189397 /            
       DATA TRIGS(110) /  0.382683432365089771728459984030 /            
       DATA TRIGS(111) /  0.896872741532688303894103936393 /            
       DATA TRIGS(112) /  0.442288690219001281995238977324 /            
       DATA TRIGS(113) /  0.866025403784438646763723170753 /            
       DATA TRIGS(114) /  0.500000000000000000000000000000 /            
       DATA TRIGS(115) /  0.831469612302545237078788377618 /            
       DATA TRIGS(116) /  0.555570233019602224742830813949 /            
       DATA TRIGS(117) /  0.793353340291235164579776961501 /            
       DATA TRIGS(118) /  0.608761429008720639416097542898 /            
       DATA TRIGS(119) /  0.751839807478977396407519406377 /            
       DATA TRIGS(120) /  0.659345815100068868425124612055 /            
       DATA TRIGS(121) /  0.707106781186547524400844362105 /            
       DATA TRIGS(122) /  0.707106781186547524400844362105 /            
       DATA TRIGS(123) /  0.659345815100068868425124612055 /            
       DATA TRIGS(124) /  0.751839807478977396407519406377 /            
       DATA TRIGS(125) /  0.608761429008720639416097542898 /            
       DATA TRIGS(126) /  0.793353340291235164579776961501 /            
       DATA TRIGS(127) /  0.555570233019602224742830813949 /            
       DATA TRIGS(128) /  0.831469612302545237078788377618 /            
       DATA TRIGS(129) /  0.500000000000000000000000000000 /            
       DATA TRIGS(130) /  0.866025403784438646763723170753 /            
       DATA TRIGS(131) /  0.442288690219001281995238977324 /            
       DATA TRIGS(132) /  0.896872741532688303894103936393 /            
       DATA TRIGS(133) /  0.382683432365089771728459984030 /            
       DATA TRIGS(134) /  0.923879532511286756128183189397 /            
       DATA TRIGS(135) /  0.321439465303161580701057624079 /            
       DATA TRIGS(136) /  0.946930129495105664255804274854 /            
       DATA TRIGS(137) /  0.258819045102520762348898837624 /            
       DATA TRIGS(138) /  0.965925826289068286749743199729 /            
       DATA TRIGS(139) /  0.195090322016128267848284868477 /            
       DATA TRIGS(140) /  0.980785280403230449126182236134 /            
       DATA TRIGS(141) /  0.130526192220051591548406227895 /            
       DATA TRIGS(142) /  0.991444861373810411144557526929 /            
       DATA TRIGS(143) /  0.065403129230143066815315558775 /            
       DATA TRIGS(144) /  0.997858923238603506738069791273 /            
      END                                                               
      SUBROUTINE BMCM (TOV,P1,P2,H1,H2,DEL,CI,BM,CM,DT,SV,AM)           
      DIMENSION BM( 12 , 12 ), CM( 12 , 12 ),                           
     1  TOV( 12 ), P1( 12 ), P2( 12 ), H1( 12 ), H2( 12 ),              
     2  DEL( 12 ), CI( 13 ), SV( 12 ), AM( 12 , 12 )                    
CX30003 CALL SBINX (3)                                                  
CY30003 CALL SBINY (3)                                                  
      RK =  287.05  /  1005.                                            
CX10515 CALL DOINX( 515 )                                               
CY10515 CALL DOINY( 515 )                                               
      DO 1 LE=1, 11                                                     
1     H1(LE) = P1(LE) * TOV(LE+1) - TOV(LE)                             
CX20516 CALL DOOUTX( 516 )                                              
CY20516 CALL DOOUTY( 516 )                                              
      H1( 12 ) = 0.                                                     
      H2( 1 ) = 0.                                                      
CX10517 CALL DOINX( 517 )                                               
CY10517 CALL DOINY( 517 )                                               
      DO 2 LE=2, 12                                                     
2     H2(LE) = TOV(LE) - P2(LE) * TOV(LE-1)                             
CX20518 CALL DOOUTX( 518 )                                              
CY20518 CALL DOOUTY( 518 )                                              
CX10519 CALL DOINX( 519 )                                               
CY10519 CALL DOINY( 519 )                                               
      DO 3 K=1, 12                                                      
      KK = K                                                            
      SAV1 = RK*TOV(K) + (CI(K+1)*H1(K)+CI(K)*H2(K)) / (2.*DEL(K))      
      SAV2 = (H1(K)+H2(K)) / (2.*DEL(K))                                
      DO 33 J=1, 12                                                     
33    BM(K,J) = -SAV1*DEL(J)                                            
      DO 34 J=1,KK                                                      
34    BM(K,J) = BM(K,J) + SAV2*DEL(J)                                   
      BM(K,K) = BM(K,K) - H2(K) / 2.                                    
3     CONTINUE                                                          
CX20520 CALL DOOUTX( 520 )                                              
CY20520 CALL DOOUTY( 520 )                                              
      RAA =  287.05  / ( 6370000. * 6370000. )                          
CX10521 CALL DOINX( 521 )                                               
CY10521 CALL DOINY( 521 )                                               
      DO 4 I=1, 12                                                      
      DO 4 J=1, 12                                                      
      CM(I,J) = 0.                                                      
      DO 44 K=1, 12                                                     
44    CM(I,J) = CM(I,J) + AM(I,K) * BM(K,J)                             
      CM(I,J) = (CM(I,J)+RAA*TOV(I)*SV(J))*DT*DT                        
4     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20522 CALL DOOUTX( 522 )                                              
CY20522 CALL DOOUTY( 522 )                                              
CX40004 CALL SBOUTX (4)                                                 
CY40004 CALL SBOUTY (4)                                                 
      RETURN                                                            
      END                                                               
      SUBROUTINE CONKUO(DT,SATCRI,ACUMT,NKUO,PS,DEL,SL,QN,QN1,TN1,      
     1 ACUM,DQ,GESHEM,RELKUO,MSTA)                                      
      DIMENSION GESHEM( 96 )                                            
      DIMENSION PS( 96 ),DEL( 12 ),SL( 12 ),                            
     1QN( 96 , 7 ),QN1( 96 , 7 ),DQ( 96 , 7 ),                          
     2TN1( 96 , 12 )                                                    
      DIMENSION ACUM( 96 ),PRESS( 7 ),TIN( 7 ),QIN( 7 ),                
     1TMST( 7 ),QMST( 7 ),DTKUO( 7 ),DQKUO( 7 ),ESAT( 7 )               
CX30005 CALL SBINX (5)                                                  
CY30005 CALL SBINY (5)                                                  
      IF(DT.EQ.0.) THEN                                                 
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
      RETURN                                                            
      ENDIF                                                             
C***                                                                    
C*** DEL=DEL SIGMA (DEL P OVR PSFC)                                     
C*** PS=SFC PRES (CB),  RELKUO(=.80) IS THE SATURATION CRITERION        
C*** ACMT = TEN TO MINUS SEVEN PER SEC                                  
      RDT=1./DT                                                         
      CPOVL= 1005. /2.5E+6                                              
CX10523 CALL DOINX( 523 )                                               
CY10523 CALL DOINY( 523 )                                               
      DO 5 K=1, 7                                                       
      DO 5 I=1, 96                                                      
5     DQ(I,K)=QN1(I,K)-QN(I,K)                                          
C*** NET TIME RATE MOISTURE INFLOW ..LWST KACUM.. LYRS                  
CX20524 CALL DOOUTX( 524 )                                              
CY20524 CALL DOOUTY( 524 )                                              
      KACUM=4                                                           
CX10525 CALL DOINX( 525 )                                               
CY10525 CALL DOINY( 525 )                                               
      DO  10 I=1, 96                                                    
10    ACUM(I)=ACUMT*0.1                                                 
CX20526 CALL DOOUTX( 526 )                                              
CY20526 CALL DOOUTY( 526 )                                              
CX10527 CALL DOINX( 527 )                                               
CY10527 CALL DOINY( 527 )                                               
      DO 15 K=1,KACUM                                                   
      DO 15 I=1, 96                                                     
15    ACUM(I)=ACUM(I)+RDT*DQ(I,K)*DEL(K)                                
CX20528 CALL DOOUTX( 528 )                                              
CY20528 CALL DOOUTY( 528 )                                              
CX10529 CALL DOINX( 529 )                                               
CY10529 CALL DOINY( 529 )                                               
C***                                                                    
C*** ENTER LON LOOP........                                             
      DO 125 I=1, 96                                                    
C***                                                                    
C*** TEST 1....SUFFICIENT WATER ACCUM? NET INFLO INTO LWST LYR?         
      IF(ACUM(I).LT.0.) GO TO 125                                       
      IF(DQ(I,1) .LE. .0) GO TO 125                                     
      DO 12345 K=1, 7                                                   
12345 PRESS(K)=SL(K)*PS(I)                                              
C*** TIN, QIN ARE PRELIM TMP AND Q AT N+1, ZERO NEGATIVE Q              
      DO 103 K=1, 7                                                     
      TIN(K)=TN1(I,K)                                                   
      QIN(K)=QN1(I,K)                                                   
      IF(QN1(I,K) .LT. .0) QIN(K) = .0                                  
103   CONTINUE                                                          
       IF(TIN(1) .LT. 278.16) GO TO 125                                 
       IF(TIN(2) .GT. TIN(1)) GO TO 125                                 
       CALL SATVAP(TIN,ESAT,2)                                          
       QSATK = .65*.622*ESAT(1)/(PRESS(1)-0.378*ESAT(1))                
       QDEF = QIN(1) - QSATK                                            
C*** LWST LYR MIN  TMP AND RH CONDS PLUS NO TMP INVERSION EXISTS        
       IF(QDEF .LT. .0) GO TO 125                                       
C*** MODIFY TMP FOR MSTADB (SUBR PERMITS AN INTERVENING STBL LYR)       
      TIN(1) =  TIN(1) + .25* (TIN(2) - TIN(1))                         
      HNEW=PS(I)/100.                                                   
      CALL MSTADB(HNEW, SL  ,TIN,1, 7 ,THE,TMST,QMST,KTOP,IER)          
C***  TMST AND QMST,K=1...KTOP CONTAIN CLOUD TEMPS AND SPECIFIC HUMID.  
C**** -------------------------------                                   
      IF(IER)105,106,105                                                
105   MSTA=MSTA+1                                                       
      GO TO 125                                                         
106   CONTINUE                                                          
C**** -------------------------------                                   
C***                                                                    
C*** TEST 2...INSTABILITY.  KTOP .GE. 4,  NO STBL LYRS PERMITTED        
       IF(KTOP.LT.4) GO TO 125                                          
       DO 150 K = 2,KTOP                                                
       X = TMST(K) - TIN(K)                                             
       IF(X .LT. 0.) GO TO 125                                          
150    CONTINUE                                                         
C*** HERE  KUO REQUIREMENTS ARE MET, FIRST COMPUTE WATER=TONS WATER     
C*** SUBSTANCE  ACCUMULATED IN ONE TIME STEP PER SQ M IN LAYERS 2-KTOP  
      NKUO=NKUO+1                                                       
C*** AVL  WATER  (WATER=WATER*G/PS) .GT. .0, COMES FM TWO LWST LYRS     
      WATER=0.                                                          
      DO 107 K=1,2                                                      
      WATER=WATER+   DQ(I,K)*DEL(K)                                     
107   CONTINUE                                                          
      IF(WATER .LE. .0) GO TO 125                                       
C***                                                                    
      Q1=0.                                                             
      Q2=0.                                                             
        DO 111 K=2,KTOP                                                 
        X=QMST(K)-QIN(K)                                                
        IF(X)109,109,108                                                
108     Q1=Q1+X*DEL(K)                                                  
109     X=TMST(K)-TIN(K)                                                
        IF(X)111,111,110                                                
110     Q2=Q2+X*DEL(K)                                                  
111     CONTINUE                                                        
C*** CONTRIB. FACTOR (FRAC CVRG?) .GE. MIN, NOT .GT. 1.0                
      Q2=Q2*CPOVL                                                       
      QEFF=WATER/(Q1+Q2)                                                
      IF(QEFF .LT. .001) GO TO 125                                      
      IF(QEFF .GT. 1.00) QEFF = 1.00                                    
C***                                                                    
C*** CONVECTIVE  CONTRIBUTION TO HEATING AND MOISTENING                 
        DO 304 K=2,KTOP                                                 
        X=QEFF*(TMST(K)-TIN(K))                                         
C***  CHANGE T ONLY IF X GT.0                                           
        IF(X)300,300,301                                                
300     X=0.                                                            
301     DTKUO(K)=X                                                      
        X=QEFF*(QMST(K)-QIN(K))                                         
C***  CHANGE Q ONLY IF X GT.0                                           
        IF(X)302,302,303                                                
302     X=0.                                                            
303     DQKUO(K)=X                                                      
304     CONTINUE                                                        
      DTKUO(1)=0.                                                       
      DQKUO(1)=0.                                                       
      CEVAP=0.                                                          
        DO 305 K=2,KTOP                                                 
        TIN(K)=TIN(K)+DTKUO(K)                                          
        QIN(K)=DQKUO(K)   + QIN(K)                                      
        CEVAP=CEVAP+DQKUO(K)*DEL(K)                                     
305     CONTINUE                                                        
C***                                                                    
C*** CONVECTIVE PCPN  = AVL WATER, REMOVE WATER FM TWO LWST LYRS        
      FALL = WATER - CEVAP                                              
      GESHEM(I)=GESHEM(I)+PS(I)*FALL/9.8                                
       QIN(1) = QIN(1) - DQ(I,1)                                        
       IF(DQ(I,2) .GT. 0.) QIN(2) = QIN(2)-DQ(I,2)                      
C***  RESTORE BASIC FIELDS                                              
      TIN(1) = TN1(I,1)                                                 
      DO 310 K=1,KTOP                                                   
      QN1(I,K)=QIN(K)                                                   
      TN1(I,K)=TIN(K)                                                   
310   CONTINUE                                                          
125   CONTINUE                                                          
CX20530 CALL DOOUTX( 530 )                                              
CY20530 CALL DOOUTY( 530 )                                              
CX40006 CALL SBOUTX (6)                                                 
CY40006 CALL SBOUTY (6)                                                 
C***  END LON LOOP                                                      
C***                                                                    
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
      RETURN                                                            
      END                                                               
      SUBROUTINE DAMP1(X,S,DIM)                                         
      COMPLEX   X( 961 ),DIM( 961 )                                     
      DIMENSION S( 961 )                                                
CX30007 CALL SBINX (7)                                                  
CY30007 CALL SBINY (7)                                                  
CX10531 CALL DOINX( 531 )                                               
CY10531 CALL DOINY( 531 )                                               
      DO 45 I=1, 961                                                    
      X(I)=X(I)+CMPLX(S(I)*REAL(DIM(I)),S(I)*AIMAG(DIM(I)))             
 45   CONTINUE                                                          
CX20532 CALL DOOUTX( 532 )                                              
CY20532 CALL DOOUTY( 532 )                                              
CX40008 CALL SBOUTX (8)                                                 
CY40008 CALL SBOUTY (8)                                                 
      RETURN                                                            
      END                                                               
      SUBROUTINE DAMP2(Y,S,TEM)                                         
      COMPLEX   Y( 31 , 31 ),TEM( 31 , 31 )                             
      DIMENSION S( 31 , 31 )                                            
CX30009 CALL SBINX (9)                                                  
CY30009 CALL SBINY (9)                                                  
CX10533 CALL DOINX( 533 )                                               
CY10533 CALL DOINY( 533 )                                               
      DO 45 LL=1, 31                                                    
      N=LL-2                                                            
      DO 43 I=1, 31                                                     
      N=N+1                                                             
      IF(N.LT. 31 ) GO TO 43                                            
      Y(I,LL)=Y(I,LL)+CMPLX(S(I,LL)* REAL(TEM(I,LL)),                   
     1                      S(I,LL)*AIMAG(TEM(I,LL)))                   
 43   CONTINUE                                                          
 45   CONTINUE                                                          
CX20534 CALL DOOUTX( 534 )                                              
CY20534 CALL DOOUTY( 534 )                                              
CX40010 CALL SBOUTX (10)                                                
CY40010 CALL SBOUTY (10)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE DEL4(DK,TK,S,RT,RM,X,DIM,Y,TEM,W,ZEM,Z,QM)             
      COMPLEX   RT( 31 , 31 , 7 ),RM( 31 , 31 , 7 ),                    
     1X( 31 , 31 , 12 ),DIM( 31 , 31 , 12 ),                            
     2Y( 31 , 31 , 12 ),TEM( 31 , 31 , 12 ),                            
     3W( 31 , 31 , 12 ),ZEM( 31 , 31 , 12 ),                            
     4Z( 31 , 31 ),QM( 31 , 31 )                                        
      DIMENSION S( 31 , 31 )                                            
CX30011 CALL SBINX (11)                                                 
CY30011 CALL SBINY (11)                                                 
      RDK=-DK/( 6370000. * 6370000. )                                   
      RDK=RDK/( 6370000. * 6370000. )                                   
CX10535 CALL DOINX( 535 )                                               
CY10535 CALL DOINY( 535 )                                               
      DO 10 LL=1, 31                                                    
      N=LL-2                                                            
      DO 10 I=1, 31                                                     
      N=N+1                                                             
      S(I,LL)=RDK*N*N*(N+1.)*(N+1.)                                     
10    CONTINUE                                                          
CX20536 CALL DOOUTX( 536 )                                              
CY20536 CALL DOOUTY( 536 )                                              
CX10537 CALL DOINX( 537 )                                               
CY10537 CALL DOINY( 537 )                                               
      DO 30 K=1, 12                                                     
      CALL DAMP1(X(1,1,K),S,DIM(1,1,K))                                 
30    CONTINUE                                                          
CX20538 CALL DOOUTX( 538 )                                              
CY20538 CALL DOOUTY( 538 )                                              
      K= 12                                                             
      CALL DAMP2(Y(1,1,K),S,TEM(1,1,K))                                 
      CALL DAMP2(W(1,1,K),S,ZEM(1,1,K))                                 
      IF(DK.EQ.TK) GO TO 41                                             
      RTK=-TK/( 6370000. * 6370000. )                                   
      RTK=RTK/( 6370000. * 6370000. )                                   
CX10539 CALL DOINX( 539 )                                               
CY10539 CALL DOINY( 539 )                                               
      DO 40 LL=1, 31                                                    
      N=LL-2                                                            
      DO 40 I=1, 31                                                     
      N=N+1                                                             
      S(I,LL)=RTK*N*N*(N+1.)*(N+1.)                                     
40    CONTINUE                                                          
CX20540 CALL DOOUTX( 540 )                                              
CY20540 CALL DOOUTY( 540 )                                              
41    CONTINUE                                                          
CX10541 CALL DOINX( 541 )                                               
CY10541 CALL DOINY( 541 )                                               
C***  DAMP TEMP AND VORTICITY TENDENCIES                                
      DO 50 K=1, 11                                                     
      CALL DAMP2(Y(1,1,K),S,TEM(1,1,K))                                 
      CALL DAMP2(W(1,1,K),S,ZEM(1,1,K))                                 
50    CONTINUE                                                          
CX20542 CALL DOOUTX( 542 )                                              
CY20542 CALL DOOUTY( 542 )                                              
CX10543 CALL DOINX( 543 )                                               
CY10543 CALL DOINY( 543 )                                               
      DO 70 K=1, 7                                                      
      CALL DAMP2(RT(1,1,K),S,RM(1,1,K))                                 
70    CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20544 CALL DOOUTX( 544 )                                              
CY20544 CALL DOOUTY( 544 )                                              
CX40012 CALL SBOUTX (12)                                                
CY40012 CALL SBOUTY (12)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE DELLNP (QLN, DPDPHI, DPDLAM, EPS)                      
C      SUBROUTINE COMPUTES SPHERICAL HARMONIC COEFFICIENTS OF           
C         COS(LAT) * GRAD(LN(PI))                                       
C      FROM SPHERICAL HARMONIC COEFFICIENTS OF LN(PI)                   
C      USED IN NONLIN TERMS A AND B                                     
C      WHEN GRAD(LN(PI)) IS NEEDED IN ZLN AND SIGMADOT                  
C      DIVIDE BY COS(LAT) AFTER EXPANSION ON LAT. CIRCLE,               
      DIMENSION QLN( 31 , 31 ), DPDPHI( 32 , 31 ),                      
     1  DPDLAM( 31 , 31 ), EPS( 32 , 31 )                               
C      NOTE DIFFERENT TRUNCATION IN GRADIENT COMPONENTS.                
      COMPLEX QLN, DPDPHI, DPDLAM                                       
      DOUBLE PRECISION EPS                                              
CX30013 CALL SBINX (13)                                                 
CY30013 CALL SBINY (13)                                                 
CX10545 CALL DOINX( 545 )                                               
CY10545 CALL DOINY( 545 )                                               
C      DPDPHI = COS(LAT) D(LN(PI)) / A*D(LAT)                           
C      DPDLAM =          D(LN(PI)) / A*D(LON) = I*L*Q(L,N)              
      DO 2 LL=1, 31                                                     
      FL=LL-1                                                           
      DO 2 IND=1, 31                                                    
C      RI=CMPLX(0.,1.)                                                  
C      DPDLAM(IND,LL)=RI*FL*QLN(IND,LL)                                 
      DPDLAM(IND,LL)=CMPLX(-AIMAG(QLN(IND,LL))*FL,REAL(QLN(IND,LL))*FL) 
    2 CONTINUE                                                          
CX20546 CALL DOOUTX( 546 )                                              
CY20546 CALL DOOUTY( 546 )                                              
CX10547 CALL DOINX( 547 )                                               
CY10547 CALL DOINY( 547 )                                               
      DO 4 LL=1, 31                                                     
      L=LL-1                                                            
      N=L                                                               
C      DPDPHI(1,LL)=(L+2.) * EPS(2,LL) * QLN(2,LL)                      
      DPDPHI(1,LL)=CMPLX(    (REAL(QLN(2,LL))*((L+2.)*EPS(2,LL))),      
     1        (AIMAG(QLN(2,LL))*((L+2.)*EPS(2,LL))))                    
      DO 3 IND=2, 30                                                    
      N=N+1                                                             
C      DPDPHI(IND,LL)=(N+2.)*EPS(IND+1,LL)*QLN(IND+1,LL)                
C    1   +(1.-N)*EPS(IND,LL)*QLN(IND-1,LL)                              
      DPDPHI(IND,LL)=CMPLX(    ((N+2.)*EPS(IND+1,LL)*                   
     1   REAL(QLN(IND+1,LL)) + (1.-N)*EPS(IND,LL)*REAL(QLN(IND-1,LL))), 
     2       ((N+2.)*EPS(IND+1,LL)*AIMAG(QLN(IND+1,LL)) +               
     3   (1.-N)*EPS(IND,LL)*AIMAG(QLN(IND-1,LL))))                      
    3 CONTINUE                                                          
C     N=N+1                                                             
C      DPDPHI( 31 ,LL)=-(L+ 30 -1) * EPS( 31 ,LL) * QLN( 30 ,LL)        
      DPDPHI( 31 ,LL)=CMPLX(-    ((L+ 30 -1)*EPS( 31 ,LL)*              
     1    REAL(QLN( 30 ,LL))), -    ((L+ 30 -1)*EPS( 31 ,LL)*           
     2    AIMAG(QLN( 30 ,LL))))                                         
C     N=N+1                                                             
C      DPDPHI( 32 ,LL)=-(L+ 30 )*EPS( 32 ,LL)*QLN( 31 ,LL)              
      DPDPHI( 32 ,LL)=CMPLX(-    ((L+ 30 )*EPS( 32 ,LL)*                
     1 REAL(QLN( 31 ,LL))), -    ((L+ 30 )*EPS( 32 ,LL)*                
     2 AIMAG(QLN( 31 ,LL))))                                            
    4 CONTINUE                                                          
CX20548 CALL DOOUTX( 548 )                                              
CY20548 CALL DOOUTY( 548 )                                              
      AA=1./ 6370000.                                                   
CX10549 CALL DOINX( 549 )                                               
CY10549 CALL DOINY( 549 )                                               
      DO 5 J=1, 31                                                      
      DO 5 I=1, 31                                                      
      DPDLAM(I,J)=CMPLX( REAL(DPDLAM(I,J))*AA  ,                        
     1                  AIMAG(DPDLAM(I,J))*AA  )                        
    5 CONTINUE                                                          
CX20550 CALL DOOUTX( 550 )                                              
CY20550 CALL DOOUTY( 550 )                                              
CX10551 CALL DOINX( 551 )                                               
CY10551 CALL DOINY( 551 )                                               
      DO 6 J=1, 31                                                      
      DO 6 I=1, 32                                                      
      DPDPHI(I,J)=CMPLX( REAL(DPDPHI(I,J))*AA  ,                        
     1                  AIMAG(DPDPHI(I,J))*AA  )                        
    6 CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20552 CALL DOOUTX( 552 )                                              
CY20552 CALL DOOUTY( 552 )                                              
CX40014 CALL SBOUTX (14)                                                
CY40014 CALL SBOUTY (14)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE DRUK(EPS,PLN,Q,PS,QF,KDT)                              
      DOUBLE PRECISION EPS,COL                                          
      DIMENSION EPS( 992 ),PLN( 992 )                                   
      DIMENSION PS( 96 ),QF( 96 ),PRS(12)                               
      COMPLEX Q( 961 )                                                  
      DATA IFIR/0/                                                      
CX30015 CALL SBINX (15)                                                 
CY30015 CALL SBINY (15)                                                 
      LON1=1                                                            
      LON2= 96 /4                                                       
      LON3=2*LON2                                                       
      LON4=3*LON2                                                       
C                                                                       
      PI6=0.523598                                                      
      TENDEG=0.174533                                                   
C                                                                       
      LAT=1                                                             
      I=-4                                                              
C     IF(IFIR.EQ.0) PRINT 99,LON1,LON2,LON3,LON4                        
      IF(IFIR.EQ.0) WRITE(6,99)LON1,LON2,LON3,LON4                      
99    FORMAT(1H0,' LON1...LON4 ',4(I4,2X))                              
CX10553 CALL DOINX( 553 )                                               
CY10553 CALL DOINY( 553 )                                               
      DO 1000 LA=1,3                                                    
      COLR=TENDEG+(LA-1)*PI6                                            
      COL=COLR                                                          
C     IF(IFIR.EQ.0) PRINT 98,COLR                                       
      IF(IFIR.EQ.0) WRITE(6,98)COLR                                     
98    FORMAT(1H0,'COLR=',E12.4)                                         
      CALL PLN2(PLN,COL,LAT,EPS)                                        
      CALL SUMPLS(Q,PS,QF,PLN)                                          
      CALL FFS99(PS,QF,1)                                               
      I=I+4                                                             
      PRS(I+1)=10.*EXP(PS(LON1))                                        
      PRS(I+2)=10.*EXP(PS(LON2))                                        
      PRS(I+3)=10.*EXP(PS(LON3))                                        
      PRS(I+4)=10.*EXP(PS(LON4))                                        
1000  CONTINUE                                                          
CX20554 CALL DOOUTX( 554 )                                              
CY20554 CALL DOOUTY( 554 )                                              
C     PRINT 100,KDT,(PRS(I),I=1,12)                                     
      WRITE(6,100)KDT,(PRS(I),I=1,12)                                   
100   FORMAT(1H ,'KDT=',I3,2X,12(F7.1,1X))                              
      IFIR=1                                                            
CX40016 CALL SBOUTX (16)                                                
CY40016 CALL SBOUTY (16)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE EPSLON (EPS,JCAP)                                      
      DIMENSION EPS(1)                                                  
      DOUBLE PRECISION EPS, A                                           
CX30017 CALL SBINX (17)                                                 
CY30017 CALL SBINY (17)                                                 
      JCAP1 = JCAP + 1                                                  
      JCAP2 = JCAP + 2                                                  
CX10555 CALL DOINX( 555 )                                               
CY10555 CALL DOINY( 555 )                                               
      DO 1 LL=1,JCAP1                                                   
      L = LL - 1                                                        
      JLE = (LL-1) * JCAP2                                              
      DO 1 INDE=2,JCAP2                                                 
      N = L + INDE - 1                                                  
      A=(N*N-L*L)/(4.D+0*N*N-1.D+0)                                     
      EPS(JLE+INDE) = SQRT(A)                                           
1     CONTINUE                                                          
CX20556 CALL DOOUTX( 556 )                                              
CY20556 CALL DOOUTY( 556 )                                              
CX10557 CALL DOINX( 557 )                                               
CY10557 CALL DOINY( 557 )                                               
      DO 2 LL=1,JCAP1                                                   
      JLE = (LL-1) * JCAP2                                              
      EPS(JLE+1) = 0.D+0                                                
2     CONTINUE                                                          
CX20558 CALL DOOUTX( 558 )                                              
CY20558 CALL DOOUTY( 558 )                                              
CX40018 CALL SBOUTX (18)                                                
CY40018 CALL SBOUTY (18)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE FFA99 (G,A,W,LOT)                                      
      DIMENSION G(96), A(96), W(96)                                     
C                                                                       
      COMMON /FFTS/ TRIGS(144)                                          
C                                                                       
C                                                                       
C  MULTIPLE FAST FOURIER TRANSFORM - ANALYSIS.                          
C  GOOD FOR ZONAL WAVE NUMBER  30 .                                     
C                                                                       
C     DIMENSION G(96,LOT), A(96,LOT), W(96,LOT)                         
C                                                                       
C   INPUT - LOT SETS OF GRID VALUES IN                                  
C           G(1,J), G(2,J), ..., G(96,J), J=1,...,LOT.                  
C           G ARRAY IS NOT CHANGED BY SUBROUTINE FFA99.                 
C                                                                       
C  OUTPUT - LOT SETS OF COMPLEX COEFFICIENTS IN                         
C           A(1,J), A(2,J), ..., A( 31 *2,J), J=1,...,LOT.              
C           A( 31 *2+1,J), ..., A(96,J), J=1,...,LOT ARE NOT SET.       
C                                                                       
C  W IS A WORK ARRAY.                                                   
C  TRIGS ARRAY IS PRESET WITH SIN AND COS VALUES.                       
C                                                                       
C                                                                       
      DATA HM /-0.5/                                                    
      DATA SIN60/0.866025403784438646763723170753/                      
CX30019 CALL SBINX (19)                                                 
CY30019 CALL SBINY (19)                                                 
C                                                                       
      I=1                                                               
CX10559 CALL DOINX( 559 )                                               
CY10559 CALL DOINY( 559 )                                               
      DO 210 L=1,LOT                                                    
      W(I)=G(I)+(G(I+32)+G(I+64))                                       
      W(I+1)=G(I+1)+(G(I+33)+G(I+65))                                   
      W(I+2)=(G(I)+HM*(G(I+32)+G(I+64)))-(SIN60*(G(I+33)-G(I+65)))      
      W(I+4)=(G(I)+HM*(G(I+32)+G(I+64)))+(SIN60*(G(I+33)-G(I+65)))      
      W(I+3)=(G(I+1)+HM*(G(I+33)+G(I+65)))+(SIN60*(G(I+32)-G(I+64)))    
      W(I+5)=(G(I+1)+HM*(G(I+33)+G(I+65)))-(SIN60*(G(I+32)-G(I+64)))    
      I=I+96                                                            
  210 CONTINUE                                                          
CX20560 CALL DOOUTX( 560 )                                              
CY20560 CALL DOOUTY( 560 )                                              
      IBASE=3                                                           
      JBASE=7                                                           
      KC=5                                                              
CX10561 CALL DOINX( 561 )                                               
CY10561 CALL DOINY( 561 )                                               
      DO 230 KB=3,31,2                                                  
      C1=TRIGS(KB)                                                      
      S1=TRIGS(KB+1)                                                    
      C2=TRIGS(KC)                                                      
      S2=TRIGS(KC+1)                                                    
      I=IBASE                                                           
      J=JBASE                                                           
      DO 220 L=1,LOT                                                    
      W(J)=G(I)+(G(I+32)+G(I+64))                                       
      W(J+1)=G(I+1)+(G(I+33)+G(I+65))                                   
      W(J+2)=                                                           
     *    C1*((G(I)+HM*(G(I+32)+G(I+64)))-(SIN60*(G(I+33)-G(I+65))))    
     *   -S1*((G(I+1)+HM*(G(I+33)+G(I+65)))+(SIN60*(G(I+32)-G(I+64))))  
      W(J+3)=                                                           
     *    S1*((G(I)+HM*(G(I+32)+G(I+64)))-(SIN60*(G(I+33)-G(I+65))))    
     *   +C1*((G(I+1)+HM*(G(I+33)+G(I+65)))+(SIN60*(G(I+32)-G(I+64))))  
      W(J+4)=                                                           
     *    C2*((G(I)+HM*(G(I+32)+G(I+64)))+(SIN60*(G(I+33)-G(I+65))))    
     *   -S2*((G(I+1)+HM*(G(I+33)+G(I+65)))-(SIN60*(G(I+32)-G(I+64))))  
      W(J+5)=                                                           
     *    S2*((G(I)+HM*(G(I+32)+G(I+64)))+(SIN60*(G(I+33)-G(I+65))))    
     *   +C2*((G(I+1)+HM*(G(I+33)+G(I+65)))-(SIN60*(G(I+32)-G(I+64))))  
      I=I+96                                                            
      J=J+96                                                            
  220 CONTINUE                                                          
      IBASE=IBASE+2                                                     
      JBASE=JBASE+6                                                     
      KC=KC+4                                                           
  230 CONTINUE                                                          
CX20562 CALL DOOUTX( 562 )                                              
CY20562 CALL DOOUTY( 562 )                                              
C                                                                       
      IBASE=1                                                           
CX10563 CALL DOINX( 563 )                                               
CY10563 CALL DOINY( 563 )                                               
      DO 320 K=1,3                                                      
      I=IBASE                                                           
      DO 310 L=1,LOT                                                    
      A(I)=(W(I)+W(I+48))+(W(I+24)+W(I+72))                             
      A(I+12)=(W(I)+W(I+48))-(W(I+24)+W(I+72))                          
      A(I+1)=(W(I+1)+W(I+49))+(W(I+25)+W(I+73))                         
      A(I+13)=(W(I+1)+W(I+49))-(W(I+25)+W(I+73))                        
      A(I+6)=(W(I)-W(I+48))-(W(I+25)-W(I+73))                           
      A(I+18)=(W(I)-W(I+48))+(W(I+25)-W(I+73))                          
      A(I+7)=(W(I+1)-W(I+49))+(W(I+24)-W(I+72))                         
      A(I+19)=(W(I+1)-W(I+49))-(W(I+24)-W(I+72))                        
      I=I+96                                                            
  310 CONTINUE                                                          
      IBASE=IBASE+2                                                     
  320 CONTINUE                                                          
CX20564 CALL DOOUTX( 564 )                                              
CY20564 CALL DOOUTY( 564 )                                              
      JBASE=25                                                          
      KC=13                                                             
      KD=19                                                             
CX10565 CALL DOINX( 565 )                                               
CY10565 CALL DOINY( 565 )                                               
      DO 350 KB=7,19,6                                                  
      C1=TRIGS(KB)                                                      
      S1=TRIGS(KB+1)                                                    
      C2=TRIGS(KC)                                                      
      S2=TRIGS(KC+1)                                                    
      C3=TRIGS(KD)                                                      
      S3=TRIGS(KD+1)                                                    
      DO 340 K=1,3                                                      
      I=IBASE                                                           
      J=JBASE                                                           
      DO 330 L=1,LOT                                                    
      A(J)=(W(I)+W(I+48))+(W(I+24)+W(I+72))                             
      A(J+1)=(W(I+1)+W(I+49))+(W(I+25)+W(I+73))                         
      A(J+12)=                                                          
     *    C2*((W(I)+W(I+48))-(W(I+24)+W(I+72)))                         
     *   -S2*((W(I+1)+W(I+49))-(W(I+25)+W(I+73)))                       
      A(J+13)=                                                          
     *    S2*((W(I)+W(I+48))-(W(I+24)+W(I+72)))                         
     *   +C2*((W(I+1)+W(I+49))-(W(I+25)+W(I+73)))                       
      A(J+6)=                                                           
     *    C1*((W(I)-W(I+48))-(W(I+25)-W(I+73)))                         
     *   -S1*((W(I+1)-W(I+49))+(W(I+24)-W(I+72)))                       
      A(J+7)=                                                           
     *    S1*((W(I)-W(I+48))-(W(I+25)-W(I+73)))                         
     *   +C1*((W(I+1)-W(I+49))+(W(I+24)-W(I+72)))                       
      A(J+18)=                                                          
     *    C3*((W(I)-W(I+48))+(W(I+25)-W(I+73)))                         
     *   -S3*((W(I+1)-W(I+49))-(W(I+24)-W(I+72)))                       
      A(J+19)=                                                          
     *    S3*((W(I)-W(I+48))+(W(I+25)-W(I+73)))                         
     *   +C3*((W(I+1)-W(I+49))-(W(I+24)-W(I+72)))                       
      I=I+96                                                            
      J=J+96                                                            
  330 CONTINUE                                                          
      IBASE=IBASE+2                                                     
      JBASE=JBASE+2                                                     
  340 CONTINUE                                                          
      JBASE=JBASE+18                                                    
      KC=KC+12                                                          
      KD=KD+18                                                          
  350 CONTINUE                                                          
CX20566 CALL DOOUTX( 566 )                                              
CY20566 CALL DOOUTY( 566 )                                              
C                                                                       
      IBASE=1                                                           
CX10567 CALL DOINX( 567 )                                               
CY10567 CALL DOINY( 567 )                                               
      DO 420 K=1,12                                                     
      I=IBASE                                                           
      DO 410 L=1,LOT                                                    
      W(I)=(A(I)+A(I+48))+(A(I+24)+A(I+72))                             
      W(I+48)=(A(I)+A(I+48))-(A(I+24)+A(I+72))                          
      W(I+1)=(A(I+1)+A(I+49))+(A(I+25)+A(I+73))                         
      W(I+49)=(A(I+1)+A(I+49))-(A(I+25)+A(I+73))                        
      W(I+24)=(A(I)-A(I+48))-(A(I+25)-A(I+73))                          
      W(I+72)=(A(I)-A(I+48))+(A(I+25)-A(I+73))                          
      W(I+25)=(A(I+1)-A(I+49))+(A(I+24)-A(I+72))                        
      W(I+73)=(A(I+1)-A(I+49))-(A(I+24)-A(I+72))                        
      I=I+96                                                            
  410 CONTINUE                                                          
      IBASE=IBASE+2                                                     
  420 CONTINUE                                                          
CX20568 CALL DOOUTX( 568 )                                              
CY20568 CALL DOOUTY( 568 )                                              
C                                                                       
      DATA SCALE /  0.00520 83333 33333 33333 33333 33333 /             
      DATA SCALF /  0.01041 66666 66666 66666 66666 66667 /             
      DATA SCALG / -0.01041 66666 66666 66666 66666 66667 /             
C                                                                       
C     A(0)                                                              
      I=1                                                               
CX10569 CALL DOINX( 569 )                                               
CY10569 CALL DOINY( 569 )                                               
      DO 510 L=1,LOT                                                    
      A(I)=SCALF*(W(I)+W(I+1))                                          
      A(I+1)=0.0                                                        
      I=I+96                                                            
  510 CONTINUE                                                          
CX20570 CALL DOOUTX( 570 )                                              
CY20570 CALL DOOUTY( 570 )                                              
C     REMAINING WAVENUMBERS                                             
      IBASE=3                                                           
      JBASE=95                                                          
CX10571 CALL DOINX( 571 )                                               
CY10571 CALL DOINY( 571 )                                               
C                                                                       
C                  131 = -2* 30 +191                                    
      DO 530 K=99, 131 ,2                                               
C                          25 GE ' 30 ' LE 46  FOR PLI PREPROCESSOR.    
      I=IBASE                                                           
      J=JBASE                                                           
      C1=TRIGS(K)                                                       
      S1=TRIGS(K+1)                                                     
      DO 520 L=1,LOT                                                    
      A(I)=SCALE*((W(I)+W(J))                                           
     *   +(C1*(W(I+1)+W(J+1))+S1*(W(I)-W(J))))                          
      A(I+1)=SCALE*((C1*(W(I)-W(J))-S1*(W(I+1)+W(J+1)))                 
     *   +(W(J+1)-W(I+1)))                                              
      I=I+96                                                            
      J=J+96                                                            
  520 CONTINUE                                                          
      IBASE=IBASE+2                                                     
      JBASE=JBASE-2                                                     
  530 CONTINUE                                                          
CX20572 CALL DOOUTX( 572 )                                              
CY20572 CALL DOOUTY( 572 )                                              
CX10573 CALL DOINX( 573 )                                               
CY10573 CALL DOINY( 573 )                                               
C                                                                       
C               133 = -2* 30 +193                                       
      DO 550 K= 133 ,143,2                                              
      I=IBASE                                                           
      J=JBASE                                                           
      C1=TRIGS(K)                                                       
      S1=TRIGS(K+1)                                                     
      DO 540 L=1,LOT                                                    
      A(I)=SCALE*((W(I)+W(J))                                           
     *   +(C1*(W(I+1)+W(J+1))+S1*(W(I)-W(J))))                          
      A(J)=SCALE*((W(I)+W(J))                                           
     *   -(C1*(W(I+1)+W(J+1))+S1*(W(I)-W(J))))                          
      A(I+1)=SCALE*((C1*(W(I)-W(J))-S1*(W(I+1)+W(J+1)))                 
     *    +(W(J+1)-W(I+1)))                                             
      A(J+1)=SCALE*((C1*(W(I)-W(J))-S1*(W(I+1)+W(J+1)))                 
     *    -(W(J+1)-W(I+1)))                                             
      I=I+96                                                            
      J=J+96                                                            
  540 CONTINUE                                                          
      IBASE=IBASE+2                                                     
      JBASE=JBASE-2                                                     
  550 CONTINUE                                                          
CX20574 CALL DOOUTX( 574 )                                              
CY20574 CALL DOOUTY( 574 )                                              
C     WAVENUMBER N/4                                                    
      I=49                                                              
CX10575 CALL DOINX( 575 )                                               
CY10575 CALL DOINY( 575 )                                               
      DO 560 L=1,LOT                                                    
      A(I)=SCALF*W(I)                                                   
      A(I+1)=SCALG*W(I+1)                                               
      I=I+96                                                            
  560 CONTINUE                                                          
CX20576 CALL DOOUTX( 576 )                                              
CY20576 CALL DOOUTY( 576 )                                              
CX40020 CALL SBOUTX (20)                                                
CY40020 CALL SBOUTY (20)                                                
C                                                                       
      RETURN                                                            
      END                                                               
      SUBROUTINE FFS99 (A,W,LOT)                                        
      DIMENSION A(96),W(96)                                             
C                                                                       
      COMMON /FFTS/ TRIGS(144)                                          
CX30021 CALL SBINX (21)                                                 
CY30021 CALL SBINY (21)                                                 
C                                                                       
C                                                                       
C  MULTIPLE FAST FOURIER TRANSFORM - SYNTHESIS.                         
C  GOOD FOR ZONAL WAVE NUMBER  30 .                                     
C                                                                       
C     DIMENSION A(96,LOT), W(96,LOT)                                    
C                                                                       
C   INPUT - LOT SETS OF COMPLEX COEFFICIENTS IN                         
C           A(1,J), A(2,J), ..., A( 31 *2,J), J=1,...,LOT.              
C           A( 31 *2+1,J), ..., A(96,J), J=1,...,LOT ARE NOT SET.       
C                                                                       
C  OUTPUT - LOT SETS OF GRID VALUES IN                                  
C           A(1,J), A(2,J), ..., A(96,J), J=1,...,LOT.                  
C                                                                       
C  W IS A WORK ARRAY.                                                   
C  TRIGS ARRAY IS PRESET WITH SIN AND COS VALUES.                       
C                                                                       
C                                                                       
C     A(0)                                                              
      I=1                                                               
CX10577 CALL DOINX( 577 )                                               
CY10577 CALL DOINY( 577 )                                               
      DO 110 L=1,LOT                                                    
      W(I)=A(I)                                                         
      W(I+1)=A(I)                                                       
      I=I+96                                                            
  110 CONTINUE                                                          
CX20578 CALL DOOUTX( 578 )                                              
CY20578 CALL DOOUTY( 578 )                                              
C     REMAINING WAVENUMBERS                                             
      IBASE=3                                                           
      JBASE=95                                                          
CX10579 CALL DOINX( 579 )                                               
CY10579 CALL DOINY( 579 )                                               
C                                                                       
C                  131 = -2* 30 +191                                    
      DO 130 K=99, 131 ,2                                               
C                          25 GE ' 30 ' LE 46  FOR PLI PREPROCESSOR.    
      I=IBASE                                                           
      J=JBASE                                                           
      C1=TRIGS(K)                                                       
      S1=TRIGS(K+1)                                                     
      DO 120 L=1,LOT                                                    
      W(I)=A(I)-(S1*A(I)+C1*A(I+1))                                     
      W(J)=A(I)+(S1*A(I)+C1*A(I+1))                                     
      W(I+1)=(C1*A(I)-S1*A(I+1))+A(I+1)                                 
      W(J+1)=(C1*A(I)-S1*A(I+1))-A(I+1)                                 
      I=I+96                                                            
      J=J+96                                                            
  120 CONTINUE                                                          
      IBASE=IBASE+2                                                     
      JBASE=JBASE-2                                                     
  130 CONTINUE                                                          
CX20580 CALL DOOUTX( 580 )                                              
CY20580 CALL DOOUTY( 580 )                                              
CX10581 CALL DOINX( 581 )                                               
CY10581 CALL DOINY( 581 )                                               
C                                                                       
C               133 = -2* 30 +193                                       
      DO 150 K= 133 ,143,2                                              
      I=IBASE                                                           
      J=JBASE                                                           
      C1=TRIGS(K)                                                       
      S1=TRIGS(K+1)                                                     
      DO 140 L=1,LOT                                                    
      W(I)=(A(I)+A(J))-                                                 
     *    (S1*(A(I)-A(J))+C1*(A(I+1)+A(J+1)))                           
      W(J)=(A(I)+A(J))+                                                 
     *    (S1*(A(I)-A(J))+C1*(A(I+1)+A(J+1)))                           
      W(I+1)=(C1*(A(I)-A(J))-S1*(A(I+1)+A(J+1)))+                       
     *    (A(I+1)-A(J+1))                                               
      W(J+1)=(C1*(A(I)-A(J))-S1*(A(I+1)+A(J+1)))-                       
     *    (A(I+1)-A(J+1))                                               
      I=I+96                                                            
      J=J+96                                                            
  140 CONTINUE                                                          
      IBASE=IBASE+2                                                     
      JBASE=JBASE-2                                                     
  150 CONTINUE                                                          
CX20582 CALL DOOUTX( 582 )                                              
CY20582 CALL DOOUTY( 582 )                                              
C     WAVENUMBER N/4                                                    
      I=49                                                              
CX10583 CALL DOINX( 583 )                                               
CY10583 CALL DOINY( 583 )                                               
      DO 160 L=1,LOT                                                    
      W(I)=A(I)+A(I)                                                    
      W(I+1)=-(A(I+1)+A(I+1))                                           
      I=I+96                                                            
  160 CONTINUE                                                          
CX20584 CALL DOOUTX( 584 )                                              
CY20584 CALL DOOUTY( 584 )                                              
C                                                                       
      DATA HM /-0.5/                                                    
      DATA SIN60/0.866025403784438646763723170753/                      
C                                                                       
      I=1                                                               
CX10585 CALL DOINX( 585 )                                               
CY10585 CALL DOINY( 585 )                                               
      DO 210 L=1,LOT                                                    
      A(I)=W(I)+(W(I+32)+W(I+64))                                       
      A(I+1)=W(I+1)+(W(I+33)+W(I+65))                                   
      A(I+2)=(W(I)+HM*(W(I+32)+W(I+64)))-(SIN60*(W(I+33)-W(I+65)))      
      A(I+4)=(W(I)+HM*(W(I+32)+W(I+64)))+(SIN60*(W(I+33)-W(I+65)))      
      A(I+3)=(W(I+1)+HM*(W(I+33)+W(I+65)))+(SIN60*(W(I+32)-W(I+64)))    
      A(I+5)=(W(I+1)+HM*(W(I+33)+W(I+65)))-(SIN60*(W(I+32)-W(I+64)))    
      I=I+96                                                            
  210 CONTINUE                                                          
CX20586 CALL DOOUTX( 586 )                                              
CY20586 CALL DOOUTY( 586 )                                              
      IBASE=3                                                           
      JBASE=7                                                           
      KC=5                                                              
CX10587 CALL DOINX( 587 )                                               
CY10587 CALL DOINY( 587 )                                               
      DO 230 KB=3,31,2                                                  
      C1=TRIGS(KB)                                                      
      S1=TRIGS(KB+1)                                                    
      C2=TRIGS(KC)                                                      
      S2=TRIGS(KC+1)                                                    
      I=IBASE                                                           
      J=JBASE                                                           
      DO 220 L=1,LOT                                                    
      A(J)=W(I)+(W(I+32)+W(I+64))                                       
      A(J+1)=W(I+1)+(W(I+33)+W(I+65))                                   
      A(J+2)=                                                           
     *    C1*((W(I)+HM*(W(I+32)+W(I+64)))-(SIN60*(W(I+33)-W(I+65))))    
     *   -S1*((W(I+1)+HM*(W(I+33)+W(I+65)))+(SIN60*(W(I+32)-W(I+64))))  
      A(J+3)=                                                           
     *    S1*((W(I)+HM*(W(I+32)+W(I+64)))-(SIN60*(W(I+33)-W(I+65))))    
     *   +C1*((W(I+1)+HM*(W(I+33)+W(I+65)))+(SIN60*(W(I+32)-W(I+64))))  
      A(J+4)=                                                           
     *    C2*((W(I)+HM*(W(I+32)+W(I+64)))+(SIN60*(W(I+33)-W(I+65))))    
     *   -S2*((W(I+1)+HM*(W(I+33)+W(I+65)))-(SIN60*(W(I+32)-W(I+64))))  
      A(J+5)=                                                           
     *    S2*((W(I)+HM*(W(I+32)+W(I+64)))+(SIN60*(W(I+33)-W(I+65))))    
     *   +C2*((W(I+1)+HM*(W(I+33)+W(I+65)))-(SIN60*(W(I+32)-W(I+64))))  
      I=I+96                                                            
      J=J+96                                                            
  220 CONTINUE                                                          
      IBASE=IBASE+2                                                     
      JBASE=JBASE+6                                                     
      KC=KC+4                                                           
  230 CONTINUE                                                          
CX20588 CALL DOOUTX( 588 )                                              
CY20588 CALL DOOUTY( 588 )                                              
C                                                                       
      IBASE=1                                                           
CX10589 CALL DOINX( 589 )                                               
CY10589 CALL DOINY( 589 )                                               
      DO 320 K=1,3                                                      
      I=IBASE                                                           
      DO 310 L=1,LOT                                                    
      W(I)=(A(I)+A(I+48))+(A(I+24)+A(I+72))                             
      W(I+12)=(A(I)+A(I+48))-(A(I+24)+A(I+72))                          
      W(I+1)=(A(I+1)+A(I+49))+(A(I+25)+A(I+73))                         
      W(I+13)=(A(I+1)+A(I+49))-(A(I+25)+A(I+73))                        
      W(I+6)=(A(I)-A(I+48))-(A(I+25)-A(I+73))                           
      W(I+18)=(A(I)-A(I+48))+(A(I+25)-A(I+73))                          
      W(I+7)=(A(I+1)-A(I+49))+(A(I+24)-A(I+72))                         
      W(I+19)=(A(I+1)-A(I+49))-(A(I+24)-A(I+72))                        
      I=I+96                                                            
  310 CONTINUE                                                          
      IBASE=IBASE+2                                                     
  320 CONTINUE                                                          
CX20590 CALL DOOUTX( 590 )                                              
CY20590 CALL DOOUTY( 590 )                                              
      JBASE=25                                                          
      KC=13                                                             
      KD=19                                                             
CX10591 CALL DOINX( 591 )                                               
CY10591 CALL DOINY( 591 )                                               
      DO 350 KB=7,19,6                                                  
      C1=TRIGS(KB)                                                      
      S1=TRIGS(KB+1)                                                    
      C2=TRIGS(KC)                                                      
      S2=TRIGS(KC+1)                                                    
      C3=TRIGS(KD)                                                      
      S3=TRIGS(KD+1)                                                    
      DO 340 K=1,3                                                      
      I=IBASE                                                           
      J=JBASE                                                           
      DO 330 L=1,LOT                                                    
      W(J)=(A(I)+A(I+48))+(A(I+24)+A(I+72))                             
      W(J+1)=(A(I+1)+A(I+49))+(A(I+25)+A(I+73))                         
      W(J+12)=                                                          
     *    C2*((A(I)+A(I+48))-(A(I+24)+A(I+72)))                         
     *   -S2*((A(I+1)+A(I+49))-(A(I+25)+A(I+73)))                       
      W(J+13)=                                                          
     *    S2*((A(I)+A(I+48))-(A(I+24)+A(I+72)))                         
     *   +C2*((A(I+1)+A(I+49))-(A(I+25)+A(I+73)))                       
      W(J+6)=                                                           
     *    C1*((A(I)-A(I+48))-(A(I+25)-A(I+73)))                         
     *   -S1*((A(I+1)-A(I+49))+(A(I+24)-A(I+72)))                       
      W(J+7)=                                                           
     *    S1*((A(I)-A(I+48))-(A(I+25)-A(I+73)))                         
     *   +C1*((A(I+1)-A(I+49))+(A(I+24)-A(I+72)))                       
      W(J+18)=                                                          
     *    C3*((A(I)-A(I+48))+(A(I+25)-A(I+73)))                         
     *   -S3*((A(I+1)-A(I+49))-(A(I+24)-A(I+72)))                       
      W(J+19)=                                                          
     *    S3*((A(I)-A(I+48))+(A(I+25)-A(I+73)))                         
     *   +C3*((A(I+1)-A(I+49))-(A(I+24)-A(I+72)))                       
      I=I+96                                                            
      J=J+96                                                            
  330 CONTINUE                                                          
      IBASE=IBASE+2                                                     
      JBASE=JBASE+2                                                     
  340 CONTINUE                                                          
      JBASE=JBASE+18                                                    
      KC=KC+12                                                          
      KD=KD+18                                                          
  350 CONTINUE                                                          
CX20592 CALL DOOUTX( 592 )                                              
CY20592 CALL DOOUTY( 592 )                                              
C                                                                       
      IBASE=1                                                           
CX10593 CALL DOINX( 593 )                                               
CY10593 CALL DOINY( 593 )                                               
      DO 420 K=1,12                                                     
      I=IBASE                                                           
      DO 410 L=1,LOT                                                    
      A(I)=(W(I)+W(I+48))+(W(I+24)+W(I+72))                             
      A(I+48)=(W(I)+W(I+48))-(W(I+24)+W(I+72))                          
      A(I+1)=(W(I+1)+W(I+49))+(W(I+25)+W(I+73))                         
      A(I+49)=(W(I+1)+W(I+49))-(W(I+25)+W(I+73))                        
      A(I+24)=(W(I)-W(I+48))-(W(I+25)-W(I+73))                          
      A(I+72)=(W(I)-W(I+48))+(W(I+25)-W(I+73))                          
      A(I+25)=(W(I+1)-W(I+49))+(W(I+24)-W(I+72))                        
      A(I+73)=(W(I+1)-W(I+49))-(W(I+24)-W(I+72))                        
      I=I+96                                                            
  410 CONTINUE                                                          
      IBASE=IBASE+2                                                     
  420 CONTINUE                                                          
CX20594 CALL DOOUTX( 594 )                                              
CY20594 CALL DOOUTY( 594 )                                              
CX40022 CALL SBOUTX (22)                                                
CY40022 CALL SBOUTY (22)                                                
C                                                                       
      RETURN                                                            
      END                                                               
      SUBROUTINE FILTER (TEM,TE,Y,DIM,DI,X,ZEM,ZE,W,QM,Q,Z,FILTA)       
      DIMENSION                                                         
     1  Q( 961 ), TE( 961 , 12 ), DI( 961 , 12 ), ZE( 961 , 12 ),       
     2  Z( 961 ),  Y( 961 , 12 ),  X( 961 , 12 ),  W( 961 , 12 ),       
     3 QM( 961 ),TEM( 961 , 12 ),DIM( 961 , 12 ),ZEM( 961 , 12 )        
      COMPLEX TEM, TE, Y, DIM, DI, X, ZEM, ZE, W, QM, Q, Z              
CX30023 CALL SBINX (23)                                                 
CY30023 CALL SBINY (23)                                                 
      FILTB = (1. - FILTA) * 0.5                                        
CX10595 CALL DOINX( 595 )                                               
CY10595 CALL DOINY( 595 )                                               
C     QM CONTAINS N-1,  Q CONTAINS N,  Z CONTAINS N+1  ETC.             
      DO 21 LE=1, 12                                                    
      DO 21 I=1, 961                                                    
      TEM(I,LE) = CMPLX(                                                
     1 FILTA* REAL(TE(I,LE)) + FILTB* REAL(TEM(I,LE)),                  
     2 FILTA*AIMAG(TE(I,LE)) + FILTB*AIMAG(TEM(I,LE)))                  
21    CONTINUE                                                          
CX20596 CALL DOOUTX( 596 )                                              
CY20596 CALL DOOUTY( 596 )                                              
CX10597 CALL DOINX( 597 )                                               
CY10597 CALL DOINY( 597 )                                               
      DO 22 LE=1, 12                                                    
      DO 22  I=1, 961                                                   
      DIM(I,LE) = CMPLX(                                                
     1 FILTA* REAL(DI(I,LE)) + FILTB*( REAL(X(I,LE))+ REAL(DIM(I,LE))), 
     2 FILTA*AIMAG(DI(I,LE)) + FILTB*(AIMAG(X(I,LE))+AIMAG(DIM(I,LE)))) 
22    CONTINUE                                                          
CX20598 CALL DOOUTX( 598 )                                              
CY20598 CALL DOOUTY( 598 )                                              
CX10599 CALL DOINX( 599 )                                               
CY10599 CALL DOINY( 599 )                                               
      DO 23 LE=1, 12                                                    
      DO 23  I=1, 961                                                   
      ZEM(I,LE) = CMPLX(                                                
     1 FILTA* REAL(ZE(I,LE)) + FILTB*( REAL(W(I,LE))+ REAL(ZEM(I,LE))), 
     2 FILTA*AIMAG(ZE(I,LE)) + FILTB*(AIMAG(W(I,LE))+AIMAG(ZEM(I,LE)))) 
23    CONTINUE                                                          
CX20600 CALL DOOUTX( 600 )                                              
CY20600 CALL DOOUTY( 600 )                                              
CX10601 CALL DOINX( 601 )                                               
CY10601 CALL DOINY( 601 )                                               
      DO 24  I=1, 961                                                   
      QM(I) = CMPLX(                                                    
     1 FILTA* REAL(Q(I)) + FILTB*( REAL(Z(I))+ REAL(QM(I))),            
     2 FILTA*AIMAG(Q(I)) + FILTB*(AIMAG(Z(I))+AIMAG(QM(I))))            
24    CONTINUE                                                          
CX20602 CALL DOOUTX( 602 )                                              
CY20602 CALL DOOUTY( 602 )                                              
CX10603 CALL DOINX( 603 )                                               
CY10603 CALL DOINY( 603 )                                               
      DO 1003 LEV=1, 12                                                 
      DO 1003 I=1, 961                                                  
      DI(I,LEV)=X(I,LEV)                                                
      TE(I,LEV)=Y(I,LEV)                                                
      ZE(I,LEV)=W(I,LEV)                                                
1003  CONTINUE                                                          
CX20604 CALL DOOUTX( 604 )                                              
CY20604 CALL DOOUTY( 604 )                                              
CX10605 CALL DOINX( 605 )                                               
CY10605 CALL DOINY( 605 )                                               
      DO 1004 I=1, 961                                                  
      Q(I) = Z(I)                                                       
1004  CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20606 CALL DOOUTX( 606 )                                              
CY20606 CALL DOOUTY( 606 )                                              
CX40024 CALL SBOUTX (24)                                                
CY40024 CALL SBOUTY (24)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE FIXTET(PS,WORK,DEL,RPI,RPIREC,TG,QG,TETA,NUN,SL100K)   
      DIMENSION DEL( 12 ),RPI( 11 ),RPIREC( 11 ),                       
     1    PS( 96 ),                                                     
     2  WORK( 96 ),                                                     
     3    TG( 96 , 12 ),                                                
     4   TETA( 96 , 12 ),                                               
     5   QG( 96 , 7 )                                                   
      DIMENSION TET( 12 ),Q( 7 )                                        
CX30025 CALL SBINX (25)                                                 
CY30025 CALL SBINY (25)                                                 
C***  PS IS ASSUMED TRUE SURF. PRESS. IN CB.----------------------      
      RK =  287.05 / 1005.                                              
CX10607 CALL DOINX( 607 )                                               
CY10607 CALL DOINY( 607 )                                               
      DO 1 I=1, 96                                                      
      WORK(I) = SL100K*PS(I)**RK                                        
1     CONTINUE                                                          
CX20608 CALL DOOUTX( 608 )                                              
CY20608 CALL DOOUTY( 608 )                                              
CX10609 CALL DOINX( 609 )                                               
CY10609 CALL DOINY( 609 )                                               
C     COMPUTE BOTTOM POT. TEMP. IN TETA( 96 ,1)                         
      DO 2 I=1, 96                                                      
      TETA(I,1) = TG(I,1)/WORK(I)                                       
2     CONTINUE                                                          
CX20610 CALL DOOUTX( 610 )                                              
CY20610 CALL DOOUTY( 610 )                                              
CX10611 CALL DOINX( 611 )                                               
CY10611 CALL DOINY( 611 )                                               
C     COMPUTE REST OF POT TEMPS FROM BOTTOM ONE UP                      
      DO 4 K=1, 11                                                      
      DO 3 I=1, 96                                                      
      TETA(I,K+1) = TG(I,K+1)*TETA(I,K)*RPIREC(K)/TG(I,K)               
3     CONTINUE                                                          
4     CONTINUE                                                          
CX20612 CALL DOOUTX( 612 )                                              
CY20612 CALL DOOUTY( 612 )                                              
CX10613 CALL DOINX( 613 )                                               
CY10613 CALL DOINY( 613 )                                               
C     NOW HAVE POT. TEMP. IN TETA, ADJUST.                              
C     .................................                                 
C***   QG IS MIXING RATIO ---------------------                         
C     .................................                                 
      DO 66666 I=1, 96                                                  
      DO 5 K=1, 12                                                      
5     TET(K) = TETA(I,K)                                                
      DO 6 K=1, 7                                                       
6     Q(K)=QG(I,K)                                                      
      ITEST=0                                                           
C     TET(K) = POT. TEMP. AT LON.                                       
55555  CONTINUE                                                         
      DO    13 K=2, 12                                                  
      KK = K                                                            
      IF (TET(K)-TET(K-1)) 14, 13, 13                                   
   13  CONTINUE                                                         
      GO TO    26                                                       
   14 KBOT = KK - 1                                                     
      KTOP = KK                                                         
   15 SUMTET = 0.                                                       
      SUMDEL = 0.                                                       
      SUMCOL = 0.                                                       
      SUMQ = 0.                                                         
      ITEST=1                                                           
      DO    16 K=KBOT,KTOP                                              
      SUMDEL = SUMDEL + DEL(K)                                          
   16 SUMTET = SUMTET + DEL(K)*TET(K)                                   
      TET(KBOT) = SUMTET/SUMDEL                                         
      IF (KTOP.GT. 7 )  GO TO 19                                        
      DO    17 K=KBOT,KTOP                                              
   17 SUMQ = SUMQ + DEL(K)*Q(K)                                         
      Q(KBOT) = SUMQ/SUMDEL                                             
      DO    18 K=KK,KTOP                                                
   18 Q(K) = Q(KBOT)                                                    
      GO TO    22                                                       
   19 IF (KBOT.GE. 7 )  GO TO 22                                        
      DO    20 K=KBOT, 7                                                
      SUMCOL = SUMCOL + DEL(K)                                          
      SUMQ = SUMQ + DEL(K)*Q(K)                                         
   20 CONTINUE                                                          
      Q(KBOT) = SUMQ/SUMCOL                                             
      DO    21 K=KK, 7                                                  
      Q(K) = Q(KBOT)                                                    
   21 CONTINUE                                                          
   22 CONTINUE                                                          
      DO      23      K=KK,KTOP                                         
   23 TET(K) = TET(KBOT)                                                
      IF (KTOP- 12 ) 24, 26, 26                                         
   24 IF (TET(KTOP+1)-TET(KTOP)) 25, 55555, 55555                       
   25 KTOP = KTOP + 1                                                   
      GO TO    15                                                       
   26 CONTINUE                                                          
      IF(ITEST)27,66666,27                                              
27    NUN=NUN+1                                                         
C     RETURN ADJUSTED POTENTIAL TEMPS TO TETA                           
      DO 7 K=1, 12                                                      
7     TETA(I,K) = TET(K)                                                
      DO 8 K=1, 7                                                       
8     QG(I,K)=Q(K)                                                      
66666   CONTINUE                                                        
CX20614 CALL DOOUTX( 614 )                                              
CY20614 CALL DOOUTY( 614 )                                              
CX10615 CALL DOINX( 615 )                                               
CY10615 CALL DOINY( 615 )                                               
C     .................................                                 
C     .................................                                 
C     CONVERT POT. TEMPS. TO TEMP. IN TG( 96 , 12 )                     
C     START BY COMPUTING BOTTOM TEMP. AND BUILD UP.                     
      DO 9 I=1, 96                                                      
      TG(I,1) = TETA(I,1)*WORK(I)                                       
9     CONTINUE                                                          
CX20616 CALL DOOUTX( 616 )                                              
CY20616 CALL DOOUTY( 616 )                                              
CX10617 CALL DOINX( 617 )                                               
CY10617 CALL DOINY( 617 )                                               
      DO 11 K=1, 11                                                     
      DO 10 I=1, 96                                                     
      TG(I,K+1) = TG(I,K)*TETA(I,K+1)*RPI(K)/TETA(I,K)                  
10    CONTINUE                                                          
11    CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20618 CALL DOOUTX( 618 )                                              
CY20618 CALL DOOUTY( 618 )                                              
CX40026 CALL SBOUTX (26)                                                
CY40026 CALL SBOUTY (26)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE FL22(FP,FM,FLN,PLN)                                    
C     SR. COMPUTES CONTRIBUTION TO FLN FROM 2 LATITUDES                 
C     FP CONTAINS COEFS. OF SYMMETRIC PLN                               
C     FM CONTAINS COEFS. OF ANTISYM. PLN                                
C     WG IS THE GAUSSIAN WEIGHT                                         
      DIMENSION FP( 31 ), FM( 31 ),                                     
     1FLN( 31 , 31 ),PLN( 32 , 31 )                                     
      COMPLEX FP,FM,FLN                                                 
CX30027 CALL SBINX (27)                                                 
CY30027 CALL SBINY (27)                                                 
CX10619 CALL DOINX( 619 )                                               
CY10619 CALL DOINY( 619 )                                               
      DO 4 LL=1, 31                                                     
      DO 2 I=1, 31 ,2                                                   
      FLN(I,LL) = FLN(I,LL)+CMPLX(                                      
     1REAL(FP(LL))   *PLN(I,LL),                                        
     2AIMAG(FP(LL))   *PLN(I,LL))                                       
2     CONTINUE                                                          
      DO 3 I=2, 30  ,2                                                  
      FLN(I,LL)=FLN(I,LL)+CMPLX(                                        
     1REAL(FM(LL))   *PLN(I,LL),                                        
     2AIMAG(FM(LL))   *PLN(I,LL))                                       
3     CONTINUE                                                          
4     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20620 CALL DOOUTX( 620 )                                              
CY20620 CALL DOOUTY( 620 )                                              
CX40028 CALL SBOUTX (28)                                                
CY40028 CALL SBOUTY (28)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GFIDI(DG,TG,UG,VG,ZG,RQG,PS,DPHI,DLAM,CG,              
     * A,B,F,G,EG,RTG,TAU,UQ,VQ,RCL,SINLAT,CD,TSEA,U1,V1,T1,Q1)         
      DIMENSION PS( 96 ),DPHI( 96 ),DLAM( 96 ),                         
     * DG( 96 , 12 ), ZG( 96 , 12 ), TG( 96 , 12 ),                     
     * UG( 96 , 12 ), VG( 96 , 12 ), CG( 96 , 12 ),                     
     * EG( 96 , 12 ),TAU( 96 , 12 ),                                    
     *  A( 96 , 12 ),  B( 96 , 12 ),                                    
     *  F( 96 , 12 ),  G( 96 , 12 ),                                    
     *RQG( 96 , 7 ),RTG( 96 , 7 ),                                      
     * UQ( 96 , 7 ), VQ( 96 , 7 )                                       
      COMMON/VERCOM/AM( 12 , 12 ),HM( 12 , 12 ),TM( 12 , 12 ),          
     O              BM( 12 , 12 ),CM( 12 , 12 ),EKIN( 12 ),             
     1 SI( 13 ),SL( 12 ),DEL( 12 ),RDEL2( 12 ),RMSDOT( 11 ),            
     2 CI( 13 ),CL( 12 ),TOV( 12 ),   SV( 12 ),   RPI( 11 ),            
     3 P1( 12 ),P2( 12 ), H1( 12 ),   H2( 12 ),RPIREC( 11 ),            
     4 ROTSIN,ROTCOS,SEADRY,SL1KAP,C1,C2,RLRV,SL100K                    
      COMMON/MXSTP/ MSTP,ITERS                                          
      DIMENSION CD( 96 ),TSEA( 96 ),                                    
     * U1( 96 ),V1( 96 ),T1( 96 ),Q1( 96 )                              
      DIMENSION DB( 12 ),CB( 12 ),DOT( 13 )                             
      DIMENSION DUP( 12 ),DVP( 12 ),DUM( 12 ),DVM( 12 )                 
CX30029 CALL SBINX (29)                                                 
CY30029 CALL SBINY (29)                                                 
CX10621 CALL DOINX( 621 )                                               
CY10621 CALL DOINY( 621 )                                               
C      DATA DOT/ 13 *0./,                                               
C     1DUM/ 12 *0./,DVM/ 12 *0./,DUP/ 12 *0./,DVP/ 12 *0./              
      DO 77 K=1,12                                                      
      DOT(K)=0.                                                         
      DUM(K)=0.                                                         
      DVM(K)=0.                                                         
      DUP(K)=0.                                                         
      DVP(K)=0.                                                         
77    CONTINUE                                                          
CX20622 CALL DOOUTX( 622 )                                              
CY20622 CALL DOOUTY( 622 )                                              
      DOT(13)=0.                                                        
      CALL FFS99(DG,CG, 12 )                                            
      CALL FFS99(TG,CG, 12 )                                            
      CALL FFS99(UG,CG, 12 )                                            
      CALL FFS99(VG,CG, 12 )                                            
      CALL FFS99(ZG,CG, 12 )                                            
      CALL FFS99(RQG,CG, 7 )                                            
      CALL FFS99(PS,CG,1)                                               
      CALL FFS99(DPHI,CG,1)                                             
      CALL FFS99(DLAM,CG,1)                                             
C*** ABOVE ASSUMES PS DPHI DLAM ARE CONSECUTIVE IN CORE                 
      RCL2=0.5*RCL                                                      
CX10623 CALL DOINX( 623 )                                               
CY10623 CALL DOINY( 623 )                                               
      DO 2 LE=1, 12                                                     
      DO 1 LO=1, 96                                                     
      A(LO,LE)=ZG(LO,LE)*UG(LO,LE)+ 287.05 *TG(LO,LE)*DPHI(LO)          
      B(LO,LE)=ZG(LO,LE)*VG(LO,LE)- 287.05 *TG(LO,LE)*DLAM(LO)          
1     CONTINUE                                                          
      DO 11 LO=1, 96                                                    
      F(LO,LE)=TG(LO,LE)*UG(LO,LE)                                      
      G(LO,LE)=TG(LO,LE)*VG(LO,LE)                                      
      ZG(LO,LE)=(UG(LO,LE)*UG(LO,LE)+VG(LO,LE)*VG(LO,LE))*RCL2          
11    CONTINUE                                                          
2     CONTINUE                                                          
CX20624 CALL DOOUTX( 624 )                                              
CY20624 CALL DOOUTY( 624 )                                              
      CALL FFA99(ZG,EG,TAU, 12 )                                        
CX10625 CALL DOINX( 625 )                                               
CY10625 CALL DOINY( 625 )                                               
      DO 201 K=1, 7                                                     
      DO 200 LO=1, 96                                                   
      RTG(LO,K)=UG(LO,K)*RQG(LO,K)                                      
200   CONTINUE                                                          
201   CONTINUE                                                          
CX20626 CALL DOOUTX( 626 )                                              
CY20626 CALL DOOUTY( 626 )                                              
      CALL FFA99(RTG,UQ,CG, 7 )                                         
CX10627 CALL DOINX( 627 )                                               
CY10627 CALL DOINY( 627 )                                               
      DO 203 K=1, 7                                                     
      DO 202 LO=1, 96                                                   
      RTG(LO,K)=VG(LO,K)*RQG(LO,K)                                      
202   CONTINUE                                                          
203   CONTINUE                                                          
CX20628 CALL DOOUTX( 628 )                                              
CY20628 CALL DOOUTY( 628 )                                              
      CALL FFA99(RTG,VQ,CG, 7 )                                         
CX10629 CALL DOINX( 629 )                                               
CY10629 CALL DOINY( 629 )                                               
C     COMPUTE C=V(TRUE)*DEL(LN(PS)).DIVIDE BY COS FOR DEL COS FOR V     
      DO 3 LO=1, 96                                                     
      DPHI(LO)=DPHI(LO)*RCL                                             
      DLAM(LO)=DLAM(LO)*RCL                                             
3     CONTINUE                                                          
CX20630 CALL DOOUTX( 630 )                                              
CY20630 CALL DOOUTY( 630 )                                              
CX10631 CALL DOINX( 631 )                                               
CY10631 CALL DOINY( 631 )                                               
      DO 5 LE=1, 12                                                     
      DO 4 LO=1, 96                                                     
      CG(LO,LE)=UG(LO,LE)*DLAM(LO)+VG(LO,LE)*DPHI(LO)                   
4     CONTINUE                                                          
5     CONTINUE                                                          
C----- LOWER LAYER DRAG ----------------------------                    
CX20632 CALL DOOUTX( 632 )                                              
CY20632 CALL DOOUTY( 632 )                                              
      CALL FFS99(U1,TAU,1)                                              
      CALL FFS99(V1,TAU,1)                                              
      CALL FFS99(T1,TAU,1)                                              
      CALL FFS99(Q1,TAU,1)                                              
      ROTCS=ROTSIN*SINLAT                                               
CX10633 CALL DOINX( 633 )                                               
CY10633 CALL DOINY( 633 )                                               
      DO 50 I=1, 96                                                     
      DLAM(I)=CD(I)*SQRT(ZG(I,1))                                       
50    CONTINUE                                                          
CX20634 CALL DOOUTX( 634 )                                              
CY20634 CALL DOOUTY( 634 )                                              
CX10635 CALL DOINX( 635 )                                               
CY10635 CALL DOINY( 635 )                                               
      DO 52 I=1, 96                                                     
      A(I,1)=A(I,1)+DLAM(I)*(V1(I)*ROTCOS+U1(I)*ROTCS)                  
      B(I,1)=B(I,1)-DLAM(I)*(U1(I)*ROTCOS-V1(I)*ROTCS)                  
52    CONTINUE                                                          
CX20636 CALL DOOUTX( 636 )                                              
CY20636 CALL DOOUTY( 636 )                                              
CX10637 CALL DOINX( 637 )                                               
CY10637 CALL DOINY( 637 )                                               
      DO 54 I=1, 96                                                     
      DLAM(I)=DLAM(I)+SEADRY*ZG(I,1)                                    
      IF(TSEA(I).LT.0.) DLAM(I)=0.                                      
54    CONTINUE                                                          
CX20638 CALL DOOUTX( 638 )                                              
CY20638 CALL DOOUTY( 638 )                                              
C***  DLAM NOW CONTAINS LAND SEA TERM. USE IN EVAP. AND SURF. HEATING   
C---------------- END FRICTION ---------                                
      RK= 287.05 / 1005.                                                
CX10639 CALL DOINX( 639 )                                               
CY10639 CALL DOINY( 639 )                                               
      DO 1000 LO=1, 96                                                  
      DB(1)=DEL(1)*DG(LO,1)                                             
      CB(1)=DEL(1)*CG(LO,1)                                             
      DO 6 LE=1, 11                                                     
      DB(LE+1)=DB(LE)+DEL(LE+1)*DG(LO,LE+1)                             
      CB(LE+1)=CB(LE)+DEL(LE+1)*CG(LO,LE+1)                             
6     CONTINUE                                                          
C---  STORE INTEGRAL OF CG IN DPHI ----------                           
      DPHI(LO)=-CB( 12 )                                                
C                                                                       
C     SIGMA DOT COMPUTED ONLY AT INTERIOR INTERFACES                    
      DO 7 K=1, 11                                                      
      DOT(K+1)=DOT(K)+DEL(K)*(DB( 12 )+CB( 12 )-DG(LO,K)-CG(LO,K))      
7     CONTINUE                                                          
      DO 8 K=1, 12                                                      
      TAU(LO,K)=TG(LO,K)*(DG(LO,K)-RK*DB( 12 ))+                        
     1 RK*(TOV(K)+TG(LO,K))*(CG(LO,K)-CB( 12 ))                         
8      CONTINUE                                                         
      DO 12 K=1, 11                                                     
      DVP(K)=VG(LO,K+1)-VG(LO,K)                                        
      DUP(K)=UG(LO,K+1)-UG(LO,K)                                        
12    CONTINUE                                                          
      DO 13 K=1, 11                                                     
      DVM(K+1)=VG(LO,K+1)-VG(LO,K)                                      
      DUM(K+1)=UG(LO,K+1)-UG(LO,K)                                      
13    CONTINUE                                                          
      DO 14 K=1, 12                                                     
      A(LO,K)=A(LO,K)+ RDEL2(K)*(DOT(K+1)*DVP(K)+DOT(K)*DVM(K))         
14    CONTINUE                                                          
      DO 15 K=1, 12                                                     
      B(LO,K)=B(LO,K)- RDEL2(K)*(DOT(K+1)*DUP(K)+DOT(K)*DUM(K))         
15    CONTINUE                                                          
      DO 16 K=1, 11                                                     
      DUP(K)=P1(K)*TG(LO,K+1)-TG(LO,K)                                  
      DVP(K)=CI(K+1)*CB( 12 )-CB(K)                                     
16    CONTINUE                                                          
      DO 17 K=1, 11                                                     
      DUM(K+1)=TG(LO,K+1)-P2(K+1)*TG(LO,K)                              
      DVM(K+1)=CI(K+1)*CB( 12 )-CB(K)                                   
17    CONTINUE                                                          
      DO 18 K=1, 12                                                     
C     H1 H2 COME FROM BMCM                                              
      TAU(LO,K)=TAU(LO,K)-RDEL2(K)*(DOT(K+1)*DUP(K)+DOT(K)*DUM(K) +     
     1H1(K)*DVP(K) + H2(K)*DVM(K)  )                                    
18    CONTINUE                                                          
      DO 30 K=1, 6                                                      
      DUP(K)=RQG(LO,K+1)-RQG(LO,K)                                      
30    CONTINUE                                                          
      DUP( 7 )=0.                                                       
      DO 32 K=1, 6                                                      
      DUM(K+1)=RQG(LO,K+1)-RQG(LO,K)                                    
32    CONTINUE                                                          
      DO 34 K=1, 7                                                      
      RTG(LO,K)=DG(LO,K)*RQG(LO,K)   -                                  
     1     RDEL2(K)*(DOT(K+1)*DUP(K)+DOT(K)*DUM(K))                     
34    CONTINUE                                                          
1000  CONTINUE                                                          
C*** EVAPORATION INTO LOWEST LAYER -----                                
CX20640 CALL DOOUTX( 640 )                                              
CY20640 CALL DOOUTY( 640 )                                              
CX10641 CALL DOINX( 641 )                                               
CY10641 CALL DOINY( 641 )                                               
      DO 56 I=1, 96                                                     
      U1(I)=PS(I)-RLRV/TSEA(I)                                          
56    CONTINUE                                                          
CX20642 CALL DOOUTX( 642 )                                              
CY20642 CALL DOOUTY( 642 )                                              
CX10643 CALL DOINX( 643 )                                               
CY10643 CALL DOINY( 643 )                                               
      DO 58 I=1, 96                                                     
      U1(I)=EXP(U1(I))                                                  
58    CONTINUE                                                          
CX20644 CALL DOOUTX( 644 )                                              
CY20644 CALL DOOUTY( 644 )                                              
CX10645 CALL DOINX( 645 )                                               
CY10645 CALL DOINY( 645 )                                               
      DO 60 I=1, 96                                                     
      V1(I)=C1*0.9/(U1(I)+C2)                                           
60    CONTINUE                                                          
C*** V1 IS QSAT AT SEA TEMP SCALED BY 0.9                               
CX20646 CALL DOOUTX( 646 )                                              
CY20646 CALL DOOUTY( 646 )                                              
CX10647 CALL DOINX( 647 )                                               
CY10647 CALL DOINY( 647 )                                               
C**** QS AT LAYER 1 IS USED TO SHUT OFF EVAP. INTO SATURATED LAYER      
      DO 61 I=1, 96                                                     
      U1(I)=PS(I)-RLRV/T1(I)                                            
61    CONTINUE                                                          
CX20648 CALL DOOUTX( 648 )                                              
CY20648 CALL DOOUTY( 648 )                                              
CX10649 CALL DOINX( 649 )                                               
CY10649 CALL DOINY( 649 )                                               
      DO 63 I=1, 96                                                     
      U1(I)=EXP(U1(I))                                                  
63    CONTINUE                                                          
CX20650 CALL DOOUTX( 650 )                                              
CY20650 CALL DOOUTY( 650 )                                              
CX10651 CALL DOINX( 651 )                                               
CY10651 CALL DOINY( 651 )                                               
      DO 64 I=1, 96                                                     
      U1(I)=C1*0.9/(SL(1)*U1(I)+C2)                                     
64    CONTINUE                                                          
C*** U1 IS SAT Q AT LAYER 1 SCALED TO .9                                
CX20652 CALL DOOUTX( 652 )                                              
CY20652 CALL DOOUTY( 652 )                                              
CX10653 CALL DOINX( 653 )                                               
CY10653 CALL DOINY( 653 )                                               
      DO 62 I=1, 96                                                     
      IF( (Q1(I).GT.U1(I))  .OR. (Q1(I).GT.V1(I)) ) GO TO 62            
      RTG(I,1)=RTG(I,1)+DLAM(I)*(V1(I)-Q1(I))*.5                        
62    CONTINUE                                                          
C*** HEATING OF LOWEST LAYER FROM UNDERLYING WATER.T1*SL1KAP=T SURF     
CX20654 CALL DOOUTX( 654 )                                              
CY20654 CALL DOOUTY( 654 )                                              
CX10655 CALL DOINX( 655 )                                               
CY10655 CALL DOINY( 655 )                                               
      DO 70 I=1, 96                                                     
      TDIFS1=TSEA(I)-T1(I)*SL1KAP                                       
      IF(TDIFS1.LT.0.) GO TO 70                                         
      TAU(I,1)=TAU(I,1)+DLAM(I)*TDIFS1*.5                               
70    CONTINUE                                                          
CX20656 CALL DOOUTX( 656 )                                              
CY20656 CALL DOOUTY( 656 )                                              
      IF (ITERS.EQ.MSTP) CALL VERIFY(TAU,ITERS,SINLAT)                  
C*** RETURN TO COEFS.                                                   
      CALL FFA99(DPHI,DLAM,CG,1)                                        
      CALL FFA99(A,ZG,CG, 12 )                                          
      CALL FFA99(B,DG,CG, 12 )                                          
      CALL FFA99(F,UG,CG, 12 )                                          
      CALL FFA99(G,VG,CG, 12 )                                          
      CALL FFA99(TAU,TG,CG, 12 )                                        
      CALL FFA99(RTG,RQG,CG, 7 )                                        
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX40030 CALL SBOUTX (30)                                                
CY40030 CALL SBOUTY (30)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GLATS(KHALF,COLRAD,WGT,WGTCS,RCS2)                     
      DIMENSION COLRAD( 38 ),WGT( 38 ),WGTCS( 38 ),RCS2( 38 )           
      DOUBLE PRECISION COLRAD, EPS, SI, RK2, SCALE, PI, DRADZ, RAD,     
     1 DRAD, P1, P2, PHI, W, X, SN, RC                                  
CX30031 CALL SBINX (31)                                                 
CY30031 CALL SBINY (31)                                                 
      EPS=1.D-12                                                        
C     PRINT 101                                                         
101   FORMAT ('0 I   COLAT   COLRAD     WGT', 12X, 'WGTCS',             
     1 10X, 'ITER  RES')                                                
      SI=1.D+00                                                         
      K2=2*KHALF                                                        
      RK2=K2                                                            
      SCALE=2.D+00/(RK2**2)                                             
      K1=K2-1                                                           
      PI=ATAN(SI)*4.D+00                                                
      DRADZ=PI/360.D+00                                                 
      RAD=0.D+00                                                        
CX10657 CALL DOINX( 657 )                                               
CY10657 CALL DOINY( 657 )                                               
      DO 1000 K=1,KHALF                                                 
      ITER=0                                                            
      DRAD=DRADZ                                                        
1     CALL POLY(K2,RAD,P2)                                              
2     P1 =P2                                                            
      ITER=ITER+1                                                       
      RAD=RAD+DRAD                                                      
      CALL POLY(K2,RAD,P2)                                              
      IF(SIGN(SI,P1).EQ.SIGN(SI,P2))GO TO 2                             
      IF(DRAD.LT.EPS)GO TO 3                                            
      RAD=RAD-DRAD                                                      
      DRAD=DRAD*0.25D+00                                                
      GO TO 1                                                           
3     CONTINUE                                                          
      COLRAD(K)=RAD                                                     
      PHI=RAD   *180.D+00/PI                                            
      CALL POLY(K1,RAD,P1)                                              
      X=COS(RAD)                                                        
      W=SCALE*(1.D+00-X*X)/(P1*P1)                                      
      WGT(K)=REAL(W)                                                    
      SN=SIN(RAD)                                                       
      W=W/(SN*SN)                                                       
      WGTCS(K)=REAL(W)                                                  
      RC=1./(SN*SN)                                                     
      RCS2(K)=REAL(RC)                                                  
      CALL POLY(K2,RAD,P1)                                              
102   FORMAT(1H ,I2,2X,F6.2,2X,F10.7,2X,E13.7,2X,E13.7,2X,I4,2X,D13.7)  
1000  CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20658 CALL DOOUTX( 658 )                                              
CY20658 CALL DOOUTY( 658 )                                              
CX40032 CALL SBOUTX (32)                                                
CY40032 CALL SBOUTY (32)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GLOOP(TPS)                                             
      COMMON IDATE(4),                                                  
     * GESHEM( 96 , 76 ),TSEA( 96 , 76 ),GRGRID( 96 , 76 )              
      COMMON DPDPHI( 992 ),DPDLAM( 961 ),RELVOR( 12 ),ABSVOR( 12 ),     
     1  Q( 961 ), TE( 961 , 12 ), DI( 961 , 12 ), ZE( 961 , 12 ),       
     2  Z( 961 ),  Y( 961 , 12 ),  X( 961 , 12 ),  W( 961 , 12 ),       
     3 QM( 961 ),TEM( 961 , 12 ),DIM( 961 , 12 ),ZEM( 961 , 12 ),       
     5 GZ( 961 ), RQ( 961 , 7  ), RM( 961 , 7  ), RT( 961 , 7  )        
      COMPLEX GZ,Q,TE,DI,ZE,Z,Y,X,W,QM,TEM,DIM,ZEM,RELVOR,              
     1DPDLAM,DPDPHI,ABSVOR,RQ,RM,RT                                     
      COMMON                                                            
     * PS( 96 ), DPHI( 96 ), DLAM( 96 ),                                
     * QF( 96 ),DPHIF( 96 ),DLAMF( 96 ),                                
     *  EF( 96 , 12 ), EEF( 96 , 12 ), CG( 96 , 12 ),                   
     *  TF( 96 , 12 ), TEF( 96 , 12 ),TAU( 96 , 12 ),                   
     *  AF( 96 , 12 ), ZEF( 96 , 12 ),  A( 96 , 12 ),                   
     *  BF( 96 , 12 ), DIF( 96 , 12 ),  B( 96 , 12 ),                   
     *  FF( 96 , 12 ),  UF( 96 , 12 ),  F( 96 , 12 ),                   
     *  GF( 96 , 12 ),  VF( 96 , 12 ),  G( 96 , 12 ),                   
     * UQF( 96 , 7 ), UQFF( 96 , 7 ),                                   
     * VQF( 96 , 7 ), VQFF( 96 , 7 ),                                   
     * RQF( 96 , 7 ), RQFF( 96 , 7 ), RTG( 96 , 7 )                     
      COMMON                                                            
     * U1( 96 ),V1( 96 ),T1( 96 ),Q1( 96 ),                             
     * U2( 96 ),V2( 96 ),T2( 96 ),Q2( 96 )                              
      COMMON/PLNCOM/PLN( 992 ),DER( 961 ),PLNWCS( 961 ),EPS( 992 ),     
     1 COLRAD( 38 ),WGT( 38 ),WGTCS( 38 ),RCS2( 38 )                    
      DOUBLE PRECISION COLRAD,EPS                                       
      COMMON/VERCOM/AM( 12 , 12 ),HM( 12 , 12 ),TM( 12 , 12 ),          
     O              BM( 12 , 12 ),CM( 12 , 12 ),EKIN( 12 ),             
     1 SI( 13 ),SL( 12 ),DEL( 12 ),RDEL2( 12 ),RMSDOT( 11 ),            
     2 CI( 13 ),CL( 12 ),TOV( 12 ),   SV( 12 ),   RPI( 11 ),            
     3 P1( 12 ),P2( 12 ), H1( 12 ),   H2( 12 ),RPIREC( 11 ),            
     4 ROTSIN,ROTCOS,SEADRY,SL1KAP,C1,C2,RLRV,SL100K                    
C ----------------------MEMBER GCOM---------------------------------    
CX30033 CALL SBINX (33)                                                 
CY30033 CALL SBINY (33)                                                 
      CALL STIMER                                                       
      CALL XSTORE(EKIN(1),0., 12 )                                      
      CALL XSTORE(RMSDOT(1),0., 11 )                                    
      CALL DELLNP (Q, DPDPHI, DPDLAM, EPS)                              
CX10659 CALL DOINX( 659 )                                               
CY10659 CALL DOINY( 659 )                                               
C     REMOVE MEAN FROM TEMP.                                            
      DO 2 LEV=1, 12                                                    
      RELVOR(LEV) = ZE(2,LEV)                                           
      ABSVOR(LEV) = RELVOR(LEV) + CMPLX(1.19045E-4,0.)                  
      TE(1,LEV)=TE(1,LEV)-CMPLX(TOV(LEV)*1.4142135,0.)                  
2     CONTINUE                                                          
CX20660 CALL DOOUTX( 660 )                                              
CY20660 CALL DOOUTY( 660 )                                              
      CALL XSTORE( Y(1,1),0.,2* 961 * 12 )                              
      CALL XSTORE( W(1,1),0.,2* 961 * 12 )                              
      CALL XSTORE(Z(1),0.,2* 961 )                                      
      CALL XSTORE(RT(1,1),0.,2* 961 * 7 )                               
CX10661 CALL DOINX( 661 )                                               
CY10661 CALL DOINY( 661 )                                               
      DO 5  LEV=1, 12                                                   
      CALL XMOVEX(X(1,LEV),GZ(1),8* 961 )                               
5     CONTINUE                                                          
CX20662 CALL DOOUTX( 662 )                                              
CY20662 CALL DOOUTY( 662 )                                              
CX10663 CALL DOINX( 663 )                                               
CY10663 CALL DOINY( 663 )                                               
C     X NOW CONTAINS CONTRIB. FROM TOPOG.                               
C************************************************                       
C*****        LAT LOOP                                                  
C************************************************                       
      DO 1000 LAT = 1, 38                                               
      CALL PLN2(PLN,COLRAD,LAT,EPS)                                     
      CALL SUMPLS(Q,PS,QF,PLN)                                          
      CALL SUMPLS(DPDLAM,DLAM,DLAMF,PLN)                                
      CALL SUMPLV(DPDPHI,DPHI,DPHIF,PLN)                                
            CALL LEGUV(PLN,PLNWCS,DER,EPS)                              
      DO 9 K=1, 12                                                      
      ZE(2,K) = RELVOR(K)                                               
      CALL UVGLOB(PLNWCS,DER,DI(1,K),ZE(1,K),FF(1,K),UF(1,K),           
     1                                       GF(1,K),VF(1,K))           
      ZE(2,K) = ABSVOR(K)                                               
9     CONTINUE                                                          
      CALL UVGLOB(PLNWCS,DER,DIM,ZEM,U1,U2,V1,V2)                       
      DO 10 K=1, 12                                                     
      CALL SUMPLS (DI(1,K),BF(1,K),DIF(1,K),PLN)                        
      CALL SUMPLS (TE(1,K),TF(1,K),TEF(1,K),PLN)                        
      CALL SUMPLS (ZE(1,K),AF(1,K),ZEF(1,K),PLN)                        
10    CONTINUE                                                          
      DO 8 K=1, 7                                                       
      CALL SUMPLS(RQ(1,K),RQF(1,K),RQFF(1,K),PLN)                       
8     CONTINUE                                                          
      CALL SUMPLS(TEM,T1,T2,PLN)                                        
      CALL SUMPLS( RM,Q1,Q2,PLN)                                        
      SINLAT=REAL(COS(COLRAD(LAT)))                                     
      LATCO=LAT                                                         
      CALL GFIDI(BF,TF,FF,GF,AF,RQF,PS,DPHI,DLAM,CG,                    
     *A,B,F,G,EF,RTG,TAU,UQF,VQF,RCS2(LAT),SINLAT,                      
     * GRGRID(1, LATCO),TSEA(1,LATCO),U1,V1,T1,Q1)                      
      SINLAT=-SINLAT                                                    
      LATCO= 76 +1-LAT                                                  
      CALL GFIDI(DIF,TEF,UF,VF,ZEF,RQFF,QF,DPHIF,DLAMF,CG,              
     *A,B,F,G,EEF,RTG,TAU,UQFF,VQFF,RCS2(LAT),SINLAT,                   
     * GRGRID(1, LATCO),TSEA(1,LATCO),U2,V2,T2,Q2)                      
C                                                                       
      DO 88 K=1, 7                                                      
      CALL SYMASY(UQF(1,K),UQFF(1,K))                                   
      CALL SYMASY(VQF(1,K),VQFF(1,K))                                   
      CALL SYMASY(RQF(1,K),RQFF(1,K))                                   
88    CONTINUE                                                          
      CALL SYMASY(DLAM,DLAMF)                                           
      DO 15 K=1, 12                                                     
      CALL SYMASY (EF(1,K),EEF(1,K))                                    
      CALL SYMASY (TF(1,K),TEF(1,K))                                    
      CALL SYMASY (AF(1,K),ZEF(1,K))                                    
      CALL SYMASY (BF(1,K),DIF(1,K))                                    
      CALL SYMASY (FF(1,K), UF(1,K))                                    
      CALL SYMASY (GF(1,K), VF(1,K))                                    
15    CONTINUE                                                          
      DO 90 I=1, 992                                                    
      PLN(I)=PLN(I)*WGT(LAT)                                            
90    CONTINUE                                                          
C     COMPUTE NON,LIN,CONTRIB. TO V DEL LNP                             
      CALL FL22(DLAM,DLAMF,Z,PLN)                                       
C     COMPUTE NON,LIN, CONTRIB. TO Y FROM TAU                           
      DO 16 K=1, 12                                                     
      CALL FL22(TF(1,K),TEF(1,K),Y(1,K),PLN)                            
16    CONTINUE                                                          
      DO 99 K=1, 7                                                      
      CALL FL22(RQF(1,K),RQFF(1,K),RT(1,K),PLN)                         
99    CONTINUE                                                          
            CALL GOZRIM (PLN, DER, EPS, LAT, PLNWCS, RCS2)              
      DO 109 K=1, 7                                                     
      CALL MSU22(UQFF(1,K),UQF(1,K),VQFF(1,K),VQF(1,K),RT(1,K),PLNWCS,  
     1 DER )                                                            
109   CONTINUE                                                          
      DO 17 K=1, 12                                                     
      CALL PSU22(ZEF(1,K),AF(1,K),DIF(1,K),BF(1,K),X(1,K),PLNWCS,DER)   
      CALL MSU22(ZEF(1,K),AF(1,K),DIF(1,K),BF(1,K),W(1,K),PLNWCS,DER)   
      CALL MSU22( UF(1,K),FF(1,K), VF(1,K),GF(1,K),Y(1,K),PLNWCS,DER)   
C***  FL22 ACTS AS LAPLACIAN ACCUMULATOR . SEE PLN FROM GOZRIM          
      CALL FL22(EF(1,K),EEF(1,K),X(1,K),PLN)                            
17    CONTINUE                                                          
1000   CONTINUE                                                         
CX20664 CALL DOOUTX( 664 )                                              
CY20664 CALL DOOUTY( 664 )                                              
CX10665 CALL DOINX( 665 )                                               
CY10665 CALL DOINY( 665 )                                               
C     RESTORE TEMP. REL. VORTICITY                                      
      DO 18  LEV=1, 12                                                  
      TE(1,LEV)=TE(1,LEV)+CMPLX(TOV(LEV)*1.4142135,0.)                  
      ZE(2,LEV) = RELVOR(LEV)                                           
18    CONTINUE                                                          
CX20666 CALL DOOUTX( 666 )                                              
CY20666 CALL DOOUTY( 666 )                                              
      CALL TTIMER(TPS)                                                  
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX40034 CALL SBOUTX (34)                                                
CY40034 CALL SBOUTY (34)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GNMINI(DK,TK,SATCRI,DISPK,PERCUT,MODS,NITER,NF)        
      COMMON IDATE(4),                                                  
     * GESHEM( 96 , 76 ),TSEA( 96 , 76 ),GRGRID( 96 , 76 )              
      COMMON DPDPHI( 992 ),DPDLAM( 961 ),RELVOR( 12 ),ABSVOR( 12 ),     
     1  Q( 961 ), TE( 961 , 12 ), DI( 961 , 12 ), ZE( 961 , 12 ),       
     2  Z( 961 ),  Y( 961 , 12 ),  X( 961 , 12 ),  W( 961 , 12 ),       
     3 QM( 961 ),TEM( 961 , 12 ),DIM( 961 , 12 ),ZEM( 961 , 12 ),       
     5 GZ( 961 ), RQ( 961 , 7  ), RM( 961 , 7  ), RT( 961 , 7  )        
      COMPLEX GZ,Q,TE,DI,ZE,Z,Y,X,W,QM,TEM,DIM,ZEM,RELVOR,              
     1DPDLAM,DPDPHI,ABSVOR,RQ,RM,RT                                     
      COMMON                                                            
     * PS( 96 ), DPHI( 96 ), DLAM( 96 ),                                
     * QF( 96 ),DPHIF( 96 ),DLAMF( 96 ),                                
     *  EF( 96 , 12 ), EEF( 96 , 12 ), CG( 96 , 12 ),                   
     *  TF( 96 , 12 ), TEF( 96 , 12 ),TAU( 96 , 12 ),                   
     *  AF( 96 , 12 ), ZEF( 96 , 12 ),  A( 96 , 12 ),                   
     *  BF( 96 , 12 ), DIF( 96 , 12 ),  B( 96 , 12 ),                   
     *  FF( 96 , 12 ),  UF( 96 , 12 ),  F( 96 , 12 ),                   
     *  GF( 96 , 12 ),  VF( 96 , 12 ),  G( 96 , 12 ),                   
     * UQF( 96 , 7 ), UQFF( 96 , 7 ),                                   
     * VQF( 96 , 7 ), VQFF( 96 , 7 ),                                   
     * RQF( 96 , 7 ), RQFF( 96 , 7 ), RTG( 96 , 7 )                     
      COMMON                                                            
     * U1( 96 ),V1( 96 ),T1( 96 ),Q1( 96 ),                             
     * U2( 96 ),V2( 96 ),T2( 96 ),Q2( 96 )                              
      COMMON/PLNCOM/PLN( 992 ),DER( 961 ),PLNWCS( 961 ),EPS( 992 ),     
     1 COLRAD( 38 ),WGT( 38 ),WGTCS( 38 ),RCS2( 38 )                    
      DOUBLE PRECISION COLRAD,EPS                                       
      COMMON/VERCOM/AM( 12 , 12 ),HM( 12 , 12 ),TM( 12 , 12 ),          
     O              BM( 12 , 12 ),CM( 12 , 12 ),EKIN( 12 ),             
     1 SI( 13 ),SL( 12 ),DEL( 12 ),RDEL2( 12 ),RMSDOT( 11 ),            
     2 CI( 13 ),CL( 12 ),TOV( 12 ),   SV( 12 ),   RPI( 11 ),            
     3 P1( 12 ),P2( 12 ), H1( 12 ),   H2( 12 ),RPIREC( 11 ),            
     4 ROTSIN,ROTCOS,SEADRY,SL1KAP,C1,C2,RLRV,SL100K                    
C ----------------------MEMBER GCOM---------------------------------    
      DIMENSION EIGG( 12 , 12 ),EIGGT( 12 , 12 ),                       
     1DOTPRO( 12 ),GH( 12 ),VERIN( 12 ),TO( 12 )                        
      COMMON/NMODF/NFILES                                               
C**** IF NFILES.EQ.1 , ONE LARGE DATASET FOR NORMAL MODES               
      DATA DT/1./,TPS/0./                                               
CX30035 CALL SBINX (35)                                                 
CY30035 CALL SBINY (35)                                                 
C     PRINT 100                                                         
      WRITE(6,100)                                                      
100   FORMAT(1H0,'INPUT VALUES')                                        
      CALL RMS(QM,DIM,ZEM,TEM)                                          
      REWIND NF                                                         
CX10667 CALL DOINX( 667 )                                               
CY10667 CALL DOINY( 667 )                                               
      DO 2000 ITER=1,NITER                                              
C     READ(NF)EIGG,EIGGT,GH,DOTPRO,TO                                   
      READ(NF,8000) EIGG,EIGGT                                          
      READ(NF,8000) GH,DOTPRO,TO                                        
8000  FORMAT(6E13.7)                                                    
      DO 10 K=1, 12                                                     
      VERIN(K)=0.                                                       
      DO 10 J=1, 12                                                     
10    VERIN(K)=VERIN(K)+SV(J)*EIGG(J,K)                                 
      IF(NFILES.NE.1)REWIND NF                                          
      CALL BMCM(TOV,P1,P2,H1,H2,DEL,CI,BM,CM,DT,SV,AM)                  
      CALL GLOOP(TPS)                                                   
      CALL TEND(Z,X,Y,QM,DIM,TEM,AM,TOV,BM,DEL)                         
C     PRINT 101,ITER                                                    
      WRITE(6,101)ITER                                                  
101   FORMAT(1H0,'TENDENCIES AT ITERATION ',I2)                         
      CALL RMS(Z,X,W,Y)                                                 
C***  CREATE COMPOSITE VARIABLE=(PHI+R*TO*Q)(DOT)                       
      DO 20 LEV=1, 12                                                   
      DO 20 I=1, 961                                                    
      ZE(I,LEV)=CMPLX( 287.05 * TO(LEV)* REAL(Z(I)),                    
     1                 287.05 * TO(LEV)*AIMAG(Z(I)))                    
      DO 15 K=1, 12                                                     
      ZE(I,LEV)=ZE(I,LEV)+CMPLX(HM(LEV,K)* REAL(Y(I,K)),                
     1                          HM(LEV,K)*AIMAG(Y(I,K)))                
15    CONTINUE                                                          
20    CONTINUE                                                          
      CALL VERTIC(X ,EIGG,EIGGT,DOTPRO,MODS,-1)                         
      CALL VERTIC(W ,EIGG,EIGGT,DOTPRO,MODS,-1)                         
      CALL VERTIC(ZE,EIGG,EIGGT,DOTPRO,MODS,-1)                         
      CALL PRIMES(W,X,ZE,GH,MODS,-1)                                    
C********                                                               
      CALL HORIZ1(W,X,ZE,EEF,ZEF,PS,RT,RTG,DK,TK,PERCUT,NF,MODS)        
C********                                                               
      CALL PRIMES(W,X,ZE,GH,MODS,+1)                                    
C***  COMPUTE DELTA(Q) FROM COMPOSITE VARIABLE                          
      DO 30 I=1, 961                                                    
      Z(I)=(0.,0.)                                                      
      DO 30 MOD=1,MODS                                                  
      Z(I)=Z(I)+CMPLX(                                                  
     1 -VERIN(MOD)* REAL(ZE(I,MOD))/GH(MOD),                            
     2 -VERIN(MOD)*AIMAG(ZE(I,MOD))/GH(MOD))                            
30    CONTINUE                                                          
      CALL VERTIC(X ,EIGG,EIGGT,DOTPRO,MODS,+1)                         
      CALL VERTIC(W ,EIGG,EIGGT,DOTPRO,MODS,+1)                         
      CALL VERTIC(ZE,EIGG,EIGGT,DOTPRO,MODS,+1)                         
C***  COMPUTE DELTA(PHI) FROM COMPOSITE VARIABLE                        
      DO 40 LEV=1, 12                                                   
      DO 40 I=1, 961                                                    
      ZE(I,LEV)=ZE(I,LEV)-CMPLX( 287.05 *TO(LEV)* REAL(Z(I)),           
     1                           287.05 *TO(LEV)*AIMAG(Z(I)))           
40    CONTINUE                                                          
C***  COMPUTE DELTA(T) FROM PHI                                         
      DO 50 LEV=1, 12                                                   
      DO 50 I=1, 961                                                    
      Y(I,LEV)=(0.,0.)                                                  
      DO 45 K=1, 12                                                     
      Y(I,LEV)=Y(I,LEV)+CMPLX(TM(LEV,K)* REAL(ZE(I,K)),                 
     1                        TM(LEV,K)*AIMAG(ZE(I,K)))                 
45    CONTINUE                                                          
50    CONTINUE                                                          
      DO 60 I=1, 961                                                    
      Q(I)=QM(I)-Z(I)                                                   
      QM(I)=Q(I)                                                        
60    CONTINUE                                                          
      DO 70 K=1, 12                                                     
      DO 70 I=1, 961                                                    
      TE(I,K)=TEM(I,K)-Y(I,K)                                           
      TEM(I,K)=TE(I,K)                                                  
      ZE(I,K)=ZEM(I,K)-W(I,K)                                           
      ZEM(I,K)=ZE(I,K)                                                  
      DI(I,K)=DIM(I,K)-X(I,K)                                           
      DIM(I,K)=DI(I,K)                                                  
70    CONTINUE                                                          
C     PRINT 105,ITER                                                    
      WRITE(6,105)ITER                                                  
105   FORMAT(1H0,'CHANGES AT ITERATION ',I2)                            
      CALL RMS(Z,X,W,Y)                                                 
C     PRINT 106,ITER                                                    
      WRITE(6,106)ITER                                                  
106   FORMAT(1H0,'INITIALIZED VALUES AT ITERATION ',I2)                 
      CALL RMS( Q,DI,ZE,TE)                                             
2000  CONTINUE                                                          
CX20668 CALL DOOUTX( 668 )                                              
CY20668 CALL DOOUTY( 668 )                                              
CX40036 CALL SBOUTX (36)                                                
CY40036 CALL SBOUTY (36)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GOZRIM (PLN, DER, EPS, LAT, PLNWCS, RCS2)              
      DIMENSION PLN( 992 ), DER( 961 ), EPS( 992 ),                     
     1  PLNWCS( 961 ),  RCS2( 38 )                                      
      DOUBLE PRECISION EPS                                              
CX30037 CALL SBINX (37)                                                 
CY30037 CALL SBINY (37)                                                 
CX10669 CALL DOINX( 669 )                                               
CY10669 CALL DOINY( 669 )                                               
C     COMPUTE PLN DERIVATIVES                                           
      DO 3 LL=1, 31                                                     
      L = LL - 1                                                        
      JL = L *  31                                                      
      JLE = L *  32                                                     
      N=L                                                               
      DER(JL+1)=-N*EPS(JLE+2)*PLN(JLE+2)                                
      DO 2 IND=2, 31                                                    
      N = N + 1                                                         
      DER(JL+IND) = (N+1.) * EPS(JLE+IND) * PLN(JLE+IND-1)              
     1  - N * EPS(JLE+IND+1) * PLN(JLE+IND+1)                           
2     CONTINUE                                                          
3     CONTINUE                                                          
CX20670 CALL DOOUTX( 670 )                                              
CY20670 CALL DOOUTY( 670 )                                              
      WCSA =  RCS2(LAT) /  6370000.                                     
      RAA=1./( 6370000. * 6370000. )                                    
CX10671 CALL DOINX( 671 )                                               
CY10671 CALL DOINY( 671 )                                               
      DO 4 I=1, 961                                                     
      DER(I) = DER(I) * WCSA                                            
4     CONTINUE                                                          
CX20672 CALL DOOUTX( 672 )                                              
CY20672 CALL DOOUTY( 672 )                                              
CX10673 CALL DOINX( 673 )                                               
CY10673 CALL DOINY( 673 )                                               
      DO 6 LL=1, 31                                                     
      L=LL-1                                                            
      N=L-1                                                             
      JL = (LL-1) *  31                                                 
      JLE = (LL-1) *  32                                                
      DO 5 IND=1, 31                                                    
      N=N+1                                                             
      PLNWCS(JL+IND) = PLN(JLE+IND) * WCSA * L                          
      PLN(JLE+IND)=(N*N+N)*PLN(JLE+IND)*RAA                             
5     CONTINUE                                                          
6     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20674 CALL DOOUTX( 674 )                                              
CY20674 CALL DOOUTY( 674 )                                              
CX40038 CALL SBOUTX (38)                                                
CY40038 CALL SBOUTY (38)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINEGRDLNF(N,FHOUR,IDATE,DPHI,DLAM,GZ,Q,TE,DI,ZE,RQ,SL,SI)  
      COMPLEX                                                           
     1 GZ( 961 ),                                                       
     2 Q( 961 ),                                                        
     3 TE( 961   , 12 ),DI( 961   , 12 ),ZE( 961   , 12 ),              
     4 RQ( 961   , 7 )                                                  
      DIMENSION IDATE(4),SI( 13 ),SL( 12 )                              
      DIMENSION DPHI( 13 ),DLAM( 12 )                                   
      COMMON /LAB8/LAB(8),KDATE(4),IFINAL                               
CX30039 CALL SBINX (39)                                                 
CY30039 CALL SBINY (39)                                                 
C     SPECTRAL DATA FILE FORMAT                                         
C     HOUR,IDATE(4),SI( 13 ),SL( 12 )                                   
C***  COMPLEX ZLN Q TE DI ZE                                            
C***  IF DATA IS GENERATED BY PREPROECESSOR, USE UNFORMAT READING       
C     READ(N)LAB                                                        
      READ(N,6000) LAB                                                  
      CALL W3FS03(LAB(3),KDATE(1),KDATE(4),KDATE(2),KDATE(3))           
C     PRINT 3000,KDATE,LAB                                              
      WRITE(6,3000)KDATE,LAB                                            
3000  FORMAT(1H0,'GRDLNF KDATE ',4I4,2X,2A4,2X,3A8,2X,3A4)              
C     READ(N)FHOUR,IDATE,(DPHI(K),K=1, 13 ),(DLAM(K),K=1, 12 )          
      READ(N,7000) FHOUR,IDATE,DPHI,DLAM                                
6000  FORMAT(1X,2A4,2X,3A8,2X,3A4)                                      
7000  FORMAT(F5.1,4(I4,1X)/(5E13.7))                                    
8000  FORMAT(6E13.7)                                                    
C     READ (N) GZ                                                       
C     READ (N) Q                                                        
      READ(N,8000) GZ                                                   
      READ(N,8000) Q                                                    
CX10675 CALL DOINX( 675 )                                               
CY10675 CALL DOINY( 675 )                                               
      DO 400 K=1, 12                                                    
C     READ(N) (TE(I,K),I=1, 961   )                                     
      READ(N,8000)(TE(I,K),I=1,961 )                                    
400   CONTINUE                                                          
CX20676 CALL DOOUTX( 676 )                                              
CY20676 CALL DOOUTY( 676 )                                              
CX10677 CALL DOINX( 677 )                                               
CY10677 CALL DOINY( 677 )                                               
      DO 401 K=1, 12                                                    
C     READ(N)(DI(I,K),I=1, 961   )                                      
C     READ(N)(ZE(I,K),I=1, 961   )                                      
      READ(N,8000)(DI(I,K),I=1,961 )                                    
      READ(N,8000)(ZE(I,K),I=1,961 )                                    
401   CONTINUE                                                          
CX20678 CALL DOOUTX( 678 )                                              
CY20678 CALL DOOUTY( 678 )                                              
       XRH =0.9                                                         
      K=1                                                               
C     READ (N) (RQ(I,K),I=1, 961   )                                    
      READ(N,8000)(RQ(I,K),I=1,961 )                                    
CX10679 CALL DOINX( 679 )                                               
CY10679 CALL DOINY( 679 )                                               
      DO 190 I=1, 961                                                   
      RQ(I,K)=CMPLX(XRH*REAL(RQ(I,K)),XRH*AIMAG(RQ(I,K)))               
190   CONTINUE                                                          
CX20680 CALL DOOUTX( 680 )                                              
CY20680 CALL DOOUTY( 680 )                                              
CX10681 CALL DOINX( 681 )                                               
CY10681 CALL DOINY( 681 )                                               
      DO 200 K=2, 7                                                     
C     READ (N) (RQ(I,K),I=1, 961   )                                    
      READ(N,8000)(RQ(I,K),I=1,961 )                                    
200   CONTINUE                                                          
CX20682 CALL DOOUTX( 682 )                                              
CY20682 CALL DOOUTY( 682 )                                              
CX10683 CALL DOINX( 683 )                                               
CY10683 CALL DOINY( 683 )                                               
      DO 3 K=1, 12                                                      
3     DLAM(K)=DLAM(K)-SL(K)                                             
CX20684 CALL DOOUTX( 684 )                                              
CY20684 CALL DOOUTY( 684 )                                              
C     PRINT 100,(DLAM(K),K=1, 12 )                                      
      WRITE(6,100)(DLAM(K),K=1, 12 )                                    
CX10685 CALL DOINX( 685 )                                               
CY10685 CALL DOINY( 685 )                                               
      DO 4 K=1, 13                                                      
4     DPHI(K)=DPHI(K)-SI(K)                                             
CX20686 CALL DOOUTX( 686 )                                              
CY20686 CALL DOOUTY( 686 )                                              
C     PRINT 100,(DPHI(K),K=1, 13 )                                      
      WRITE(6,100)(DPHI(K),K=1, 13 )                                    
100   FORMAT(1H0, 13 (E9.3))                                            
C     PRINT 101,N,FHOUR,IDATE                                           
      WRITE(6,101)N,FHOUR,IDATE                                         
101   FORMAT (1H0, 'IF ABOVE TWO ROWS NOT ZERO, ',                      
     1 'INCONSISTENCY IN SIGMA DEFINITION ON N=',I2,2X,F6.1,2X,4(I4))   
CX40040 CALL SBOUTX (40)                                                
CY40040 CALL SBOUTY (40)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GSSTCD(Z00,SATCRI,ISST)                                
      COMMON IDATE(4),                                                  
     * GESHEM( 96 , 76 ),TSEA( 96 , 76 ),GRGRID( 96 , 76 )              
      COMMON DPDPHI( 992 ),DPDLAM( 961 ),RELVOR( 12 ),ABSVOR( 12 ),     
     1  Q( 961 ), TE( 961 , 12 ), DI( 961 , 12 ), ZE( 961 , 12 ),       
     2  Z( 961 ),  Y( 961 , 12 ),  X( 961 , 12 ),  W( 961 , 12 ),       
     3 QM( 961 ),TEM( 961 , 12 ),DIM( 961 , 12 ),ZEM( 961 , 12 ),       
     5 GZ( 961 ), RQ( 961 , 7  ), RM( 961 , 7  ), RT( 961 , 7  )        
      COMPLEX GZ,Q,TE,DI,ZE,Z,Y,X,W,QM,TEM,DIM,ZEM,RELVOR,              
     1DPDLAM,DPDPHI,ABSVOR,RQ,RM,RT                                     
      COMMON                                                            
     * PS( 96 ), DPHI( 96 ), DLAM( 96 ),                                
     * QF( 96 ),DPHIF( 96 ),DLAMF( 96 ),                                
     *  EF( 96 , 12 ), EEF( 96 , 12 ), CG( 96 , 12 ),                   
     *  TF( 96 , 12 ), TEF( 96 , 12 ),TAU( 96 , 12 ),                   
     *  AF( 96 , 12 ), ZEF( 96 , 12 ),  A( 96 , 12 ),                   
     *  BF( 96 , 12 ), DIF( 96 , 12 ),  B( 96 , 12 ),                   
     *  FF( 96 , 12 ),  UF( 96 , 12 ),  F( 96 , 12 ),                   
     *  GF( 96 , 12 ),  VF( 96 , 12 ),  G( 96 , 12 ),                   
     * UQF( 96 , 7 ), UQFF( 96 , 7 ),                                   
     * VQF( 96 , 7 ), VQFF( 96 , 7 ),                                   
     * RQF( 96 , 7 ), RQFF( 96 , 7 ), RTG( 96 , 7 )                     
      COMMON                                                            
     * U1( 96 ),V1( 96 ),T1( 96 ),Q1( 96 ),                             
     * U2( 96 ),V2( 96 ),T2( 96 ),Q2( 96 )                              
      COMMON/PLNCOM/PLN( 992 ),DER( 961 ),PLNWCS( 961 ),EPS( 992 ),     
     1 COLRAD( 38 ),WGT( 38 ),WGTCS( 38 ),RCS2( 38 )                    
      DOUBLE PRECISION COLRAD,EPS                                       
      COMMON/VERCOM/AM( 12 , 12 ),HM( 12 , 12 ),TM( 12 , 12 ),          
     O              BM( 12 , 12 ),CM( 12 , 12 ),EKIN( 12 ),             
     1 SI( 13 ),SL( 12 ),DEL( 12 ),RDEL2( 12 ),RMSDOT( 11 ),            
     2 CI( 13 ),CL( 12 ),TOV( 12 ),   SV( 12 ),   RPI( 11 ),            
     3 P1( 12 ),P2( 12 ), H1( 12 ),   H2( 12 ),RPIREC( 11 ),            
     4 ROTSIN,ROTCOS,SEADRY,SL1KAP,C1,C2,RLRV,SL100K                    
C ----------------------MEMBER GCOM---------------------------------    
CX30041 CALL SBINX (41)                                                 
CY30041 CALL SBINY (41)                                                 
CX10687 CALL DOINX( 687 )                                               
CY10687 CALL DOINY( 687 )                                               
C     RPI(K) = (SL(K+1)/SL(K))**RK  FROM SETSIG  K=1... 11              
      DO 1 K=1, 11                                                      
1     RPIREC(K) = 1./RPI(K)                                             
CX20688 CALL DOOUTX( 688 )                                              
CY20688 CALL DOOUTY( 688 )                                              
CX10689 CALL DOINX( 689 )                                               
CY10689 CALL DOINY( 689 )                                               
      DO 2 K=1, 12                                                      
2     RDEL2(K)=1./(2.*DEL(K))                                           
CX20690 CALL DOOUTX( 690 )                                              
CY20690 CALL DOOUTY( 690 )                                              
      REWIND ISST                                                       
      READ(ISST,777) GRGRID                                             
777   FORMAT(6E13.7)                                                    
      GA2= 9.8 /( 6370000. * 6370000. )                                 
      Z00=REAL(GZ(1))                                                   
CX10691 CALL DOINX( 691 )                                               
CY10691 CALL DOINY( 691 )                                               
      DO 4 LL=1, 31                                                     
      N=LL-2                                                            
      DO 4 I=1, 31                                                      
      N=N+1                                                             
      FACT=N*(N+1.)*GA2                                                 
      GZ((LL-1)* 31 +I)=CMPLX(                                          
     1 FACT* REAL(GZ((LL-1)* 31 +I)),                                   
     1 FACT*AIMAG(GZ((LL-1)* 31 +I)))                                   
4     CONTINUE                                                          
CX20692 CALL DOOUTX( 692 )                                              
CY20692 CALL DOOUTY( 692 )                                              
C     NOW READ TSEA                                                     
55    CONTINUE                                                          
      READ(ISST,777)TSEA                                                
      REWIND ISST                                                       
C-------                                                                
      RLRV=-2.51E+10/(1.61*2.87E+6)                                     
      CONST=6.11*EXP(-RLRV/273.16)/10.                                  
      C1=CONST* 0.622                                                   
      C2=CONST*(-0.378)                                                 
C-------                                                                
      HBOUND= 287.05 *280.*DEL(1)/( 9.8 *SL(1))                         
      SEADRY=(2.*0.07*0.001)/HBOUND                                     
      RK= 287.05 / 1005.                                                
      SL1KAP=1./SL(1)**RK                                               
      SL100K = (SL(1)/100.)**RK                                         
      ROTSIN=SIN(20./57.295779)                                         
      ROTCOS=COS(20./57.295779)                                         
C     PRINT 9877,HBOUND                                                 
      WRITE(6,9877)HBOUND                                               
9877  FORMAT(1H0,' HBOUND=',E12.3)                                      
      SQR2=SQRT(2.)                                                     
CX10693 CALL DOINX( 693 )                                               
CY10693 CALL DOINY( 693 )                                               
      DO 444 LA=1, 76                                                   
      DO 444 LO=1, 96                                                   
444   GRGRID(LO,LA)=GRGRID(LO,LA)*SQR2/HBOUND                           
CX20694 CALL DOOUTX( 694 )                                              
CY20694 CALL DOOUTY( 694 )                                              
102   FORMAT(1H0,'GSSTCD COMPLETED FOR ',4I4)                           
C     PRINT 102, IDATE                                                  
      WRITE(6,102) IDATE                                                
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX40042 CALL SBOUTX (42)                                                
CY40042 CALL SBOUTY (42)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GWATER(FILTA,DT,IFILT,SATCRI,ACUMT,RELKUO,NUN,NKUO)    
      COMMON IDATE(4),                                                  
     * GESHEM( 96 , 76 ),TSEA( 96 , 76 ),GRGRID( 96 , 76 )              
      COMMON DPDPHI( 992 ),DPDLAM( 961 ),RELVOR( 12 ),ABSVOR( 12 ),     
     1  Q( 961 ), TE( 961 , 12 ), DI( 961 , 12 ), ZE( 961 , 12 ),       
     2  Z( 961 ),  Y( 961 , 12 ),  X( 961 , 12 ),  W( 961 , 12 ),       
     3 QM( 961 ),TEM( 961 , 12 ),DIM( 961 , 12 ),ZEM( 961 , 12 ),       
     5 GZ( 961 ), RQ( 961 , 7  ), RM( 961 , 7  ), RT( 961 , 7  )        
      COMPLEX GZ,Q,TE,DI,ZE,Z,Y,X,W,QM,TEM,DIM,ZEM,RELVOR,              
     1DPDLAM,DPDPHI,ABSVOR,RQ,RM,RT                                     
      COMMON                                                            
     * PS( 96 ), DPHI( 96 ), DLAM( 96 ),                                
     * QF( 96 ),DPHIF( 96 ),DLAMF( 96 ),                                
     *  EF( 96 , 12 ), EEF( 96 , 12 ), CG( 96 , 12 ),                   
     *  TF( 96 , 12 ), TEF( 96 , 12 ),TAU( 96 , 12 ),                   
     *  AF( 96 , 12 ), ZEF( 96 , 12 ),  A( 96 , 12 ),                   
     *  BF( 96 , 12 ), DIF( 96 , 12 ),  B( 96 , 12 ),                   
     *  FF( 96 , 12 ),  UF( 96 , 12 ),  F( 96 , 12 ),                   
     *  GF( 96 , 12 ),  VF( 96 , 12 ),  G( 96 , 12 ),                   
     * UQF( 96 , 7 ), UQFF( 96 , 7 ),                                   
     * VQF( 96 , 7 ), VQFF( 96 , 7 ),                                   
     * RQF( 96 , 7 ), RQFF( 96 , 7 ), RTG( 96 , 7 )                     
      COMMON                                                            
     * U1( 96 ),V1( 96 ),T1( 96 ),Q1( 96 ),                             
     * U2( 96 ),V2( 96 ),T2( 96 ),Q2( 96 )                              
      COMMON/PLNCOM/PLN( 992 ),DER( 961 ),PLNWCS( 961 ),EPS( 992 ),     
     1 COLRAD( 38 ),WGT( 38 ),WGTCS( 38 ),RCS2( 38 )                    
      DOUBLE PRECISION COLRAD,EPS                                       
      COMMON/VERCOM/AM( 12 , 12 ),HM( 12 , 12 ),TM( 12 , 12 ),          
     O              BM( 12 , 12 ),CM( 12 , 12 ),EKIN( 12 ),             
     1 SI( 13 ),SL( 12 ),DEL( 12 ),RDEL2( 12 ),RMSDOT( 11 ),            
     2 CI( 13 ),CL( 12 ),TOV( 12 ),   SV( 12 ),   RPI( 11 ),            
     3 P1( 12 ),P2( 12 ), H1( 12 ),   H2( 12 ),RPIREC( 11 ),            
     4 ROTSIN,ROTCOS,SEADRY,SL1KAP,C1,C2,RLRV,SL100K                    
C ----------------------MEMBER GCOM---------------------------------    
CX30043 CALL SBINX (43)                                                 
CY30043 CALL SBINY (43)                                                 
C     RT CONTAINS SPECTRAL Q TENDENCIES FROM DYNAMICS COMPUTED IN       
C     NSLOOP AND WALLFD. SEA SURF EVAPORATION SHOULD BE IN WALLFD       
C***                                                                    
      DT2=DT+DT                                                         
      NUN=0                                                             
      NKUO=0                                                            
      MSTA=0                                                            
CX10695 CALL DOINX( 695 )                                               
CY10695 CALL DOINY( 695 )                                               
C***                                                                    
C     FIRST,STEP Q IN TIME TO OBTAIN PRELIM. FORECAST. PUT IN X         
      DO 49 K=1, 7                                                      
      DO 49 I=1, 961                                                    
      X(I,K)=CMPLX(REAL(RM(I,K))+DT2* REAL(RT(I,K))  ,                  
     1            AIMAG(RM(I,K))+DT2*AIMAG(RT(I,K))  )                  
49    CONTINUE                                                          
CX20696 CALL DOOUTX( 696 )                                              
CY20696 CALL DOOUTY( 696 )                                              
C     RT WILL CONTAIN ADJUESTED Q AT END OF LAT LOOP                    
C     RT WAS USED IN X=ESTIMATED NEW Q                                  
C     ZERO ADJUSTED TEMP. COEFS.                                        
      CALL XSTORE(RT(1,1),0.,2* 961 * 7 )                               
      CALL XSTORE( Y(1,1),0.,2* 961 * 12 )                              
CX10697 CALL DOINX( 697 )                                               
CY10697 CALL DOINY( 697 )                                               
C                                                                       
      DO 1000 LAT=1, 38                                                 
      CALL PLN2 (PLN, COLRAD, LAT, EPS)                                 
      DO 10 K=1, 12                                                     
      CALL SUMPLS (TE(1,K), CG(1,K),  UF(1,K), PLN)                     
10    CONTINUE                                                          
      CALL SUMPLS (Q, PS, QF, PLN)                                      
      DO 11 K=1, 7                                                      
      CALL SUMPLS(RQ(1,K),A(1,K),AF(1,K),PLN)                           
C     A(I,K) CONTAIN Q AT N                                             
      CALL SUMPLS( X(1,K),RTG(1,K), EEF(1,K),PLN )                      
11    CONTINUE                                                          
      LATCO=LAT                                                         
      CALL FFS99( CG,G, 12 )                                            
      CALL FFS99(RTG,G, 7 )                                             
      CALL FFS99(  A,G, 7 )                                             
      CALL FFS99( PS,G,1)                                               
      DO 2 I=1, 96                                                      
2     PS(I)=EXP(PS(I))                                                  
      CALL CONKUO(DT,SATCRI,ACUMT,NKUO,PS,DEL,SL,A,RTG,CG,B,F,          
     1 GESHEM(1,LATCO),RELKUO,MSTA)                                     
C                                                                       
      CALL LRGSCL(SATCRI,GESHEM(1,LATCO),CG,G,RTG,PS,B,DEL,             
     1 DLAM,SL,DLAMF,DPHIF,F)                                           
C                                                                       
      CALL FIXTET(PS,DLAM,DEL,RPI,RPIREC,CG,RTG,TAU,NUN,SL100K)         
      CALL FFA99( CG, TF,G, 12 )                                        
      CALL FFA99(RTG,RQF,G, 7 )                                         
      LATCO= 76 +1-LAT                                                  
      CALL FFS99( UF,G, 12 )                                            
      CALL FFS99(EEF,G, 7 )                                             
      CALL FFS99( AF,G, 7 )                                             
      CALL FFS99( QF,G,1)                                               
      DO 3 I=1, 96                                                      
3     QF(I)=EXP(QF(I))                                                  
      CALL CONKUO(DT,SATCRI,ACUMT,NKUO,QF,DEL,SL,AF,EEF,UF,B,F,         
     1 GESHEM(1,LATCO),RELKUO,MSTA)                                     
C                                                                       
      CALL LRGSCL(SATCRI,GESHEM(1,LATCO),UF,G,EEF,QF,B,DEL,             
     1 DLAM,SL,DLAMF,DPHIF,F)                                           
C                                                                       
            CALL FIXTET(QF,DLAM,DEL,RPI,RPIREC,UF,EEF,TAU,NUN,SL100K)   
C                                                                       
      CALL FFA99( UF, TEF,G, 12 )                                       
      CALL FFA99(EEF,RQFF,G, 7 )                                        
      DO 90 I=1, 992                                                    
      PLN(I)=PLN(I)*WGT(LAT)                                            
90    CONTINUE                                                          
      DO 16 K=1, 12                                                     
      CALL SYMASY (TF(1,K), TEF(1,K))                                   
      CALL FL22 (TF(1,K), TEF(1,K), Y(1,K), PLN)                        
16    CONTINUE                                                          
      DO 17 K=1, 7                                                      
      CALL SYMASY(RQF(1,K),RQFF(1,K))                                   
      CALL   FL22(RQF(1,K),RQFF(1,K),RT(1,K),PLN)                       
17    CONTINUE                                                          
C     RAITOT=RAITOT+WGT(LAT)*RAILAT/(2.* 96 .* 9.8 )                    
1000  CONTINUE                                                          
CX20698 CALL DOOUTX( 698 )                                              
CY20698 CALL DOOUTY( 698 )                                              
CX10699 CALL DOINX( 699 )                                               
CY10699 CALL DOINY( 699 )                                               
C***  MOVE CORRECTED TEMP. INTO TE. COMPLETE TIME FILTER IN MAIN        
      DO 70 K=1, 12                                                     
      DO 70 I=1, 961                                                    
      TE(I,K) = Y(I,K)                                                  
70    CONTINUE                                                          
CX20700 CALL DOOUTX( 700 )                                              
CY20700 CALL DOOUTY( 700 )                                              
C     RT NOW CONTAINS ADJUSTED PREDICTED Q. NEXT TIME FILTER            
      IF(IFILT.EQ.0) GO TO 750                                          
      FILTB=(1.-FILTA)*0.5                                              
CX10701 CALL DOINX( 701 )                                               
CY10701 CALL DOINY( 701 )                                               
      DO 75 K=1, 7                                                      
      DO 75 I=1, 961                                                    
      RM(I,K)=CMPLX(                                                    
     1FILTA* REAL(RQ(I,K))+FILTB*( REAL(RT(I,K))+ REAL(RM(I,K))) ,      
     2FILTA*AIMAG(RQ(I,K))+FILTB*(AIMAG(RT(I,K))+AIMAG(RM(I,K))) )      
75    CONTINUE                                                          
CX20702 CALL DOOUTX( 702 )                                              
CY20702 CALL DOOUTY( 702 )                                              
750   CONTINUE                                                          
CX10703 CALL DOINX( 703 )                                               
CY10703 CALL DOINY( 703 )                                               
      DO 76 K=1, 7                                                      
      DO 76 I=1, 961                                                    
      RQ(I,K)=RT(I,K)                                                   
76    CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20704 CALL DOOUTX( 704 )                                              
CY20704 CALL DOOUTY( 704 )                                              
CX40044 CALL SBOUTX (44)                                                
CY40044 CALL SBOUTY (44)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE GWRITE(N,FHOUR,IDATE,Z,Q,TE,DI,ZE,RQ,SL,SI,GZ,Z00)     
      INTEGER DIGITS(3)                                                 
      COMPLEX                                                           
     1  Z( 961   ),GZ( 961 ),                                           
     2 Q( 961   ),                                                      
     3 TE( 961   , 12 ),DI( 961   , 12 ),ZE( 961   , 12 ),              
     4 RQ( 961   , 7 )                                                  
      DIMENSION IDATE(4),SI( 13 ),SL( 12 ) ,LAB1(2),LAB2(4)             
      COMMON /LAB8/LAB(8),KDATE(4),IFINAL                               
      DATA        LAB1/'SIGI','SIGG'/                                   
      DATA      LAB2/'NIB2','N6B2','E6B2','ESB2'/                       
CX30045 CALL SBINX (45)                                                 
CY30045 CALL SBINY (45)                                                 
      IF(IFINAL)600,600,500                                             
500   IHOUR=FHOUR+0.5                                                   
      IF(IHOUR.GT.100)IHOUR=IHOUR-100                                   
      IF(IHOUR.GT.100)IHOUR=IHOUR-100                                   
      LAB(1)=LAB1(1)                                                    
      IF(IHOUR.EQ.6)LAB(1)=LAB1(2)                                      
      IF(IHOUR.EQ.0.AND.(IDATE(1).EQ.0.OR.IDATE(1).EQ.12))              
     1 LAB(2)=LAB2(1)                                                   
      IF(IHOUR.EQ.0.AND.(IDATE(1).EQ.6.OR.IDATE(1).EQ.18))              
     1 LAB(2)=LAB2(2)                                                   
      IF(IHOUR.EQ.6.AND.(IDATE(1).EQ.0.OR.IDATE(1).EQ.12))              
     1 LAB(2)=LAB2(3)                                                   
      IF(IHOUR.EQ.6.AND.(IDATE(1).EQ.6.OR.IDATE(1).EQ.18))              
     1 LAB(2)=LAB2(4)                                                   
      GO TO 900                                                         
600   IHOUR=FHOUR+0.5                                                   
900   CONTINUE                                                          
      WRITE(N,800)LAB                                                   
C     PRINT 3000,KDATE,LAB                                              
      WRITE(6,3000)KDATE,LAB                                            
3000  FORMAT(1H0,'GWRITE KDATE ',4I4,2X,2A4,2X,3A8,2X,3A4)              
      WRITE(N,801) FHOUR,IDATE,SI,SL                                    
C***  COMPUTE TOPOG. FROM ITS LAPLACIAN                                 
      ASQ= 6370000. * 6370000. / 9.8                                    
      K=0                                                               
CX10705 CALL DOINX( 705 )                                               
CY10705 CALL DOINY( 705 )                                               
      DO 50 I=2, 31                                                     
      K=K+1                                                             
      FACT=ASQ/(K*(K+1.))                                               
      Z(I)=CMPLX(FACT*REAL(GZ(I)),FACT*AIMAG(GZ(I)))                    
50    CONTINUE                                                          
CX20706 CALL DOOUTX( 706 )                                              
CY20706 CALL DOOUTY( 706 )                                              
CX10707 CALL DOINX( 707 )                                               
CY10707 CALL DOINY( 707 )                                               
      DO 51 LL=2, 31                                                    
      K=LL-2                                                            
      DO 51 I=1, 31                                                     
      K=K+1                                                             
      FACT=ASQ/(K*(K+1.))                                               
       Z((LL-1)* 31 +I)=CMPLX(                                          
     1 FACT* REAL(GZ((LL-1)* 31 +I)),                                   
     1 FACT*AIMAG(GZ((LL-1)* 31 +I)))                                   
51    CONTINUE                                                          
CX20708 CALL DOOUTX( 708 )                                              
CY20708 CALL DOOUTY( 708 )                                              
      Z(1)=CMPLX(Z00,0.)                                                
C***                                                                    
      WRITE(N,802)Z                                                     
      WRITE (N,802) Q                                                   
CX10709 CALL DOINX( 709 )                                               
CY10709 CALL DOINY( 709 )                                               
      DO 15 K=1, 12                                                     
      WRITE(N,802) (TE(I,K),I=1, 961   )                                
15    CONTINUE                                                          
CX20710 CALL DOOUTX( 710 )                                              
CY20710 CALL DOOUTY( 710 )                                              
CX10711 CALL DOINX( 711 )                                               
CY10711 CALL DOINY( 711 )                                               
      DO 20 K=1, 12                                                     
      WRITE(N,802)(DI(I,K),I=1, 961   )                                 
      WRITE(N,802)(ZE(I,K),I=1, 961   )                                 
20    CONTINUE                                                          
CX20712 CALL DOOUTX( 712 )                                              
CY20712 CALL DOOUTY( 712 )                                              
       XRH = 0.9                                                        
      K=1                                                               
CX10713 CALL DOINX( 713 )                                               
CY10713 CALL DOINY( 713 )                                               
      DO 23 I=1,  961                                                   
      Z(I)=CMPLX(REAL(RQ(I,K))/XRH,AIMAG(RQ(I,K))/XRH)                  
23    CONTINUE                                                          
CX20714 CALL DOOUTX( 714 )                                              
CY20714 CALL DOOUTY( 714 )                                              
      WRITE(N,802)(Z(I),I=1, 961 )                                      
CX10715 CALL DOINX( 715 )                                               
CY10715 CALL DOINY( 715 )                                               
      DO 25 K=2, 7                                                      
      WRITE (N,802) (RQ(I,K),I=1, 961 )                                 
25    CONTINUE                                                          
CX20716 CALL DOOUTX( 716 )                                              
CY20716 CALL DOOUTY( 716 )                                              
C     PRINT 3001,FHOUR,IDATE,N                                          
      WRITE(6,3001)FHOUR,IDATE,N                                        
800   FORMAT(1X,2A4,2X,3A8,2X,3A4)                                      
801   FORMAT(F5.1,4(I4,1X)/(5E13.7))                                    
802   FORMAT(6E13.7)                                                    
3001   FORMAT(1H0,'GWRITE FHOUR=',F6.2,2X,4I4,2X,'N=',I2)               
CX40046 CALL SBOUTX (46)                                                
CY40046 CALL SBOUTY (46)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE HORIZ1(VORD,DIVD,COMD,SDOT,ADOT,PER,G,WORK,DK,TK,V,    
     1 NF,MODS)                                                         
      COMMON/NMODF/NFILES                                               
      COMPLEX   SDOT( 47 ),ADOT( 46 ),WORK( 32 ),                       
     1VORD( 31 , 31 , 12 ),                                             
     2DIVD( 31 , 31 , 12 ),                                             
     2COMD( 31 , 31 , 12 )                                              
      DIMENSION PER( 32 ),G( 32 , 47 )                                  
      DATA NS/ 47 /,NA/ 46 /,                                           
     1 JSEV/ 16 /,JSOD/ 15 /                                            
CX30047 CALL SBINX (47)                                                 
CY30047 CALL SBINY (47)                                                 
      JEVPOD=JSEV+JSOD                                                  
CX10717 CALL DOINX( 717 )                                               
CY10717 CALL DOINY( 717 )                                               
      DO 2000 MOD=1,MODS                                                
C     PRINT 100,MOD                                                     
      WRITE(6,100)MOD                                                   
100   FORMAT(1H0,'DOING MODE ',I2)                                      
      NFG=NF+MOD                                                        
      IF(NFILES.EQ.1)NFG=NF                                             
      DO 1000 LL=1, 31                                                  
      K=0                                                               
      L=LL-1                                                            
      DO 10 I=2, 30 ,2                                                  
      K=K+1                                                             
      SDOT(K)=VORD(I,LL,MOD)                                            
      ADOT(K+JSEV)=DIVD(I,LL,MOD)                                       
      ADOT(K+JEVPOD) =COMD(I,LL,MOD)                                    
10    CONTINUE                                                          
      K=0                                                               
      DO 20 I=1, 31 ,2                                                  
      K=K+1                                                             
      ADOT(K)=VORD(I,LL,MOD)                                            
      SDOT(K+JSOD)=DIVD(I,LL,MOD)                                       
      SDOT(K+JEVPOD) =COMD(I,LL,MOD)                                    
20    CONTINUE                                                          
      IF(L.EQ.0) GO TO 500                                              
      CALL HORIZ2(SDOT,NS,PER,G, 32 ,WORK,V,L,NFG,DK,TK)                
500   IF(L.GT.0)GO TO 600                                               
      CALL HORIZ2(SDOT,NS,PER,G, 30 ,WORK,V,L,NFG,DK,TK)                
600   CALL HORIZ2(ADOT,NA,PER,G, 30 ,WORK,V,L,NFG,DK,TK)                
      K=0                                                               
      DO 30 I=2, 30 ,2                                                  
      K=K+1                                                             
      VORD(I,LL,MOD)=SDOT(K)                                            
      DIVD(I,LL,MOD)=ADOT(K+JSEV)                                       
      COMD(I,LL,MOD)=ADOT(K+JEVPOD)                                     
30    CONTINUE                                                          
      K=0                                                               
      DO 40 I=1, 31 ,2                                                  
      K=K+1                                                             
      VORD(I,LL,MOD)=ADOT(K)                                            
      DIVD(I,LL,MOD)=SDOT(K+JSOD)                                       
      COMD(I,LL,MOD)=SDOT(K+JEVPOD)                                     
40    CONTINUE                                                          
1000  CONTINUE                                                          
      IF(NFILES.NE.1) REWIND NFG                                        
2000  CONTINUE                                                          
CX20718 CALL DOOUTX( 718 )                                              
CY20718 CALL DOOUTY( 718 )                                              
      REWIND NFG                                                        
CX40048 CALL SBOUTX (48)                                                
CY40048 CALL SBOUTY (48)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE HORIZ2(DOT,NAS,PER,G,JG,Y,PERCUT,L,NFG,DK,TK)          
      COMMON/NMODF/NFILES                                               
      COMPLEX   DOT(NAS),Y(JG)                                          
C     NAS=VECTOR SIZE OF SYM. OR  ASY. TENDENCIES STORED IN DOT         
C     JG=NUNBER OF GRAVITY MODES. JG=JCAP FOR L=0 SYM AND ASY CASES     
C     FOR L.NE.0 JG=JCAP2 FOR SYM,JG=JCAP FOR ASY.                      
C     PER STORES PERIODS OF GRAVITY MODES,G STORES EIGENVECTORS         
      DIMENSION PER(JG),G(JG,NAS)                                       
CX30049 CALL SBINX (49)                                                 
CY30049 CALL SBINY (49)                                                 
C     READ(NFG) PER,G                                                   
      READ(NFG,8000) PER,G                                              
8000  FORMAT(6E13.7)                                                    
      DIF=DK/( 6370000. **4)                                            
      JT=0                                                              
      JC=0                                                              
      JD=0                                                              
      G1=0.                                                             
      G2=0.                                                             
      N=L+ 31                                                           
      I=0                                                               
100   I=I+1                                                             
      PERIOD=ABS(PER(I))                                                
      IF(PERIOD-PERCUT)10,10,200                                        
10    JC=JC+1                                                           
      Y(I)=(0.,0.)                                                      
      N=N-1                                                             
      DIFCUT=DIF*N*N*(N+1.)*(N+1.)  +  1.E-6                            
      DIFCUT=0.5/DIFCUT                                                 
      IF(DIFCUT-PERIOD)100,20,20                                        
20    JD=JD+1                                                           
C***  IF ARRIVED HERE,ADJUST......                                      
      JT=JT+1                                                           
CX10719 CALL DOINX( 719 )                                               
CY10719 CALL DOINY( 719 )                                               
      DO 30 J=1,NAS                                                     
      Y(I)=Y(I)+CMPLX(G(I,J)*REAL(DOT(J)),G(I,J)*AIMAG(DOT(J)))         
30    CONTINUE                                                          
CX20720 CALL DOOUTX( 720 )                                              
CY20720 CALL DOOUTY( 720 )                                              
      G1=G1+REAL(Y(I))*REAL(Y(I)) + AIMAG(Y(I))*AIMAG(Y(I))             
      Y(I)=CMPLX(PER(I)*AIMAG(Y(I)),-PER(I)*REAL(Y(I)))                 
      G2=G2+REAL(Y(I))*REAL(Y(I)) + AIMAG(Y(I))*AIMAG(Y(I))             
C                                                                       
      IF(I-JG)100,300,300                                               
C                                                                       
200   I=I-1                                                             
300   CONTINUE                                                          
CX10721 CALL DOINX( 721 )                                               
CY10721 CALL DOINY( 721 )                                               
      DO 40 J=1,NAS                                                     
   40 CONTINUE                                                          
CX20722 CALL DOOUTX( 722 )                                              
CY20722 CALL DOOUTY( 722 )                                              
      DOT(J)=(0.,0.)                                                    
      IF(I.EQ.0) THEN                                                   
CX40050 CALL SBOUTX (50)                                                
CY40050 CALL SBOUTY (50)                                                
      RETURN                                                            
      ENDIF                                                             
      G1=SQRT(G1)                                                       
      G2=SQRT(G2)                                                       
CX10723 CALL DOINX( 723 )                                               
CY10723 CALL DOINY( 723 )                                               
      DO 50 J=1,NAS                                                     
      DO 50 K=1,I                                                       
      DOT(J)=DOT(J)+CMPLX(G(K,J)*REAL(Y(K)),G(K,J)*AIMAG(Y(K)))         
50    CONTINUE                                                          
CX20724 CALL DOOUTX( 724 )                                              
CY20724 CALL DOOUTY( 724 )                                              
C     IF(NAS.EQ. 47 )PRINT 500,L,JT,G1,G2,JC,JD,I                       
C     IF(NAS.EQ. 46 )PRINT 501,L,JT,G1,G2,JC,JD,I                       
500   FORMAT(1H ,'SYM,L=',I2,2X,'TOT ADJUSTED=',I3,2X,'GDOT=',E12.4     
     1,2X,'DELTA(G)=',E12.4,2X,'JC=',I3,2X,'JD=',I3,2X,'I=',I3)         
501   FORMAT(1H ,'ASY,L=',I2,2X,'TOT ADJUSTED=',I3,2X,'GDOT=',E12.4     
     1,2X,'DELTA(G)=',E12.4,2X,'JC=',I3,2X,'JD=',I3,2X,'I=',I3)         
CX40050 CALL SBOUTX (50)                                                
CY40050 CALL SBOUTY (50)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE IMINV (A,N,D,L,M)                                      
C                                                                       
C     ..................................................................
C                                                                       
C        SUBROUTINE IMINV                                               
C                                                                       
C        PURPOSE                                                        
C           INVERT A MATRIX                                             
C                                                                       
C        USAGE                                                          
C           CALL IMINV (A,N,D,L,M)                                      
C                                                                       
C        DESCRIPTION OF PARAMETERS                                      
C           A - INPUT MATRIX, DESTROYED IN COMPUTATION AND REPLACED BY  
C               RESULTANT INVERSE.                                      
C           N - ORDER OF MATRIX A                                       
C           D - RESULTANT DETERMINANT                                   
C           L - WORK VECTOR OF LENGTH N                                 
C           M - WORK VECTOR OF LENGTH N                                 
C                                                                       
C        REMARKS                                                        
C           MATRIX A MUST BE A GENERAL MATRIX                           
C                                                                       
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  
C           NONE                                                        
C                                                                       
C        METHOD                                                         
C           THE STANDARD GAUSS-JORDAN METHOD IS USED. THE DETERMINANT   
C           IS ALSO CALCULATED. A DETERMINANT OF ZERO INDICATES THAT    
C           THE MATRIX IS SINGULAR.                                     
      DIMENSION A(1),L(1),M(1)                                          
CX30051 CALL SBINX (51)                                                 
CY30051 CALL SBINY (51)                                                 
C                                                                       
C        ...............................................................
C                                                                       
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE  
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION      
C        STATEMENT WHICH FOLLOWS.                                       
C                                                                       
C     DOUBLE PRECISION A, D, BIGA, HOLD                                 
C                                                                       
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS    
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS      
C        ROUTINE.                                                       
C                                                                       
C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO      
C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  ABS IN STATEMENT  
C        10 MUST BE CHANGED TO DABS.                                    
C        THIS IS NOT TRUE!  ABS IS A GENERIC NAME.  L. KIPP, 08-22-89   
C        ...............................................................
C                                                                       
C        SEARCH FOR LARGEST ELEMENT                                     
C                                                                       
      D=1.0                                                             
      NK=-N                                                             
CX10725 CALL DOINX( 725 )                                               
CY10725 CALL DOINY( 725 )                                               
      DO 80 K=1,N                                                       
      NK=NK+N                                                           
      L(K)=K                                                            
      M(K)=K                                                            
      KK=NK+K                                                           
      BIGA=A(KK)                                                        
      DO 20 J=K,N                                                       
      IZ=N*(J-1)                                                        
      DO 20 I=K,N                                                       
      IJ=IZ+I                                                           
C  10 IF (DABS(BIGA)-DABS(A(IJ))) 15,20,20                              
   10 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20                                 
   15 BIGA=A(IJ)                                                        
      L(K)=I                                                            
      M(K)=J                                                            
   20 CONTINUE                                                          
C                                                                       
C        INTERCHANGE ROWS                                               
C                                                                       
      J=L(K)                                                            
      IF(J-K) 35,35,25                                                  
   25 KI=K-N                                                            
      DO 30 I=1,N                                                       
      KI=KI+N                                                           
      HOLD=-A(KI)                                                       
      JI=KI-K+J                                                         
      A(KI)=A(JI)                                                       
   30 A(JI) =HOLD                                                       
C                                                                       
C        INTERCHANGE COLUMNS                                            
C                                                                       
   35 I=M(K)                                                            
      IF(I-K) 45,45,38                                                  
   38 JP=N*(I-1)                                                        
      DO 40 J=1,N                                                       
      JK=NK+J                                                           
      JI=JP+J                                                           
      HOLD=-A(JK)                                                       
      A(JK)=A(JI)                                                       
   40 A(JI) =HOLD                                                       
C                                                                       
C        DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS        
C        CONTAINED IN BIGA)                                             
C                                                                       
   45 IF(BIGA) 48,46,48                                                 
   46 D=0.0                                                             
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
   48 CONTINUE                                                          
      DO      55      I=1,N                                             
      IF(I-K) 50,55,50                                                  
   50 IK=NK+I                                                           
      A(IK)=A(IK)/(-BIGA)                                               
   55 CONTINUE                                                          
C                                                                       
C        REDUCE MATRIX                                                  
C                                                                       
      DO 65 I=1,N                                                       
      IK=NK+I                                                           
      IJ=I-N                                                            
      DO 65 J=1,N                                                       
      IJ=IJ+N                                                           
      IF(I-K) 60,65,60                                                  
   60 IF(J-K) 62,65,62                                                  
   62 KJ=IJ-I+K                                                         
      A(IJ)=A(IK)*A(KJ)+A(IJ)                                           
   65 CONTINUE                                                          
C                                                                       
C        DIVIDE ROW BY PIVOT                                            
C                                                                       
      KJ=K-N                                                            
      DO 75 J=1,N                                                       
      KJ=KJ+N                                                           
      IF(J-K) 70,75,70                                                  
   70 A(KJ)=A(KJ)/BIGA                                                  
   75 CONTINUE                                                          
C                                                                       
C        PRODUCT OF PIVOTS                                              
C                                                                       
      D=D*BIGA                                                          
C                                                                       
C        REPLACE PIVOT BY RECIPROCAL                                    
C                                                                       
      A(KK)=1.0/BIGA                                                    
   80 CONTINUE                                                          
CX20726 CALL DOOUTX( 726 )                                              
CY20726 CALL DOOUTY( 726 )                                              
C                                                                       
C        FINAL ROW AND COLUMN INTERCHANGE                               
C                                                                       
      K=N                                                               
  100 K=(K-1)                                                           
      IF(K) 150,150,105                                                 
  105 I=L(K)                                                            
      IF(I-K) 120,120,108                                               
  108 JQ=N*(K-1)                                                        
      JR=N*(I-1)                                                        
CX10727 CALL DOINX( 727 )                                               
CY10727 CALL DOINY( 727 )                                               
      DO 110 J=1,N                                                      
      JK=JQ+J                                                           
      HOLD=A(JK)                                                        
      JI=JR+J                                                           
      A(JK)=-A(JI)                                                      
  110 A(JI) =HOLD                                                       
CX20728 CALL DOOUTX( 728 )                                              
CY20728 CALL DOOUTY( 728 )                                              
  120 J=M(K)                                                            
      IF(J-K) 100,100,125                                               
  125 KI=K-N                                                            
CX10729 CALL DOINX( 729 )                                               
CY10729 CALL DOINY( 729 )                                               
      DO 130 I=1,N                                                      
      KI=KI+N                                                           
      HOLD=A(KI)                                                        
      JI=KI-K+J                                                         
      A(KI)=-A(JI)                                                      
  130 A(JI) =HOLD                                                       
CX20730 CALL DOOUTX( 730 )                                              
CY20730 CALL DOOUTY( 730 )                                              
      GO TO 100                                                         
  150 CONTINUE                                                          
CX40052 CALL SBOUTX (52)                                                
CY40052 CALL SBOUTY (52)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE LEGUV(PLN,PLNWCS,DER,EPS)                              
      DIMENSION PLNWCS( 31 , 31 ),                                      
     1             DER( 31 , 31 ),                                      
     2             PLN( 32 , 31 )                                       
      DOUBLE PRECISION EPS( 32 , 31 )                                   
CX30053 CALL SBINX (53)                                                 
CY30053 CALL SBINY (53)                                                 
CX10731 CALL DOINX( 731 )                                               
CY10731 CALL DOINY( 731 )                                               
      DO 3 LL=1, 31                                                     
      N=LL-1                                                            
      DER(1,LL)=- 6370000. *EPS(2,LL)*PLN(2,LL)/LL                      
      PLNWCS(1,LL)= 6370000. *PLN(1,LL)/LL                              
      DO 2 I=2, 31                                                      
      N=N+1                                                             
      DER(I,LL)=EPS(I,LL)*PLN(I-1,LL)/N-EPS(I+1,LL)*PLN(I+1,LL)/(N+1.)  
      DER(I,LL)= 6370000. *DER(I,LL)                                    
      PLNWCS(I,LL)=(LL-1.)*PLN(I,LL)* 6370000. /(N*N+N)                 
2     CONTINUE                                                          
3     CONTINUE                                                          
CX20732 CALL DOOUTX( 732 )                                              
CY20732 CALL DOOUTY( 732 )                                              
CX40054 CALL SBOUTX (54)                                                
CY40054 CALL SBOUTY (54)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE LRGSCL(SATCRI,GESHEM,      TF,QS,QF,PS,SATC,DEL,PREC,  
     1 SL,B,SUPER,DPOVG)                                                
C***  QF=Q(N+1),QS=SAT. Q AT T(N+1)=TF. PS=SURF. PRESS(CB)              
C***                                                                    
      DIMENSION GESHEM( 96 ),TF( 96 , 12 ),                             
     1QF( 96 , 7 ),QS( 96 , 7 ),PS( 96 ),                               
     2SATC( 96 , 7 ),DEL( 12 ),PREC( 96 ), SL( 12 ),                    
     3DPOVG( 96 , 7 ),SUPER( 96 ),B( 96 ),FACT( 7 )                     
CX30055 CALL SBINX (55)                                                 
CY30055 CALL SBINY (55)                                                 
      RLOCP=2.5E+6/ 1005.                                               
      RLRV=-2.51E+10/(1.61*2.87E+6)                                     
      CONST=6.11*EXP(-RLRV/273.16) / 10.                                
      C1=CONST*  0.622                                                  
      C2=CONST*(-0.378)                                                 
CX10733 CALL DOINX( 733 )                                               
CY10733 CALL DOINY( 733 )                                               
      DO 10 K=1, 7                                                      
      DO  9 I=1, 96                                                     
      DPOVG(I,K)=RLRV/TF(I,K)                                           
9     CONTINUE                                                          
10    CONTINUE                                                          
CX20734 CALL DOOUTX( 734 )                                              
CY20734 CALL DOOUTY( 734 )                                              
CX10735 CALL DOINX( 735 )                                               
CY10735 CALL DOINY( 735 )                                               
      DO 12 K=1, 7                                                      
      DO 11 I=1, 96                                                     
      QS(I,K) = EXP(DPOVG(I,K))                                         
11    CONTINUE                                                          
12    CONTINUE                                                          
CX20736 CALL DOOUTX( 736 )                                              
CY20736 CALL DOOUTY( 736 )                                              
CX10737 CALL DOINX( 737 )                                               
CY10737 CALL DOINY( 737 )                                               
      DO 14 K=1, 7                                                      
      DO 13 I=1, 96                                                     
      QS(I,K)=C1*QS(I,K)/(SL(K)*PS(I) + C2*QS(I,K))                     
13    CONTINUE                                                          
14    CONTINUE                                                          
CX20738 CALL DOOUTX( 738 )                                              
CY20738 CALL DOOUTY( 738 )                                              
CX10739 CALL DOINX( 739 )                                               
CY10739 CALL DOINY( 739 )                                               
      DO 2 K=1, 7                                                       
      DO 1 I=1, 96                                                      
      DPOVG(I,K)=DEL(K)*PS(I)/ 9.8                                      
1     CONTINUE                                                          
2     CONTINUE                                                          
CX20740 CALL DOOUTX( 740 )                                              
CY20740 CALL DOOUTY( 740 )                                              
      CALL XSTORE(PREC(1),0.0, 96 )                                     
      CALL XSTORE(SATC,1.0, 96 * 7 )                                    
C*** SET EVAPORATION FACTOR                                             
      FACT(1)=.25                                                       
      FACT(2)=.16                                                       
      FACT(3)=.09                                                       
CX10741 CALL DOINX( 741 )                                               
CY10741 CALL DOINY( 741 )                                               
      DO 22 K=4, 7                                                      
      FACT(K)=.04                                                       
22    CONTINUE                                                          
C*** PCPN PROCESS.....TOP LYR DOWNWARD                                  
CX20742 CALL DOOUTX( 742 )                                              
CY20742 CALL DOOUTY( 742 )                                              
      RHSAT = .80                                                       
CX10743 CALL DOINX( 743 )                                               
CY10743 CALL DOINY( 743 )                                               
      DO 2190 K=1, 7                                                    
      L= 7 +1-K                                                         
      IF(L.EQ.1)RHSAT=.90                                               
      DO 2110 I=1, 96                                                   
      SUPER(I)=(QF(I,L)-QS(I,L))*DPOVG(I,L)                             
      B(I)=PREC(I)                                                      
2110  CONTINUE                                                          
      DO 3000 I=1, 96                                                   
      IF (SUPER(I).LE.0.)  GO TO 2120                                   
C     HERE IF SUPER GT 0., COMPUTE WET-BULB                             
C     ADJUSTMENT TO T AND Q, AND AUGMENT                                
C     PRECIPITATION FALLING THROUGH COLUMN.                             
      ALPHA = -RLRV/(TF(I,L)**2)                                        
      BETA = 1./RLOCP + ALPHA*QS(I,L)                                   
      QCOND = (SUPER(I)/DPOVG(I,L)) * (1./(RLOCP*BETA))                 
      TF(I,L) = TF(I,L) + RLOCP*QCOND                                   
      QF(I,L) = QF(I,L) - QCOND                                         
      PREC(I) = PREC(I) + QCOND*DPOVG(I,L)                              
C     FINISHED WITH SUPER-SATURATED POINT                               
      GO TO 3000                                                        
2120  CONTINUE                                                          
C     HERE IF POINT IS NOT SUPER SATURATED                              
C     CHECK IF FALLING PRECIPITATION IS                                 
C     PRESENT AND IF IT SHOULD BE                                       
C     EVAPORATED.                                                       
      IF (PREC(I).LE.0.)  GO TO 3000                                    
      EVAP = (QF(I,L)-QS(I,L)*RHSAT) * DPOVG(I,L)                       
      IF (EVAP.GE.0.)  GO TO 3000                                       
      EVAP = -EVAP                                                      
C     HERE IF EVAPORATION OF FALLING                                    
C     PRECIPITATION IS TO BE COMPUTED                                   
C     COMPUTE POTENTIAL EVAPORATION                                     
      POTEVP = FACT(L)*EVAP                                             
C     ASK IF POTEVP EXCEEDS AVAILABLE                                   
C     FALLING PRECIPITATION, COMPUTE                                    
C     AMOUNT OF WATER EVAPORATED.                                       
      EXCESS = PREC(I) - POTEVP                                         
      AMTEVP = POTEVP                                                   
      IF (EXCESS.LT.0.)  AMTEVP=PREC(I)                                 
C     ADJUST T AND Q TO ACCOUNT                                         
C     FOR EVAPORATION, MODIFY PREC.                                     
      ALPHA = -RLRV/(TF(I,L)**2)                                        
      BETA = 1./RLOCP + ALPHA*QS(I,L)                                   
      QEVAP = (AMTEVP/DPOVG(I,L)) * (1./(RLOCP*BETA))                   
      TF(I,L) = TF(I,L) - RLOCP*QEVAP                                   
      QF(I,L) = QF(I,L) + QEVAP                                         
      PREC(I) = PREC(I) - QEVAP*DPOVG(I,L)                              
      PREC(I) = MAX(PREC(I),0.)                                         
C     FINISHED WITH POINT AT WHICH                                      
C     FALLING PRECIP IS EVAPORATED                                      
      GO TO 3000                                                        
3000  CONTINUE                                                          
2190  CONTINUE                                                          
C*** PCPN REACHES GROUND LEVEL....FACTOR OF .5 SINCE PCPN IS FOR        
C*** TWO (LEAPFROG) TIME-STEPS....ACCUM PCPN (GESHEM) IS NOT LEAPFROGGED
CX20744 CALL DOOUTX( 744 )                                              
CY20744 CALL DOOUTY( 744 )                                              
CX10745 CALL DOINX( 745 )                                               
CY10745 CALL DOINY( 745 )                                               
      DO 20 I=1, 96                                                     
      GESHEM(I)=GESHEM(I)+PREC(I)*.5                                    
  20   CONTINUE                                                         
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20746 CALL DOOUTX( 746 )                                              
CY20746 CALL DOOUTY( 746 )                                              
CX40056 CALL SBOUTX (56)                                                
CY40056 CALL SBOUTY (56)                                                
      RETURN                                                            
      END                                                               
      PROGRAM SPEC77                                                    
      COMMON IDATE(4),                                                  
     * GESHEM( 96 , 76 ),TSEA( 96 , 76 ),GRGRID( 96 , 76 )              
      COMMON DPDPHI( 992 ),DPDLAM( 961 ),RELVOR( 12 ),ABSVOR( 12 ),     
     1  Q( 961 ), TE( 961 , 12 ), DI( 961 , 12 ), ZE( 961 , 12 ),       
     2  Z( 961 ),  Y( 961 , 12 ),  X( 961 , 12 ),  W( 961 , 12 ),       
     3 QM( 961 ),TEM( 961 , 12 ),DIM( 961 , 12 ),ZEM( 961 , 12 ),       
     5 GZ( 961 ), RQ( 961 , 7  ), RM( 961 , 7  ), RT( 961 , 7  )        
      COMPLEX GZ,Q,TE,DI,ZE,Z,Y,X,W,QM,TEM,DIM,ZEM,RELVOR,              
     1DPDLAM,DPDPHI,ABSVOR,RQ,RM,RT                                     
      COMMON                                                            
     * PS( 96 ), DPHI( 96 ), DLAM( 96 ),                                
     * QF( 96 ),DPHIF( 96 ),DLAMF( 96 ),                                
     *  EF( 96 , 12 ), EEF( 96 , 12 ), CG( 96 , 12 ),                   
     *  TF( 96 , 12 ), TEF( 96 , 12 ),TAU( 96 , 12 ),                   
     *  AF( 96 , 12 ), ZEF( 96 , 12 ),  A( 96 , 12 ),                   
     *  BF( 96 , 12 ), DIF( 96 , 12 ),  B( 96 , 12 ),                   
     *  FF( 96 , 12 ),  UF( 96 , 12 ),  F( 96 , 12 ),                   
     *  GF( 96 , 12 ),  VF( 96 , 12 ),  G( 96 , 12 ),                   
     * UQF( 96 , 7 ), UQFF( 96 , 7 ),                                   
     * VQF( 96 , 7 ), VQFF( 96 , 7 ),                                   
     * RQF( 96 , 7 ), RQFF( 96 , 7 ), RTG( 96 , 7 )                     
      COMMON                                                            
     * U1( 96 ),V1( 96 ),T1( 96 ),Q1( 96 ),                             
     * U2( 96 ),V2( 96 ),T2( 96 ),Q2( 96 )                              
      COMMON/PLNCOM/PLN( 992 ),DER( 961 ),PLNWCS( 961 ),EPS( 992 ),     
     1 COLRAD( 38 ),WGT( 38 ),WGTCS( 38 ),RCS2( 38 )                    
      DOUBLE PRECISION COLRAD,EPS                                       
      COMMON/VERCOM/AM( 12 , 12 ),HM( 12 , 12 ),TM( 12 , 12 ),          
     O              BM( 12 , 12 ),CM( 12 , 12 ),EKIN( 12 ),             
     1 SI( 13 ),SL( 12 ),DEL( 12 ),RDEL2( 12 ),RMSDOT( 11 ),            
     2 CI( 13 ),CL( 12 ),TOV( 12 ),   SV( 12 ),   RPI( 11 ),            
     3 P1( 12 ),P2( 12 ), H1( 12 ),   H2( 12 ),RPIREC( 11 ),            
     4 ROTSIN,ROTCOS,SEADRY,SL1KAP,C1,C2,RLRV,SL100K                    
      REAL TIME1,TIME2,TIME3,TIME4                                      
C ----------------------MEMBER GCOM---------------------------------    
      COMMON /LAB8/LAB(8),KDATE(4),IFINAL                               
      LOGICAL W3AQ09,ISOP                                               
      COMMON/NMODF/NFILES                                               
      COMMON/MXSTP/ MSTP,ITERS                                          
      DIMENSION NUM(28)                                                 
      DATA JDT/0/, SHOUR/0./, JCAP/ 30 /, LEVS/ 12 /,                   
     1FILTA/ 0.92 /,DT/ 1029. /,LIMLOW/1/                               
      DATA PERCUT/2.7502E+4/,MODS/4/,NITER/2/,DISPK/0./                 
CX00001 CALL MAINIX (1)                                                 
CY00001 CALL MAINIY (1)                                                 
CX60000 CALL CNTINX ()                                                  
CY60000 CALL CNTINY ()                                                  
C                                                                       
      CALL CPUTIM(TIME1)                                                
      CALL ELAPSE(TIME2)                                                
C                                                                       
      OPEN (UNIT=6,FILE= 'WSO6')                                        
      CALL W3STRT                                                       
C     ISOP=W3AQ09(0)                                                    
      ISOP=.FALSE.                                                      
C     PRINT 100, JCAP, LEVS                                             
      WRITE(6,100) JCAP, LEVS                                           
100   FORMAT (1H0, 'SMF', I2, I2, 'G  CREATED JULY 8 1982')             
      SATCRI=0.8                                                        
      RELKUO=0.8                                                        
      ACUMT=-1.E-7                                                      
      FILTB = (1. - FILTA) * 0.5                                        
C     CALL ERRSET (208, 0, -1, 1)                                       
      CALL SETSIG(CI,SI,DEL,SL,CL,RPI)                                  
      CALL AMHMTM(DEL,RPI,SV,P1,P2,AM,HM,TM)                            
      CALL GLATS ( 38 , COLRAD, WGT, WGTCS, RCS2)                       
      CALL EPSLON (EPS,  30 )                                           
CX10747 CALL DOINX( 747 )                                               
CY10747 CALL DOINY( 747 )                                               
      DO 4 LEV=1, 12                                                    
4     TOV(LEV)=300.                                                     
CX20748 CALL DOOUTX( 748 )                                              
CY20748 CALL DOOUTY( 748 )                                              
C***                                                                    
      OPEN (UNIT=5,FILE= 'WSI5')                                        
      OPEN (UNIT=4,FILE= 'WSV')                                         
      CALL PRMFLD(NUM)                                                  
      IF(NUM(18).NE.0)MODS=NUM(18)                                      
      IF(NUM(19).NE.0)NITER=NUM(19)                                     
      N1=NUM(1)                                                         
      N2=NUM(2)                                                         
      N3=NUM(3)                                                         
      N4=NUM(4)                                                         
      IBRAD=NUM(5)                                                      
      ISEMI=NUM(6)                                                      
      NFILES=NUM(11)                                                    
      NF=NUM(12)                                                        
C***                                                                    
      OPEN (UNIT=N3,FILE= 'WSO14')                                      
      OPEN (UNIT=N4,FILE= 'WSO15')                                      
      OPEN (UNIT=N1,FILE= 'WSI18')                                      
      OPEN (UNIT=NF,FILE= 'WSI17')                                      
      CALL                                                              
     1    GRDLNF(N1,FHOUR,IDATE,DPHI,DLAM,GZ,QM,TEM,DIM,ZEM,RM,SL,SI)   
      REWIND N1                                                         
      IF(ISOP) CALL SHOWHR(FHOUR,1)                                     
      CALL                                                              
     1    GRDLNF(N2,FHOUR,IDATE,DPHI,DLAM,GZ,Q,TE,DI,ZE,RQ,SL,SI)       
      IF(N1.EQ.N2) GO TO 801                                            
C N1.NE.N2 *** IF FHOUR=MOD12 ZERO RAIN,IF NOT READ FROM SIG B FILE     
      IHOUR=FHOUR+0.5                                                   
      MOD12=MOD(IHOUR,12)                                               
      IF(MOD12.EQ.0)CALL XSTORE(GESHEM(1,1),0., 96 * 76 )               
      IF(MOD12.NE.0)READ(N2) GESHEM                                     
777   CONTINUE                                                          
      GO TO 802                                                         
801   CONTINUE                                                          
      CALL XSTORE(GESHEM(1,1),0., 96 * 76 )                             
802   CONTINUE                                                          
      REWIND N2                                                         
      IFINAL=NUM(15)                                                    
      ISST=NUM(16)                                                      
C***                                                                    
      OPEN (UNIT=ISST,FILE= 'WSI16')                                    
      CALL GSSTCD(Z00,SATCRI,ISST)                                      
      DK=NUM(9)                                                         
      DK=DK*10.**NUM(10)                                                
      TK=NUM(20)                                                        
      TK=TK*10.**NUM(21)                                                
      IF(NUM(20).EQ.0)TK=DK                                             
      WRITE(6,105)SATCRI,RELKUO,ACUMT,DT,FILTA,DK,TK                    
      ITERS=0                                                           
      MSTP=1                                                            
105   FORMAT(1H ,F4.2,1X,F4.2,1X,E8.2,1X,F5.0,1X,F4.2,1X,E8.2,1X,E8.2)  
C--------------------------------------------------------------         
      IF(IBRAD.EQ.1)                                                    
     1      CALL GNMINI(DK,TK,SATCRI,DISPK,PERCUT,MODS,NITER,NF)        
      IF(NUM(13).EQ.0) GO TO 5                                          
      NI=NUM(14)                                                        
      OPEN (UNIT=NI,FILE= 'WSO19')                                      
      CALL GWRITE(NI,SHOUR,IDATE,Z,Q,TE,DI,ZE,RQ,SL,SI,GZ,Z00)          
      REWIND NI                                                         
5     CONTINUE                                                          
C--------------------------------------------------------------         
      IF(N1.NE.N2) GO TO 6000                                           
      LIMLOW=2                                                          
C****************************************************************       
C***  FIRST STEP IS FORWARD. CORRESPONDS TO DT/2 IN SEMI-IMPLICIT       
C****************************************************************       
      DT= 1029. /4.                                                     
CX10749 CALL DOINX( 749 )                                               
CY10749 CALL DOINY( 749 )                                               
      DO 5000 JDT=1,2                                                   
      KDT=JDT                                                           
      CALL DRUK(EPS,PLN,Q,PS,QF,KDT)                                    
      CALL BMCM (TOV,P1,P2,H1,H2,DEL,CI,BM,CM,DT,SV,AM)                 
      CALL GLOOP(TPS)                                                   
      CALL DEL4(DK,TK,DER,RT,RM,X,DIM,Y,TEM,W,ZEM,Z,QM)                 
      IF(ISEMI.EQ.1)CALL SICDKD                                         
     1            (DIM,TEM,QM,X,Y,Z,AM,BM,CM,TOV,SV,DT,ZEM,W,GZ,DISPK)  
      IF(ISEMI.NE.1)CALL SIBCKD                                         
     1            (DIM,TEM,QM,X,Y,Z,AM,BM,CM,TOV,SV,DT,ZEM,W,GZ,DISPK)  
      DO 1003 LEV=1, 12                                                 
      DO 1003 I=1, 961                                                  
      DI(I,LEV)=X(I,LEV)                                                
      TE(I,LEV)=Y(I,LEV)                                                
      ZE(I,LEV)=W(I,LEV)                                                
1003  CONTINUE                                                          
      DO 1004 I=1, 961                                                  
      Q(I) = Z(I)                                                       
1004  CONTINUE                                                          
      CALL GWATER(FILTA,DT,0,SATCRI,ACUMT,RELKUO,NUN,NKUO)              
      DT=2.*DT                                                          
5000  CONTINUE                                                          
CX20750 CALL DOOUTX( 750 )                                              
CY20750 CALL DOOUTY( 750 )                                              
      SHOUR=DT                                                          
C****************************************************************       
C***  END SMOOTH START                                                  
C****************************************************************       
6000  CONTINUE                                                          
      MAXSTP=NUM(7)                                                     
      ITERS=MAXSTP                                                      
      ISTEPS=NUM(8)                                                     
      XHOUR=SHOUR                                                       
CX10751 CALL DOINX( 751 )                                               
CY10751 CALL DOINY( 751 )                                               
      DO 20000 ISTEP=1,ISTEPS                                           
C*************************************************                      
C*****        TIME LOOP                                                 
C*************************************************                      
C                                                                       
      CALL BMCM (TOV,P1,P2,H1,H2,DEL,CI,BM,CM,DT,SV,AM)                 
      DO 10000 JDT=LIMLOW,MAXSTP                                        
      KDT=JDT                                                           
      MSTP=KDT                                                          
      CALL DRUK(EPS,PLN,Q,PS,QF,KDT)                                    
      SHOUR=SHOUR+DT                                                    
      XHOUR=XHOUR+DT                                                    
      CHOUR=SHOUR/3600.                                                 
      IHOUR=CHOUR+0.5                                                   
      CHOUR=IHOUR                                                       
      THOUR=FHOUR+CHOUR                                                 
      CALL GLOOP(TPS)                                                   
      CALL DEL4(DK,TK,DER,RT,RM,X,DIM,Y,TEM,W,ZEM,Z,QM)                 
      IF(ISEMI.EQ.1)CALL SICDKD                                         
     1            (DIM,TEM,QM,X,Y,Z,AM,BM,CM,TOV,SV,DT,ZEM,W,GZ,DISPK)  
      IF(ISEMI.NE.1)CALL SIBCKD                                         
     1            (DIM,TEM,QM,X,Y,Z,AM,BM,CM,TOV,SV,DT,ZEM,W,GZ,DISPK)  
      CALL FILTER (TEM,TE,Y,DIM,DI,X,ZEM,ZE,W,QM,Q,Z,FILTA)             
      CALL GWATER(FILTA,DT,1,SATCRI,ACUMT,RELKUO,NUN,NKUO)              
C***  COMPLETE TIME FILTER OF T(N-1). TE IS NOW ADJUSTED T(N+1)         
      DO 6 K=1, 12                                                      
      DO 6 I=1, 961                                                     
      TEM(I,K)=CMPLX( REAL(TEM(I,K))+FILTB* REAL(TE(I,K))  ,            
     1               AIMAG(TEM(I,K))+FILTB*AIMAG(TE(I,K))  )            
6     CONTINUE                                                          
      IF(XHOUR.LT.3000.) GO TO 10000                                    
      IF(ISOP) CALL SHOWHR(THOUR,4)                                     
      XHOUR=0.                                                          
10000 CONTINUE                                                          
C************************************************                       
C*****     END TIME LOOP                                                
C************************************************                       
C     PRINT 102,DT,CHOUR                                                
      WRITE(6,102)DT,CHOUR                                              
102   FORMAT(1H0,'TIME STEP ',F8.1,2X,'FCST SEGMENT OF ',F7.2,' H')     
      IF(ISTEP.EQ.2)N4=NUM(17)                                          
      CALL GWRITE(N4,THOUR,IDATE,Z,Q,TE,DI,ZE,RQ,SL,SI,GZ,Z00)          
C     WRITE (N4) GESHEM                                                 
      WRITE(N4,8000)GESHEM                                              
8000  FORMAT(6E13.7)                                                    
      REWIND N4                                                         
      IHOUR=THOUR+0.5                                                   
      MOD12=MOD(IHOUR,12)                                               
      IF(MOD12.EQ.0)CALL XSTORE(GESHEM(1,1),0., 96 * 76 )               
      LIMLOW=1                                                          
20000 CONTINUE                                                          
CX20752 CALL DOOUTX( 752 )                                              
CY20752 CALL DOOUTY( 752 )                                              
      CALL GWRITE(N3,THOUR,IDATE,Z,QM,TEM,DIM,ZEM,RM,SL,SI,GZ,Z00)      
      REWIND N3                                                         
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
2200  CALL W3STOP                                                       
      CALL CPUTIM(TIME3)                                                
      CALL ELAPSE(TIME4)                                                
C                                                                       
      WRITE(4,*)                                                        
      WRITE(4,*) 'ELAPSED CPU TIME IN SECONDS: ',TIME3-TIME1            
C MFLOPS = MFLOP ON 1 CPU CRAY X-MP/CPU TIME                            
      WRITE(4,*) 'MFLOPS RATE: ', 1852.767/(TIME3-TIME1)                
      WRITE(4,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ',                 
     *            MOD(TIME4-TIME2+1000000.0,1000000.0)                  
C                                                                       
      CLOSE(UNIT=5)                                                     
      CLOSE(UNIT=6)                                                     
      CLOSE(UNIT=N1)                                                    
      CLOSE(UNIT=N3)                                                    
      CLOSE(UNIT=N4)                                                    
      CLOSE(UNIT=NF)                                                    
      CLOSE(UNIT=NI)                                                    
      CLOSE(UNIT=ISST)                                                  
      CLOSE(UNIT=4)                                                     
CX70000 CALL CNTOTX ()                                                  
CY70000 CALL CNTOTY ()                                                  
CX50000 CALL MAINOX ()                                                  
CY50000 CALL MAINOY ()                                                  
      STOP                                                              
      END                                                               
      SUBROUTINE MATVEC (A, B, C)                                       
      DIMENSION A( 12 , 12 ),B( 12 ),C( 12 )                            
      COMPLEX B, C, SUM                                                 
CX30059 CALL SBINX (59)                                                 
CY30059 CALL SBINY (59)                                                 
CX10753 CALL DOINX( 753 )                                               
CY10753 CALL DOINY( 753 )                                               
      DO 2 I=1, 12                                                      
      SUM = CMPLX(0.,0.)                                                
      DO 1 J=1, 12                                                      
C     SUM = SUM + A(I,J) * B(J)                                         
      SUM=SUM+CMPLX(A(I,J)*REAL(B(J)),A(I,J)*AIMAG(B(J)))               
1     CONTINUE                                                          
      C(I) = SUM                                                        
2     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20754 CALL DOOUTX( 754 )                                              
CY20754 CALL DOOUTY( 754 )                                              
CX40060 CALL SBOUTX (60)                                                
CY40060 CALL SBOUTY (60)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE MINMAX (A, LGTH, CNST, SMIN, SMAX, IFCN, N)            
      DIMENSION A(LGTH), CNST(4)                                        
CX30061 CALL SBINX (61)                                                 
CY30061 CALL SBINY (61)                                                 
      SMIN = A(1)                                                       
      SMAX = A(1)                                                       
CX10755 CALL DOINX( 755 )                                               
CY10755 CALL DOINY( 755 )                                               
      DO 25 I=2,LGTH                                                    
      IF (A(I).GT.SMAX)  SMAX = A(I)                                    
      IF (A(I).LT.SMIN)  SMIN = A(I)                                    
25    CONTINUE                                                          
CX20756 CALL DOOUTX( 756 )                                              
CY20756 CALL DOOUTY( 756 )                                              
C                                                                       
C     IF IFCN IS EQUAL TO 1, CNST SUPPLIED                              
C     IF IFCN IS NOT EQUAL TO 1, CNST CALCULATED WITH SMIN AND SMAX     
C                                                                       
      IF (IFCN.EQ.1) THEN                                               
CX40062 CALL SBOUTX (62)                                                
CY40062 CALL SBOUTY (62)                                                
      RETURN                                                            
      ENDIF                                                             
      SCALE = ABS(SMAX)                                                 
      IF (ABS(SMIN).GT.SCALE)  SCALE = ABS(SMIN)                        
      IF (SCALE.NE.0.0)  GO TO 1                                        
      CNST(1) = 0.0                                                     
      CNST(2) = 0.0                                                     
      CNST(3) = 0.0                                                     
      CNST(4) = 0.0                                                     
C     PRINT 100, SCALE                                                  
      WRITE(6,100) SCALE                                                
100   FORMAT (1H0, 'SCALE = ', E10.3)                                   
CX40062 CALL SBOUTX (62)                                                
CY40062 CALL SBOUTY (62)                                                
      RETURN                                                            
1     CONTINUE                                                          
      CNST(1) = 0.0                                                     
      CNST(2) = 9999./SCALE                                             
      CNST(3) = (SMAX-SMIN)/N * CNST(2)                                 
      CNST(4) = 0.                                                      
CX40062 CALL SBOUTX (62)                                                
CY40062 CALL SBOUTY (62)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE MOVEFF(A,B)                                            
      DIMENSION A( 31 ), B( 31 )                                        
      COMPLEX A,B                                                       
CX30063 CALL SBINX (63)                                                 
CY30063 CALL SBINY (63)                                                 
CX10757 CALL DOINX( 757 )                                               
CY10757 CALL DOINY( 757 )                                               
C     NOTE. A IN CALL IS REAL LATITUDE CIRCLE                           
      DO 1 I=1, 31                                                      
1     A(I) =B(I)                                                        
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20758 CALL DOOUTX( 758 )                                              
CY20758 CALL DOOUTY( 758 )                                              
CX40064 CALL SBOUTX (64)                                                
CY40064 CALL SBOUTY (64)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE MSTADB(H,PRESS,TIN,IT,KM,THE,TMST,QMST,KTOP,IER)       
C     THIS ROUTINE ACCEPTS INPUT DATA OF                                
C         H=SURFACE PRESSURE DVDED BY 100 CBS                           
C         PRESS(K=1,--,KM) = FACTOR SUCH THAT P AT LVL K = H*PRESS(K)   
C         TIN(1 + (K-1)*IT ) = INPUT TEMPS IN A COLUMN                  
C         KM=NUMBER OF LEVELS                                           
C     AND RETURNS                                                       
C         THE IS EQUIV POT TEMP OF LAYER 1 (ASSUMING SATN)              
C         KTOP (G.T.E 1 AND L.T.E. KM) IS HIGHEST LVL FOR WHICH         
C              TIN IS COLDER THAN MOIST ADIABAT GIVEN BY THE            
C              (KTOP=1 DENOTES TIN(K=2) IS ALREADY G.T.E. MOIST ADB.)   
C                   ALLOWANCE IS MADE FOR PERHAPS ONE LEVEL BELOW KTOP  
C                  AT WHICH TIN WAS WARMER THAN TMST.                   
C         TMST(K) AND QMST(K) FOR K=2,--,KTOP, ARE TEMP AND SAT SPEC    
C              HUMS ON MOIST ADB THE.                                   
C                                                                       
C     OTHER SUBROUTINES USED ARE SATVAP AND SATVP1                      
C                                                                       
C     PRESSURE IN LAYER ONE MUST LIE BETWEEN 50 AND 110 CBS.            
C     TEMP IN LAYER ONE MUST LIE BETWEEN 220 AND 330 DEGREES            
C     THE RESULTING THE MUST LIE BETWEEN 220 AND 500 DEGREES.           
C         ( THE IS TESTED FOR THIS POSSIBILITY, AN ERROR RETURN OF      
C         IER=0 DENOTING OKEH CONDITIONS, A RETURN OF IER=1 DENOTING    
C         VIOLATION OF THIS RANGE.)                                     
C                                                                       
      REAL   T3,CP,L3,RD,RV,CL,THETAE(13,12),T,ES(13),KAPPA             
      REAL   TFMTHE(29,12),ESS                                          
      REAL   QFMTHE(29,12)                                              
      DIMENSION PRESS(1),TIN(1),TMST(1),QMST(1)                         
      DATA T3/273.16/,CP/1005./,RD/287.05/,RV/461.5/,L3/2.501E6/,       
     C  CL/4187./                                                       
      DATA ICALL /0/                                                    
CX30065 CALL SBINX (65)                                                 
CY30065 CALL SBINY (65)                                                 
      IER=0                                                             
      IF (ICALL) 1,1,100                                                
C     SET UP TABLES OF THETA E OF P AND T, AND OF T OF THETA E AND P    
    1 ICALL = 1                                                         
      KAPPA=RD/CP                                                       
      EPS=RD/RV                                                         
      T=220.                                                            
CX10759 CALL DOINX( 759 )                                               
CY10759 CALL DOINY( 759 )                                               
      DO 3 I=1,13                                                       
      EL=L3-(CL-CP)*(T-T3)                                              
      CALL SATVAP(T,ES(I),1)                                            
      P=.5                                                              
          DO 2 K=1,12                                                   
          RATIO=EPS*ES(I)/(100.*P-ES(I) )                               
          POWER=EL*RATIO/(CP*T)                                         
          PD=P-.01*ES(I)                                                
          THETAE(I,K)=T*EXP(POWER)/(PD**KAPPA)                          
          P=P+.05                                                       
    2     CONTINUE                                                      
      T=T+10.                                                           
    3 CONTINUE                                                          
CX20760 CALL DOOUTX( 760 )                                              
CY20760 CALL DOOUTY( 760 )                                              
C     CONSTRUCT TABLE TO GET TEMPERATURE FROM THETAE AND P              
C         THE TABLE IS TFMTHE(I,K), WHERE                               
C              THETAE=220,230,--,500 = 210+10*I, AND                    
C                  P=PRESS/100 CB = 0,.1,.2,--,1.1 =  0.1*(K-1)         
C     FOR P=0, TFMTHE(I,1) IS IDENTICALLY ZERO                          
C                                                                       
      CRIT=.01                                                          
      THEE=220.                                                         
CX10761 CALL DOINX( 761 )                                               
CY10761 CALL DOINY( 761 )                                               
      DO 16 I=1,29                                                      
      TFMTHE(I,1)=0.                                                    
      QFMTHE(I,1) = 0.                                                  
      P=.1                                                              
      DO 12 K=2,12                                                      
C     FIRST GUESS                                                       
      T=THEE*(P**KAPPA)                                                 
      T=310.                                                            
C     NEWTON ITERATION METHOD                                           
          DO 10 L=1,20                                                  
          CALL SATVAP(T,ESS,1)                                          
          RATIO=EPS*ESS/(100.*P-ESS)                                    
          EL=L3-(CL-CP)*(T-T3)                                          
          POWER=RATIO*EL/(CP*T)                                         
          FUN=T*EXP(POWER)/((P-.01*ESS)**KAPPA)                         
          DFUN=(FUN/T)*(1.+(RATIO/(CP*T) )*((CP-CL)*T                   
     C        +(P/(P-.01*ESS))*(EL*EL/(RV*T))))                         
          CHG=(THEE-FUN)/DFUN                                           
          T=T+CHG                                                       
          IF(ABS(CHG)-CRIT) 11,10,10                                    
   10     CONTINUE                                                      
   11 TFMTHE(I,K)=T                                                     
      CALL SATVAP(T,ESS,1)                                              
      QFMTHE(I,K)=EPS*ESS/(100.*P-(1.-EPS)*ESS)                         
      P=P+.1                                                            
   12 CONTINUE                                                          
C                                                                       
      THEE=THEE+10.                                                     
   16 CONTINUE                                                          
CX20762 CALL DOOUTX( 762 )                                              
CY20762 CALL DOOUTY( 762 )                                              
C                                                                       
  100 CONTINUE                                                          
C     SURFACE PRESS IS H, PRESSURE AT OTHER LVLS IS H TIMES PRESS(K)    
C     COMPUTE THETA E  = THE  , FOR LAYER ONE                           
      P=H*PRESS(1)                                                      
      TI=.1*TIN(1)-21.                                                  
      JT=TI                                                             
      X=TI-JT                                                           
      PK=20.*P-9.                                                       
      KP=PK                                                             
      Y=PK-KP                                                           
      YY=1.-Y                                                           
      XX=1.-X                                                           
      THE=XX*(YY*THETAE(JT,KP)+Y*THETAE(JT,KP+1))                       
     C     + X*(YY*THETAE(JT+1,KP)+Y*THETAE(JT+1,KP+1))                 
C                                                                       
C     GET T AND Q ON MST ADIABAT FOR LAYERS WHICH ARE COLDER THAN       
C         THE MST ADIABAT                                               
      KTOP=1                                                            
      TK=.1*THE-21.                                                     
      KT=TK                                                             
      Y=TK-KT                                                           
      YY=1.-Y                                                           
      KADD=IT+1                                                         
C                                                                       
C     TEST LIMITS ON THE                                                
      IF(KT-1) 105,106,106                                              
C     THE IS L.T. 220  OR  G.T. 500.                                    
  104 FORMAT(1H1,4X,17HEQUIV POT TEMP = ,E14.6,23H OUT OF RANGE IN MSTAD
     CB)                                                                
C*******************************************                            
105   CONTINUE                                                          
C 105 PRINT 104,THE                                                     
C     CALL PDUMP(TIN(1),TIN(1),5,H,H,5,IT,IT,4)                         
C*******************************************                            
      IER=1                                                             
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
      RETURN                                                            
  106 IF(KT-28) 107,107,105                                             
  107 CONTINUE                                                          
C                                                                       
C     FIRST PUT TMST=TIN AND QMST=QSAT(TIN) FOR K=1                     
      TMST(1)=TIN(1)                                                    
      P=H*PRESS(1)                                                      
      PI=10.*P+1.                                                       
      IP=PI                                                             
      X=PI-IP                                                           
      QMST(1)=(1.-X)*(YY*QFMTHE(KT,IP)+Y*QFMTHE(KT+1,IP) )              
     C        +X*(YY*QFMTHE(KT,IP+1)+Y*QFMTHE(KT+1,IP+1) )              
C                                                                       
C     WE WILL ALLOW ONE STABLE LAYER (WITH TIN G.T. TMST) TO            
C         INTERRUPT A SEQUENCE OF UNSTABLE LAYERS.                      
      LSTB=0                                                            
CX10763 CALL DOINX( 763 )                                               
CY10763 CALL DOINY( 763 )                                               
      DO 102 K=2,KM                                                     
      P=H*PRESS(K)                                                      
      PI=10.*P+1.                                                       
      IP=PI                                                             
      X=PI-IP                                                           
      TMST(K)=(1.-X)*(YY*TFMTHE(KT,IP)+Y*TFMTHE(KT+1,IP))               
     C          +X*(YY*TFMTHE(KT,IP+1)+Y*TFMTHE(KT+1,IP+1))             
      QMST(K)=(1.-X)*(YY*QFMTHE(KT,IP)+Y*QFMTHE(KT+1,IP))               
     C          +X*(YY*QFMTHE(KT,IP+1)+Y*QFMTHE(KT+1,IP+1))             
      IF(TMST(K   )-TIN(KADD) ) 103,103,101                             
C     WE HAVE REACHED A STABLE LAYER.  IS IT THE FIRST SUCH,,,,         
  103 IF(LSTB) 108,108,109                                              
C     YES, IT IS THE FIRST TIME.  RECORD KTOP AND, IF K IS L.T. KM,     
C         REJOIN THE K-LOOP                                             
  108 LSTB=1                                                            
      KTOP1=KTOP                                                        
      IF(K-KM) 101,111,111                                              
C                                                                       
C     TIN IS L.T. T OF MST ADIAB,,CONTINUE                              
  101 KTOP=KTOP+1                                                       
      KADD=KADD+IT                                                      
  102 CONTINUE                                                          
CX20764 CALL DOOUTX( 764 )                                              
CY20764 CALL DOOUTY( 764 )                                              
C                                                                       
C                                                                       
C     WE HAVE GONE THRU THE LOOP WITH ALL LAYERS UNSTABLE,,WOW,,,       
      GO TO 111                                                         
C                                                                       
C     FOUND THE SECOND STBL LAYER.  WAS IT JUST ABOVE THE FIRST STBL LYR
  109 IF(KTOP-KTOP1-1) 110,110,111                                      
C     SECOND STBL LYR WAS JUST ABV FIRST ONE.  RESTORE KTOP             
  110 KTOP=KTOP1                                                        
C                                                                       
  111 CONTINUE                                                          
CX40066 CALL SBOUTX (66)                                                
CY40066 CALL SBOUTY (66)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE MSU22(AM,AP,BM,BP,FLN,PLNWCS,DER)                      
      DIMENSION AM( 31 ),AP( 31 ),BM( 31 ),BP( 31 ),                    
     1FLN( 31 , 31 ),PLNWCS( 31 , 31 ),DER( 31 , 31 )                   
      COMPLEX AM,AP,BM,BP,FLN                                           
CX30067 CALL SBINX (67)                                                 
CY30067 CALL SBINY (67)                                                 
CX10765 CALL DOINX( 765 )                                               
CY10765 CALL DOINY( 765 )                                               
      DO 3 LL=1, 31                                                     
      DO 1 I=1, 31 ,2                                                   
      FLN(I,LL) = FLN(I,LL) +                                           
     1 CMPLX(REAL(BM(LL))*DER(I,LL)        +AIMAG(AP(LL))*PLNWCS(I,LL), 
     2      AIMAG(BM(LL))*DER(I,LL)        - REAL(AP(LL))*PLNWCS(I,LL)) 
1     CONTINUE                                                          
      DO 2 I=2, 30  ,2                                                  
      FLN(I,LL) = FLN(I,LL) +                                           
     1 CMPLX(REAL(BP(LL))*DER(I,LL)        +AIMAG(AM(LL))*PLNWCS(I,LL), 
     2      AIMAG(BP(LL))*DER(I,LL)        - REAL(AM(LL))*PLNWCS(I,LL)) 
2     CONTINUE                                                          
3     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20766 CALL DOOUTX( 766 )                                              
CY20766 CALL DOOUTY( 766 )                                              
CX40068 CALL SBOUTX (68)                                                
CY40068 CALL SBOUTY (68)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE NEGPLN(PLN)                                            
      DIMENSION PLN( 992 )                                              
C SUB. CHANGES SIGN OF ASYM. PLN IN S. HEMI.                            
CX30069 CALL SBINX (69)                                                 
CY30069 CALL SBINY (69)                                                 
CX10767 CALL DOINX( 767 )                                               
CY10767 CALL DOINY( 767 )                                               
      DO 2 LL=1, 31                                                     
      JLE=(LL-1)* 32                                                    
      DO 1 IND=2, 32 ,2                                                 
      PLN(JLE+IND)=-PLN(JLE+IND)                                        
1     CONTINUE                                                          
2     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20768 CALL DOOUTX( 768 )                                              
CY20768 CALL DOOUTY( 768 )                                              
CX40070 CALL SBOUTX (70)                                                
CY40070 CALL SBOUTY (70)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE PLN2 (PLN, COLRAD, LAT, EPS)                           
      DIMENSION PLN( 992 ), COLRAD( 38 ), EPS( 992 )                    
      DOUBLE PRECISION COLRAD, EPS, SINLAT, COS2, PROD,                 
     1 A, B, FL, ALP1, ALP2, ALP3                                       
CX30071 CALL SBINX (71)                                                 
CY30071 CALL SBINY (71)                                                 
      SINLAT = COS(COLRAD(LAT))                                         
      COS2 = 1.D+0 - SINLAT * SINLAT                                    
      PROD = 1.D+0                                                      
      A = 1.D+0                                                         
      B = 0.D+0                                                         
CX10769 CALL DOINX( 769 )                                               
CY10769 CALL DOINY( 769 )                                               
      DO 30 LL=1, 31                                                    
      L = LL - 1                                                        
      FL = L                                                            
      JLE = L *  32                                                     
      IF (L.EQ.0) GO TO 12                                              
      A = A + 2.D+0                                                     
      B = B + 2.D+0                                                     
      PROD = PROD * COS2 * A / B                                        
12    ALP1 = SQRT (0.5D+0 * PROD)                                       
      PLN(JLE+1) = REAL(ALP1)                                           
      ALP2 = SQRT(2.D+0 * FL + 3.D+0) * SINLAT * ALP1                   
      PLN(JLE+2) = REAL(ALP2)                                           
      DO 20 N=3, 32                                                     
      MN = JLE + N                                                      
      ALP3 = (SINLAT*ALP2 - EPS(MN-1)*ALP1)/EPS(MN)                     
      PLN(MN) = REAL(ALP3)                                              
      ALP1 = ALP2                                                       
      ALP2 = ALP3                                                       
20    CONTINUE                                                          
30    CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20770 CALL DOOUTX( 770 )                                              
CY20770 CALL DOOUTY( 770 )                                              
CX40072 CALL SBOUTX (72)                                                
CY40072 CALL SBOUTY (72)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE POLY(N,RAD,P)                                          
      DOUBLE PRECISION X,RAD,P,Y1,Y2,Y3,G                               
CX30073 CALL SBINX (73)                                                 
CY30073 CALL SBINY (73)                                                 
      X=COS(RAD)                                                        
      Y1=1.D+0                                                          
      Y2=X                                                              
CX10771 CALL DOINX( 771 )                                               
CY10771 CALL DOINY( 771 )                                               
      DO 1 I=2,N                                                        
      G=X*Y2                                                            
      Y3=G-Y1+G-(G-Y1)/FLOAT(I)                                         
      Y1=Y2                                                             
      Y2=Y3                                                             
1     CONTINUE                                                          
CX20772 CALL DOOUTX( 772 )                                              
CY20772 CALL DOOUTY( 772 )                                              
      P=Y3                                                              
CX40074 CALL SBOUTX (74)                                                
CY40074 CALL SBOUTY (74)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE PRIMES(VORD,DIVD,COMD,GH,MODS,INOUT)                   
      COMPLEX                                                           
                                                                        
     1VORD( 31 , 31 , 12 ),                                             
     2DIVD( 31 , 31 , 12 ),                                             
     2COMD( 31 , 31 , 12 )                                              
      DIMENSION GH( 12 )                                                
CX30075 CALL SBINX (75)                                                 
CY30075 CALL SBINY (75)                                                 
C***  INOUT=-1 SCALE TO PRIMES,=+1 RETURN FROM PRIMES                   
      IF(INOUT)10,20,30                                                 
10    CONTINUE                                                          
CX10773 CALL DOINX( 773 )                                               
CY10773 CALL DOINY( 773 )                                               
      DO 4 K=1,MODS                                                     
      T=1./( 6370000. *SQRT(GH(K)))                                     
      DO 3 LL=1, 31                                                     
      N=LL-2                                                            
      DO 2 I=1, 31                                                      
      N=N+1.                                                            
      F=1.                                                              
      IF(N.EQ.0)GO TO 1                                                 
      F=1./SQRT(N*(N+1.))                                               
1     CONTINUE                                                          
      VORD(I,LL,K)=CMPLX(F*REAL(VORD(I,LL,K)),F*AIMAG(VORD(I,LL,K)))    
      DIVD(I,LL,K)=CMPLX(-F*AIMAG(DIVD(I,LL,K)), F*REAL(DIVD(I,LL,K)))  
      COMD(I,LL,K)=CMPLX(T*REAL(COMD(I,LL,K)),T*AIMAG(COMD(I,LL,K)))    
2     CONTINUE                                                          
3     CONTINUE                                                          
4     CONTINUE                                                          
CX20774 CALL DOOUTX( 774 )                                              
CY20774 CALL DOOUTY( 774 )                                              
   20 CONTINUE                                                          
CX40076 CALL SBOUTX (76)                                                
CY40076 CALL SBOUTY (76)                                                
      RETURN                                                            
30    CONTINUE                                                          
CX10775 CALL DOINX( 775 )                                               
CY10775 CALL DOINY( 775 )                                               
      DO 8 K=1,MODS                                                     
      T=     6370000. *SQRT(GH(K))                                      
      DO 7 LL=1, 31                                                     
      N=LL-2                                                            
      DO 6 I=1, 31                                                      
      N=N+1.                                                            
      F=   SQRT(N*(N+1.))                                               
      VORD(I,LL,K)=CMPLX(F*REAL(VORD(I,LL,K)),F*AIMAG(VORD(I,LL,K)))    
      DIVD(I,LL,K)=CMPLX( F*AIMAG(DIVD(I,LL,K)),-F*REAL(DIVD(I,LL,K)))  
      COMD(I,LL,K)=CMPLX(T*REAL(COMD(I,LL,K)),T*AIMAG(COMD(I,LL,K)))    
6     CONTINUE                                                          
7     CONTINUE                                                          
8     CONTINUE                                                          
CX20776 CALL DOOUTX( 776 )                                              
CY20776 CALL DOOUTY( 776 )                                              
CX40076 CALL SBOUTX (76)                                                
CY40076 CALL SBOUTY (76)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE PRINTA (A, JCAP1)                                      
      DIMENSION A(JCAP1,JCAP1)                                          
CX30077 CALL SBINX (77)                                                 
CY30077 CALL SBINY (77)                                                 
      R = -1.E+30                                                       
CX10777 CALL DOINX( 777 )                                               
CY10777 CALL DOINY( 777 )                                               
      DO 1 I=1,JCAP1                                                    
      DO 1 J=1,JCAP1                                                    
      IF (ABS(A(I,J)).GT.R) R=ABS(A(I,J))                               
1     CONTINUE                                                          
CX20778 CALL DOOUTX( 778 )                                              
CY20778 CALL DOOUTY( 778 )                                              
C     PRINT 100, R                                                      
      WRITE(6,100) R                                                    
100   FORMAT (1H0, 'SCALE OF MATRIX  =', E12.5)                         
CX10779 CALL DOINX( 779 )                                               
CY10779 CALL DOINY( 779 )                                               
      DO 2 I=1,JCAP1                                                    
      DO 2 J=1,JCAP1                                                    
      A(I,J) = A(I,J) / R                                               
2     CONTINUE                                                          
CX20780 CALL DOOUTX( 780 )                                              
CY20780 CALL DOOUTY( 780 )                                              
CX10781 CALL DOINX( 781 )                                               
CY10781 CALL DOINY( 781 )                                               
      DO 3 I=1,JCAP1                                                    
C     PRINT 101, (A(I,J), J=1,JCAP1)                                    
      WRITE(6,101) (A(I,J), J=1,JCAP1)                                  
C 101 FORMAT (1X, 21(F5.2, 1X))                                         
101   FORMAT (1X, 16(F6.3, 1X))                                         
3     CONTINUE                                                          
CX20782 CALL DOOUTX( 782 )                                              
CY20782 CALL DOOUTY( 782 )                                              
CX10783 CALL DOINX( 783 )                                               
CY10783 CALL DOINY( 783 )                                               
      DO 4 I=1,JCAP1                                                    
      DO 4 J=1,JCAP1                                                    
      A(I,J) = A(I,J) * R                                               
4     CONTINUE                                                          
CX20784 CALL DOOUTX( 784 )                                              
CY20784 CALL DOOUTY( 784 )                                              
CX40078 CALL SBOUTX (78)                                                
CY40078 CALL SBOUTY (78)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE PRMFLD(NUM)                                            
      DIMENSION NUM(28)                                                 
CX30079 CALL SBINX (79)                                                 
CY30079 CALL SBINY (79)                                                 
CX10785 CALL DOINX( 785 )                                               
CY10785 CALL DOINY( 785 )                                               
      DO 1 I=1,28                                                       
1     NUM(I)=0                                                          
CX20786 CALL DOOUTX( 786 )                                              
CY20786 CALL DOOUTY( 786 )                                              
      READ(5,100)NUM                                                    
100   FORMAT(28I2)                                                      
C     PRINT 101,NUM                                                     
      WRITE(6,101)NUM                                                   
101   FORMAT(1H0,'NUM=',28(1X,I2))                                      
CX40080 CALL SBOUTX (80)                                                
CY40080 CALL SBOUTY (80)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE PSU22(AM,AP,BM,BP,FLN,PLNWCS,DER)                      
      DIMENSION AM( 31 ),AP( 31 ),BM( 31 ),BP( 31 ),                    
     1FLN( 31 , 31 ),PLNWCS( 31 , 31 ),DER( 31 , 31 )                   
      COMPLEX AM,AP,BM,BP,FLN                                           
CX30081 CALL SBINX (81)                                                 
CY30081 CALL SBINY (81)                                                 
CX10787 CALL DOINX( 787 )                                               
CY10787 CALL DOINY( 787 )                                               
      DO 3 LL=1, 31                                                     
      DO 1 I=1, 31 ,2                                                   
      FLN(I,LL) = FLN(I,LL) +                                           
     1 CMPLX(REAL(AM(LL))*DER(I,LL)        -AIMAG(BP(LL))*PLNWCS(I,LL), 
     2      AIMAG(AM(LL))*DER(I,LL)        + REAL(BP(LL))*PLNWCS(I,LL)) 
1     CONTINUE                                                          
      DO 2 I=2, 30  ,2                                                  
      FLN(I,LL) = FLN(I,LL) +                                           
     1 CMPLX(REAL(AP(LL))*DER(I,LL)        -AIMAG(BM(LL))*PLNWCS(I,LL), 
     2      AIMAG(AP(LL))*DER(I,LL)        + REAL(BM(LL))*PLNWCS(I,LL)) 
2     CONTINUE                                                          
3     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20788 CALL DOOUTX( 788 )                                              
CY20788 CALL DOOUTY( 788 )                                              
CX40082 CALL SBOUTX (82)                                                
CY40082 CALL SBOUTY (82)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE RMS(Q,X,Y,Z)                                           
      COMPLEX X( 31 , 31 , 12 ),Y( 31 , 31 , 12 ),                      
     1Z( 31 , 31 , 12 ),Q( 31 , 31 )                                    
      DIMENSION D( 12 ),V( 12 ),T( 12 )                                 
CX30083 CALL SBINX (83)                                                 
CY30083 CALL SBINY (83)                                                 
      DT=0.                                                             
      QT=0.                                                             
      VT=0.                                                             
      TT=0.                                                             
CX10789 CALL DOINX( 789 )                                               
CY10789 CALL DOINY( 789 )                                               
      DO 10 LEV=1, 12                                                   
      D(LEV)=0.                                                         
      V(LEV)=0.                                                         
      T(LEV)=0.                                                         
      DO 15 LL=1, 31                                                    
      A=2.                                                              
      IF(LL .EQ. 1) A=1.                                                
      DO 15 NN=1, 31                                                    
      F=REAL(X(NN,LL,LEV))                                              
      G=AIMAG(X(NN,LL,LEV))                                             
      D(LEV)=D(LEV)+A*(F*F+G*G)                                         
      F=REAL(Y(NN,LL,LEV))                                              
      G=AIMAG(Y(NN,LL,LEV))                                             
      V(LEV)=V(LEV)+A*(F*F+G*G)                                         
      F=REAL(Z(NN,LL,LEV))                                              
      G=AIMAG(Z(NN,LL,LEV))                                             
   15 T(LEV)=T(LEV)+A*(F*F+G*G)                                         
      DT=DT+D(LEV)                                                      
      VT=VT+V(LEV)                                                      
      TT=TT+T(LEV)                                                      
      D(LEV)=SQRT(D(LEV)/2.)                                            
      V(LEV)=SQRT(V(LEV)/2.)                                            
   10 T(LEV)=SQRT(T(LEV)/2.)                                            
CX20790 CALL DOOUTX( 790 )                                              
CY20790 CALL DOOUTY( 790 )                                              
CX10791 CALL DOINX( 791 )                                               
CY10791 CALL DOINY( 791 )                                               
      DO 20 LL=1, 31                                                    
      A=2.                                                              
      IF(LL .EQ. 1) A=1.                                                
      DO 20 NN=1, 31                                                    
      F=REAL(Q(NN,LL))                                                  
      G=AIMAG(Q(NN,LL))                                                 
   20 QT=QT+A*(F*F+G*G)                                                 
CX20792 CALL DOOUTX( 792 )                                              
CY20792 CALL DOOUTY( 792 )                                              
      QT=SQRT(QT/2.)                                                    
      F=2.* 12                                                          
      DT=DT/F                                                           
      VT=VT/F                                                           
      TT=TT/F                                                           
      DT=SQRT(DT)                                                       
      VT=SQRT(VT)                                                       
      TT=SQRT(TT)                                                       
      WRITE(6,100) DT,(D(L),L=1, 12 )                                   
  100 FORMAT(' DIV',9E12.4)                                             
      WRITE(6,110) VT,(V(L),L=1, 12 )                                   
  110 FORMAT(' VORT',9E12.4)                                            
      WRITE(6,120) TT,(T(L),L=1, 12 )                                   
  120 FORMAT(' TEMP',9E12.4)                                            
      WRITE(6,130) QT                                                   
  130 FORMAT(' Q',E15.6)                                                
CX40084 CALL SBOUTX (84)                                                
CY40084 CALL SBOUTY (84)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE SATVAP(T,E,JMAX)                                       
      REAL   A(25),B(25),C(25),D(25),TT(25),AI,AL,BI,BL,E3,T3,          
     C  T(1),E(1)                                                       
      DATA E3/.611/,T3/273.16/                                          
      DATA AI/.0985921/,BI/22.488695/,AL/5.006501/,BL/19.839233/        
      DATA ICALL /0/                                                    
      SAVE TLAST,LM                                                     
CX30085 CALL SBINX (85)                                                 
CY30085 CALL SBINY (85)                                                 
C                                                                       
C     ********  THE FLWG TWO CARDS MAKE ICE SAME AS WATER***************
      AI=AL                                                             
      BI=BL                                                             
      IF(ICALL) 1,1,2                                                   
C     SET UP THE TABLES                                                 
    1 LM=25                                                             
      TLO=200.                                                          
      THI=320.                                                          
      CALL SATVP1(A,B,C,D,TT,TLO,THI,LM)                                
C                                                                       
C     E=SAT VAP PRESS IS NOW GIVEN BY                                   
C                                                                       
C         E(X) = A(L) + X*( B(L) + X*( C(L) + X*D(L) ) )                
C                  WHERE                                                
C     X = T - TT(L-1),   L=2,3,--,LM                                    
C                  AND                                                  
C              TT(L-1) IS L.T.E.  T  IS L.T.E. TT(L)                    
C                                                                       
C     FOR T L.T. TT(1) = TLO,  OR,  T G.T. TT(LM) = THI, WE USE         
C         E = E3 * (TAU**A) * EXP( (A+B)*(1.-TAU) ), WITH TAU=T3/T.     
C         AND A AND B DIFFER FOR ICE AND LIQUID.                        
C     LAST AND TLAST REMEMBER THE LAST USES TO SPEED UP THE TABLE USE   
      LAST=2                                                            
      TLAST=TLO+.5*(THI-TLO)/(LM-1.)                                    
      ICALL=1                                                           
C                                                                       
C     REJOIN THE JOB AT HAND                                            
C                                                                       
    2 CONTINUE                                                          
CX10793 CALL DOINX( 793 )                                               
CY10793 CALL DOINY( 793 )                                               
      DO      9      J=1,JMAX                                           
      TE=T(J)                                                           
      IF(TE-TLAST) 3,3,6                                                
C     TE IS L.T.E. TLAST.  SEARCH IN DECREASING T                       
    3 LL=LAST                                                           
          DO 4 M=2,LL                                                   
          LAST=LL+2-M                                                   
          IF(TE-TT(LAST-1) ) 4,5,5                                      
    4     CONTINUE                                                      
C     TE IS L.T. TT(1) = TLO                                            
      TLAST=.5*(TT(1)+TT(2) )                                           
      LAST=2                                                            
      TAU=T3/TE                                                         
      E(J)=E3*(TAU**AI)*EXP((AI+BI)*(1.-TAU) )                          
      GO TO 9                                                           
C       TE IS G.T.E. TT(LAST-1)                                         
    5 X=TE-TT(LAST-1)                                                   
      TLAST=TE                                                          
      E(J)=A(LAST)+X*(B(LAST)+X*(C(LAST)+X*D(LAST) ) )                  
      GO TO 9                                                           
C                                                                       
C     TE IS G.T. TLAST.  SEARCH IN INCREASING T                         
    6 LL=LAST                                                           
          DO 7 L=LL,LM                                                  
          LAST=L                                                        
          IF(TE-TT(L) ) 8,8,7                                           
    7     CONTINUE                                                      
C       T=TE IS G.T. TT(LM) = THI                                       
      LAST = LM                                                         
      TLAST=.5*(TT(LM)+TT(LM-1) )                                       
      TAU=T3/TE                                                         
      E(J)=E3*(TAU**AL)*EXP( (AL+BL)*(1.-TAU) )                         
      GO TO 9                                                           
C     T=TE IS L.T.E. TT(L)                                              
    8 X=TE-TT(LAST-1)                                                   
      E(J)=A(LAST)+X*(B(LAST)+X*(C(LAST)+X*D(LAST) ) )                  
      TLAST=TE                                                          
C         END OF J LOOP                                                 
    9 CONTINUE                                                          
CX20794 CALL DOOUTX( 794 )                                              
CY20794 CALL DOOUTY( 794 )                                              
CX40086 CALL SBOUTX (86)                                                
CY40086 CALL SBOUTY (86)                                                
C                                                                       
      RETURN                                                            
      END                                                               
      SUBROUTINE SATVP1 ( A,B,C,D,T,TLO,THI,LM)                         
C     ROUTINE TO COMPUTE INTERPOLATION TABLES FOR SAT VAP PRESS         
C     THAT WILL BE USED BY SUBROUTINE SATVAP                            
C     PHYSICAL CONSTANTS USED ARE FOR LIQUID WHEN T IS G.T.E. T3=       
C         TRIPLE POINT = 273.16, AND FOR ICE WHEN T IS L.T.E.           
C     TF = -30 DEG CELSIUS.  BETWEEN T3 AND TF A SMOOTH MIXTURE         
C         IS USED.                                                      
      REAL   A(LM),B(LM),C(LM),D(LM),T(LM)                              
      DOUBLE PRECISION R,CP,CI,CL,LATI,LATL,T3,TF,E3,DT,AI,AL,BI,BL,    
     C  TD,TAU,PL,PI,PDL,PDI,TEE,AA,AD,P1,PP1,P2,PP2                    
      DATA R/461.5/,CP/1876.5/,E3/.611/,T3/273.16/,TF/243.16/,          
     C  CI/1922./,CL/4187./,LATI/2.835D6/,LATL/2.501D6/                 
CX30087 CALL SBINX (87)                                                 
CY30087 CALL SBINY (87)                                                 
C                                                                       
C     ********  THE FLWG TWO CARDS MAKE ICE SAME AS WATER***************
      CI=CL                                                             
      LATI=LATL                                                         
      DT=(THI-TLO)/(LM-1.)                                              
      AI=(CI-CP)/R                                                      
      AL=(CL-CP)/R                                                      
      BI=LATI/(R*T3)                                                    
      BL=LATL/(R*T3)                                                    
      T(1)=TLO                                                          
      TD=TLO                                                            
CX10795 CALL DOINX( 795 )                                               
CY10795 CALL DOINY( 795 )                                               
      DO 18 L=2,LM                                                      
      TAU=T3/TD                                                         
      PL=E3*(TAU**AL)*EXP( (AL+BL)*(1.-TAU) )                           
      PI=E3*(TAU**AI)*EXP( (AI+BI)*(1.-TAU) )                           
      PDL=(PL/(R*TD*TD) )*(LATL-(CL-CP)*(TD-T3) )                       
      PDI=(PI/(R*TD*TD) )*(LATI-(CI-CP)*(TD-T3) )                       
      TEE=(TD-TF)/(T3-TF)                                               
      AA=TEE*TEE*(3.-2.*TEE)                                            
      AD=6.*TEE*(1.-TEE)/(T3-TF)                                        
      IF(TEE) 10,10,11                                                  
   10 AA=0.                                                             
      AD=0.                                                             
      GO TO 13                                                          
   11 IF(TEE-1.) 13,12,12                                               
   12 AA=1.                                                             
      AD=0.                                                             
   13 P1=AA*PL+(1.-AA)*PI                                               
      PP1=AD*(PL-PI)+AA*PDL+(1.-AA)*PDI                                 
      TD=TD+DT                                                          
      T(L)=TD                                                           
      TAU=T3/TD                                                         
      PL=E3*(TAU**AL)*EXP( (AL+BL)*(1.-TAU) )                           
      PI=E3*(TAU**AI)*EXP( (AI+BI)*(1.-TAU) )                           
      PDL=(PL/(R*TD*TD) )*(LATL-(CL-CP)*(TD-T3) )                       
      PDI=(PI/(R*TD*TD) )*(LATI-(CI-CP)*(TD-T3) )                       
      TEE=(TD-TF)/(T3-TF)                                               
      AA=TEE*TEE*(3.-2.*TEE)                                            
      AD=6.*TEE*(1.-TEE)/(T3-TF)                                        
      IF(TEE)  14,14,15                                                 
   14 AA=0.                                                             
      AD=0.                                                             
      GO TO 17                                                          
   15 IF(TEE-1.) 17,16,16                                               
   16 AA=1.                                                             
      AD=0.                                                             
   17 P2=AA*PL+(1.-AA)*PI                                               
      PP2=AD*(PL-PI)+AA*PDL+(1.-AA)*PDI                                 
      A(L)=P1                                                           
      B(L)=PP1                                                          
      C(L)=(3.*(P2-P1)-DT*(2.*PP1+PP2) )/(DT*DT)                        
      D(L)=(2.*(P1-P2)+DT*(PP1+PP2) )/(DT*DT*DT)                        
   18 CONTINUE                                                          
CX20796 CALL DOOUTX( 796 )                                              
CY20796 CALL DOOUTY( 796 )                                              
      A(1)=0.                                                           
      B(1)=0.                                                           
      C(1)=0.                                                           
      D(1)=0.                                                           
CX40088 CALL SBOUTX (88)                                                
CY40088 CALL SBOUTY (88)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE SETSIG (CI, SI, DEL, SL, CL, RPI)                      
      DIMENSION CI( 13 ), SI( 13 ),                                     
     1 DEL( 12 ), SL( 12 ), CL( 12 ), RPI( 11 )                         
CX30089 CALL SBINX (89)                                                 
CY30089 CALL SBINY (89)                                                 
C     PRINT 98                                                          
      WRITE(6,98)                                                       
98    FORMAT (1H0, 'BEGIN SETSIG')                                      
      DEL(1)=0.075                                                      
      DEL(2)=0.125                                                      
      DEL(3)=0.15                                                       
      DEL(4)=0.15                                                       
      DEL(5)=0.125                                                      
      DEL(6)=0.075                                                      
CX10797 CALL DOINX( 797 )                                               
CY10797 CALL DOINY( 797 )                                               
      DO 9889 K=7, 12                                                   
9889  DEL(K)=0.05                                                       
CX20798 CALL DOOUTX( 798 )                                              
CY20798 CALL DOOUTY( 798 )                                              
      CI(1) = 0.                                                        
CX10799 CALL DOINX( 799 )                                               
CY10799 CALL DOINY( 799 )                                               
      DO 54 K=1, 12                                                     
54    CI(K+1)=CI(K)+DEL(K)                                              
CX20800 CALL DOOUTX( 800 )                                              
CY20800 CALL DOOUTY( 800 )                                              
      CI( 13 )=1.                                                       
      RK =  287.05  /  1005.                                            
      RK1 = RK + 1.                                                     
CX10801 CALL DOINX( 801 )                                               
CY10801 CALL DOINY( 801 )                                               
      DO 1 LI=1, 13                                                     
1     SI(LI) = 1. - CI(LI)                                              
CX20802 CALL DOOUTX( 802 )                                              
CY20802 CALL DOOUTY( 802 )                                              
CX10803 CALL DOINX( 803 )                                               
CY10803 CALL DOINY( 803 )                                               
      DO 3 LE=1, 12                                                     
      DIF = SI(LE)**RK1 - SI(LE+1)**RK1                                 
      DIF = DIF / (RK1*(SI(LE)-SI(LE+1)))                               
      SL(LE) = DIF**(1./RK)                                             
      CL(LE) = 1. - SL(LE)                                              
3     CONTINUE                                                          
CX20804 CALL DOOUTX( 804 )                                              
CY20804 CALL DOOUTY( 804 )                                              
CX10805 CALL DOINX( 805 )                                               
CY10805 CALL DOINY( 805 )                                               
C     COMPUTE PI RATIOS FOR TEMP. MATRIX.                               
      DO 4 LE=1, 11                                                     
4     RPI(LE) = (SL(LE+1)/SL(LE))**RK                                   
CX20806 CALL DOOUTX( 806 )                                              
CY20806 CALL DOOUTY( 806 )                                              
CX10807 CALL DOINX( 807 )                                               
CY10807 CALL DOINY( 807 )                                               
      DO 5 LE=1, 13                                                     
C     PRINT 100, LE, CI(LE), SI(LE)                                     
      WRITE(6,100) LE, CI(LE), SI(LE)                                   
100   FORMAT (1H , 'LEVEL=', I2, 2X, 'CI=', F6.3, 2X, 'SI=', F6.3)      
5     CONTINUE                                                          
CX20808 CALL DOOUTX( 808 )                                              
CY20808 CALL DOOUTY( 808 )                                              
C     PRINT 97                                                          
      WRITE(6,97)                                                       
97    FORMAT (1H0)                                                      
CX10809 CALL DOINX( 809 )                                               
CY10809 CALL DOINY( 809 )                                               
      DO 6 LE=1, 12                                                     
C     PRINT 101, LE, CL(LE), SL(LE), DEL(LE)                            
      WRITE(6,101) LE, CL(LE), SL(LE), DEL(LE)                          
101   FORMAT (1H , 'LAYER=', I2, 2X, 'CL=', F6.3, 2X, 'SL=', F6.3, 2X,  
     1 'DEL=', F6.3)                                                    
6     CONTINUE                                                          
CX20810 CALL DOOUTX( 810 )                                              
CY20810 CALL DOOUTY( 810 )                                              
C     PRINT 102, (RPI(LE), LE=1, 11 )                                   
      WRITE(6,102) (RPI(LE), LE=1, 11 )                                 
102   FORMAT (1H0, 'RPI=',  11 (F6.3,2X))                               
C     PRINT 99                                                          
      WRITE(6,99)                                                       
99    FORMAT (1H0, 'END SETSIG')                                        
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX40090 CALL SBOUTX (90)                                                
CY40090 CALL SBOUTY (90)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE SHOWHR(THOUR,K)                                        
CX30091 CALL SBINX (91)                                                 
CY30091 CALL SBINY (91)                                                 
C     PRINT 100,THOUR                                                   
      WRITE(6,100)THOUR                                                 
100   FORMAT(1H ,'FORECAST REACHED ',F6.2,' HOURS')                     
CX40092 CALL SBOUTX (92)                                                
CY40092 CALL SBOUTY (92)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE SIBCKD (D,T,Q,X,Y,Z,AM,BM,CM,TOV,SV,DT,ZE,W,GZ,DISPK)  
C     D, T, Q ARE SPECTRAL COMPONENTS AT (T-1)                          
C     X, Y, Z ARE SPECTRAL COMPONENTS OF NONLIN. TERMS                  
      DIMENSION D( 961 , 12 ), T( 961 , 12 ),                           
     1  Q( 961 ), X( 961 , 12 ),                                        
     2  Y( 961 , 12 ), Z( 961 ), AM( 12 , 12 ), BM( 12 , 12 ),          
     3  CM( 12 , 12 ), TOV( 12 ), SV( 12 )                              
      DIMENSION ZE( 961 , 12 ), W( 961 , 12 ), GZ( 961 )                
C     VARIABLES USED ONLY IN SUB.                                       
      DIMENSION RIM( 12 , 12 ), DM( 12 , 12 ),                          
     1  EM( 12 , 12 ), LV( 12 ), MV( 12 ), XV( 12 ),                    
     2  YV( 12 ), SAV( 12 ), SBV( 12 ),                                 
     3  RHSV( 12 ), DV( 12 ), TV( 12 ), TOVRA( 12 )                     
      COMPLEX XV, YV, SAV, SBV, RHSV, DV, TV, SUM                       
      COMPLEX T, D, Q, X, Y, Z, ZE, W, GZ                               
CX30093 CALL SBINX (93)                                                 
CY30093 CALL SBINY (93)                                                 
C       DT2 = 2. * DT                                                   
      DT2 = DT + DT                                                     
CX10811 CALL DOINX( 811 )                                               
CY10811 CALL DOINY( 811 )                                               
C     TOV = R*TO/A**2                                                   
C     SET UNIT MATRIX RIM                                               
      DO 1 J=1, 12                                                      
C       TOVRA(J)= 287.05 *TOV(J)/( 6370000. * 6370000. )                
      TOVRA(J)=( 287.05 /( 6370000. * 6370000. ))*TOV(J)                
      DO 1 I=1, 12                                                      
      RIM(I,J) = 0.                                                     
1     CONTINUE                                                          
CX20812 CALL DOOUTX( 812 )                                              
CY20812 CALL DOOUTY( 812 )                                              
CX10813 CALL DOINX( 813 )                                               
CY10813 CALL DOINY( 813 )                                               
      DO 2 I=1, 12                                                      
      RIM(I,I) = 1. + DT2*DISPK                                         
2     CONTINUE                                                          
CX20814 CALL DOOUTX( 814 )                                              
CY20814 CALL DOOUTY( 814 )                                              
CCC   1000 IS (L,N) SPECTRAL COMPONENT LOOP                             
      LL1=1                                                             
      JMN=2* 30 +1                                                      
CX10815 CALL DOINX( 815 )                                               
CY10815 CALL DOINY( 815 )                                               
      DO 1000 NN=1,JMN                                                  
      N=NN-1                                                            
C***  CALCULATE MATRIX                                                  
C     COMPUTE MATRICES DM, EM FROM CM                                   
      DO 4 J=1, 12                                                      
      DO 4 I=1, 12                                                      
      DM(I,J) = RIM(I,J) - 4.*N*(N+1)*CM(I,J)                           
4     CONTINUE                                                          
      CALL IMINV (DM,  12 , DET, LV, MV)                                
C     DM NOW CONTAINS ITS INVERSE.                                      
      LL2=NN                                                            
      IF(N- 30 )900,900,800                                             
800   LL1=LL1+1                                                         
      LL2= 31                                                           
900   CONTINUE                                                          
      IUP=0                                                             
      DO 1000 LL=LL1,LL2                                                
      L=LL2-LL                                                          
      IUP=IUP+1                                                         
      IF(N.GT. 30 )L=LL2-IUP                                            
      IN =L* 31 +LL                                                     
C     SET COLUMN VECTORS XV, YV, DV, TV                                 
      DO 5 LEV=1, 12                                                    
      DV(LEV) = D(IN,LEV)                                               
      TV(LEV) = T(IN,LEV)                                               
      XV(LEV) = X(IN,LEV)                                               
      YV(LEV) = Y(IN,LEV)                                               
5     CONTINUE                                                          
C     BEGIN CONSTRUCTING RHSV  (FROM R TO L)                            
      DO 6 LEV=1, 12                                                    
C       RHSV(LEV) = (DT2*Z(IN) + Q(IN))*TOVRA(LEV)                      
      RHSV(LEV) = CMPLX(( REAL(Z(IN))*DT2 +  REAL(Q(IN)))*TOVRA(LEV),   
     1                  (AIMAG(Z(IN))*DT2 + AIMAG(Q(IN)))*TOVRA(LEV))   
6     CONTINUE                                                          
      DO 7 LEV=1, 12                                                    
C       SAV(LEV) = DT2 * YV(LEV) + TV(LEV)                              
      SAV(LEV) = CMPLX(DT2* REAL(YV(LEV)) +  REAL(TV(LEV)),             
     1                 DT2*AIMAG(YV(LEV)) + AIMAG(TV(LEV)))             
7     CONTINUE                                                          
C     MULTIPLY MATRIX AM BY SAV   PUT IN SBV                            
      CALL MATVEC (AM, SAV, SBV)                                        
      DO 10 LEV=1, 12                                                   
C       RHSV(LEV) = N * (N+1.) * (RHSV(LEV) + SBV(LEV))                 
      RHSV(LEV)= CMPLX(( REAL(RHSV(LEV)) +  REAL(SBV(LEV)))*(N*(N+1.)), 
     1                 (AIMAG(RHSV(LEV)) + AIMAG(SBV(LEV)))*(N*(N+1.))) 
10    CONTINUE                                                          
      DO 11 LEV=1, 12                                                   
C       RHSV(LEV) = DT2 * (XV(LEV) + RHSV(LEV))  +  DV(LEV)             
      RHSV(LEV) = CMPLX (                                               
     1   ( REAL(XV(LEV)) +  REAL(RHSV(LEV))) * DT2 +  REAL(DV(LEV)) ,   
     2   (AIMAG(XV(LEV)) + AIMAG(RHSV(LEV))) * DT2 + AIMAG(DV(LEV)) )   
11    CONTINUE                                                          
C     MULTIPLY RHSV BY DM INVERSE  (NOW IN DM)                          
       CALL MATVEC (DM, RHSV, SAV)                                      
C     SAV IS NOW D VECTOR AT T+1                                        
C     STORE D AT T+1 IN X                                               
      DO 15 LEV=1, 12                                                   
      X(IN,LEV) = SAV(LEV)                                              
15    CONTINUE                                                          
C     STEP T, Q IN TIME.                                                
C     MULTIPLY BM BY SAV   PUT IN SBV                                   
      CALL MATVEC (BM, SAV, SBV)                                        
      DO 20 LEV=1, 12                                                   
C       Y(IN,LEV) = TV(LEV) + DT2 * (YV(LEV) + SBV(LEV))                
      Y(IN,LEV) = CMPLX (                                               
     1   ( REAL(YV(LEV)) +  REAL(SBV(LEV))) * DT2 +  REAL(TV(LEV)) ,    
     2   (AIMAG(YV(LEV)) + AIMAG(SBV(LEV))) * DT2 + AIMAG(TV(LEV)) )    
20    CONTINUE                                                          
C     MULTIPLY S TRANSPOSE BY SAV                                       
      SUM = CMPLX (0., 0.)                                              
      DO 21 LEV=1, 12                                                   
C       SUM = SUM + SV(LEV) * SAV(LEV)                                  
      SUM = CMPLX ( REAL(SUM) + SV(LEV)* REAL(SAV(LEV)),                
     1             AIMAG(SUM) + SV(LEV)*AIMAG(SAV(LEV)))                
21    CONTINUE                                                          
C       Z(IN) = Q(IN) + DT2 * (Z(IN) + SUM)                             
      Z(IN) = CMPLX (                                                   
     1    REAL(Q(IN)) + ( REAL(Z(IN)) +  REAL(SUM)) * DT2 ,             
     2   AIMAG(Q(IN)) + (AIMAG(Z(IN)) + AIMAG(SUM)) * DT2 )             
1000  CONTINUE                                                          
CX20816 CALL DOOUTX( 816 )                                              
CY20816 CALL DOOUTY( 816 )                                              
CX10817 CALL DOINX( 817 )                                               
CY10817 CALL DOINY( 817 )                                               
C     COMPUTE VORTICITY EXPLICITLY                                      
      DO 30 LEV=1, 12                                                   
      DO 30 I=1, 961                                                    
C       W(I,LEV) = ZE(I,LEV) + W(I,LEV)*DT2                             
      W(I,LEV) = CMPLX( REAL(ZE(I,LEV)) +  REAL(W(I,LEV))*DT2,          
     1                 AIMAG(ZE(I,LEV)) + AIMAG(W(I,LEV))*DT2)          
30    CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20818 CALL DOOUTX( 818 )                                              
CY20818 CALL DOOUTY( 818 )                                              
CX40094 CALL SBOUTX (94)                                                
CY40094 CALL SBOUTY (94)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE SICDKD (D,T,Q,X,Y,Z,AM,BM,CM,TOV,SV,DT,ZE,W,GZ,DISPK)  
C     D, T, Q, ARE SPECTRAL COMPONENTS AT (T-1)                         
C     X, Y, Z ARE SPECTRAL COMPONENTS OF NONLIN. TERMS                  
      DIMENSION D( 961 , 12 ), T( 961 , 12 ),                           
     1  Q( 961 ), X( 961 , 12 ),                                        
     2  Y( 961 , 12 ), Z( 961 ), AM( 12 , 12 ), BM( 12 , 12 ),          
     3  CM( 12 , 12 ), TOV( 12 ), SV( 12 )                              
      DIMENSION ZE( 961 , 12 ), W( 961 , 12 ) ,GZ( 961 )                
      COMPLEX ZE, W                                                     
C     VARIABLES USED ONLY IN SUB.                                       
      DIMENSION RIMDIS( 12 , 12 )                                       
      DIMENSION RIM( 12 , 12 ), DM( 12 , 12 ),                          
     1  EM( 12 , 12 ), LV( 12 ), MV( 12 ), XV( 12 ),                    
     2  YV( 12 ), SAV( 12 ), SBV( 12 ),                                 
     3  RHSV( 12 ), DV( 12 ), TV( 12 )                                  
      DIMENSION TOVRA( 12 )                                             
      COMPLEX T, D, Q, X, Y, Z,GZ                                       
      COMPLEX XV, YV, SAV, SBV, RHSV, DV, TV, SUM                       
CX30095 CALL SBINX (95)                                                 
CY30095 CALL SBINY (95)                                                 
      JCAP =  30                                                        
      LEVS =  12                                                        
C       DT2 = 2. * DT                                                   
      DT2 = DT + DT                                                     
CX10819 CALL DOINX( 819 )                                               
CY10819 CALL DOINY( 819 )                                               
C     TOV = R*TO/Q**2                                                   
C     SET UNIT MATRIX RIM                                               
      DO 1 J=1,LEVS                                                     
C       TOVRA(J)= 287.05 *TOV(J)/( 6370000. * 6370000. )                
      TOVRA(J)=( 287.05 /( 6370000. * 6370000. ))*TOV(J)                
      DO 1 I=1,LEVS                                                     
      RIMDIS(I,J)=0.                                                    
      RIM(I,J) = 0.                                                     
1     CONTINUE                                                          
CX20820 CALL DOOUTX( 820 )                                              
CY20820 CALL DOOUTY( 820 )                                              
CX10821 CALL DOINX( 821 )                                               
CY10821 CALL DOINY( 821 )                                               
      DO 2 I=1,LEVS                                                     
      RIM(I,I) = 1.                                                     
      RIMDIS(I,I)=DT2*DISPK +1.                                         
2     CONTINUE                                                          
CX20822 CALL DOOUTX( 822 )                                              
CY20822 CALL DOOUTY( 822 )                                              
CCC   1000 IS (L,N) SPECTRAL COMPONENT LOOP                             
      LL1=1                                                             
      JMN=2* 30 +1                                                      
CX10823 CALL DOINX( 823 )                                               
CY10823 CALL DOINY( 823 )                                               
      DO 1000 NN=1,JMN                                                  
      N=NN-1                                                            
C     COMPUTE MATRICES DM, EM FROM CM                                   
      DO 4 J=1,LEVS                                                     
      DO 4 I=1,LEVS                                                     
      TEMP = N * (N+1.) * CM(I,J)                                       
C     DM(I,J) = RIM(I,J) - TEMP                                         
      DM(I,J)=RIMDIS(I,J)-TEMP                                          
      EM(I,J) = RIM(I,J) + TEMP                                         
4     CONTINUE                                                          
      CALL IMINV (DM, LEVS, DET, LV, MV)                                
C     DM NOW CONTAINS ITS INVERSE.                                      
      LL2=NN                                                            
      IF(N- 30 )900,900,800                                             
800   LL1=LL1+1                                                         
      LL2= 31                                                           
900   CONTINUE                                                          
      IUP=0                                                             
      DO 1000 LL=LL1,LL2                                                
      L=LL2-LL                                                          
      IUP=IUP+1                                                         
      IF(N.GT. 30 )L=LL2-IUP                                            
      IN =L* 31 +LL                                                     
C     SET COLUMN VECTORS XV, YV, DV, TV                                 
      DO 5 LEV=1,LEVS                                                   
      DV(LEV) = D(IN,LEV)                                               
      TV(LEV) = T(IN,LEV)                                               
      XV(LEV) = X(IN,LEV)                                               
      YV(LEV) = Y(IN,LEV)                                               
5     CONTINUE                                                          
C     BEGIN CONSTRUCTING RHSV  (FROM R TO L)                            
      DO 6 LEV=1,LEVS                                                   
C       RHSV(LEV)=(DT*Z(IN)+Q(IN))*TOVRA(LEV)                           
      RHSV(LEV) = CMPLX((REAL(Z(IN))*DT + REAL(Q(IN)))*TOVRA(LEV),      
     1                  (AIMAG(Z(IN))*DT + AIMAG(Q(IN)))*TOVRA(LEV))    
6     CONTINUE                                                          
      DO 7 LEV=1,LEVS                                                   
C      SAV(LEV) = DT * YV(LEV) + TV(LEV)                                
      SAV(LEV) = CMPLX(DT*REAL(YV(LEV)) + REAL(TV(LEV)),                
     1                 DT*AIMAG(YV(LEV)) + AIMAG(TV(LEV)))              
7     CONTINUE                                                          
C     MULTIPLY MATRIX AM BY SAV   PUT IN SBV                            
      CALL MATVEC (AM, SAV, SBV)                                        
      DO 10 LEV=1,LEVS                                                  
C       RHSV(LEV) = N * (N+1.) * (RHSV(LEV) + SBV(LEV))                 
      RHSV(LEV) = CMPLX((REAL(RHSV(LEV)) + REAL(SBV(LEV)))*(N*(N+1)),   
     1                  (AIMAG(RHSV(LEV)) + AIMAG(SBV(LEV)))*(N*(N+1.)))
10    CONTINUE                                                          
      DO 11 LEV=1,LEVS                                                  
C       RHSV(LEV) = DT2 * (XV(LEV) + RHSV(LEV))                         
      RHSV(LEV) = CMPLX((REAL(XV(LEV)) + REAL(RHSV(LEV)))*DT2,          
     1                   (AIMAG(XV(LEV)) + AIMAG(RHSV(LEV)))*DT2)       
11    CONTINUE                                                          
C     MULTIPLY EM MATRIX BY DV(T-1)  PUT IN SAV                         
      CALL MATVEC (EM, DV, SAV)                                         
      DO 14 LEV=1,LEVS                                                  
      RHSV(LEV) = RHSV(LEV) + SAV(LEV)                                  
14    CONTINUE                                                          
C     MULTIPLY RHSV BY DM INVERSE  (NOW IN DM)                          
      CALL MATVEC (DM, RHSV, SAV)                                       
C     SAV IS NOW D VECTOR AT T+1                                        
C     STORE D AT T+1 IN X                                               
      DO 15 LEV=1,LEVS                                                  
      X(IN,LEV) = SAV(LEV)                                              
15    CONTINUE                                                          
C     STEP T, Q IN TIME.                                                
C     SET  SAV = D(T+1) + D(T-1)                                        
      DO 16 LEV=1,LEVS                                                  
      SAV(LEV) = SAV(LEV) + DV(LEV)                                     
16    CONTINUE                                                          
C     MULTIPLY BM BY SAV  PUT IN SBV                                    
      CALL MATVEC (BM, SAV, SBV)                                        
      DO 20 LEV=1,LEVS                                                  
C       Y(IN,LEV) = TV(LEV) + DT * (2. * YV(LEV) + SBV(LEV))            
      Y(IN,LEV) = CMPLX((REAL(YV(LEV))+REAL(YV(LEV))+REAL(SBV(LEV)))*DT+
     1                        REAL(TV(LEV)),                            
     2                  (AIMAG(YV(LEV))+AIMAG(YV(LEV))+AIMAG(SBV(LEV)))*
     3                         DT + AIMAG(TV(LEV)))                     
20    CONTINUE                                                          
C     MULTIPLY S TRANSPOSE BY SAV                                       
      SUM = CMPLX(0.,0.)                                                
      DO 21 LEV=1,LEVS                                                  
C       SUM = SUM + SV(LEV) * SAV(LEV)                                  
      SUM = CMPLX(REAL(SUM) + SV(LEV)*REAL(SAV(LEV)),                   
     1            AIMAG(SUM) + SV(LEV)*AIMAG(SAV(LEV)))                 
21    CONTINUE                                                          
C       Z(IN) = Q(IN) + DT * (2. * Z(IN) + SUM)                         
      Z(IN) = CMPLX(REAL(Q(IN))+(REAL(Z(IN))+REAL(Z(IN))+REAL(SUM))*DT, 
     1          AIMAG(Q(IN))+(AIMAG(Z(IN))+AIMAG(Z(IN))+AIMAG(SUM))*DT) 
1000  CONTINUE                                                          
CX20824 CALL DOOUTX( 824 )                                              
CY20824 CALL DOOUTY( 824 )                                              
CX10825 CALL DOINX( 825 )                                               
CY10825 CALL DOINY( 825 )                                               
C     COMPUTE VORTICITY EXPLICITLY                                      
      DO 30 LEV=1, 12                                                   
      DO 30 I=1, 961                                                    
C       W(I,LEV) = ZE(I,LEV) + 2.*DT*W(I,LEV)                           
      W(I,LEV) = CMPLX(REAL(ZE(I,LEV)) + REAL(W(I,LEV))*DT2,            
     1 AIMAG(ZE(I,LEV)) + AIMAG(W(I,LEV))*DT2)                          
30    CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20826 CALL DOOUTX( 826 )                                              
CY20826 CALL DOOUTY( 826 )                                              
CX40096 CALL SBOUTX (96)                                                
CY40096 CALL SBOUTY (96)                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE STIMER                                                 
CX30097 CALL SBINX (97)                                                 
CY30097 CALL SBINY (97)                                                 
CX40098 CALL SBOUTX (98)                                                
CY40098 CALL SBOUTY (98)                                                
C***   FAKE I.B.M. TIMER                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE SUMPLS(FLN,AP,AM,PLN)                                  
C     SR COMPUTES FOURIER MULTIPLIERS                                   
C     IN AP COEFFS. OF SYMMETRIC LEGENDRES                              
C     IN AM COEFFS. OF ANTISYM.  LEGENDRES                              
      DIMENSION AP( 31 ),AM( 31 ),                                      
     1PLN( 32 , 31 ),FLN( 31 , 31 )                                     
      COMPLEX AP,AM,FLN                                                 
CX30099 CALL SBINX (99)                                                 
CY30099 CALL SBINY (99)                                                 
CX10827 CALL DOINX( 827 )                                               
CY10827 CALL DOINY( 827 )                                               
      DO 3 L=1, 31                                                      
      CRP=0.                                                            
      CIP=0.                                                            
      DO 1 I=1, 31 ,2                                                   
      CRP=CRP+PLN(I,L)* REAL(FLN(I,L))                                  
      CIP=CIP+PLN(I,L)*AIMAG(FLN(I,L))                                  
1     CONTINUE                                                          
      CRM=0.                                                            
      CIM=0.                                                            
      DO 2 I=2, 30  ,2                                                  
      CRM=CRM+PLN(I,L)* REAL(FLN(I,L))                                  
      CIM=CIM+PLN(I,L)*AIMAG(FLN(I,L))                                  
2     CONTINUE                                                          
      AP(L)=CMPLX(CRP+CRM,CIP+CIM)                                      
      AM(L)=CMPLX(CRP-CRM,CIP-CIM)                                      
3     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20828 CALL DOOUTX( 828 )                                              
CY20828 CALL DOOUTY( 828 )                                              
CX40100 CALL SBOUTX (100)                                               
CY40100 CALL SBOUTY (100)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE SUMPLV(VLN,AP,AM,PLN)                                  
C     SR COMPUTES FOURIER MULTIPLIERS                                   
C     IN AP COEFFS. OF SYMMETRIC LEGENDRES                              
C     IN AM COEFFS. OF ANTISYM.  LEGENDRES                              
      DIMENSION AP( 31 ),AM( 31 ),                                      
     1PLN( 32 , 31 ),VLN( 32 , 31 )                                     
      COMPLEX AP,AM,VLN                                                 
CX30101 CALL SBINX (101)                                                
CY30101 CALL SBINY (101)                                                
CX10829 CALL DOINX( 829 )                                               
CY10829 CALL DOINY( 829 )                                               
      DO 3 L=1, 31                                                      
      CRP=0.                                                            
      CIP=0.                                                            
      DO 1 I=1, 31 ,2                                                   
      CRP=CRP+PLN(I,L)* REAL(VLN(I,L))                                  
      CIP=CIP+PLN(I,L)*AIMAG(VLN(I,L))                                  
1     CONTINUE                                                          
      CRM=0.                                                            
      CIM=0.                                                            
      DO 2 I=2, 32 ,2                                                   
      CRM=CRM+PLN(I,L)* REAL(VLN(I,L))                                  
      CIM=CIM+PLN(I,L)*AIMAG(VLN(I,L))                                  
2     CONTINUE                                                          
      AP(L)=CMPLX(CRP+CRM,CIP+CIM)                                      
      AM(L)=CMPLX(CRP-CRM,CIP-CIM)                                      
3     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20830 CALL DOOUTX( 830 )                                              
CY20830 CALL DOOUTY( 830 )                                              
CX40102 CALL SBOUTX (102)                                               
CY40102 CALL SBOUTY (102)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE SYMASY (A,B)                                           
      DIMENSION A( 31 ),B( 31 )                                         
      COMPLEX A,B                                                       
CX30103 CALL SBINX (103)                                                
CY30103 CALL SBINY (103)                                                
CX10831 CALL DOINX( 831 )                                               
CY10831 CALL DOINY( 831 )                                               
C     ON INPUT A CORRESP. TO N. HEMIS. FOURIERS                         
C     ON INPUT B CORRESP. TO S. HEMIS. FOURIERS                         
C     ON OUTPUT A CORRESP. TO MULTIPLIERS OF SYM. PLN                   
C     ON OUTPUT B CORRESP. TO MULTIPLIERS OF ANTISYM. PLN               
      DO 1 I = 1, 31                                                    
      CR= REAL(A(I))                                                    
      CI=AIMAG(A(I))                                                    
      A(I)=CMPLX(CR+REAL(B(I)),CI+AIMAG(B(I)))                          
      B(I)=CMPLX(CR-REAL(B(I)),CI-AIMAG(B(I)))                          
1     CONTINUE                                                          
C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC.                 
CX20832 CALL DOOUTX( 832 )                                              
CY20832 CALL DOOUTY( 832 )                                              
CX40104 CALL SBOUTX (104)                                               
CY40104 CALL SBOUTY (104)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE TEND(Z,X,Y,QM,DIM,TEM,A,TO,B,DEL)                      
      COMPLEX Z( 31 , 31 ),QM( 31 , 31 ),DD,                            
     1X( 31 , 31 , 12 ),Y( 31 , 31 , 12 ),                              
     2DIM( 31 , 31 , 12 ),TEM( 31 , 31 , 12 )                           
      DIMENSION A( 12 , 12 ),B( 12 , 12 ),                              
     1TO( 12 ),DEL( 12 )                                                
CX30105 CALL SBINX (105)                                                
CY30105 CALL SBINY (105)                                                
CX10833 CALL DOINX( 833 )                                               
CY10833 CALL DOINY( 833 )                                               
      DO 10 LL=1, 31                                                    
      N=LL-2                                                            
      DO 10 I=1, 31                                                     
      N=N+1                                                             
      DO 15 K=1, 12                                                     
      F=-DEL(K)                                                         
      DD=DIM(LL,I,K)                                                    
      Z(LL,I)=Z(LL,I)+CMPLX(F*REAL(DD),F*AIMAG(DD))                     
      F=N*(N+1.)*TO(K)* 287.05 /( 6370000. * 6370000. )                 
      DD=QM(LL,I)                                                       
      X(LL,I,K)=X(LL,I,K)+CMPLX(F*REAL(DD),F*AIMAG(DD))                 
      DO 20 KV=1, 12                                                    
      F=N*(N+1.)*A(K,KV)                                                
      DD=TEM(LL,I,KV)                                                   
      X(LL,I,K)=X(LL,I,K)+CMPLX(F*REAL(DD),F*AIMAG(DD))                 
      F=B(K,KV)                                                         
      DD=DIM(LL,I,KV)                                                   
   20 Y(LL,I,K)=Y(LL,I,K)+CMPLX(F*REAL(DD),F*AIMAG(DD))                 
   15 CONTINUE                                                          
   10 CONTINUE                                                          
CX20834 CALL DOOUTX( 834 )                                              
CY20834 CALL DOOUTY( 834 )                                              
CX40106 CALL SBOUTX (106)                                               
CY40106 CALL SBOUTY (106)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE TTIMER(T)                                              
CX30107 CALL SBINX (107)                                                
CY30107 CALL SBINY (107)                                                
CX40108 CALL SBOUTX (108)                                               
CY40108 CALL SBOUTY (108)                                               
C***   FAKE I.B.M. TIMER                                                
      RETURN                                                            
      END                                                               
      SUBROUTINE UVGLOB(PLNWCS,DER,DI,ZE,UP,UM,VP,VM)                   
      DIMENSION PLNWCS( 31 , 31 ),                                      
     1             DER( 31 , 31 )                                       
      COMPLEX   DI( 31 , 31 )                                           
      COMPLEX   ZE( 31 , 31 )                                           
      COMPLEX   UP( 31 ),UM( 31 )                                       
      COMPLEX   VP( 31 ),VM( 31 )                                       
CX30109 CALL SBINX (109)                                                
CY30109 CALL SBINY (109)                                                
CX10835 CALL DOINX( 835 )                                               
CY10835 CALL DOINY( 835 )                                               
      DO 50 LL=1, 31                                                    
      CRP=0.                                                            
      CIP=0.                                                            
      DO 10 I=1, 31 ,2                                                  
      CRP=CRP+PLNWCS(I,LL)*AIMAG(DI(I,LL))                              
      CIP=CIP+PLNWCS(I,LL)* REAL(DI(I,LL))                              
10    CONTINUE                                                          
      CIP=-CIP                                                          
      CRM=0.                                                            
      CIM=0.                                                            
      DO 20 I=2, 30 ,2                                                  
      CRM=CRM+PLNWCS(I,LL)*AIMAG(DI(I,LL))                              
      CIM=CIM+PLNWCS(I,LL)* REAL(DI(I,LL))                              
20    CONTINUE                                                          
      CIM=-CIM                                                          
      DO 30 I=1, 31 ,2                                                  
      CRM=CRM+DER(I,LL) *REAL(ZE(I,LL))                                 
      CIM=CIM+DER(I,LL)*AIMAG(ZE(I,LL))                                 
30    CONTINUE                                                          
      DO 40 I=2, 30 ,2                                                  
      CRP=CRP+DER(I,LL) *REAL(ZE(I,LL))                                 
      CIP=CIP+DER(I,LL)*AIMAG(ZE(I,LL))                                 
40    CONTINUE                                                          
      UP(LL)=CMPLX(CRP+CRM,CIP+CIM)                                     
      UM(LL)=CMPLX(CRP-CRM,CIP-CIM)                                     
      CRP=0.                                                            
      CIP=0.                                                            
      DO 11 I=1, 31 ,2                                                  
      CRP=CRP+PLNWCS(I,LL)*AIMAG(ZE(I,LL))                              
      CIP=CIP+PLNWCS(I,LL)* REAL(ZE(I,LL))                              
11    CONTINUE                                                          
      CIP=-CIP                                                          
      CRM=0.                                                            
      CIM=0.                                                            
      DO 21 I=2, 30 ,2                                                  
      CRM=CRM+PLNWCS(I,LL)*AIMAG(ZE(I,LL))                              
      CIM=CIM+PLNWCS(I,LL)* REAL(ZE(I,LL))                              
21    CONTINUE                                                          
      CIM=-CIM                                                          
      DO 31 I=1, 31 ,2                                                  
      CRM=CRM-DER(I,LL) *REAL(DI(I,LL))                                 
      CIM=CIM-DER(I,LL)*AIMAG(DI(I,LL))                                 
31    CONTINUE                                                          
      DO 41 I=2, 30 ,2                                                  
      CRP=CRP-DER(I,LL) *REAL(DI(I,LL))                                 
      CIP=CIP-DER(I,LL)*AIMAG(DI(I,LL))                                 
41    CONTINUE                                                          
      VP(LL)=CMPLX(CRP+CRM,CIP+CIM)                                     
      VM(LL)=CMPLX(CRP-CRM,CIP-CIM)                                     
50    CONTINUE                                                          
CX20836 CALL DOOUTX( 836 )                                              
CY20836 CALL DOOUTY( 836 )                                              
      UP(1)=CMPLX(REAL(UP(1)),0.)                                       
      UM(1)=CMPLX(REAL(UM(1)),0.)                                       
      VP(1)=CMPLX(REAL(VP(1)),0.)                                       
      VM(1)=CMPLX(REAL(VM(1)),0.)                                       
CX40110 CALL SBOUTX (110)                                               
CY40110 CALL SBOUTY (110)                                               
C***  UP ARE FOURIERS FOR N. HEMI.  UM FOR S. HEMI                      
      RETURN                                                            
      END                                                               
      SUBROUTINE VERIFY(U,ITERS,SLAT)                                   
C**   THIS SR WILL PRINT THE MAXIMUM VALUES OF U ARRAY AT EVERY LEVEL   
      DIMENSION U(96,12),UEXACT(12)                                     
CX30111 CALL SBINX (111)                                                
CY30111 CALL SBINY (111)                                                
      LAT=INT(ACOS(SLAT)*180./3.141592653589793)                        
      IF (LAT.EQ.37) THEN                                               
                                                                        
      IFLAG= 0                                                          
      UEXACT(1) =   0.32882E-03                                         
      UEXACT(2) =   0.30024E-03                                         
      UEXACT(3) =   0.23677E-03                                         
      UEXACT(4) =   0.17480E-03                                         
      UEXACT(5) =   0.14639E-03                                         
      UEXACT(6) =  -0.42280E-03                                         
      UEXACT(7) =  -0.88541E-03                                         
      UEXACT(8) =  -0.72669E-03                                         
      UEXACT(9) =  -0.34500E-03                                         
      UEXACT(10) =  0.49664E-04                                         
      UEXACT(11) =  0.31366E-03                                         
      UEXACT(12) =  0.49537E-04                                         
                                                                        
      WRITE(4,1)                                                        
1     FORMAT(1X,
     * 'OUTPUT FOR PERFECT CLUB BENCHMARK: SPEC77 INSTRUMENTED VERSION'/ 
     *       1X,'$Revision: 1.A.1 $ $Author: psharma $'/)     
      WRITE(4,2)                                                        
2     FORMAT(1X,'VALIDATION PARAMETERS:'/)                              
      WRITE(4,5) (U(48,I),I=1,12)                                       
CX10837 CALL DOINX( 837 )                                               
CY10837 CALL DOINY( 837 )                                               
      DO 599 I=1,12                                                     
C     IF (ABS((APPROX-EXACT)/EXACT).GE.1E-2)                            
      IF (ABS((U(48,I)-UEXACT(I))/UEXACT(I)).GT.1E-2) IFLAG = 1         
599   CONTINUE                                                          
CX20838 CALL DOOUTX( 838 )                                              
CY20838 CALL DOOUTY( 838 )                                              
C                                                                       
      IF( IFLAG .EQ. 0) THEN                                            
        WRITE(4,3)                                                      
3       FORMAT(//1X,'RESULTS FOR THIS RUN ARE:   VALID')                
      ELSE                                                              
        WRITE(4,4)                                                      
4       FORMAT(//1X,'RESULTS FOR THIS RUN ARE: INVALID')                
      ENDIF                                                             
      ENDIF                                                             
  5   FORMAT(5E15.5)                                                    
CX40112 CALL SBOUTX (112)                                               
CY40112 CALL SBOUTY (112)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINEVERTIC(F,EIGG,EIGGT,DOTPRO,MODS,INOUT)                  
      DIMENSION EIGG( 12 , 12 ),EIGGT( 12 , 12 ),DOTPRO( 12 )           
      COMPLEX   F( 961 , 12 ),COL( 12 ),SUM                             
CX30113 CALL SBINX (113)                                                
CY30113 CALL SBINY (113)                                                
C     INOUT=-1 TO OBTAIN VERTICAL MODE EXPANSION                        
C     INOUT=+1 TO OBTAIN SPET. COEFS. FROM VERTICAL EXPANSION           
      IF(INOUT)10,20,30                                                 
10    CONTINUE                                                          
CX10839 CALL DOINX( 839 )                                               
CY10839 CALL DOINY( 839 )                                               
      DO 4 I=1, 961                                                     
      DO 2 MOD=1,MODS                                                   
      SUM=(0.,0.)                                                       
      DO 1 LEV=1, 12                                                    
      SUM=SUM+CMPLX(EIGGT(LEV,MOD)* REAL(F(I,LEV)),                     
     1              EIGGT(LEV,MOD)*AIMAG(F(I,LEV)))                     
1     CONTINUE                                                          
2     COL(MOD)=CMPLX(DOTPRO(MOD)*REAL(SUM),DOTPRO(MOD)*AIMAG(SUM))      
      DO 3 MOD=1,MODS                                                   
3     F(I,MOD)=COL(MOD)                                                 
4     CONTINUE                                                          
CX20840 CALL DOOUTX( 840 )                                              
CY20840 CALL DOOUTY( 840 )                                              
   20 CONTINUE                                                          
CX40114 CALL SBOUTX (114)                                               
CY40114 CALL SBOUTY (114)                                               
      RETURN                                                            
30    CONTINUE                                                          
CX10841 CALL DOINX( 841 )                                               
CY10841 CALL DOINY( 841 )                                               
      DO 8 I=1, 961                                                     
      DO 5 MOD=1,MODS                                                   
5     COL(MOD)=F(I,MOD)                                                 
      DO 7 LEV=1, 12                                                    
      F(I,LEV)=(0.,0.)                                                  
      DO 6 MOD=1,MODS                                                   
      F(I,LEV)=F(I,LEV)+CMPLX(                                          
     1EIGG(LEV,MOD)*REAL(COL(MOD)),EIGG(LEV,MOD)*AIMAG(COL(MOD)))       
6     CONTINUE                                                          
7     CONTINUE                                                          
8     CONTINUE                                                          
CX20842 CALL DOOUTX( 842 )                                              
CY20842 CALL DOOUTY( 842 )                                              
CX40114 CALL SBOUTX (114)                                               
CY40114 CALL SBOUTY (114)                                               
      RETURN                                                            
      END                                                               
      LOGICAL FUNCTION W3AQ09(N)                                        
CX30115 CALL SBINX (115)                                                
CY30115 CALL SBINY (115)                                                
      N=0                                                               
      W3AQ09=.FALSE.                                                    
CX40116 CALL SBOUTX (116)                                               
CY40116 CALL SBOUTY (116)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE W3EXIT(R,I)                                            
CX30117 CALL SBINX (117)                                                
CY30117 CALL SBINY (117)                                                
C     PRINT 100                                                         
      WRITE(6,100)                                                      
100   FORMAT(1H0,'FAKE OPERATIONAL W3EXIT')                             
CX40118 CALL SBOUTX (118)                                               
CY40118 CALL SBOUTY (118)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE W3FS03(I,J,K,L,M)                                      
CX30119 CALL SBINX (119)                                                
CY30119 CALL SBINY (119)                                                
C     PRINT 100                                                         
      WRITE(6,100)                                                      
100   FORMAT(1H0,'FAKE OPERATIONAL W3FS03')                             
CX40120 CALL SBOUTX (120)                                               
CY40120 CALL SBOUTY (120)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE W3STOP                                                 
CX30121 CALL SBINX (121)                                                
CY30121 CALL SBINY (121)                                                
C     PRINT 100                                                         
      WRITE(6,100)                                                      
100   FORMAT(1H0,'FAKE OPERATIONAL W3STOP')                             
CX40122 CALL SBOUTX (122)                                               
CY40122 CALL SBOUTY (122)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE W3STRT                                                 
CX30123 CALL SBINX (123)                                                
CY30123 CALL SBINY (123)                                                
C     PRINT 100                                                         
      WRITE(6,100)                                                      
100   FORMAT(1H0,'FAKE OPERATIONAL W3STRT')                             
CX40124 CALL SBOUTX (124)                                               
CY40124 CALL SBOUTY (124)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE XMOVEX(A,B,N)                                          
C***   FAKE 360/195 BIT MOVING ROUTINE                                  
      COMPLEX A( 961 ),B( 961 )                                         
CX30125 CALL SBINX (125)                                                
CY30125 CALL SBINY (125)                                                
CX10843 CALL DOINX( 843 )                                               
CY10843 CALL DOINY( 843 )                                               
      DO 1 I=1, 961                                                     
1     A(I)=B(I)                                                         
CX20844 CALL DOOUTX( 844 )                                              
CY20844 CALL DOOUTY( 844 )                                              
CX40126 CALL SBOUTX (126)                                               
CY40126 CALL SBOUTY (126)                                               
      RETURN                                                            
      END                                                               
      SUBROUTINE XSTORE(A,X,N)                                          
C***   FAKE 360/195 BYTE MOVING ROUTINE                                 
      DIMENSION A(1)                                                    
CX30127 CALL SBINX (127)                                                
CY30127 CALL SBINY (127)                                                
CX10845 CALL DOINX( 845 )                                               
CY10845 CALL DOINY( 845 )                                               
      DO 1 I=1,N                                                        
1     A(I)=X                                                            
CX20846 CALL DOOUTX( 846 )                                              
CY20846 CALL DOOUTY( 846 )                                              
CX40128 CALL SBOUTX (128)                                               
CY40128 CALL SBOUTY (128)                                               
      RETURN                                                            
      END                                                               
