C**********************************************************                  TSTEP.2     
C             SUBROUTINE TSTEP                                               TSTEP.3     
C**********************************************************                  TSTEP.4     
      SUBROUTINE TSTEP                                                       TSTEP.5     
C                                                                            TSTEP.6     
C     Takes an adiabatic timestep in spectral space                          TSTEP.7     
C     either a centred semi-implicit timestep for KOUNT>KITS                 TSTEP.8     
C     or an initial short timestep for KOUNT<KITS                            TSTEP.9     
C                                                                            TSTEP.10    
C     A time filter used for centred timestep.                               TSTEP.11    
C                                                                            TSTEP.12    
C                                                                            TSTEP.13    
C                                                                            PARAM1.2     
C     Determines model resolution                                            PARAM1.3     
C                                                                            PARAM1.4     
      PARAMETER(NN=21,MM=21,NHEM=2,NL=5,MOCT=1,MG=64,JG=16,NWJ2=121          NETCDF.35    
     P         ,NCRAY=64,JGL=JG,NTRAC=1,NLEVRF=1)                            NETCDF.36    
                                                                             PARAM1.7     
C                                                                            PARAM1.8     
C                                                                            PARAM2.2     
C     Sets basic constants, especially those needed for array dimensions     PARAM2.3     
C                                                                            PARAM2.4     
      PARAMETER(MH=2,PI=3.14159265359,PI2=2.0*PI                             PARAM2.5     
     +,NNP=NN+1,MGPP=MG+2,JGP=JG+1,JGG=JG*NHEM,JGGP=JGG+1,MJP=NWJ2+NWJ2      PARAM2.6     
     +,NLM=NL-1,NLP=NL+1,NLPP=NL+2,NLA=NL+3,NLB=NL+4,NL2=NL*NL               PARAM2.7     
     +,IDA=(MG+MG+MG)/2+1,IDB=NWJ2*NL,IDC=IDB+IDB,IDD=MGPP*NL                PARAM2.8     
     +,IDE=NL2*NN,IDF=NCRAY*(MG+1),IDG=JG*NL,IDH=JG*MG                       PARAM2.9     
     +,IDI=NNP/2,IDJ=IDI*IDI,IDK=NL*IDI,IDL=MGPP/2,IDM=NNP/2,IDN=IDM*NL      PARAM2.10    
     +,NWW=1+(MM-1)/MOCT)                                                    PARAM2.11    
      PARAMETER(IGA=NWJ2*NHEM,IGB=IDB*NHEM,IGC=MGPP*NHEM,IGD=IDD*NHEM        PARAM2.12    
     +,IGG=IDG*NHEM,IGL=IDL*NHEM,IGM=IDM*NHEM,IGN=IDN*NHEM                   PARAM2.13    
     +,IGO=IGA+IGA,IGP=IGB+IGB,NFTWG=(5+NTRAC)*NL+3                          PARAM2.14    
     +,NFTGW=(6+3*NTRAC)*NL+2,NFTGD=(3+NTRAC)*NL,NLTR=NL*NTRAC)              PARAM2.15    
C     Number of 2D (surface) output fields. This value is                    PARAM2.16    
C     Doubled due to averaged and instantaneous fields.                      PARAM2.17    
      PARAMETER (N2DFLD=21,NGRPAD=N2DFLD*2*IGC)                              PARAM2.18    
C                                                                            PARAM2.19    
      PARAMETER(NLT=NL+NL)                                                   TSTEP.16    
C                                                                            BLANK.2     
C     Basic planetary parameters for run plus information about              BLANK.3     
C     vertical grid structure                                                BLANK.4     
C                                                                            BLANK.5     
C     Note that RD and GASCON are identical and CPD is set from RD,AKAP.     BLANK.6     
      COMMON        SQ(NNP),RSQ(NNP),SIGMAH(NLM),SIGMA(NL)                   BLANK.7     
     +              ,T01S2(NLM),T0(NL),ALPHA(NL),DSIGMA(NL),RDSIG(NL)        BLANK.8     
     +              ,TKP(NL),C(NL2),SQH(NNP)                                 BLANK.9     
     +              ,MF,MFP,JZF,NF,NFP                                       BLANK.10    
     +              ,AKAP,GA,GASCON,RADEA,WW,P0,PFAC,EZ,AIOCT                BLANK.11    
     +              ,RD,RV,CPD,CLATNT                                        BLANK.12    
     +              ,LRSTRT,LSHORT,LTVEC,LSTRETCH                            BLANK.13    
     +              ,LFLUX                                                   BLANK.14    
     +              ,LBALAN,LRESTIJ                                          BLANK.15    
     +              ,LCLIM, LPERPET, L22L,LOROG ,LCSFCT                      BLANK.16    
     +              ,LNOISE                                                  BLANK.17    
      COMPLEX EZ,AIOCT                                                       BLANK.18    
      LOGICAL LRSTRT,LSHORT,LTVEC,LSTRETCH,LBALAN,LRESTIJ                    BLANK.19    
     +       ,LFLUX,LNOISE                                                   BLANK.20    
     +       ,LCLIM, LPERPET, L22L,LOROG,LCSFCT                              BLANK.21    
