      SUBROUTINE IMPINT(DT,ALPH)                                                                                                    
                                                                                                                                    
      INCLUDE 'param.h'                                                                                                             
      INCLUDE 'param1.h'                                                                                                            
      INCLUDE 'const1.h'                                                                                                            
      INCLUDE 'const2.h'                                                                                                            
                                                                                                                                    
      DIMENSION DSUM(KX),INDX(KX),YA(KX,KX)                                                                                         
                                                                                                                                    
      XI=DT*ALPH                                                                                                                    
      XXI = XI/(A*A)                                                                                                                
                                                                                                                                    
      DO 201 K=1,KX                                                                                                                 
        DHSX(K)=XI*DHS(K)                                                                                                           
  201 CONTINUE                                                                                                                      
                                                                                                                                    
      do 60 n=1,jx                                                                                                                  
      do 60 m=1,mx                                                                                                                  
        MN=M+N                                                                                                                      
        MM=ISC*(M-1)+1                                                                                                              
        LL=MM+N-2                                                                                                                   
        elz(m,n)=FLOAT(LL)*FLOAT(LL+1)*XXI                                                                                          
   60 continue                                                                                                                      
                                                                                                                                    
c  define reference atmosphere, function of height only                                                                             
                                                                                                                                    
      DO 101 K=1,KX                                                                                                                 
        TREF(K)=255.                                                                                                                
        TREF1(K)=RGAS*TREF(K)                                                                                                       
  101 CONTINUE                                                                                                                      
                                                                                                                                    
c   T(K) = TEX(K)+YA(K,K')*D(K') + XA(K,K')*SIG(K')                                                                                 
                                                                                                                                    
      DO 1 K=1,KX                                                                                                                   
      DO 1 K1=1,KXM                                                                                                                 
        XA(K,K1)=0.                                                                                                                 
    1 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 2 K=1,KX                                                                                                                   
      DO 2 K1=1,KX                                                                                                                  
        YA(K,K1)=-AKAP*TREF(K)*DHS(K1)                                                                                              
    2 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 3 K=2,KX                                                                                                                   
        XA(K,K-1)=0.5*(AKAP*TREF(K)/FSG(K)                                                                                          
     *      -(TREF(K)-TREF(K-1))/DHS(K))                                                                                            
    3 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 4 K=1,KXM                                                                                                                  
        XA(K,K)=0.5*(AKAP*TREF(K)/FSG(K)                                                                                            
     *    -(TREF(K+1)-TREF(K))/DHS(K))                                                                                              
    4 CONTINUE                                                                                                                      
                                                                                                                                    
c     SIG(K)=XB(K,K')*D(K')                                                                                                         
                                                                                                                                    
      DSUM(1)=DHS(1)                                                                                                                
      DO 6 K=2,KX                                                                                                                   
        DSUM(K)=DSUM(K-1)+DHS(K)                                                                                                    
    6 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 5 K=1,KXM                                                                                                                  
      DO 5 K1=1,KX                                                                                                                  
        XB(K,K1)=DHS(K1)*DSUM(K)                                                                                                    
        IF(K1.LE.K) XB(K,K1)=XB(K,K1)-DHS(K1)                                                                                       
    5 CONTINUE                                                                                                                      
                                                                                                                                    
c     T(K)=TEX(K)+XC(K,K')*D(K')                                                                                                    
                                                                                                                                    
      DO 7 K=1,KX                                                                                                                   
      DO 7 K1=1,KX                                                                                                                  
        XC(K,K1)=YA(K,K1)                                                                                                           
        DO 77 K2=1,KXM                                                                                                              
          XC(K,K1)=XC(K,K1)+XA(K,K2)*XB(K2,K1)                                                                                      
   77   CONTINUE                                                                                                                    
    7 CONTINUE                                                                                                                      
                                                                                                                                    
c      P(K)=XD(K,K')*T(K')                                                                                                          
                                                                                                                                    
      DO 8 K=1,KX                                                                                                                   
      DO 8 K1=1,KX                                                                                                                  
        XD(K,K1)=0.                                                                                                                 
    8 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 91 K=1,KX                                                                                                                  
      DO 91 K1=K+1,KX                                                                                                               
        XD(K,K1)=RGAS*ALOG(HSG(K1+1)/HSG(K1))                                                                                       
   91 CONTINUE                                                                                                                      
      DO 92 K=1,KX                                                                                                                  
        XD(K,K)=RGAS*ALOG(HSG(K+1)/FSG(K))                                                                                          
   92 CONTINUE                                                                                                                      
                                                                                                                                    
c      P(K)=YE(K)+XE(K,K')*D(K')                                                                                                    
                                                                                                                                    
      DO 10 K=1,KX                                                                                                                  
      DO 10 K1=1,KX                                                                                                                 
        XE(K,K1)=0.                                                                                                                 
      DO 10 K2=1,KX                                                                                                                 
        XE(K,K1)=XE(K,K1)+XD(K,K2)*XC(K2,K1)                                                                                        
   10 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 11 L=1,LMAX                                                                                                                
        XXX=(FLOAT(L)*FLOAT(L+1))/(A*A)                                                                                             
        DO 12 K=1,KX                                                                                                                
        DO 12 K1=1,KX                                                                                                               
          XF(K,K1,L)=XI*XI*XXX*(RGAS*TREF(K)*DHS(K1)-XE(K,K1))                                                                      
   12   CONTINUE                                                                                                                    
        DO 13 K=1,KX                                                                                                                
          XF(K,K,L)=XF(K,K,L)+1.                                                                                                    
   13   CONTINUE                                                                                                                    
   11 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 30 L=1,LMAX                                                                                                                
        CALL INV(XF(1,1,L),XJ(1,1,L),INDX,KX)                                                                                       
   30 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 50 K=1,KX                                                                                                                  
      DO 50 K1=1,KX                                                                                                                 
        XC(K,K1)=XC(K,K1)*XI                                                                                                        
   50 CONTINUE                                                                                                                      
                                                                                                                                    
      RETURN                                                                                                                        
      END                                                                                                                           
