SUBROUTINE COEFF

      INCLUDE 'PARAMETERS.INC'

      INCLUDE 'PLANET.INC'

      DIMENSION D(MB,NALT) , H(MB,NALT) , DELTAZ(NALT)

      COMMON

     &      /ONED7/   A(MB,NALT) , C(MB,NALT)

     &      /ONED17/  NTRANS

     &      /BLOCK18/ CONC(MB,NALT)

     &      /ONED24/  VMAX1(MB) , VMAXN(MB) , R0 , RN

     &      /ONED26/  XMASS(MB)

     &      /ONED28/  DEN(NALT)

     &      /BLOCK10/ Z(NALT)      ! FROM CHEM

     &      /BLOCK11/ T(NALT)      ! FROM CHEM

     &      /ONED29/  EDDY(NALT)

     &      /ONED30/  B(MB,NALT)

     &      /ONED40/  MTRANS(MB)

     &      /ONED61/  G1 , GN      ! USED IN BOUNDRY FLUX CALCULATIONS



C THIS SUBROUTINE CALCULATES COEFFICIENTS FOR TRANSPORT FOR ALL
C LEVELS. ONLY NEEDS TO BE CALLED ONCE UNLESS MAJOR PERTURBATION
C OCCURS

      DATA ADIFH2 , SDIFH2 / 3.64d-5 , 1.75 /

      R = RADIUS*1.0d5 ! CONVERT TO CM
      DO K = 1 , NALT-1
        DELTAZ(K) = (1d5)*(Z(K+1) - Z(K))
        SCALE = (DELTAZ(K)/DLOG(DEN(K)/DEN(K+1)))

C COMPUTE MEAN MASS. ONLY THE FIRST 40 ARE NEEDED. THE CONCENTRATIONS OF THE REST
C ARE SMALL ENOUGH THAT THEY WILL NOT SIGNIFICANTLY CONTRIBUTE TO THE AVERAGE MASS

        DENSITY = 0.D0
        AVGMASS = 0.D0
        DO M = 1 , 40
          DENSITY = DENSITY + CONC(M,K)
          AVGMASS = AVGMASS + CONC(M,K)*XMASS(M)
        ENDDO

        AVGMASS = AVGMASS/DENSITY

C DIFFUSION COEFFICIENT  SCALE HEIGHT WIND

        TAVG = 0.5D0*(T(K) + T(K+1))
        DENAVG = 0.5D0*(DEN(K) + DEN(K+1))
        EAVG = 0.5D0*(EDDY(K) + EDDY(K+1))

        DO M = 1 , NTRANS
          MVARY = MTRANS(M)
          DIFF = ADIFH2*(TAVG**(SDIFH2-1.0d0))*(7.3439D+21/DENAVG)*
     &                    SQRT(2.01594D0/XMASS(MVARY))

          D(M,K) = (DIFF + EAVG)

C THERMAL DIFFUSION OPTION. DEFAULT ALPHA = 0

          ALPHA = 0.0d0
C          IF( AMOL(M) .EQ. 'H       ') ALPHA = -0.39
C          IF( AMOL(M) .EQ. 'H2      ') ALPHA = -0.31
C          IF( AMOL(M) .EQ. 'HE      ') ALPHA = -0.36
C          IF( AMOL(M) .EQ. 'O       ') ALPHA = -0.27


          ALPHA = ALPHA*((DEN(K) - CONC(MVARY,K))/DEN(K))
          ZMM = XMASS(MVARY)/AVGMASS
          FA = 1.0D0 + ((1.0D0 + ALPHA - ZMM)*(SCALE/ZMM)*
     &               (2.0D0/(T(K+1) + T(K)))*
     &               ((T(K+1) - T(K))/DELTAZ(K)))
          H(M,K) = D(M,K)/((EAVG + (DIFF*ZMM*FA))/SCALE)
        ENDDO
      ENDDO

C FILL IN LAST LEVEL

      DO M = 1 , NTRANS
        D(M,NALT) = D(M,NALT-1)
        H(M,NALT) = H(M,NALT-1)
      ENDDO


