C**********************************************************                  HANALV.2     
C             SUBROUTINE HANALV                                              HANALV.3     
C**********************************************************                  HANALV.4     
      SUBROUTINE HANALV(SPA,VP,DTE,TT,QT,DT,ZT,SPG,VPG,EG                    HANALV.5     
     *                 ,TNLG,QNLG,FUG,FVG,VTG,VQG,FVGT,FUGT)                 HANALV.6     
C                                                                            HANALV.7     
C     Perform all the direct Legendre transforms for the adiabatic           HANALV.8     
C     part of the timestep at the current latitude (pair), in place          HANALV.9     
C     of separate calls to HANAL for individual transform types.             HANALV.10    
C                                                                            HANALV.11    
C     The following types of Legendre function and thence types of           HANALV.12    
C     transform may be used:                                                 HANALV.13    
C        ITYPE=1,2  :  ALP   :  ALPN(,,,1)   :  normal transform.            HANALV.14    
C        ITYPE=3,4  :  DALP  :  ALPN(,,,2)   :  y-derivative.                HANALV.15    
C        ITYPE=5,6  :  RLP   :  ALPN(,,,3)   :  del(-2).                     HANALV.16    
C        ITYPE=7,8  :  RDLP  :  ALPN(,,,4)   :  y-derivative of del(-2).     HANALV.17    
C     An even/odd value of ITYPE denotes a spectral field of even/odd        HANALV.18    
C     symmetry.                                                              HANALV.19    
C                                                                            HANALV.20    
C     Maximum vector efficiency is achieved by chaining all multi-level      HANALV.21    
C     transforms in one loop and by chaining all single-level transforms     HANALV.22    
C     in a second loop.                                                      HANALV.23    
C                                                                            HANALV.24    
C     All dummy argument arrays are declared complex.                        HANALV.25    
C     All array dimensions are parameters.                                   HANALV.26    
C     Multi-level arrays are 3-dimensional.                                  HANALV.27    
C                                                                            HANALV.28    
C     Tracers (QNLG,QT) are included conditionally inside inner loops.       HANALV.29    
C     This does *not* affect vectorisation, since the number of tracers,     HANALV.30    
C     NTRAC, is a parameter and the relevant IF constructs are removed       HANALV.31    
C     at compilation.  The same applies to code conditional on NHEM.         HANALV.32    
C     (Only tested on Cray J90 with cf77).                                   HANALV.33    
C                                                                            HANALV.34    
C     NOTE: The y-derivative transforms use integration by parts and         HANALV.35    
C           are valid only if the input field has zero zonal mean at         HANALV.36    
C           both poles.                                                      HANALV.37    
C                                                                            HANALV.38    
C     NOTE: *** THE INPUT FOURIER FIELDS ARE MODIFIED IF GLOBAL ***          HANALV.39    
C                                                                            HANALV.40    
C     Version for RSGUP3.                     Mike Blackburn,  12.01.95.     HANALV.41    
C     Include tracers for IGCM2.              Piers Forster,   ??.04.97.     HANALV.42    
C     Conditional separate loops for tracers. Mike Blackburn,  04.09.98.     HANALV.43    
C                                                                            HANALV.44    
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                                                                            LEGAU.2     
C     Legendre polynomials and information about gaussian latitudes          LEGAU.3     
C                                                                            LEGAU.4     
      COMMON/LEGAU/ ALPJ(MJP),DALPJ(MJP)                                     LEGAU.5     
     +              ,ALP(NWJ2,2,JGL),DALP(NWJ2,2,JGL)                        LEGAU.6     
     +              ,RLP(NWJ2,2,JGL),RDLP(NWJ2,2,JGL)                        LEGAU.7     
     +              ,SI(JGG),CS(JGG),SISQ(JGG),CSSQ(JGG),SECSQ(JGG)          LEGAU.8     
     +              ,ALAT(JGG),GWT(JGG),AW(JGG),JH,JL,JINC                   LEGAU.9     
C                                                                            LEGAU.10    
C                                                                            POLYNO.2     
C     Polynomial used to aid vectorization of Legendre transforms            POLYNO.3     
C                                                                            POLYNO.4     
      COMMON/POLYNO/POLY(NWJ2,2,4),CMPA(IGL)                                 POLYNO.5     
      COMPLEX CMPA                                                           POLYNO.6     