****************************************************************                                                                    
      SUBROUTINE IMPLIC                                                                                                             
                                                                                                                                    
      INCLUDE 'param.h'                                                                                                             
      INCLUDE 'param1.h'                                                                                                            
      PARAMETER (MXJXKX=MX*JX*KX)                                                                                                   
      INCLUDE 'const1.h'                                                                                                            
      INCLUDE 'const2.h'                                                                                                            
                                                                                                                                    
      COMPLEX YE(MX,JX,KX),YF(MX,JX,KX),ZERO                                                                                        
                                                                                                                                    
      ZERO=(0.,0.)                                                                                                                  
                                                                                                                                    
      DO 1 K=1,KX                                                                                                                   
      DO 1 N=1,JX                                                                                                                   
      DO 1 M=1,MX                                                                                                                   
        YE(M,N,K)=ZERO                                                                                                              
    1 CONTINUE                                                                                                                      
                                                                                                                                    
c      DO 1 M=1,MXJXKX                                                                                                              
c        YE(M,1,1)=ZERO                                                                                                             
c    1 CONTINUE                                                                                                                     
                                                                                                                                    
      DO 2 K1=1,KX                                                                                                                  
      DO 2 K=1,KX                                                                                                                   
        DO 11 N=1,JX                                                                                                                
        DO 11 M=1,MX                                                                                                                
          YE(M,N,K)=YE(M,N,K)+XD(K,K1)*TDT(M,N,K1)                                                                                  
   11   CONTINUE                                                                                                                    
c        DO 11 M=1,MXJX                                                                                                             
c          YE(M,1,K)=YE(M,1,K)+XD(K,K1)*TDT(M,1,K1)                                                                                 
c   11   CONTINUE                                                                                                                   
    2 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 21 K=1,KX                                                                                                                  
      DO 22 N=1,JX                                                                                                                  
      DO 22 M=1,MX                                                                                                                  
        YE(M,N,K)=YE(M,N,K)+TREF1(K)*PSDT(M,N)                                                                                      
   22 CONTINUE                                                                                                                      
c      DO 22 M=1,MXJX                                                                                                               
c        YE(M,1,K)=YE(M,1,K)+TREF1(K)*PSDT(M,1)                                                                                     
c   22 CONTINUE                                                                                                                     
   21 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 4 K=1,KX                                                                                                                   
        DO 44 N=1,JX                                                                                                                
        DO 44 M=1,MX                                                                                                                
          YF(M,N,K)=DIVDT(M,N,K)+ELZ(M,N)*YE(M,N,K)                                                                                 
   44   CONTINUE                                                                                                                    
c        DO 44 M=1,MXJX                                                                                                             
c          YF(M,1,K)=DIVDT(M,1,K)+ELZ(M,1)*YE(M,1,K)                                                                                
c   44   CONTINUE                                                                                                                   
    4 CONTINUE                                                                                                                      
                                                                                                                                    
c      DO 5 K=1,KX                                                                                                                  
c      DO 5 N=1,JX                                                                                                                  
c      DO 5 M=1,MX                                                                                                                  
c        DIVDT(M,N,K)=ZERO                                                                                                          
c    5 CONTINUE                                                                                                                     
                                                                                                                                    
      DO 5 M=1,MXJXKX                                                                                                               
        DIVDT(M,1,1)=ZERO                                                                                                           
    5 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 6 N=1,JX                                                                                                                   
      DO 6 M=1,MX                                                                                                                   
        MM=ISC*(M-1)+1                                                                                                              
        LL=MM+N-2                                                                                                                   
        IF(LL.NE.0) THEN                                                                                                            
          DO 66 K1=1,KX                                                                                                             
          DO 66 K=1,KX                                                                                                              
            DIVDT(M,N,K)=DIVDT(M,N,K)+XJ(K,K1,LL)*YF(M,N,K1)                                                                        
   66     CONTINUE                                                                                                                  
        ENDIF                                                                                                                       
    6 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 46 K=1,KX                                                                                                                  
c        DO 47 N=1,JX                                                                                                               
c        DO 47 M=1,MX                                                                                                               
c          PSDT(M,N)=PSDT(M,N)-DIVDT(M,N,K)*DHSX(K)                                                                                 
c   47   CONTINUE                                                                                                                   
        DO 47 M=1,MXJX                                                                                                              
          PSDT(M,1)=PSDT(M,1)-DIVDT(M,1,K)*DHSX(K)                                                                                  
   47   CONTINUE                                                                                                                    
   46 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 7 K=1,KX                                                                                                                   
      DO 7 K1=1,KX                                                                                                                  
C        DO 48 N=1,JX                                                                                                               
C        DO 48 M=1,MX                                                                                                               
C          TDT(M,N,K)=TDT(M,N,K)+XC(K,K1)*DIVDT(M,N,K1)                                                                             
C   48   CONTINUE                                                                                                                   
        DO 48 M=1,MXJX                                                                                                              
          TDT(M,1,K)=TDT(M,1,K)+XC(K,K1)*DIVDT(M,1,K1)                                                                              
   48   CONTINUE                                                                                                                    
    7 CONTINUE                                                                                                                      
                                                                                                                                    
      RETURN                                                                                                                        
      END                                                                                                                           
*******************************************************************                                                                 
      SUBROUTINE TRSUN(JJ)                                                                                                          
                                                                                                                                    
      INCLUDE 'param.h'                                                                                                             
      INCLUDE 'param1.h'                                                                                                            
      INCLUDE 'const1.h'                                                                                                            
                                                                                                                                    
c   writes out some formatted data for graphics on SUN                                                                              
                                                                                                                                    
      DIMENSION SAVE(IX,IL,2)                                                                                                       
      DIMENSION DUMZ(IL,9,10)                                                                                                       
                                                                                                                                    
      DO 1 L=1,10                                                                                                                   
      DO 1 J=1,IL                                                                                                                   
      DO 1 I=1,IX                                                                                                                   
        SAVE(I,J,L)=0.                                                                                                              
    1 CONTINUE                                                                                                                      
                                                                                                                                    
      CALL GRID(T(1,1,KX,JJ),SAVE(1,1,1),1)                                                                                         
      CALL GRID(VOR(1,1,4,JJ),SAVE(1,1,2),1)                                                                                        
      OPEN(90,FILE='sun',FORM='FORMATTED')                                                                                          
      WRITE(90,3000) SAVE                                                                                                           
      CLOSE (90)                                                                                                                    
 3000 FORMAT(10E12.5)                                                                                                               
                                                                                                                                    
      RETURN                                                                                                                        
      END                                                                                                                           