C CALCULATE COEFFICIENTS  A  B  C

      DO K = 1 , NALT-2
        TOP15 = R + 0.5D0*(1E5)*(Z(K+1) + Z(K+2))
        TOP05 = R + 0.5D0*(1E5)*(Z(K) + Z(K+1))
        BOTTOM = R + (1d5)*Z(K+1)
        R15 = (TOP15/BOTTOM)**2
        R05 = (TOP05/BOTTOM)**2
        DEL05 = DELTAZ(K)
        DEL15 = DELTAZ(K+1)
        BIGDEL = 0.5D0*(1d5)*(Z(K+2)-Z(K))

        DO M = 1 , NTRANS
          D05P = DEXP(0.5D0*DEL05/H(M,K))/DEL05
          D05M = DEXP(-0.5D0*DEL05/H(M,K))/DEL05
          D15P = DEXP(0.5D0*DEL15/H(M,K+1))/DEL15
          D15M = DEXP(-0.5D0*DEL15/H(M,K+1))/DEL15

          A(M,K+1) = -R05*D(M,K)*D05M/BIGDEL
          B(M,K+1) = R15*D(M,K+1)*D15M/BIGDEL +
     &               R05*D(M,K)*D05P/BIGDEL
          C(M,K+1) = -R15*D(M,K+1)*D15P/BIGDEL
        ENDDO

      ENDDO

C LOWER BOUNDARY

      R0 = ((R + (1d5)*Z(1))/R)**2
      DEL15 = DELTAZ(1)
      TOP = R + (1d5)*Z(1) + 0.5D0*DEL15
      BOTTOM = R + (1d5)*Z(1) + 0.25D0*DEL15
      R15 = (TOP/BOTTOM)**2
      R1 = (R/BOTTOM)**2
      G1= 2.0D0*R1/DEL15

      DO M = 1 , NTRANS
        D15P = DEXP(0.5D0*DEL15/H(M,1))/DEL15
        D15M = DEXP(-0.5D0*DEL15/H(M,1))/DEL15
        VMAX1(M) = -D(M,1)/H(M,1)/R0

        A(M,1) = R15*2.0D0*D(M,1)*D15P/DEL15
        B(M,1) = R15*2.0D0*D(M,1)*D15M/DEL15
        C(M,1) = 2.0D0*R1/DEL15
      ENDDO

C UPPER BOUNDARY

      RN = ((R + (1d5)*Z(NALT))/R)**2
      DELN05 = DELTAZ(NALT-1)
      TOP = R + (1d5)*Z(NALT-1) + 0.5D0*DELN05
      BOTTOM = R + (1d5)*Z(NALT-1) + 0.75D0*DELN05
      RN05 = (TOP/BOTTOM)**2
      R = (R/BOTTOM)**2
      GN = 2.0D0*R/DELN05

      DO M = 1 , NTRANS
        DN05P = DEXP( 0.5D0*DELN05/H(M,NALT-1))/DELN05
        DN05M = DEXP(-0.5D0*DELN05/H(M,NALT-1))/DELN05
        VMAXN(M) = D(M,NALT-1)/H(M,NALT-1)/RN

        A(M,NALT) = RN05*2.0D0*D(M,NALT-1)*DN05P/DELN05
        B(M,NALT) = RN05*2.0D0*D(M,NALT-1)*DN05M/DELN05
        C(M,NALT) = 2.0D0*R/DELN05
      ENDDO
 


      RETURN
      END


C-------------------------------------------------------------------------------
C-------------------------------------------------------------------------------
C-------------------------------------------------------------------------------


      SUBROUTINE BMATRIX

      INCLUDE 'PARAMETERS.INC'

      COMMON

     &      /ONED14/ DELTNEW

     &      /ONED16/ DJNET(MB,MB,NALT)

     &      /ONED17/ NTRANS

     &      /ONED30/ B(MB,NALT)

     &      /ONED31/ BDELTA(MB,MB,NALT)

     &      /ONED40/ MTRANS(MB)


C CALCULATES JACOBIAN AND BDELTA MATRIX FOR LEVEL 2 TO NALT-1
C BOUNDARY TERMS ARE DONE SEPARATELY IN SUBROUTINE BNDRY

 

C COMPUTE JACOBIAN FOR SOURCE TERM

      CALL JACOB  ! RETURNS DJNET

C ZERO OUT ARRAY

c$doacross
c$& local(k,m2,m1) ,
c$& share(bdelta)

      DO K = 1 , NALT
        DO M2 = 1 , NTRANS
          DO M1 = 1 , NTRANS
            BDELTA(M1,M2,K) = 0.0D0
          ENDDO
        ENDDO
      ENDDO

