subroutine moldiffcoeff_red(nq,dij,gcmind,ncompdiff2)

       use diffusion_h
       use tracer_h, only: noms, mmol
       IMPLICIT NONE
c=======================================================================
c   subject:
c   --------
c   Computing molecular diffusion coefficients
c   following Nair 94 (pg 131)
c   author:  MAC 2002
c   ------
c
c=======================================================================
#include "dimensions.h"

c-----------------------------------------------------------------------
c    Input/Output
c    ------------
c       integer,parameter :: ncompmoldiff = 12
	integer ncompdiff2
	integer gcmind(ncompdiff2)
      real dij(ncompdiff2,ncompdiff2)

c    Local variables:
c    ---------------
      INTEGER nq, n, nn, i,iq
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     tracer numbering in the molecular diffusion
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!      real :: dijh2co,dijh2n2,dijh2co2,dijh2o2,dijho,dijref,dijch4n2
       real :: dijref,dijch4n2,dijch4co  !TB16
!       integer :: g_h2,g_h,g_o
       integer :: g_ch4,g_n2   ! TB16
!	integer :: i_h2,i_h,i_o
!      integer,parameter :: i_ch4 = 1

      real dnh
      logical,save :: firstcall=.true.
      logical,parameter :: outputcoeffs=.false. ! to output 'coeffs.dat' file,
                                                ! set outputcoeffs=.true.

! Initializations at first call (and some sanity checks)
      if (firstcall) then
        ! identify the indexes of the tracers we'll need
!        g_n2=igcm_n2
!        if (g_n2.eq.0) then
!          write(*,*) "moldiffcoeff: Error; no N2 tracer !!!"
!          stop
!        endif

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c    fill array to relate local indexes to gcm indexes
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!        gcmind(i_n2)  =   g_n2
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        firstcall= .false.
      endif ! of if (firstcall)

       dijch4n2 = 0.000018524
       dijch4co = 0.000018524

! find h2, h and o index in gcm
! these species are used to define the diffusion coefficients

	do n=1,nq
	if (noms(n) .eq. 'ch4_gas') g_ch4=n
	if (noms(n) .eq. 'n2') g_n2=n
	enddo

       print*,'moldiffcoeff_red: COEFF CALC'


      do n=1,ncompdiff2
	dijref=0.
	if (noms(gcmind(n)) .eq. 'n2') dijref=dijch4n2
	if (noms(gcmind(n)) .eq. 'co_gas') dijref=dijch4co

        if (dijref .gt. 0.0) then
          do nn=n,ncompdiff2
            dij(nn,n)=dijref
     &                  *sqrt(mmol(g_ch4)/mmol(gcmind(nn)))
            if(n.eq.nn) dij(nn,n)=1.0
            if(mmol(gcmind(n)).eq.mmol(gcmind(nn))) dij(nn,n)=1.0
            dij(n,nn)=dij(nn,n)
          enddo
        endif

        if (dijref .eq. 0.0) then
	dijref=dijch4n2
          dnh=dijref*sqrt(mmol(g_n2)/mmol(gcmind(n)))
          do nn=n,ncompdiff2
            dij(nn,n)=dnh*sqrt(mmol(g_ch4)/mmol(gcmind(nn)))
            if(n.eq.nn) dij(nn,n)=1.0
            dij(n,nn)=dij(nn,n)
          enddo
        endif

      enddo
      print*, 'TB16: dij=',dij(1,:)
      print*, 'TB16: dij=',dij(2,:)
      print*, 'TB16: dij=',dij(3,:)


      return
      end