*******************************************************************                                                                 
      SUBROUTINE ANAL(JJ)                                                                                                           
                                                                                                                                    
      INCLUDE 'param.h'                                                                                                             
      INCLUDE 'param1.h'                                                                                                            
      INCLUDE 'const1.h'                                                                                                            
                                                                                                                                    
c  monitors a few selected variables                                                                                                
                                                                                                                                    
      COMPLEX TEMP(MX,JX)                                                                                                           
      DIMENSION EK(KX)                                                                                                              
                                                                                                                                    
      PRINT*,ISTEP                                                                                                                  
                                                                                                                                    
      DO 1 K=1,KX                                                                                                                   
        CALL INVLAP(VOR(1,1,K,JJ),TEMP)                                                                                             
        EK(K)=0.                                                                                                                    
        DO 2 M=2,MX                                                                                                                 
        DO 2 N=1,JX                                                                                                                 
          EK(K)=EK(K)-TEMP(M,N)*CONJG(VOR(M,N,K,JJ))                                                                                
    2   CONTINUE                                                                                                                    
    1 CONTINUE                                                                                                                      
                                                                                                                                    
      WRITE(6,2000) (EK(K),K=1,KX)                                                                                                  
 2000 FORMAT(10(1X,E12.5))                                                                                                          
                                                                                                                                    
      CALL UVSPEC(VOR(1,1,3,JJ),DIV(1,1,3,JJ),                                                                                      
     * DUMC(1,1,1),DUMC(1,1,2))                                                                                                     
                                                                                                                                    
      CALL GRID(DUMC(1,1,1),DUMR(1,1,1),2)                                                                                          
      CALL GRID(DUMC(1,1,2),DUMR(1,1,2),2)                                                                                          
      CALL GRID(T(1,1,8,JJ),DUMR(1,1,3),1)                                                                                          
                                                                                                                                    
      WRITE(6,1000) (DUMR(1,J,1),J=1,40,4)                                                                                          
      WRITE(6,1000) (DUMR(1,J,2),J=1,40,4)                                                                                          
      WRITE(6,1000) (DUMR(1,J,3),J=1,40,4)                                                                                          
 1000 FORMAT(1X,10(1X,E12.5))                                                                                                       
                                                                                                                                    
      RETURN                                                                                                                        
      END                                                                                                                           
*******************************************************************                                                                 
                                                             
      SUBROUTINE CONST                                                                                                              
                                                                                                                                    
      INCLUDE 'param.h'                                                                                                             
      INCLUDE 'param1.h'                                                                                                            
      INCLUDE 'const1.h'                                                                                                            
                                                                                                                                    
      COMPLEX ZERO                                                                                                                  
                                                                                                                                    
      ZERO=(0.,0.)                                                                                                                  
      PI=4.*ATAN(1.)                                                                                                                
      OMEGA=2.*PI/86400.                                                                                                            
      AKAP=2./7.                                                                                                                    
      G=9.8                                                                                                                         
      PREF=1.0E05                                                                                                                   
      RGAS=AKAP*1004.                                                                                                               
                                                                                                                                    
      RTAU1=20.                                                                                                                     
      DELH=60.                                                                                                                      
      DRAG1=0.2                                                                                                                     
                                                                                                                                    
      RTAU=1./(86400.*RTAU1)                                                                                                        
      DRAG=1./(86400.*DRAG1)                                                                                                        
                                                                                                                                    
      DO 50  K=1,KXP                                                                                                                
        HSG(K)=FLOAT(K-1)/FLOAT(KX)                                                                                                 
   50 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 6 K=1,KX                                                                                                                   
        DHS(K)=HSG(K+1)-HSG(K)                                                                                                      
        FSG(K)=0.5*(HSG(K+1)+HSG(K))                                                                                                
    6 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 166 K=1,KX                                                                                                                 
        DHSR(K)=0.5/DHS(K)                                                                                                          
        FSGR(K)=AKAP/(2.*FSG(K))                                                                                                    
  166 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 5 J=1,IY                                                                                                                   
        JJ=IL+1-J                                                                                                                   
        RADANG(J)=-ASIN(SIA(J))                                                                                                     
        RADANG(JJ)=ASIN(SIA(J))                                                                                                     
    5 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 66 J=1,IL                                                                                                                  
        XXX=RADANG(J)                                                                                                               
        THETA(J)=XXX*180./PI                                                                                                        
        GSIN(J)=SIN(XXX)                                                                                                            
        CORIOL(J)=2.*OMEGA*GSIN(J)                                                                                                  
   66 CONTINUE                                                                                                                      
                                                                                                                                    
c  radiative equilibrtium temperature                                                                                               
                                                                                                                                    
      DO 116 J=1,IL                                                                                                                 
      DO 116 K=1,KX                                                                                                                 
        RADEQ(J,K)=288.-60.*(1.-FSG(K))                                                                                             
     *  -DELH*(3.*GSIN(J)*GSIN(J)-1.)/3.                                                                                            
  116 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 2 J=1,IL                                                                                                                   
      DO 2 I=1,IX                                                                                                                   
        ZEROR(I,J)=0.                                                                                                               
    2 CONTINUE                                                                                                                      
      DO 3 N=1,JX                                                                                                                   
      DO 3 M=1,MX                                                                                                                   
        ZEROC(M,N)=ZERO                                                                                                             
    3 CONTINUE                                                                                                                      
                                                                                                                                    