C CALCULATE BDELTA MATRIX; LEVEL 2 TO NALT-1
C THE COMPUTATION OF BDELTA FOR ALTITUDE LEVELS 1 AND NALT IS DONE IN
C SUBROUTINE BNDRY
 
      IF( IMARCH .EQ. 0 ) THEN

c$doacross
c$& local(k,m,m2,m1) ,
c$& share(bdelta,b,deltnew,DJNET,mtrans)
 
        DO K = 2 , NALT-1

          DO M = 1 , NTRANS
            BDELTA(M,M,K) = B(M,K) + 1.0D0/DELTNEW ! DIAGONAL TERMS
          ENDDO

          DO M2 = 1 , NTRANS
            DO M1 = 1 , NTRANS
              BDELTA(M1,M2,K) = BDELTA(M1,M2,K) -
     &                          DJNET(M1,MTRANS(M2),K)
            ENDDO
          ENDDO

        ENDDO

      ELSE   ! NO TRANSPORT - COEFFICIENT MATRICES ARE NOT NEEDED

        DO K = 1 , NALT

          DO M = 1 , NTRANS
            BDELTA(M,M,K) = 1.0D0/DELTNEW ! DIAGONAL TERMS
          ENDDO

          DO M2 = 1 , NTRANS
            DO M1 = 1 , NTRANS
              BDELTA(M1,M2,K) = BDELTA(M1,M2,K) -
     &                          DJNET(M1,MTRANS(M2),K)
            ENDDO
          ENDDO

        ENDDO

      ENDIF


      RETURN
      END



C-------------------------------------------------------------------------------
C-------------------------------------------------------------------------------
C-------------------------------------------------------------------------------



      SUBROUTINE DIAG

      INCLUDE 'PARAMETERS.INC'

      COMMON

     &      /ONED9/  DCONC(MB,NALT)

     &      /ONED17/ NTRANS

     &      /ONED31/ BDELTA(MB,MB,NALT)

     &      /ONED39/ RESID(MB,NALT)

     &      /ONED40/ MTRANS(MB)



      DOUBLE PRECISION BTEMP(MB,MB) , RTEMP(MB)
      DIMENSION INDEX(MB)

      DO K = 1 , NALT
 
        DO M2 = 1 , NTRANS
          DO M1 = 1 , NTRANS
             BTEMP(M1,M2) = BDELTA(M1,M2,K)
          ENDDO
          RTEMP(M2) = RESID(M2,K)
        ENDDO

        CALL LU(BTEMP,NTRANS,INDEX)           ! LU DECOMPOSITION
        CALL CROUT(BTEMP,NTRANS,INDEX,RTEMP)  ! SOLVE MATRIX SYSTEM
                                              ! SOLUTION IS RETURNED AS RTEMP

        DO M = 1 , NTRANS
          DCONC(MTRANS(M),K) = RTEMP(M)
        ENDDO

      ENDDO


 

      RETURN
      END



C-------------------------------------------------------------------------------
C-------------------------------------------------------------------------------
C-------------------------------------------------------------------------------


      SUBROUTINE BNDRY

      INCLUDE 'PARAMETERS.INC'

      COMMON

     &      /ONED1/   CONCO(MB,NALT)

     &      /ONED7/   A(MB,NALT) , C(MB,NALT)

     &      /ONED22/  RNET(MB,NALT)

     &      /ONED12/  ILOWER(MB) , IUPPER(MB) ,
     &                XLOWER(MB) , XUPPER(MB)

     &      /ONED13/  ADELTA(MB,MB,NALT) , CDELTA(MB,MB,NALT)

     &      /ONED14/  DELTNEW

     &      /ONED16/  DJNET(MB,MB,NALT)

     &      /ONED17/  NTRANS

     &      /BLOCK18/ CONC(MB,NALT)

     &      /ONED24/  VMAX1(MB) , VMAXN(MB) , R0 , RN

     &      /ONED28/  DEN(NALT)

     &      /ONED30/  B(MB,NALT)

     &      /ONED31/  BDELTA(MB,MB,NALT)

     &      /ONED39/  RESID(MB,NALT)

     &      /ONED40/  MTRANS(MB)

     &      /ONED55/  MH2O , MH2OTRAN   ! H2O INDEX IN MOLECULE LISTING AND