C                                                                            BLANK.22    
C                                                                            SPECTR.2     
C     Array ordering in SPECTR must correspond to that in GRIDP.             SPECTR.3     
C                                                                            SPECTR.4     
      COMMON/SPECTR/Z(IGB),D(IGB),T(IGB),TRA(IGB,NTRAC),SP(IGA),GS(IGA)      SPECTR.5     
     :              ,SPA(IGA),VP(IGA),DTE(IGB),TT(IGB)                       SPECTR.6     
     :              ,TRAT(IGB,NTRAC),DT(IGB),ZT(IGB)                         SPECTR.7     
     :              ,ZMI(IGB),DMI(IGB),TMI(IGB)                              SPECTR.8     
     :              ,TRAMI(IGB,NTRAC),SPMI(IGA)                              SPECTR.9     
      COMPLEX Z,D,T,TRA,SP,GS,SPA,VP,DTE,TT,TRAT,DT,ZT                       SPECTR.10    
     :       ,ZMI,DMI,TMI,TRAMI,SPMI                                         SPECTR.11    
C                                                                            SPECTR.12    
C                                                                            BATS.2     
C     Constant arrays and variables associated with time and vertical        BATS.3     
C     differencing. Also counters.                                           BATS.4     
C                                                                            BATS.5     
      COMMON/BATS/  BM1(IDE),AK(NNP),AQ(NL2),G(NL2),TAU(NL2)                 BATS.6     
     +              ,KOUNT,KITS,KSTART,KTOTAL,KRUN,BEGDAY,ITSPD              BATS.7     
     +              ,DELT,DELT2,CV,CG,CT,CQ,PNU,PNU2,PNU21                   BATS.8     
     +              ,NTRACO,CTRA(NTRAC),KOLOUR(NTRAC),RGG(NL2)               BATS.9     
     +              ,BEGDOY,DOY                                              BATS.10    
C                                                                            BATS.11    
C                                                                            RESTOR.2     
C     Restoration fields and timescale                                       RESTOR.3     
C                                                                            RESTOR.4     
      COMMON/RESTOR/ZRES(IGN),DRES(IGN),TRES(IGN),SPRES(IGM),DAMP            RESTOR.5     
C                                                                            RESTOR.6     
      DIMENSION TBM1(NL2),D1(NL),WA(NL)                                      TSTEP.21    
      INTEGER IWA(NL)                                                        TSTEP.22    
      DIMENSION TBME(NL,NL)                                                  TSTEP.23    
      EQUIVALENCE (TBM1(1),TBME(1,1))                                        TSTEP.24    
      COMPLEX TMPA(NL),TMPB(NL)                                              TSTEP.25    
      REAL RMPA(NLT),RMPB(NLT)                                               TSTEP.26    
      EQUIVALENCE (TMPA(1),RMPA(1)),(TMPB(1),RMPB(1))                        TSTEP.27    
      DIMENSION TTR(IGP),DTR(IGP),TMIR(IGP),D1R(NLT)                         TSTEP.28    
      EQUIVALENCE (TTR(1),TT(1)),(DTR(1),DT(1)),(TMIR(1),TMI(1))             TSTEP.29    
     +,(D1R(1),D1(1))                                                        TSTEP.30    
      COMPLEX ZPAV,TPAV,DPAV,DTI,SPAV,TRAPAV                                 TSTEP.31    
      COMPLEX D1,SPPA,VPS,GSI1                                               TSTEP.32    
C                                                                            TSTEP.33    
      IF (KOUNT.GT.KITS) THEN                                                TSTEP.34    
