C**********************************************************                  DIFUSE.2     
C             SUBROUTINE DIFUSE                                              DIFUSE.3     
C**********************************************************                  DIFUSE.4     
      SUBROUTINE DIFUSE                                                      DIFUSE.5     
C                                                                            DIFUSE.6     
C     Calculates spectral tendencies from restoration (if included)          DIFUSE.7     
C     and biharmonic diffusion.                                              DIFUSE.8     
C                                                                            DIFUSE.9     
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    
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     
C                                                                            RESTIJ.2     
C     Restoration temperature field and constants which determine it,        RESTIJ.3     
C     also contains timescales                                               RESTIJ.4     
C                                                                            RESTIJ.5     
      COMMON/RESTIJ/TTRES(IGB)                                               RESTIJ.6     
     + ,DTNS,DTEP,DTTRP,FAC(NL),DDAMP(NL),TFRC(NL),YRLEN,TRS(NL)             RESTIJ.7     
     +  ,ALR,ZTROP,TGR                                                       RESTIJ.8     
      COMPLEX TTRES                                                          RESTIJ.9     
      real rfcoeff(nlevrf)                                                   DIFUSE.17    
      logical ifirstrf                                                       DIFUSE.18    
      save ifirstrf, rfcoeff                                                 DIFUSE.19    
      data ifirstrf /.true./                                                 DIFUSE.20    
C                                                                            DIFUSE.21    
C Add newtonian cooling and surface Rayleigh friction.                       DIFUSE.22    
C                                                                            DIFUSE.23    
      IF (LRESTIJ) THEN                                                      DIFUSE.24    
      DO 21 IHEM=1,NHEM                                                      DIFUSE.25    
         DO 22 L=1,NL                                                        DIFUSE.26    
            I=(IHEM-1)*NWJ2+(L-1)*IGA                                        DIFUSE.27    
            DO 23 J=1,NWJ2                                                   DIFUSE.28    
               TT(I+J)=TT(I+J)-DDAMP(L)*(T(I+J)-TTRES(I+J))                  DIFUSE.29    
               ZT(I+J)=ZT(I+J)-TFRC(L)*Z(I+J)                                DIFUSE.30    
               DT(I+J)=DT(I+J)-TFRC(L)*D(I+J)                                DIFUSE.31    
23         CONTINUE                                                          DIFUSE.32    
22       CONTINUE                                                            DIFUSE.33    
21    CONTINUE                                                               DIFUSE.34    
C                                                                            DIFUSE.35    
C       No friction on planetary vorticity                                   DIFUSE.36    
C                                                                            DIFUSE.37    
      DO 24 L=1,NL                                                           DIFUSE.38    
         I=(L-1)*IGA+1                                                       DIFUSE.39    
         ZT(I)=ZT(I)+TFRC(L)*EZ                                              DIFUSE.40    
 24   CONTINUE                                                               DIFUSE.41    
      ELSE                                                                   DIFUSE.42    
      IF(DAMP.GT.0.0) THEN                                                   DIFUSE.43    
         I=0                                                                 DIFUSE.44    
         IR=0                                                                DIFUSE.45    
         DO 800 IHEM=1,NHEM                                                  DIFUSE.46    
            DO 20 L=1,NL                                                     DIFUSE.47    
               DO 10 J=1,IDM                                                 DIFUSE.48    
                  I=I+1                                                      DIFUSE.49    
                  IR=IR+1                                                    DIFUSE.50    
                  ZT(I)=ZT(I)-DAMP*(Z(I)-ZRES(IR))                           DIFUSE.51    
                  DT(I)=DT(I)-DAMP*(D(I)-DRES(IR))                           DIFUSE.52    
                  TT(I)=TT(I)-DAMP*(T(I)-TRES(IR))                           DIFUSE.53    
   10          CONTINUE                                                      DIFUSE.54    
               I=I+IGA-IDM                                                   DIFUSE.55    
               IR=IR+IGM-IDM                                                 DIFUSE.56    
   20       CONTINUE                                                         DIFUSE.57    
            I=NWJ2                                                           DIFUSE.58    
            IR=IDM                                                           DIFUSE.59    
  800    CONTINUE                                                            DIFUSE.60    
      ENDIF                                                                  DIFUSE.61    
      ENDIF                                                                  DIFUSE.62    