C                                       ! TRANSPORTED SPECIES LISTING

     &      /ONED60/  BNDRYFLX(MB,NALT)

     &      /ONED61/  G1 , GN


C THIS SUBROUTINE SETS THE BOUNDARY CONDITIONS AT LEVEL 1 AND NALT

C DEFAULT OPTION IS IBND = 1  IF SECOND ORDER BOUNDARY IS NOT
C NEEDED SET IBND = 0

 

C IBND = 0 FIRST ORDER BOUNDRY = 1 SECOND ORDER

      IBND = 1

C ILOWER IUPPER  CODE FOR BOUNDARY CONDITION
C TYPE = 1   FIX NUMBER DENSITY (INCLUDES TYPE 5 MIXING RATIO)
C TYPE = 2   FIX VELOCITY (SIGN FOR FLOW IS POSITIVE FOR UP)
C TYPE = 3   MAXIMUM VELOCITY FOR LOSS ALLOWED BY SYSTEM
C TYPE = 4   FIX FLUX (SIGN FOR FLOW IS POSITIVE FOR UP)
C TYPE = 5   FIXED MIXING RATIO

C ZERO OUT BOUNDARY ARRAYS

      DO M2 = 1 , NTRANS
        DO M1 = 1 , NTRANS
          ADELTA(M1,M2,NALT) = 0.0D0  ! TOP
          CDELTA(M1,M2,1) = 0.0D0     ! BOTTOM
        ENDDO
      ENDDO

      DO M = 1 , NTRANS
        MVARY = MTRANS(M)

C SET UP LOWER BOUNDARY----------------------------------------


C Rearrange this expression to try to minimize roundoff/truncation errors
C **** THIS ASSUMES THAT IBND = 1 ! ****
C        RESID(M,1) = ( 0.75D0*RNET(M,1) + 0.25D0*RNET(M,2) +
C     &               (0.75D0*CONCO(MVARY,1) +
C     &                0.25D0*CONCO(MVARY,2))/DELTNEW )*IBND
C     &             - (B(M,1) + 0.75D0/DELTNEW*IBND)*CONC(MVARY,1) -
C     &               (-A(M,1) + 0.25D0/DELTNEW*IBND)*CONC(MVARY,2)

        DC1 = CONCO(MVARY,1) - CONC(MVARY,1)
        DC2 = CONCO(MVARY,2) - CONC(MVARY,2)
        DCFLUX = A(M,1)*CONC(MVARY,2) - B(M,1)*CONC(MVARY,1)
        DC0 = 0.25D0*( 3.0D0*RNET(M,1) + RNET(M,2) )

        SUM = 3.0D0*DC1 + DC2
        SUM = 0.25D0*SUM/DELTNEW
        SUM = SUM + DCFLUX
        SUM = SUM + DC0
        RESID(M,1) = SUM

C TYPE 1 PRESCRIBED NUMBER DENSITY

        IF (ILOWER(MVARY) .EQ. 1 ) THEN

          BNDRYFLX(M,1) = -RESID(M,1)/G1
          B(M,1) = B(M,2)
          C(M,1) = 0.0D0
          RESID(M,1) = 0.0D0
          CONC(MVARY,1) = XLOWER(MVARY)
          BDELTA(M,M,1) = B(M,1) + 0.75D0/DELTNEW*IBND

          DO M2 = 1 , NTRANS
            CDELTA(M,M2,1) = 0.0D0
          ENDDO

        ENDIF

C TYPE 2 PRESCRIBED VELOCITY

        IF (ILOWER(MVARY) .EQ. 2 ) THEN

          VSTAR1 = C(M,1)*R0*XLOWER(MVARY)
          BDELTA(M,M,1) = B(M,1) + 0.75D0/DELTNEW*IBND - VSTAR1
          CDELTA(M,M,1) = -A(M,1) + 0.25D0/DELTNEW*IBND

          DO M2 = 1 , NTRANS
            BDELTA(M,M2,1) = BDELTA(M,M2,1) -
     &                       0.75D0*DJNET(M,MTRANS(M2),1)*IBND
            CDELTA(M,M2,1) = CDELTA(M,M2,1) -
     &                       0.25D0*DJNET(M,MTRANS(M2),2)*IBND
          ENDDO

          RESID(M,1) = RESID(M,1) + VSTAR1*CONC(MVARY,1)
          BNDRYFLX(M,1) = VSTAR1*CONC(MVARY,1)/G1

        ENDIF