c  DMP is strength of biharmonic damping                                                                                            
                                                                                                                                    
      DO 10 M=1,MX                                                                                                                  
      DO 10 N=1,JX                                                                                                                  
        MM=ISC*(M-1)                                                                                                                
        ELA=FLOAT(MM+N-1)                                                                                                           
        EL2A=ELA*(ELA+1.)/(A*A)                                                                                                     
        DMP(M,N)=VNU*EL2A*EL2A                                                                                                      
   10 CONTINUE                                                                                                                      
                                                                                                                                    
c   topography                                                                                                                      
                                                                                                                                    
      DO 20 I=1,IX                                                                                                                  
      DO 20 J=1,IL                                                                                                                  
      ZSTAR(I,J)=0.                                                                                                                 
   20 CONTINUE                                                                                                                      
      CALL SPEC(ZSTAR,ZS)                                                                                                           
                                                                                                                                    
c  constants used by GEOP                                                                                                           
                                                                                                                                    
      DO 30 K=1,KX                                                                                                                  
      XGEOP1(K)=RGAS*ALOG(HSG(K+1)/FSG(K))                                                                                          
      IF(K.NE.KX) XGEOP2(K+1)=RGAS*ALOG(FSG(K+1)/HSG(K+1))                                                                          
   30 CONTINUE                                                                                                                      
                                                                                                                                    
      RETURN                                                                                                                        
      END                                                                                                                           
********************************************************************                                                                
      SUBROUTINE GEOP(JJ)                                                                                                           
                                                                                                                                    
      INCLUDE 'param.h'                                                                                                             
      INCLUDE 'param1.h'                                                                                                            
      INCLUDE 'const1.h'                                                                                                            
                                                                                                                                    
c   computes spectral geopotential from spectral temperature T and                                                                  
c      spectral topography ZS, as in GFDL Climate Group GCM                                                                         
                                                                                                                                    
c      DO 1 N=1,JX                                                                                                                  
c      DO 1 M=1,MX                                                                                                                  
c        PHI(M,N,KX)=G*ZS(M,N)+XGEOP1(KX)*T(M,N,KX,JJ)                                                                              
c    1 CONTINUE                                                                                                                     
      DO 1 M=1,MXJX                                                                                                                 
        PHI(M,1,KX)=G*ZS(M,1)+XGEOP1(KX)*T(M,1,KX,JJ)                                                                               
    1 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 2 KK=1,KXM                                                                                                                 
        K=KX-KK                                                                                                                     
C        DO 3 N=1,JX                                                                                                                
C        DO 3 M=1,MX                                                                                                                
C          PHI(M,N,K)=PHI(M,N,K+1)+XGEOP2(K+1)*T(M,N,K+1,JJ)                                                                        
C     *       +XGEOP1(K)*T(M,N,K,JJ)                                                                                                
C    3   CONTINUE                                                                                                                   
        DO 3 M=1,MXJX                                                                                                               
          PHI(M,1,K)=PHI(M,1,K+1)+XGEOP2(K+1)*T(M,1,K+1,JJ)                                                                         
     *       +XGEOP1(K)*T(M,1,K,JJ)                                                                                                 
    3   CONTINUE                                                                                                                    
                                                                                                                                    
    2 CONTINUE                                                                                                                      
                                                                                                                                    
      RETURN                                                                                                                        
      END                                                                                                                           
********************************************************************                                                                
      SUBROUTINE STEP(J1,J2,J3,DT,ALPH)                                                                                             
                                                                                                                                    
      INCLUDE 'param.h'                                                                                                             
      INCLUDE 'param1.h'                                                                                                            
      PARAMETER (IXIL=IX*IL)                                                                                                        
      INCLUDE 'const1.h'                                                                                                            
                                                                                                                                    
      DIMENSION XCON(MX,JX),DIVGM(IX,IL,KX),DM(IX,IL),                                                                              
     * SIGM(IX,IL,KXP),TGG(IX,IL,KX),PUV(IX,IL,KX)                                                                                  
      COMPLEX ZERO                                                                                                                  
                                                                                                                                    
      ZERO=(0.,0.)                                                                                                                  
                                                                                                                                    
      J4=J1                                                                                                                         
      IF(ALPH.EQ.0) J4=J2                                                                                                           
                                                                                                                                    
      DO 10 K=1,KX                                                                                                                  
        CALL UVSPEC(VOR(1,1,K,J2),DIV(1,1,K,J2),                                                                                    
     *   DUMC(1,1,1),DUMC(1,1,2))                                                                                                   
                                                                                                                                    
        CALL GRID(DUMC(1,1,1),UG(1,1,K),2)                                                                                          
        CALL GRID(DUMC(1,1,2),VG(1,1,K),2)                                                                                          
        CALL GRID(VOR(1,1,K,J2),VORG(1,1,K),1)                                                                                      
        CALL GRID(DIV(1,1,K,J2),DIVG(1,1,K),1)                                                                                      
        CALL GRID(DIV(1,1,K,J4),DIVGM(1,1,K),1)                                                                                     
        CALL GRID(T(1,1,K,J2),TG(1,1,K),1)                                                                                          
                                                                                                                                    
        DO 11 J=1,IL                                                                                                                
        DO 11 I=1,IX                                                                                                                
          VORG(I,J,K)=VORG(I,J,K)+CORIOL(J)                                                                                         
   11 CONTINUE                                                                                                                      
   10 CONTINUE                                                                                                                      
                                                                                                                                    
      CALL GEOP(J4)                                                                                                                 
                                                                                                                                    
C      DO 12 I=1,IX                                                                                                                 
C      DO 12 J=1,IL                                                                                                                 
C        UMEAN(I,J)=0.                                                                                                              
C        VMEAN(I,J)=0.                                                                                                              
C        DMEAN(I,J)=0.                                                                                                              
C        DM(I,J)=0.                                                                                                                 
C   12 CONTINUE                                                                                                                     
                                                                                                                                    
      DO 12 I=1,IXIL                                                                                                                
        UMEAN(I,1)=0.                                                                                                               
        VMEAN(I,1)=0.                                                                                                               
        DMEAN(I,1)=0.                                                                                                               
        DM(I,1)=0.                                                                                                                  
   12 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 13 K=1,KX                                                                                                                  