C                                                                            POLYNO.7     
C                                                                            HANALV.50    
      COMPLEX     SPA(NWJ2,NHEM),VP(NWJ2,NHEM),DTE(NWJ2,NHEM,NL)             HANALV.51    
     *           ,TT(NWJ2,NHEM,NL),DT(NWJ2,NHEM,NL),ZT(NWJ2,NHEM,NL)         HANALV.52    
     *           ,QT(NWJ2,NHEM,NL*NTRAC)                                     HANALV.53    
      COMPLEX     SPG(IDL,NHEM),VPG(IDL,NHEM),EG(IDL,NHEM,NL)                HANALV.54    
     *           ,TNLG(IDL,NHEM,NL),FUG(IDL,NHEM,NL),FVG(IDL,NHEM,NL)        HANALV.55    
     *           ,VTG(IDL,NHEM,NL),FVGT(IDL,NHEM,NL),FUGT(IDL,NHEM,NL)       HANALV.56    
     *           ,QNLG(IDL,NHEM,NL*NTRAC),VQG(IDL,NHEM,NL*NTRAC)             HANALV.57    
      COMPLEX     TEMP                                                       HANALV.58    
      REAL        ALPN(NWJ2,2,JGL,4)                                         HANALV.59    
      EQUIVALENCE (ALPN(1,1,1,1),ALP(1,1,1))                                 HANALV.60    
C                                                                            HANALV.61    
C     For a global run, sum and difference the complete Fourier              HANALV.62    
C     transforms at the northern and southern latitude rows to give          HANALV.63    
C     the even and odd contributions : E=(N+S)/2, O=(N-S)/2.                 HANALV.64    
C     For Fourier fields symmetric about equator  : even precedes odd.       HANALV.65    
C     For Fourier fields asymmetric about equator : odd precedes even.       HANALV.66    
C                                                                            HANALV.67    
      IF (NHEM.EQ.2) THEN                                                    HANALV.68    
C                                                                            HANALV.69    
C        Single-level fields.                                                HANALV.70    
C                                                                            HANALV.71    
         DO 10 IM=1,NWW                                                      HANALV.72    
C           Surface pressure : symmetric.                                    HANALV.73    
            TEMP=SPG(IM,1)                                                   HANALV.74    
            SPG(IM,1)=0.5*(TEMP+SPG(IM,2))                                   HANALV.75    
            SPG(IM,2)=0.5*(TEMP-SPG(IM,2))                                   HANALV.76    
C           Surface pressure tendency : symmetric.                           HANALV.77    
            TEMP=VPG(IM,1)                                                   HANALV.78    
            VPG(IM,1)=0.5*(TEMP+VPG(IM,2))                                   HANALV.79    
            VPG(IM,2)=0.5*(TEMP-VPG(IM,2))                                   HANALV.80    
   10    CONTINUE                                                            HANALV.81    
C                                                                            HANALV.82    
C        Multi-level fields.                                                 HANALV.83    
C                                                                            HANALV.84    
         DO 20 IM=1,NWW                                                      HANALV.85    
            DO 20 L=1,NL                                                     HANALV.86    
C              Divergence tendency : energy term : symmetric.                HANALV.87    
               TEMP=EG(IM,1,L)                                               HANALV.88    
               EG(IM,1,L)=0.5*(TEMP+EG(IM,2,L))                              HANALV.89    
               EG(IM,2,L)=0.5*(TEMP-EG(IM,2,L))                              HANALV.90    
C              Temperature tendency : main + d/dx part : symmetric.          HANALV.91    
               TEMP=TNLG(IM,1,L)                                             HANALV.92    
               TNLG(IM,1,L)=0.5*(TEMP+TNLG(IM,2,L))                          HANALV.93    
               TNLG(IM,2,L)=0.5*(TEMP-TNLG(IM,2,L))                          HANALV.94    
               IF (NTRAC.EQ.1) THEN                                          HANALV.95    
C                 Tracer tendency : main + d/dx part : symmetric.            HANALV.96    
                  TEMP=QNLG(IM,1,L)                                          HANALV.97    
                  QNLG(IM,1,L)=0.5*(TEMP+QNLG(IM,2,L))                       HANALV.98    
                  QNLG(IM,2,L)=0.5*(TEMP-QNLG(IM,2,L))                       HANALV.99    
               ENDIF                                                         HANALV.100   