C TYPE 3 MAXIMUM VELOCITY

        IF (ILOWER(MVARY) .EQ. 3 ) THEN

          VSTAR1 = C(M,1)*R0*VMAX1(M)
          BDELTA(M,M,1) = B(M,1) + 0.75D0/DELTNEW*IBND - VSTAR1
          CDELTA(M,M,1) = -A(M,1) + 0.25D0/DELTNEW*IBND

          DO M2 = 1 , NTRANS
            BDELTA(M,M2,1) = BDELTA(M,M2,1) -
     &                       0.75D0*DJNET(M,MTRANS(M2),1)*IBND
            CDELTA(M,M2,1) = CDELTA(M,M2,1) -
     &                       0.25D0*DJNET(M,MTRANS(M2),2)*IBND
          ENDDO

          RESID(M,1) = RESID(M,1) + VSTAR1*CONC(MVARY,1)
          BNDRYFLX(M,1) = VSTAR1*CONC(MVARY,1)/G1

        ENDIF


C TYPE 4 PRESCRIBED FLUX

        IF (ILOWER(MVARY) .EQ. 4 ) THEN

          RESID(M,1) = RESID(M,1) + C(M,1)*XLOWER(MVARY)
          BDELTA(M,M,1) = B(M,1) + 0.75D0/DELTNEW*IBND
          CDELTA(M,M,1) = -A(M,1) + 0.25D0/DELTNEW*IBND

          DO M2 = 1 , NTRANS
            BDELTA(M,M2,1) = BDELTA(M,M2,1) -
     &                       0.75D0*DJNET(M,MTRANS(M2),1)*IBND
            CDELTA(M,M2,1) = CDELTA(M,M2,1) -
     &                       0.25D0*DJNET(M,MTRANS(M2),2)*IBND
          ENDDO

          BNDRYFLX(M,1) = XLOWER(MVARY)

        ENDIF

C TYPE 5 PRESCRIBED MIXING RATIO

        IF (ILOWER(MVARY) .EQ. 5 ) THEN

          BNDRYFLX(M,1) = -RESID(M,1)/G1
          B(M,1) = B(M,2)
          C(M,1) = 0.0D0
          RESID(M,1) = 0.0D0
          CONC(MVARY,1) = XLOWER(MVARY)*DEN(1)
          BDELTA(M,M,1) = B(M,1) + 0.75D0/DELTNEW*IBND

          DO M2 = 1 , NTRANS
            CDELTA(M,M2,1) = 0.0D0
          ENDDO

        ENDIF

C SET H2O BELOW 16 KM
 
        IF( MH2OTRAN .NE. 0 ) THEN
          B(MH2OTRAN,9) = B(MH2OTRAN,10)
          C(MH2OTRAN,9) = 0.0D0
          DO K = 1 , 9
            RESID(MH2OTRAN,K) = 0.0D0
            CONC(MH2O,K) = CONCO(MH2O,K)
          ENDDO
          BDELTA(MH2OTRAN,MH2OTRAN,9) = B(MH2OTRAN,9) +
     &                                  0.75D0/DELTNEW*IBND
          DO M2 = 1 , NTRANS
            CDELTA(MH2OTRAN,M2,9) = 0.0D0
          ENDDO
        ENDIF


C SET UP UPPER BOUNDARY ---------------------------------------



        RESID(M,NALT) = ( 0.75D0*RNET(M,NALT) + 0.25D0*RNET(M,NALT-1) +
     &                (0.75D0*CONCO(MVARY,NALT) +
     &                 0.25D0*CONCO(MVARY,NALT-1))/DELTNEW )*IBND
     &                - (-B(M,NALT) +
     &                 0.25D0/DELTNEW*IBND)*CONC(MVARY,NALT-1) -
     &                (A(M,NALT) + 0.75D0/DELTNEW*IBND)*CONC(MVARY,NALT)