C                                                                            TSTEP.35    
C        Ordinary centred semi-implicit timestep                             TSTEP.36    
C                                                                            TSTEP.37    
         I1=0                                                                TSTEP.38    
         DO 800 IHEM=1,NHEM                                                  TSTEP.39    
            DO 4 MP=1,MFP,MOCT                                               TSTEP.40    
               IBM1=NL2*(MP-2 +IHEM-1)                                       TSTEP.41    
               DO 5 IN=MP,NFP,MH                                             TSTEP.42    
                  I1=I1+1                                                    TSTEP.43    
                  IF (IBM1.NE.-NL2) THEN                                     TSTEP.44    
                     SPPA=SPMI(I1)                                           TSTEP.45    
                     VPS=VP(I1)                                              TSTEP.46    
                     RCN=RSQ(IN+IHEM-1)                                      TSTEP.47    
                     GSI1=GS(I1)                                             TSTEP.48    
                     K=I1                                                    TSTEP.49    
                     CALL SGEMM('N','N',2,NL,NL,1.0,TTR(2*I1-1),IGO,         TSTEP.50    
     +                          G,NL,0.0,RMPA,2)                             TSTEP.51    
                     CALL SGEMM('N','N',2,NL,NL,1.0,TMIR(2*I1-1),IGO,        TSTEP.52    
     +                          G,NL,0.0,RMPB,2)                             TSTEP.53    
                     DO 12 L=1,NL                                            TSTEP.54    
                        D1(L)=RCN*DMI(K)+DELT*(TMPB(L)+GSI1+T0(L)*SPPA       TSTEP.55    
     1                       +RCN*DT(K)+DELT*(TMPA(L)-T0(L)*VPS))            TSTEP.56    
                        K=K+IGA                                              TSTEP.57    
   12                CONTINUE                                                TSTEP.58    
                     K=I1                                                    TSTEP.59    
                     CALL SGEMM('N','N',2,NL,NL,1.0,D1R,2,                   TSTEP.60    
     +                          BM1(IBM1+1),NL,0.0,DTR(2*I1-1),IGO)          TSTEP.61    
                     DO 14 L=1,NL                                            TSTEP.62    
                        VP(I1)=VP(I1)+DT(K)*DSIGMA(L)                        TSTEP.63    
                        K=K+IGA                                              TSTEP.64    
   14                CONTINUE                                                TSTEP.65    
                     IBM1=IBM1+NL2+NL2                                       TSTEP.66    
                     K=I1                                                    TSTEP.67    
                     CALL SGEMM('N','N',2,NL,NL,1.0,DTR(2*I1-1),IGO,         TSTEP.68    
     +                          TAU,NL,0.0,RMPA,2)                           TSTEP.69    
                     DO 16 L=1,NL                                            TSTEP.70    
                        TT(K)=TT(K)-TMPA(L)                                  TSTEP.71    
                        K=K+IGA                                              TSTEP.72    
   16                CONTINUE                                                TSTEP.73    
                  ELSE                                                       TSTEP.74    
                     IBM1=NL2                                                TSTEP.75    
                  ENDIF                                                      TSTEP.76    
    5          CONTINUE                                                      TSTEP.77    
    4       CONTINUE                                                         TSTEP.78    
            I1=NWJ2                                                          TSTEP.79    
  800    CONTINUE                                                            TSTEP.80    
         DO 40 I=1,IGB                                                       TSTEP.81    
            ZPAV=ZMI(I)                                                      TSTEP.82    
            ZMI(I)=PNU21*Z(I)+PNU*ZPAV                                       TSTEP.83    
            Z(I)=ZPAV+DELT2*ZT(I)                                            TSTEP.84    
            TPAV=TMI(I)                                                      TSTEP.85    
            TMI(I)=PNU21*T(I)+PNU*TPAV                                       TSTEP.86    
            T(I)=TPAV+DELT2*TT(I)                                            TSTEP.87    
            DPAV=DMI(I)                                                      TSTEP.88    
            DTI=DT(I)                                                        TSTEP.89    
            DMI(I)=PNU21*D(I)+PNU*DPAV                                       TSTEP.90    
            D(I)=DTI+DTI-DPAV                                                TSTEP.91    
   40    CONTINUE                                                            TSTEP.92    
         DO 141 KK=1,NTRAC                                                   TSTEP.93    
            DO 140 I=1,IGB                                                   TSTEP.94    
               TRAPAV=TRAMI(I,KK)                                            TSTEP.95    
               TRAMI(I,KK)=PNU21*TRA(I,KK)+PNU*TRAPAV                        TSTEP.96    
               TRA(I,KK)=TRAPAV+DELT2*TRAT(I,KK)                             TSTEP.97    
  140       CONTINUE                                                         TSTEP.98    
  141    CONTINUE                                                            TSTEP.99    
         DO 41 I=1,IGA                                                       TSTEP.100   
            SPAV=SPMI(I)                                                     TSTEP.101   
            SPMI(I)=PNU21*SP(I)+PNU*SPAV                                     TSTEP.102   
            SP(I)=SPAV-DELT2*VP(I)                                           TSTEP.103   
   41    CONTINUE                                                            TSTEP.104   
         RETURN                                                              TSTEP.105   
      ELSE                                                                   TSTEP.106   