C              Divergence tendency : d/dx part : symmetric.                  HANALV.101   
               TEMP=FUG(IM,1,L)                                              HANALV.102   
               FUG(IM,1,L)=0.5*(TEMP+FUG(IM,2,L))                            HANALV.103   
               FUG(IM,2,L)=0.5*(TEMP-FUG(IM,2,L))                            HANALV.104   
C              Vorticity tendency : d/dx part : anti-symmetric.              HANALV.105   
               TEMP=FVG(IM,1,L)                                              HANALV.106   
               FVG(IM,1,L)=0.5*(TEMP-FVG(IM,2,L))                            HANALV.107   
               FVG(IM,2,L)=0.5*(TEMP+FVG(IM,2,L))                            HANALV.108   
C              Temperature tendency : d/dy part : anti-symmetric.            HANALV.109   
               TEMP=VTG(IM,1,L)                                              HANALV.110   
               VTG(IM,1,L)=0.5*(TEMP-VTG(IM,2,L))                            HANALV.111   
               VTG(IM,2,L)=0.5*(TEMP+VTG(IM,2,L))                            HANALV.112   
               IF (NTRAC.EQ.1) THEN                                          HANALV.113   
C                 Tracer tendency : d/dy part : anti-symmetric.              HANALV.114   
                  TEMP=VQG(IM,1,L)                                           HANALV.115   
                  VQG(IM,1,L)=0.5*(TEMP-VQG(IM,2,L))                         HANALV.116   
                  VQG(IM,2,L)=0.5*(TEMP+VQG(IM,2,L))                         HANALV.117   
               ENDIF                                                         HANALV.118   
C              Divergence tendency : d/dy part : anti-symmetric.             HANALV.119   
               TEMP=FVGT(IM,1,L)                                             HANALV.120   
               FVGT(IM,1,L)=0.5*(TEMP-FVGT(IM,2,L))                          HANALV.121   
               FVGT(IM,2,L)=0.5*(TEMP+FVGT(IM,2,L))                          HANALV.122   
C              Vorticity tendency : d/dy part : symmetric.                   HANALV.123   
               TEMP=FUGT(IM,1,L)                                             HANALV.124   
               FUGT(IM,1,L)=0.5*(TEMP+FUGT(IM,2,L))                          HANALV.125   
               FUGT(IM,2,L)=0.5*(TEMP-FUGT(IM,2,L))                          HANALV.126   
   20    CONTINUE                                                            HANALV.127   
C                                                                            HANALV.128   
C        Tracers, treated as a single (NL*NTRAC)-level field.                HANALV.129   
C                                                                            HANALV.130   
         IF (NTRAC.GT.1) THEN                                                HANALV.131   
            DO 22 IM=1,NWW                                                   HANALV.132   
               DO 22 L=1,NL*NTRAC                                            HANALV.133   
C                 Tracer tendency : main + d/dx part : symmetric.            HANALV.134   
                  TEMP=QNLG(IM,1,L)                                          HANALV.135   
                  QNLG(IM,1,L)=0.5*(TEMP+QNLG(IM,2,L))                       HANALV.136   
                  QNLG(IM,2,L)=0.5*(TEMP-QNLG(IM,2,L))                       HANALV.137   
C                 Tracer tendency : d/dy part : anti-symmetric.              HANALV.138   
                  TEMP=VQG(IM,1,L)                                           HANALV.139   
                  VQG(IM,1,L)=0.5*(TEMP-VQG(IM,2,L))                         HANALV.140   
                  VQG(IM,2,L)=0.5*(TEMP+VQG(IM,2,L))                         HANALV.141   
   22       CONTINUE                                                         HANALV.142   
         ENDIF                                                               HANALV.143   
C                                                                            HANALV.144   
      ENDIF                                                                  HANALV.145   
C                                                                            HANALV.146   
C     Set up the appropriate Gaussian weight for the current latitude,       HANALV.147   
C     dependent on transform type.                                           HANALV.148   
C     Assumes JH in /LEGAU/ contains latitude counter from calling loop.     HANALV.149   
C                                                                            HANALV.150   
      AW1256=AW(JH)*CSSQ(JH)                                                 HANALV.151   
      AW3478=-AW(JH)                                                         HANALV.152   