C TYPE 1 PRESCRIBED NUMBER DENSITY

        IF( IUPPER(MVARY) .EQ. 1 ) THEN

          BNDRYFLX(M,NALT) = RESID(M,NALT)/GN
          A(M,NALT) = 0.0D0
          B(M,NALT) = B(M,NALT-1)
          RESID(M,NALT) = 0.0D0
          CONC(MVARY,NALT) = XUPPER(MVARY)
          BDELTA(M,M,NALT) = A(M,NALT) + 0.75D0/DELTNEW*IBND
          DO M2 = 1 , NTRANS
            ADELTA(M,M2,NALT) = 0.0D0
          ENDDO

        ENDIF

C TYPE 2 PRESCRIBED VELOCITY

        IF( IUPPER(MVARY) .EQ. 2 ) THEN

          VSTARN = C(M,NALT)*RN*XUPPER(MVARY)
          BDELTA(M,M,NALT) = A(M,NALT) + 0.75D0/DELTNEW*IBND + VSTARN
          ADELTA(M,M,NALT) = -B(M,NALT) + 0.25D0/DELTNEW*IBND
          DO M2 = 1 , NTRANS
            ADELTA(M,M2,NALT) = ADELTA(M,M2,NALT) -
     &                          0.25D0*DJNET(M,MTRANS(M2),NALT-1)*IBND
            BDELTA(M,M2,NALT) = BDELTA(M,M2,NALT) -
     &                          0.75D0*DJNET(M,MTRANS(M2),NALT)*IBND
          ENDDO
          RESID(M,NALT) = RESID(M,NALT) - VSTARN*CONC(MVARY,NALT)
          BNDRYFLX(M,NALT) = VSTARN*CONC(MVARY,NALT)/GN

        ENDIF

C TYPE 3 MAXIMUM VELOCITY

        IF( IUPPER(MVARY) .EQ. 3 ) THEN

          VSTARN = C(M,NALT)*RN*VMAXN(M)
          BDELTA(M,M,NALT) = A(M,NALT) + 0.75D0/DELTNEW*IBND + VSTARN
          ADELTA(M,M,NALT) = -B(M,NALT) + 0.25D0/DELTNEW*IBND
          DO M2 = 1 , NTRANS
            ADELTA(M,M2,NALT) = ADELTA(M,M2,NALT) -
     &                          0.25D0*DJNET(M,MTRANS(M2),NALT-1)*IBND
            BDELTA(M,M2,NALT) = BDELTA(M,M2,NALT) -
     &                          0.75D0*DJNET(M,MTRANS(M2),NALT)*IBND
          ENDDO
          RESID(M,NALT) = RESID(M,NALT) - VSTARN*CONC(MVARY,NALT)
          BNDRYFLX(M,NALT) = VSTARN*CONC(MVARY,NALT)/GN

        ENDIF

C TYPE 4 PRESCRIBED FLUX

        IF( IUPPER(MVARY) .EQ. 4 ) THEN

          RESID(M,NALT) = RESID(M,NALT) - C(M,NALT)*XUPPER(MVARY)
          BDELTA(M,M,NALT) = A(M,NALT) + 0.75D0/DELTNEW*IBND
          ADELTA(M,M,NALT) = -B(M,NALT) + 0.25D0/DELTNEW*IBND
          DO M2 = 1 , NTRANS
            ADELTA(M,M2,NALT) = ADELTA(M,M2,NALT) -
     &                          0.25D0*DJNET(M,MTRANS(M2),NALT-1)*IBND
            BDELTA(M,M2,NALT) = BDELTA(M,M2,NALT) -
     &                          0.75D0*DJNET(M,MTRANS(M2),NALT)*IBND
          ENDDO
          BNDRYFLX(M,NALT) = XUPPER(MVARY)

        ENDIF

C TYPE 5 PRESCRIBED MIXING RATIO

        IF( IUPPER(MVARY) .EQ. 5 ) THEN

          BNDRYFLX(M,NALT) = RESID(M,NALT)/GN
          A(M,NALT) = 0.0D0
          B(M,NALT) = B(M,NALT-1)
          RESID(M,NALT) = 0.0D0
          CONC(MVARY,NALT) = XUPPER(MVARY)*DEN(NALT)
          BDELTA(M,M,NALT) = A(M,NALT) + 0.75D0/DELTNEW*IBND
          DO M2 = 1 , NTRANS
            ADELTA(M,M2,NALT) = 0.0D0
          ENDDO

        ENDIF

      ENDDO   !  END LOOP OVER M

      RETURN
      END