C                                                                            DIFUSE.63    
C     Add in biharmonic diffusion if required                                DIFUSE.64    
C                                                                            DIFUSE.65    
      IF (AK(2).GT.0.0) THEN                                                 DIFUSE.66    
         DO 820 IHEM=1,NHEM                                                  DIFUSE.67    
            DO 821 L=1,NL                                                    DIFUSE.68    
               JOFF=(L-1)*IGA+(IHEM-1)*NWJ2                                  DIFUSE.69    
               J=JOFF                                                        DIFUSE.70    
               DO 822 MP=1,MFP,MOCT                                          DIFUSE.71    
                  DO 823 IN=MP,NFP,MH                                        DIFUSE.72    
                     J=J+1                                                   DIFUSE.73    
                     AKZ =AK(IN+2-IHEM)                                      DIFUSE.74    
                     AKDT=AK(IN-1+IHEM)                                      DIFUSE.75    
                     ZT(J)=ZT(J)-AKZ *Z(J)                                   DIFUSE.76    
                     DT(J)=DT(J)-AKDT*D(J)                                   DIFUSE.77    
                     TT(J)=TT(J)-AKDT*T(J)                                   DIFUSE.78    
823               CONTINUE                                                   DIFUSE.79    
822            CONTINUE                                                      DIFUSE.80    
               DO 824 KK=1,NTRAC                                             DIFUSE.81    
                  J=JOFF                                                     DIFUSE.82    
                  DO 825 MP=1,MFP,MOCT                                       DIFUSE.83    
                     DO 826 IN=MP,NFP,MH                                     DIFUSE.84    
                        J=J+1                                                DIFUSE.85    
                        AKDT=AK(IN-1+IHEM)                                   DIFUSE.86    
                        TRAT(J,KK)=TRAT(J,KK)-AKDT*TRA(J,KK)                 DIFUSE.87    
826                  CONTINUE                                                DIFUSE.88    
825               CONTINUE                                                   DIFUSE.89    
824            CONTINUE                                                      DIFUSE.90    
821         CONTINUE                                                         DIFUSE.91    
820      CONTINUE                                                            DIFUSE.92    
C                                                                            DIFUSE.93    
C        No diffusion on EZ (planetary vorticity)                            DIFUSE.94    
C                                                                            DIFUSE.95    
         I=1                                                                 DIFUSE.96    
         DO 30 L=1,NL                                                        DIFUSE.97    
            ZT(I)=ZT(I)+AK(2)*EZ                                             DIFUSE.98    
            I=I+IGA                                                          DIFUSE.99    
30       CONTINUE                                                            DIFUSE.100   
      ENDIF                                                                  DIFUSE.101   
C                                                                            DIFUSE.102   
c                                                                            DIFUSE.103   
c 11-08-97 Very crude Rayleigh Friction in the top NLEVRF levels.            DIFUSE.104   
c          NLEVRF is set in the parameter line.                              DIFUSE.105   
c 29-06-98 Changed to use RFCOEFF instead of TFRC.  SMR                      DIFUSE.106   
C 11-05-00 NB RFCOEFF AND TFRC ARE BOTH RAYLEIGH FRICTION                    DIFUSE.107   
C          COEFFICIENT ARRAYS.                      SH                       DIFUSE.108   
c                                                                            DIFUSE.109   
      if (ifirstrf) then                                                     DIFUSE.110   
        ifirstrf=.false.                                                     DIFUSE.111   
        do l=1,nlevrf                                                        DIFUSE.112   
! The 0.0625 is a tunable parameter. 0.1875 recommended for L26.             DIFUSE.113   
          rfcoeff(l)=((1.0-((l-1)*(1.0/float(nlevrf))))/pi2)*0.0625          DIFUSE.114   
        enddo                                                                DIFUSE.115   
      endif                                                                  DIFUSE.116   
      do ihem=1,nhem                                                         DIFUSE.117   
        do l=1,nlevrf                                                        DIFUSE.118   
          i=(ihem-1)*nwj2+(l-1)*iga                                          DIFUSE.119   
          do j=1,nwj2                                                        DIFUSE.120   
            zt(i+j)=zt(i+j)-rfcoeff(l)*z(i+j)                                DIFUSE.121   
            dt(i+j)=dt(i+j)-rfcoeff(l)*d(i+j)                                DIFUSE.122   
          enddo                                                              DIFUSE.123   
        enddo                                                                DIFUSE.124   
      enddo                                                                  DIFUSE.125   
      do l=1,nlevrf                                                          DIFUSE.126   
        i=(l-1)*iga+1                                                        DIFUSE.127   
        zt(i)=zt(i)+rfcoeff(l)*ez                                            DIFUSE.128   
      enddo                                                                  DIFUSE.129   
      RETURN                                                                 DIFUSE.130   
      END                                                                    DIFUSE.131   