C        DO 113 I=1,IX                                                                                                              
C        DO 113 J=1,IL                                                                                                              
C          UMEAN(I,J)=UMEAN(I,J)+UG(I,J,K)*DHS(K)                                                                                   
C          VMEAN(I,J)=VMEAN(I,J)+VG(I,J,K)*DHS(K)                                                                                   
C          DMEAN(I,J)=DMEAN(I,J)+DIVG(I,J,K)*DHS(K)                                                                                 
C          DM(I,J)=DM(I,J)+DIVGM(I,J,K)*DHS(K)                                                                                      
C  113   CONTINUE                                                                                                                   
                                                                                                                                    
        DO 113 I=1,IXIL                                                                                                             
          UMEAN(I,1)=UMEAN(I,1)+UG(I,1,K)*DHS(K)                                                                                    
          VMEAN(I,1)=VMEAN(I,1)+VG(I,1,K)*DHS(K)                                                                                    
          DMEAN(I,1)=DMEAN(I,1)+DIVG(I,1,K)*DHS(K)                                                                                  
          DM(I,1)=DM(I,1)+DIVGM(I,1,K)*DHS(K)                                                                                       
  113   CONTINUE                                                                                                                    
                                                                                                                                    
   13 CONTINUE                                                                                                                      
                                                                                                                                    
c  compute tendency of log(surface pressure)                                                                                        
                                                                                                                                    
      PS(1,1,J2)=ZERO                                                                                                               
      CALL GRAD(PS(1,1,J2),DUMC(1,1,2),DUMC(1,1,3))                                                                                 
      CALL GRID(DUMC(1,1,2),PX,2)                                                                                                   
      CALL GRID(DUMC(1,1,3),PY,2)                                                                                                   
                                                                                                                                    
c      DO 14 I=1,IX                                                                                                                 
c      DO 14 J=1,IL                                                                                                                 
c        DUMR(I,J,1)=UMEAN(I,J)*PX(I,J)+VMEAN(I,J)*PY(I,J)                                                                          
c        DUMR(I,J,2)=-DUMR(I,J,1)-DM(I,J)                                                                                           
c   14 CONTINUE                                                                                                                     
                                                                                                                                    
      DO 14 I=1,IXIL                                                                                                                
        DUMR(I,1,1)=UMEAN(I,1)*PX(I,1)+VMEAN(I,1)*PY(I,1)                                                                           
        DUMR(I,1,2)=-DUMR(I,1,1)-DM(I,1)                                                                                            
   14 CONTINUE                                                                                                                      
                                                                                                                                    
      CALL SPEC(DUMR(1,1,2),PSDT)                                                                                                   
      PSDT(1,1)=ZERO                                                                                                                
                                                                                                                                    
c   compute sigma-dot "vertical" velocity in two ways,                                                                              
c        using divergence from j2 and j1 time steps                                                                                 
                                                                                                                                    
C      DO 16 I=1,IX                                                                                                                 
C      DO 16 J=1,IL                                                                                                                 
C        SIGDT(I,J,1)=0.                                                                                                            
C        SIGM(I,J,1)=0.                                                                                                             
C   16 CONTINUE                                                                                                                     
                                                                                                                                    
      DO 16 I=1,IXIL                                                                                                                
        SIGDT(I,1,1)=0.                                                                                                             
        SIGM(I,1,1)=0.                                                                                                              
   16 CONTINUE                                                                                                                      
                                                                                                                                    
                                                                                                                                    
      DO 250 K=1,KX                                                                                                                 
C        DO 251 I=1,IX                                                                                                              
C        DO 251 J=1,IL                                                                                                              
C          PUV(I,J,K)=(UG(I,J,K)-UMEAN(I,J))*PX(I,J)+                                                                               
C     *     (VG(I,J,K)-VMEAN(I,J))*PY(I,J)                                                                                          
C  251   CONTINUE                                                                                                                   
                                                                                                                                    
        DO 251 I=1,IXIL                                                                                                             
          PUV(I,1,K)=(UG(I,1,K)-UMEAN(I,1))*PX(I,1)+                                                                                
     *     (VG(I,1,K)-VMEAN(I,1))*PY(I,1)                                                                                           
  251   CONTINUE                                                                                                                    
                                                                                                                                    
  250 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 17 K=1,KX                                                                                                                  
C        DO 117 I=1,IX                                                                                                              
C        DO 117 J=1,IL                                                                                                              
C          SIGDT(I,J,K+1)=SIGDT(I,J,K)-DHS(K)*(PUV(I,J,K)+                                                                          
C     *     DIVG(I,J,K)-DMEAN(I,J))                                                                                                 
C          SIGM(I,J,K+1)=SIGM(I,J,K)-DHS(K)*(PUV(I,J,K)+                                                                            
C     *     DIVGM(I,J,K)-DM(I,J))                                                                                                   
C  117   CONTINUE                                                                                                                   
        DO 117 I=1,IXIL                                                                                                             
          SIGDT(I,1,K+1)=SIGDT(I,1,K)-DHS(K)*(PUV(I,1,K)+                                                                           
     *     DIVG(I,1,K)-DMEAN(I,1))                                                                                                  
          SIGM(I,1,K+1)=SIGM(I,1,K)-DHS(K)*(PUV(I,1,K)+                                                                             
     *     DIVGM(I,1,K)-DM(I,1))                                                                                                    
  117   CONTINUE                                                                                                                    
                                                                                                                                    
   17 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 119 K=1,KX                                                                                                                 
c      DO 118 I=1,IX                                                                                                                
c      DO 118 J=1,IL                                                                                                                
c        TGG(I,J,K)=TG(I,J,K)-TREF(K)                                                                                               
c  118 CONTINUE                                                                                                                     
      DO 118 I=1,IXIL                                                                                                               
        TGG(I,1,K)=TG(I,1,K)-TREF(K)                                                                                                
  118 CONTINUE                                                                                                                      
  119 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 18 K=1,KX                                                                                                                  
                                                                                                                                    
