      SUBROUTINE LIMBCALC(MR,ML,MGG,MGS,MGB,MT,MDET,NCH,NTB,NPT,NTT,NST,
     &NR,NL,NGG,IGG,NG,IDG,NGB,IDGB,RE,WAVE,ISC,IRFN,IRAD,ICOR,ITAB,ZR,
     &PR,TR,TGR,Q,QGG,ZA,ZT,PT,TT,GR,QMIX,QG,QMIXB,QGB,RLT,RUT,TEARTH,
     &EMR,RMR,ESR,RSR,BSR,NDX1,NDX2,TBL,POUT,SOUT,TOUT,jsl,jel,
     & tau)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C FORWARD IS A DRIVER ROUTINE FOR PATH AND MEGA. THE USE OF FORWARD
C  IS FAIRLY SELF EXPLANATORY. MOST OF THE VARIABLES IN THIS ROUTINE ARE
C  DEFINED IN PATH AND MEGA. NOTE THAT:
C   Q AND QGG ARE THE GAS VOLUME MIXING RATIO OR ABSORPTION(/KM) AND
C   HORIZONTAL GRADIENT ARRAYS AT INPUT ATMOSPHERE ALTITUDES.
C   QMIX AND QG ARE INTERPOLATED GAS MIXING RATIO OR ABSORPTION(/KM) AND
C   HORIZONTAL GRADIENT ARRAYS AT THE TRUE TANGENT HEIGHTS OUTPUT BY
C   PATH.
C
C  COPYRIGHT 1988 BY G & A TECHNICAL SOFTWARE INC.
C                    P.O. BOX 1885
C PARAMETER LIST:
C   NCH---- NUMBER OF CHANNELS
C   ITAB--- DETECTOR NUMBERS CORRESPONDING TO EACH RAY
C   NG----- NUMBER OF GASES TO BE SIMULATED
C   IRAD--- CALCULATION TECHNICS SELECTIONS. AIR(1) OR GAS(2) EGA;
C           AIR(3) OR GAS(4) CGA; AIR(5) OR GAS(6) AVE.
C   ICOR--- ORDER OF GASES FOR RADIANCE CORRECTION.
C                    GRAFTON, VA 23692
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PARAMETER(MAXG=50)
      LOGICAL DBGIN,IBIN
      real*8 emr,esr,rmr,rsr,   tau
      DIMENSION TBL(NTB,MDET,MT),POUT(NPT),TOUT(NTT,MT),
     &SOUT(NST,MGB,MT), tau(mr,*)
      DIMENSION IRAD(MGB,MT),ICOR(MGB,MT),ITAB(ML)
      DIMENSION ZR(MR),PR(MR),TR(MR),TGR(2,MR),Q(MR,MGG),QGG(2,MR,MGG),
     &ZA(ML),ZT(ML),PT(ML),TT(ML),GR(2,ML),RLT(ML,MGB,MT),
     &RUT(ML,MGB,MT),IGG(MGG)
      DIMENSION EMR(ML,MGB,MT),RMR(ML,MGB,MT),ESR(ML,MGB,MT),
     &RSR(ML,MGB,MT),BSR(ML,MGB,MT),QG(2,ML,MGS),IDG(MGS),QMIX(ML,MGS),
     &NGB(MT),IDGB(MGB,MT),QMIXB(ML,MGB),QGB(2,ML,MGB),INDEX(MAXG)
      DATA ALBEDO/0.0/
c      JSL=1
c      JEL=NL
c        print*, "limcalc  1", igg,   re

c      call kp_inita(dbgin,npt, pout)

      DO IT=1,NCH
c         CALL SETINDX(NG,NGB(IT),IDG,IDGB(1,IT),INDEX)
c          print*, ng, ngb(it), idg ,
c     &    (idgb(ig,it),index(ig),ig=1,ngb(it) ) 
c       print*, "limcalc  2", NGG
c         CALL SETQMIXB(ML,NL,NG,NGB(IT),INDEX,QMIX,QG,QMIXB,QGB)
c       print*, "limcalc  3", NGG, mgg
c         write(*,*) 'Begin radiation transfer'
c           print*, ng, idg, ngg, igg, isc, re, wave

c       print*,'calling mp_mass'

c          isc = 1

c      print*,mr,ml,mgs,mgg,npt,nl,nr,irfn,jsl,jel,ng,idg,ngg,
c     &   igg,isc,re,wave

c       do i = 1,ngg
c         print*, igg(i)
c      enddo
c      do i=1,nl
c        print*,za(i), zr(i), pr(i), tr(i), tgr(1,i), tgr(2,i),
c     &     q(i,1), q(i,2), q(i,3), q(i,4), q(i,5) 

c       enddo
c       stop

c      do i=1,nl
c       do j=1,ngg
c        print*, zr(i),igg(j),  qgg(1,i,j), qgg(2,i,j) 
c       enddo
c       enddo

c          ngg=4
c          ng = 4
c       stop
         CALL mp_mass(mr,ml,mgs,mgg,npt,nl,nr,irfn,1,nl,ng,idg,ngg,
     &   igg,isc,re,wave,0,za,zr,pr,tr,tgr,q,qgg,zt,pt,tt,gr,tearth,
     &   qmix,qg,tbl(1,1,IT),pout)

c         do i=1,nl
c           print*, zt(i), pt(i), tt(i) 
c         enddo

         CALL SETINDX(NG,NGB(IT),IDG,IDGB(1,IT),INDEX)
c          print*, ng, ngb(it), idg ,
c     &    (idgb(ig,it),index(ig),ig=1,ngb(it) ) 
c       print*, "limcalc  2", NGG
         CALL SETQMIXB(ML,NL,NG,NGB(IT),INDEX,QMIX,QG,QMIXB,QGB)

         
c           do il = jsl,jel
c              write(*,*) zt(il),pt(il),tt(il), qmixb(il,1), 
c     &        rut(il,1,1), rlt(il,1,1)
c           enddo
         CALL MP_INITBSR(ML,MGB,BSR(1,1,IT))
c        print*, ' alling mega '
         CALL MP_MEGA(NTB,NST,NTT,MGB,ML,NL,JSL,JEL,NDX1,NDX2,
     &   ITAB,NGB(IT),IDGB(1,IT),QMIXB,QGB,RLT(1,1,IT),RUT(1,1,IT),
     &   ALBEDO,TEARTH,IRAD(1,IT),ICOR(1,IT),EMR(1,1,IT),RMR(1,1,IT),
     &   ESR(1,1,IT),RSR(1,1,IT),BSR(1,1,IT),TBL(1,1,IT),POUT,
     &   SOUT(1,1,IT),TOUT(1,IT))
c        stop 'finished mega'

c         write(*,*) 'End radiation transfer'
         do il = jsl,jel
           tau(il,1) =  1.0d0 - emr(il,1,1)
         enddo
c         print*, jsl, jel, tau(jel+1,1)
      ENDDO
c        stop
      RETURN
      END