C                                                                            HANALV.153   
C     Calculate POLY array in a vector loop before the main transforms,      HANALV.154   
C     for the required Legendre Function types.                              HANALV.155   
C     Both even and odd functions are required, irrespective of NHEM.        HANALV.156   
C     Second subscript of ALPN denotes odd or even subset of Legendre        HANALV.157   
C     functions, and depends of symmetry of spectral field.                  HANALV.158   
C     Fourth subscript of ALPN is Legendre function type, (ITYPE+1)/2.       HANALV.159   
C                                                                            HANALV.160   
      DO 30 IHEM=1,2                                                         HANALV.161   
         DO 30 IP=1,NWJ2                                                     HANALV.162   
            POLY(IP,IHEM,1)=AW1256*ALPN(IP,IHEM,JL,1)                        HANALV.163   
            POLY(IP,IHEM,2)=AW3478*ALPN(IP,IHEM,JL,2)                        HANALV.164   
   30 CONTINUE                                                               HANALV.165   
C                                                                            HANALV.166   
C     Transform single-level fields.                                         HANALV.167   
C     Vectorisation is over total wavenumber for each zonal wavenumber.      HANALV.168   
C                                                                            HANALV.169   
      IM=0                                                                   HANALV.170   
      IP=0                                                                   HANALV.171   
      DO 40 M=0,MM-1,MOCT                                                    HANALV.172   
         IM=IM+1                                                             HANALV.173   
         DO 40 N=M,NN-1,2                                                    HANALV.174   
            IP=IP+1                                                          HANALV.175   
C           Surface pressure          : type 2.                              HANALV.176   
            SPA(IP,1)=SPA(IP,1) + POLY(IP,1,1)*SPG(IM,1)                     HANALV.177   
C           Surface pressure tendency : type 2.                              HANALV.178   
            VP (IP,1)=VP (IP,1) + POLY(IP,1,1)*VPG(IM,1)                     HANALV.179   
            IF (NHEM.EQ.2) THEN                                              HANALV.180   
C              Surface pressure          : type 2.                           HANALV.181   
               SPA(IP,2)=SPA(IP,2) + POLY(IP,2,1)*SPG(IM,2)                  HANALV.182   
C              Surface pressure tendency : type 2.                           HANALV.183   
               VP (IP,2)=VP (IP,2) + POLY(IP,2,1)*VPG(IM,2)                  HANALV.184   
            ENDIF                                                            HANALV.185   
   40 CONTINUE                                                               HANALV.186   
C                                                                            HANALV.187   
C     Transform multi-level fields.                                          HANALV.188   
C     Inner loop vectorisation is over total wavenumber, to access           HANALV.189   
C     spectral memory sequentially, avoiding skip distances being a          HANALV.190   
C     multiple of 8 (which causes memory bank conflicts on Cray vector       HANALV.191   
C     machines).                                                             HANALV.192   
C                                                                            HANALV.193   
      IM=0                                                                   HANALV.194   
      IP=0                                                                   HANALV.195   
      DO 50 M=0,MM-1,MOCT                                                    HANALV.196   
         IM=IM+1                                                             HANALV.197   
         IPM=IP                                                              HANALV.198   
         DO 50 L=1,NL                                                        HANALV.199   
            IP=IPM                                                           HANALV.200   
            DO 50 N=M,NN-1,2                                                 HANALV.201   
               IP=IP+1                                                       HANALV.202   
C              Divergence tendency  : energy term      : type 2.             HANALV.203   
               DTE(IP,1,L)=DTE(IP,1,L) + POLY(IP,1,1)*EG  (IM,1,L)           HANALV.204   
C              Temperature tendency : main + d/dx part : type 2.             HANALV.205   
C              Temperature tendency : d/dy part        : type 4.             HANALV.206   
               TT (IP,1,L)=TT (IP,1,L) + POLY(IP,1,1)*TNLG(IM,1,L)           HANALV.207   
     *                                 + POLY(IP,1,2)*VTG (IM,1,L)           HANALV.208   
               IF (NTRAC.EQ.1) THEN                                          HANALV.209   
C                 Tracer tendency   : main + d/dx part : type 2.             HANALV.210   
C                 Tracer tendency   : d/dy part        : type 4.             HANALV.211   
                  QT(IP,1,L)=QT(IP,1,L) +POLY(IP,1,1)*QNLG(IM,1,L)           HANALV.212   
     *                                  +POLY(IP,1,2)*VQG (IM,1,L)           HANALV.213   
               ENDIF                                                         HANALV.214   