c    add some surface drag to lowest layer only                                                                                     
                                                                                                                                    
        DRAGG=0.                                                                                                                    
        IF(K.EQ.KX) DRAGG=DRAG                                                                                                      
                                                                                                                                    
c   verical advection                                                                                                               
                                                                                                                                    
        IF (K.EQ.1) THEN                                                                                                            
C          DO 19 J=1,IL                                                                                                             
C          DO 19 I=1,IX                                                                                                             
C            DUMR(I,J,1)=SIGDT(I,J,2)*(UG(I,J,2)-UG(I,J,1))*DHSR(K)                                                                 
C            DUMR(I,J,2)=SIGDT(I,J,2)*(VG(I,J,2)-VG(I,J,1))*DHSR(K)                                                                 
C            DUMR(I,J,3)=SIGDT(I,J,2)*(TGG(I,J,2)-TGG(I,J,1))*DHSR(K)                                                               
C     *             +SIGM(I,J,2)*(TREF(2)-TREF(1))*DHSR(K)                                                                          
C   19     CONTINUE                                                                                                                 
          DO 19 I=1,IXIL                                                                                                            
            DUMR(I,1,1)=SIGDT(I,1,2)*(UG(I,1,2)-UG(I,1,1))*DHSR(K)                                                                  
            DUMR(I,1,2)=SIGDT(I,1,2)*(VG(I,1,2)-VG(I,1,1))*DHSR(K)                                                                  
            DUMR(I,1,3)=SIGDT(I,1,2)*(TGG(I,1,2)-TGG(I,1,1))*DHSR(K)                                                                
     *             +SIGM(I,1,2)*(TREF(2)-TREF(1))*DHSR(K)                                                                           
   19     CONTINUE                                                                                                                  
        ELSE IF(K.EQ.KX) THEN                                                                                                       
C          DO 199 J=1,IL                                                                                                            
C          DO 199 I=1,IX                                                                                                            
C            DUMR(I,J,1)=SIGDT(I,J,KX)*(UG(I,J,KX)-UG(I,J,KXM))*DHSR(K)                                                             
C            DUMR(I,J,2)=SIGDT(I,J,KX)*(VG(I,J,KX)-VG(I,J,KXM))*DHSR(K)                                                             
C            DUMR(I,J,3)=SIGDT(I,J,KX)*(TGG(I,J,KX)-TGG(I,J,KXM))*DHSR(K)                                                           
C     *             +SIGM(I,J,KX)*(TREF(KX)-TREF(KXM))*DHSR(K)                                                                      
C  199     CONTINUE                                                                                                                 
          DO 199 I=1,IXIL                                                                                                           
            DUMR(I,1,1)=SIGDT(I,1,KX)*(UG(I,1,KX)-UG(I,1,KXM))*DHSR(K)                                                              
            DUMR(I,1,2)=SIGDT(I,1,KX)*(VG(I,1,KX)-VG(I,1,KXM))*DHSR(K)                                                              
            DUMR(I,1,3)=SIGDT(I,1,KX)*(TGG(I,1,KX)-TGG(I,1,KXM))*DHSR(K)                                                            
     *             +SIGM(I,1,KX)*(TREF(KX)-TREF(KXM))*DHSR(K)                                                                       
  199     CONTINUE                                                                                                                  
                                                                                                                                    
        ELSE                                                                                                                        
C          DO 1999 J=1,IL                                                                                                           
C          DO 1999 I=1,IX                                                                                                           
C            DUMR(I,J,1)=(SIGDT(I,J,K+1)*(UG(I,J,K+1)-UG(I,J,K))                                                                    
C     *           +SIGDT(I,J,K)*(UG(I,J,K)-UG(I,J,K-1)))*DHSR(K)                                                                    
C            DUMR(I,J,2)=(SIGDT(I,J,K+1)*(VG(I,J,K+1)-VG(I,J,K))                                                                    
C     *           +SIGDT(I,J,K)*(VG(I,J,K)-VG(I,J,K-1)))*DHSR(K)                                                                    
C            DUMR(I,J,3)=(SIGDT(I,J,K+1)*(TGG(I,J,K+1)-TGG(I,J,K))                                                                  
C     *           +SIGDT(I,J,K)*(TGG(I,J,K)-TGG(I,J,K-1)))*DHSR(K)                                                                  
C     *         +(SIGM(I,J,K+1)*(TREF(K+1)-TREF(K))                                                                                 
C     *          +SIGM(I,J,K)*(TREF(K)-TREF(K-1)))*DHSR(K)                                                                          
C 1999     CONTINUE                                                                                                                 
          DO 1999 I=1,IXIL                                                                                                          
            DUMR(I,1,1)=(SIGDT(I,1,K+1)*(UG(I,1,K+1)-UG(I,1,K))                                                                     
     *           +SIGDT(I,1,K)*(UG(I,1,K)-UG(I,1,K-1)))*DHSR(K)                                                                     
            DUMR(I,1,2)=(SIGDT(I,1,K+1)*(VG(I,1,K+1)-VG(I,1,K))                                                                     
     *           +SIGDT(I,1,K)*(VG(I,1,K)-VG(I,1,K-1)))*DHSR(K)                                                                     
            DUMR(I,1,3)=(SIGDT(I,1,K+1)*(TGG(I,1,K+1)-TGG(I,1,K))                                                                   
     *           +SIGDT(I,1,K)*(TGG(I,1,K)-TGG(I,1,K-1)))*DHSR(K)                                                                   
     *         +(SIGM(I,1,K+1)*(TREF(K+1)-TREF(K))                                                                                  
     *          +SIGM(I,1,K)*(TREF(K)-TREF(K-1)))*DHSR(K)                                                                           
 1999     CONTINUE                                                                                                                  
        ENDIF                                                                                                                       
                                                                                                                                    
c  other grid point manipulations for vort, div, and t tendencies                                                                   
                                                                                                                                    