C                                                                            TSTEP.107   
C        Initial short timestep                                              TSTEP.108   
C                                                                            TSTEP.109   
         I1=0                                                                TSTEP.110   
         DO 820 IHEM=1,NHEM                                                  TSTEP.111   
            DO 8 MP=1,MFP,MOCT                                               TSTEP.112   
               DO 9 IN=MP,NFP,MH                                             TSTEP.113   
                  I1=I1+1                                                    TSTEP.114   
                  IF (I1+IHEM.NE.2) THEN                                     TSTEP.115   
                     SPPA=SPMI(I1)                                           TSTEP.116   
                     VPS=VP(I1)                                              TSTEP.117   
                     RCN=RSQ(IN+IHEM-1)                                      TSTEP.118   
                     GSI1=GS(I1)                                             TSTEP.119   
                     DO 21 IL=1,NL2                                          TSTEP.120   
                        TBM1(IL)=AQ(IL)                                      TSTEP.121   
   21                CONTINUE                                                TSTEP.122   
                     DO 30 L=1,NL                                            TSTEP.123   
                        TBME(L,L)=TBME(L,L)+RCN                              TSTEP.124   
   30                CONTINUE                                                TSTEP.125   
                     CALL MATINV(TBME,NL,NL,IWA,WA)                          TSTEP.126   
                     K=I1                                                    TSTEP.127   
                     CALL SGEMM('N','N',2,NL,NL,1.0,TTR(2*I1-1),IGO,         TSTEP.128   
     +                          G,NL,0.0,RMPA,2)                             TSTEP.129   
                     CALL SGEMM('N','N',2,NL,NL,1.0,TMIR(2*I1-1),IGO,        TSTEP.130   
     +                          G,NL,0.0,RMPB,2)                             TSTEP.131   
                     DO 22 L=1,NL                                            TSTEP.132   
                        D1(L)=RCN*DMI(K)+DELT*(TMPB(L)+GSI1+T0(L)*SPPA+      TSTEP.133   
     1                        RCN*DT(K)+DELT*(TMPA(L)-T0(L)*VPS))            TSTEP.134   
                        K=K+IGA                                              TSTEP.135   
   22                CONTINUE                                                TSTEP.136   
                     K=I1                                                    TSTEP.137   
                     CALL SGEMM('N','N',2,NL,NL,1.0,D1R,2,                   TSTEP.138   
     +                          TBM1,NL,0.0,DTR(2*I1-1),IGO)                 TSTEP.139   
                     DO 24 L=1,NL                                            TSTEP.140   
                        VP(I1)=VP(I1)+DT(K)*DSIGMA(L)                        TSTEP.141   
                        K=K+IGA                                              TSTEP.142   
   24                CONTINUE                                                TSTEP.143   
                     K=I1                                                    TSTEP.144   
                     CALL SGEMM('N','N',2,NL,NL,1.0,DTR(2*I1-1),IGO,         TSTEP.145   
     +                          TAU,NL,0.0,RMPA,2)                           TSTEP.146   
                     DO 26 L=1,NL                                            TSTEP.147   
                        TT(K)=TT(K)-TMPA(L)                                  TSTEP.148   
                        K=K+IGA                                              TSTEP.149   
   26                CONTINUE                                                TSTEP.150   
                  ENDIF                                                      TSTEP.151   
    9          CONTINUE                                                      TSTEP.152   
    8       CONTINUE                                                         TSTEP.153   
            I1=NWJ2                                                          TSTEP.154   
  820    CONTINUE                                                            TSTEP.155   
         DO 42 I=1,IGB                                                       TSTEP.156   
            Z(I)=ZMI(I)+DELT2*ZT(I)                                          TSTEP.157   
            T(I)=TMI(I)+DELT2*TT(I)                                          TSTEP.158   
            D(I)=DT(I)+DT(I)-DMI(I)                                          TSTEP.159   
   42    CONTINUE                                                            TSTEP.160   
         DO 143 KK=1,NTRAC                                                   TSTEP.161   
            DO 142 I=1,IGB                                                   TSTEP.162   
               TRA(I,KK)=TRAMI(I,KK)+DELT2*TRAT(I,KK)                        TSTEP.163   
  142       CONTINUE                                                         TSTEP.164   
  143    CONTINUE                                                            TSTEP.165   
         DO 43 I=1,IGA                                                       TSTEP.166   
            SP(I)=SPMI(I)-DELT2*VP(I)                                        TSTEP.167   
   43    CONTINUE                                                            TSTEP.168   
         DO 29 L=1,NL2                                                       TSTEP.169   
            AQ(L)=AQ(L)*4.                                                   TSTEP.170   
   29    CONTINUE                                                            TSTEP.171   
         RETURN                                                              TSTEP.172   
      ENDIF                                                                  TSTEP.173   
C                                                                            TSTEP.174   
      END                                                                    TSTEP.175   