C              Divergence tendency  : d/dx part        : type 2.             HANALV.215   
C              Divergence tendency  : d/dy part        : type 4.             HANALV.216   
               DT (IP,1,L)=DT (IP,1,L) + POLY(IP,1,1)*FUG (IM,1,L)           HANALV.217   
     *                                 + POLY(IP,1,2)*FVGT(IM,1,L)           HANALV.218   
C              Vorticity tendency   : d/dx part        : type 1.             HANALV.219   
C              Vorticity tendency   : d/dy part        : type 3.             HANALV.220   
               ZT (IP,1,L)=ZT (IP,1,L) + POLY(IP,2,1)*FVG (IM,1,L)           HANALV.221   
     *                                 + POLY(IP,2,2)*FUGT(IM,1,L)           HANALV.222   
               IF (NHEM.EQ.2) THEN                                           HANALV.223   
C                 Divergence tendency  : energy term      : type 2.          HANALV.224   
                  DTE(IP,2,L)=DTE(IP,2,L) + POLY(IP,2,1)*EG  (IM,2,L)        HANALV.225   
C                 Temperature tendency : main + d/dx part : type 2.          HANALV.226   
C                 Temperature tendency : d/dy part        : type 4.          HANALV.227   
                  TT (IP,2,L)=TT (IP,2,L) + POLY(IP,2,1)*TNLG(IM,2,L)        HANALV.228   
     *                                    + POLY(IP,2,2)*VTG (IM,2,L)        HANALV.229   
                  IF (NTRAC.EQ.1) THEN                                       HANALV.230   
C                    Tracer tendency   : d/dy part        : type 4.          HANALV.231   
                     QT(IP,2,L)=QT(IP,2,L) +POLY(IP,2,1)*QNLG(IM,2,L)        HANALV.232   
     *                                     +POLY(IP,2,2)*VQG (IM,2,L)        HANALV.233   
                  ENDIF                                                      HANALV.234   
C                 Divergence tendency  : d/dx part        : type 2.          HANALV.235   
C                 Divergence tendency  : d/dy part        : type 4.          HANALV.236   
                  DT (IP,2,L)=DT (IP,2,L) + POLY(IP,2,1)*FUG (IM,2,L)        HANALV.237   
     *                                    + POLY(IP,2,2)*FVGT(IM,2,L)        HANALV.238   
C                 Vorticity tendency   : d/dx part        : type 1.          HANALV.239   
C                 Vorticity tendency   : d/dy part        : type 3.          HANALV.240   
                  ZT (IP,2,L)=ZT (IP,2,L) + POLY(IP,1,1)*FVG (IM,2,L)        HANALV.241   
     *                                    + POLY(IP,1,2)*FUGT(IM,2,L)        HANALV.242   
               ENDIF                                                         HANALV.243   
   50 CONTINUE                                                               HANALV.244   
C                                                                            HANALV.245   
C     Transform tracers, treated as a single (NL*NTRAC)-level field.         HANALV.246   
C                                                                            HANALV.247   
      IF (NTRAC.GT.1) THEN                                                   HANALV.248   
         IM=0                                                                HANALV.249   
         IP=0                                                                HANALV.250   
         DO 60 M=0,MM-1,MOCT                                                 HANALV.251   
            IM=IM+1                                                          HANALV.252   
            IPM=IP                                                           HANALV.253   
            DO 60 L=1,NL*NTRAC                                               HANALV.254   
               IP=IPM                                                        HANALV.255   
               DO 60 N=M,NN-1,2                                              HANALV.256   
                  IP=IP+1                                                    HANALV.257   
C                 Tracer tendency   : main + d/dx part : type 2.             HANALV.258   
C                 Tracer tendency   : d/dy part        : type 4.             HANALV.259   
                  QT(IP,1,L)=QT(IP,1,L)+POLY(IP,1,1)*QNLG(IM,1,L)            HANALV.260   
     *                                 +POLY(IP,1,2)*VQG (IM,1,L)            HANALV.261   
                  IF (NHEM.EQ.2) THEN                                        HANALV.262   
                     QT(IP,2,L)=QT(IP,2,L)+POLY(IP,2,1)*QNLG(IM,2,L)         HANALV.263   
     *                                    +POLY(IP,2,2)*VQG (IM,2,L)         HANALV.264   
                  ENDIF                                                      HANALV.265   
   60    CONTINUE                                                            HANALV.266   
      ENDIF                                                                  HANALV.267   
C                                                                            HANALV.268   
      RETURN                                                                 HANALV.269   
      END                                                                    HANALV.270   