C        DO 20 J=1,IL                                                                                                               
C        DO 20 I=1,IX                                                                                                               
C          DUMR(I,J,4)=VG(I,J,K)*VORG(I,J,K)                                                                                        
C     *     -RGAS*TGG(I,J,K)*PX(I,J)-DUMR(I,J,1)                                                                                    
C     *      -DRAGG*UG(I,J,K)                                                                                                       
C          DUMR(I,J,5)=-UG(I,J,K)*VORG(I,J,K)                                                                                       
C     *     -RGAS*TGG(I,J,K)*PY(I,J)-DUMR(I,J,2)                                                                                    
C     *      -DRAGG*VG(I,J,K)                                                                                                       
C          DUMR(I,J,6)=0.5*(UG(I,J,K)*UG(I,J,K)+VG(I,J,K)*VG(I,J,K))                                                                
C          DUMR(I,J,7)=UG(I,J,K)*TGG(I,J,K)                                                                                         
C          DUMR(I,J,8)=VG(I,J,K)*TGG(I,J,K)                                                                                         
C          DUMR(I,J,9)=TGG(I,J,K)*DIVG(I,J,K)-DUMR(I,J,3)                                                                           
C     *     +FSGR(K)*TGG(I,J,K)*(SIGDT(I,J,K+1)+SIGDT(I,J,K))                                                                       
C     *     +FSGR(K)*TREF(K)*(SIGM(I,J,K+1)+SIGM(I,J,K))                                                                            
C     *     +AKAP*(TG(I,J,K)*PUV(I,J,K)                                                                                             
C     *       -TGG(I,J,K)*DMEAN(I,J)-TREF(K)*DM(I,J))                                                                               
C     *     -(TG(I,J,K)-RADEQ(J,K))*RTAU                                                                                            
C   20   CONTINUE                                                                                                                   
        DO 20 I=1,IXIL                                                                                                              
          DUMR(I,1,4)=VG(I,1,K)*VORG(I,1,K)                                                                                         
     *     -RGAS*TGG(I,1,K)*PX(I,1)-DUMR(I,1,1)                                                                                     
     *      -DRAGG*UG(I,1,K)                                                                                                        
          DUMR(I,1,5)=-UG(I,1,K)*VORG(I,1,K)                                                                                        
     *     -RGAS*TGG(I,1,K)*PY(I,1)-DUMR(I,1,2)                                                                                     
     *      -DRAGG*VG(I,1,K)                                                                                                        
          DUMR(I,1,6)=0.5*(UG(I,1,K)*UG(I,1,K)+VG(I,1,K)*VG(I,1,K))                                                                 
          DUMR(I,1,7)=UG(I,1,K)*TGG(I,1,K)                                                                                          
          DUMR(I,1,8)=VG(I,1,K)*TGG(I,1,K)                                                                                          
          DUMR(I,1,9)=TGG(I,1,K)*DIVG(I,1,K)-DUMR(I,1,3)                                                                            
     *     +FSGR(K)*TGG(I,1,K)*(SIGDT(I,1,K+1)+SIGDT(I,1,K))                                                                        
     *     +FSGR(K)*TREF(K)*(SIGM(I,1,K+1)+SIGM(I,1,K))                                                                             
     *     +AKAP*(TG(I,1,K)*PUV(I,1,K)                                                                                              
     *       -TGG(I,1,K)*DMEAN(I,1)-TREF(K)*DM(I,1))                                                                                
   20   CONTINUE                                                                                                                    
        DO 220 J=1,IL                                                                                                               
        DO 220 I=1,IX                                                                                                               
          DUMR(I,J,9)=DUMR(I,J,9)-(TG(I,J,K)-RADEQ(J,K))*RTAU                                                                       
  220   CONTINUE                                                                                                                    
                                                                                                                                    
        CALL VDSPEC(DUMR(1,1,4),DUMR(1,1,5),VORDT(1,1,K),                                                                           
     *   DIVDT(1,1,K),2)                                                                                                            
                                                                                                                                    
        CALL SPEC(DUMR(1,1,6),DUMC(1,1,3))                                                                                          
c        DO 221 N=1,JX                                                                                                              
c        DO 221 M=1,MX                                                                                                              
c          DUMC(M,N,3)=DUMC(M,N,3)+PHI(M,N,K)                                                                                       
c     *     +RGAS*TREF(K)*PS(M,N,J4)                                                                                                
c  221 CONTINUE                                                                                                                     
        DO 221 M=1,MXJX                                                                                                             
          DUMC(M,1,3)=DUMC(M,1,3)+PHI(M,1,K)                                                                                        
     *     +RGAS*TREF(K)*PS(M,1,J4)                                                                                                 
  221 CONTINUE                                                                                                                      
                                                                                                                                    
        CALL LAP(DUMC(1,1,3),DUMC(1,1,4))                                                                                           
                                                                                                                                    
C        DO 222 N=1,JX                                                                                                              
C        DO 222 M=1,MX                                                                                                              
C          DIVDT(M,N,K)=DIVDT(M,N,K)-DUMC(M,N,4)                                                                                    
C  222 CONTINUE                                                                                                                     
        DO 222 M=1,MXJX                                                                                                             
          DIVDT(M,1,K)=DIVDT(M,1,K)-DUMC(M,1,4)                                                                                     
  222 CONTINUE                                                                                                                      
                                                                                                                                    
                                                                                                                                    
        CALL DSPEC(DUMR(1,1,7),DUMR(1,1,8),DUMC(1,1,2),1,2)                                                                         
        CALL SPEC(DUMR(1,1,9),DUMC(1,1,3))                                                                                          
                                                                                                                                    
C        DO 22 N=1,JX                                                                                                               
C        DO 22 M=1,MX                                                                                                               
C          TDT(M,N,K)=DUMC(M,N,3)-DUMC(M,N,2)                                                                                       
C   22   CONTINUE                                                                                                                   
        DO 22 M=1,MXJX                                                                                                              
          TDT(M,1,K)=DUMC(M,1,3)-DUMC(M,1,2)                                                                                        
   22   CONTINUE                                                                                                                    
                                                                                                                                    
   18 CONTINUE                                                                                                                      
                                                                                                                                    
c  implicit correction                                                                                                              
                                                                                                                                    
      IF(ALPH.GT.0.) CALL IMPLIC                                                                                                    
                                                                                                                                    
c  damping                                                                                                                          
                                                                                                                                    
c      DO 339 N=1,JX                                                                                                                
c      DO 339 M=1,MX                                                                                                                
c          XCON(M,N)=1./(1.+DMP(M,N)*DT)                                                                                            
c  339 CONTINUE                                                                                                                     
      DO 339 M=1,MXJX                                                                                                               
          XCON(M,1)=1./(1.+DMP(M,1)*DT)                                                                                             
  339 CONTINUE                                                                                                                      
                                                                                                                                    
      DO 318 K=1,KX                                                                                                                 
C        DO 21 N=1,JX                                                                                                               
C        DO 21 M=1,MX                                                                                                               
C          VORDT(M,N,K)=(VORDT(M,N,K)-DMP(M,N)*VOR(M,N,K,J1))                                                                       
C     *      *XCON(M,N)                                                                                                             
C          DIVDT(M,N,K)=(DIVDT(M,N,K)-DMP(M,N)*DIV(M,N,K,J1))                                                                       
C     *      *XCON(M,N)                                                                                                             
C          TDT(M,N,K)=(TDT(M,N,K)-DMP(M,N)*T(M,N,K,J1))                                                                             
C     *      *XCON(M,N)                                                                                                             
C   21   CONTINUE                                                                                                                   
        DO 21 M=1,MXJX                                                                                                              
          VORDT(M,1,K)=(VORDT(M,1,K)-DMP(M,1)*VOR(M,1,K,J1))                                                                        
     *      *XCON(M,1)                                                                                                              
          DIVDT(M,1,K)=(DIVDT(M,1,K)-DMP(M,1)*DIV(M,1,K,J1))                                                                        
     *      *XCON(M,1)                                                                                                              
          TDT(M,1,K)=(TDT(M,1,K)-DMP(M,1)*T(M,1,K,J1))                                                                              
     *      *XCON(M,1)                                                                                                              
   21   CONTINUE                                                                                                                    
                                                                                                                                    
  318 CONTINUE                                                                                                                      
                                                                                                                                    
c  Robert filter                                                                                                                    
                                                                                                                                    
      ROBB=1.-2.*ROB                                                                                                                
      DO 89 K=1,KX                                                                                                                  
C        DO 88 N=1,JX                                                                                                               
C        DO 88 M=1,MX                                                                                                               
C            DUMC(M,N,1)=VOR(M,N,K,J1)                                                                                              
C            DUMC(M,N,2)=DIV(M,N,K,J1)                                                                                              
C            DUMC(M,N,3)=T(M,N,K,J1)                                                                                                
C            VOR(M,N,K,J3)=VOR(M,N,K,J1)+DT*VORDT(M,N,K)                                                                            
C            DIV(M,N,K,J3)=DIV(M,N,K,J1)+DT*DIVDT(M,N,K)                                                                            
C            T(M,N,K,J3)=T(M,N,K,J1)+DT*TDT(M,N,K)                                                                                  
C            VOR(M,N,K,J2)=ROB*(DUMC(M,N,1)+VOR(M,N,K,J3))                                                                          
C     *       +ROBB*VOR(M,N,K,J2)                                                                                                   
C            DIV(M,N,K,J2)=ROB*(DUMC(M,N,2)+DIV(M,N,K,J3))                                                                          
C     *       +ROBB*DIV(M,N,K,J2)                                                                                                   
C            T(M,N,K,J2)=ROB*(DUMC(M,N,3)+T(M,N,K,J3))                                                                              
C     *       +ROBB*T(M,N,K,J2)                                                                                                     
C     88 CONTINUE                                                                                                                   
        DO 88 M=1,MXJX                                                                                                              
          DUMC(M,1,1)=VOR(M,1,K,J1)                                                                                                 
          DUMC(M,1,2)=DIV(M,1,K,J1)                                                                                                 
          DUMC(M,1,3)=T(M,1,K,J1)                                                                                                   
          VOR(M,1,K,J3)=VOR(M,1,K,J1)+DT*VORDT(M,1,K)                                                                               
          DIV(M,1,K,J3)=DIV(M,1,K,J1)+DT*DIVDT(M,1,K)                                                                               
          T(M,1,K,J3)=T(M,1,K,J1)+DT*TDT(M,1,K)                                                                                     
          VOR(M,1,K,J2)=ROB*(DUMC(M,1,1)+VOR(M,1,K,J3))                                                                             
     *     +ROBB*VOR(M,1,K,J2)                                                                                                      
          DIV(M,1,K,J2)=ROB*(DUMC(M,1,2)+DIV(M,1,K,J3))                                                                             
     *     +ROBB*DIV(M,1,K,J2)                                                                                                      
          T(M,1,K,J2)=ROB*(DUMC(M,1,3)+T(M,1,K,J3))                                                                                 
     *     +ROBB*T(M,1,K,J2)                                                                                                        
   88   CONTINUE                                                                                                                    
   89 CONTINUE                                                                                                                      
                                                                                                                                    
C      DO 90 N=1,JX                                                                                                                 
C      DO 90 M=1,MX                                                                                                                 
C        DUMC(M,N,1)=PS(M,N,J1)                                                                                                     
C        PS(M,N,J3)=PS(M,N,J1)+DT*PSDT(M,N)                                                                                         
C        PS(M,N,J2)=ROB*(DUMC(M,N,1)+PS(M,N,J3))                                                                                    
C     *   +ROBB*PS(M,N,J2)                                                                                                          
C   90 CONTINUE                                                                                                                     
      DO 90 M=1,MXJX                                                                                                                
        DUMC(M,1,1)=PS(M,1,J1)                                                                                                      
        PS(M,1,J3)=PS(M,1,J1)+DT*PSDT(M,1)                                                                                          
        PS(M,1,J2)=ROB*(DUMC(M,1,1)+PS(M,1,J3))                                                                                     
     *   +ROBB*PS(M,1,J2)                                                                                                           
   90 CONTINUE                                                                                                                      
                                                                                                                                    
      RETURN                                                                                                                        
      END                                                                                                                           
******************************************************************                                                                  
