do 3000 k2=k1+1,np+1

           call TABLUP(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,

           call TABLUP(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,

           call TABLUP(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,

        call H2OKDIS(ib,m,n,np,k2-1,fkw,gkw,ne,h2oexp,conexp,

          call TABLUP(k1,k2,m,n,np,nx,nc,nt,sco3,scopre,scotem,

          call CO2KDIS(m,n,np,k2-1,co2exp,tco2,trant)

          call TABLUP(k1,k2,m,n,np,nx,no,nt,sco3,scopre,scotem,

          call N2OKDIS(ib,m,n,np,k2-1,n2oexp,tn2o,trant)

          call CH4KDIS(ib,m,n,np,k2-1,ch4exp,tch4,trant)

          call COMKDIS(ib,m,n,np,k2-1,comexp,tcom,trant)

          call CFCKDIS(m,n,np,k2-1,f11exp,tf11,trant)

          call CFCKDIS(m,n,np,k2-1,f12exp,tf12,trant)

          call CFCKDIS(m,n,np,k2-1,f22exp,tf22,trant)

          call B10KDIS(m,n,np,k2-1,h2oexp,conexp,co2exp,n2oexp

         ff=0.5+(0.3739+(0.0076+0.1185*asyal(i,j,k2-1,ib))

     *      *asyal(i,j,k2-1,ib))*asyal(i,j,k2-1,ib)

         taux=taual(i,j,k2-1,ib)*(1.-ssaal(i,j,k2-1,ib)*ff)

        if( k2 .le. ict ) then

            clrhi(i,j)=min(clr(i,j,k2-1),clrhi(i,j))

        elseif( k2 .gt. ict .and. k2 .le. icb ) then

            clrmd(i,j)=min(clr(i,j,k2-1),clrmd(i,j))

        elseif( k2 .gt. icb ) then

            clrlw(i,j)=min(clr(i,j,k2-1),clrlw(i,j))

       if (k2 .eq. k1+1) then

          flcd(i,j,k2)=flcd(i,j,k2)+blayer(i,j,k1)

          flxd(i,j,k2)=flxd(i,j,k2)+blayer(i,j,k1)

         xx=trant(i,j)*dlayer(i,j,k2)

         flcd(i,j,k2) =flcd(i,j,k2)+xx

         flxd(i,j,k2) =flxd(i,j,k2)+xx*fclr(i,j)

     *                   *dlayer(i,j,k2)

           rflx(i,j,k2) = rflx(i,j,k2)+trant(i,j)*fclr(i,j)

     *                   *dlayer(i,j,k2)

           rflc(i,j,k2) = rflc(i,j,k2)+trant(i,j)

      subroutine TABLUP(k1,k2,m,n,np,nx,nh,nt,sabs,spre,stem,w1,p1,

      integer k1,k2,m,n,np,nx,nh,nt,i,j

        x1=sabs(i,j,k2)-sabs(i,j,k1)

        x2=(spre(i,j,k2)-spre(i,j,k1))/x1

        x3=(stem(i,j,k2)-stem(i,j,k1))/x1

      integer i,j,k,ip,iw,it,ib,ik,iq,isb,k1,k2

      integer k1,k2,nymd1,nhms1,nymd2,nhms2,i

      call TIME_BOUND  ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, k1,k2 )

      zoch(i) = vgz0(k2,ityp(i))*facp + vgz0(k1,ityp(i))*facm

      rdcs(i) = vgrd(k2,ityp(i))*facp + vgrd(k1,ityp(i))*facm

      rootl   = vgrt(k2,ityp(i))*facp + vgrt(k1,ityp(i))*facm

      d       = vgdd(k2,ityp(i))*facp + vgdd(k1,ityp(i))*facm

      integer i,midmon,midm,midp,id,k1,k2,kk1,kk2

      INTEGER i1, j1, k1, i2, j2, k2

      k2 = k1 + 1

      k2 = MIN( MAX( k2, 1 ), Nr )

     &   +   ddz   *( ( (1.-ddx)*(1.-ddy)*var(i1,j1,k2,bi,bj)

     &                +     ddx * ddy    *var(i2,j2,k2,bi,bj) )

     &              + (     ddx*(1.-ddy) *var(i2,j1,k2,bi,bj)

     &                + (1.-ddx)* ddy    *var(i1,j2,k2,bi,bj) ) )

      INTEGER km, k1, k2

        k2 = myNr

        k2 = kLoc

        DO k=k1,k2

                mda(k2) = (md(k2)*dz(k2)+md(k2+1)*dz(k2+1))/

     *                    (dz(k2)+dz(k2+1))

                wda(k2) = (wd(k2)*dz(k2)+wd(k2+1)*dz(k2+1))/

     *                    (dz(k2)+dz(k2+1))

                tda(k2) = td(k2)

                sda(k2) = sd(k2)

                taa(k2) = ttemp(k2+1)

                saa(k2) = stemp(k2+1)

                 do k2=k+1,Maxdepth-1

                   ttemp(k2) = ttemp(k2) +

     *              (mda(k2-1)*(tda(k2-1)-taa(k2-1))-

     *              mda(k2)*(tda(k2)-taa(k2)))

     *              *dt*recip_drF(k2)

                  stemp(k2) = stemp(k2) +

     *              (mda(k2-1)*(sda(k2-1)-saa(k2-1))-

     *              mda(k2)*(sda(k2)-saa(k2)))

     *              *dt*recip_drF(k2)

                do k2=1,maxdepth-1

                  taa(k2) = ttemp(k2+1)

                  saa(k2) = stemp(k2+1)

          do k2=k,maxdepth

            convadj(i,j,k2) = convadj(i,j,k2) + (ttemp(k2)-

     *          ta(i,k2))

            sa(i,k2) = stemp(k2)

            ta(i,k2) = ttemp(k2)

       DO k2=k,KMax-1

        D1=state1( Pd(k2,2), Pd(k2,1),i,j,k2+1,bi,bj,myThid)

        D2=state1( tracerEnv(k2+1,2), tracerEnv(k2+1,1),

     &                                i,j,k2+1,bi,bj,myThid)

        De(k2+1)=D2

        IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN

         dz2=DrF(k2+1)

         newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)*

         PlumeEntrainment(k2+1) = newflux/StartingFlux

     &          'S/R OPPS_CALC: Plume entrained to zero at level ', k2

          maxdepth = k2

          pmix=(dz1*tracerEnv(k2,ktr)+dz2*tracerEnv(k2+1,ktr))

          Pd(k2+1,ktr)=Pd(k2,ktr)

     &         - entrainrate*(pmix - Pd(k2,ktr))

         Dd(k2+1)=state1(Pd(k2+1,2),Pd(k2+1,1),i,j,k2+1,bi,bj,myThid)

     &    'S/R OPPS_CALC: Dd,De,entr,k ',Dd(k2),De(k2),entrainrate,k2

     &        (dz1*(Dd(k2)-De(k2))/De(k2)

     &        +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1))

          maxdepth = k2

     &     'S/R OPPS_CALC: Plume velocity went to zero at level ', k2

     &          (Dd(k2)-De(k2))/De(k2), (Dd(k2+1)-De(k2+1))/De(k2+1)

         Wd(k2+1)=sqrt(wsqr)

         radius=sqrt(newflux/(Wd(k2)*Dd(k2)))

         maxdepth=k2

        DO k2=k+1,maxDepth

         Md(k2)=Md(k)*PlumeEntrainment(k2)

     &         'S/R OPPS_CALC: Md, Wd, and  k are ',Md(k2),Wd(k2),k2

        do k2=k,maxDepth-1

         IF ( Wd(K2) .NE. 0. _d 0 ) dt = min(dt,drF(k2)/Wd(k2))

         mda(k2) = (md(k2)*drF(k2)+md(k2+1)*drF(k2+1))/

     *        (drF(k2)+drF(k2+1))

         wda(k2) = (wd(k2)*drF(k2)+wd(k2+1)*drF(k2+1))/

     *        (drF(k2)+drF(k2+1))

          Pda(k2,ktr) = Pd(k2,ktr)

          Paa(k2,ktr) = tracerEnv(k2+1,ktr)

         DO k2=k+1,kmx

          k2m1 = max(k,k2-1)

          k2p1 = max(k2+1,maxDepth)

            tracerEnv(k2,ktr) = tracerEnv(k2,ktr) +

     &           -mda(k2)  *(Pda(k2,ktr)  -Paa(k2,ktr))  )

     &           *dt*recip_drF(k2)

          DO k2=1,kmx

           paa(k2,ktr) = tracerEnv(k2+1,ktr)

      integer is,ie,k2

          DO k2=k,NumGridPoints-1

            D1=state1(Sd(k2),Td(k2),k2+1)

            D2=state1(stemp(k2+1),ttemp(k2+1),k2+1)

            De(k2+1)=D2

            IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN

                 dz2=GridThickness(k2+1)

                 newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)*0.50*

                 PlumeEntrainment(k2+1) = newflux/StartingFlux

                      write(0,*)"Plume entrained to zero at ",k2

                     maxdepth = k2

                 smix=(dz1*stemp(k2)+dz2*stemp(k2+1))/(dz1+dz2)

                 thmix=(dz1*ttemp(k2)+dz2*ttemp(k2+1))/(dz1+dz2)

                  sd(k2+1)=sd(k2) - entrainrate*(smix - sd(k2))

                  td(k2+1)=td(k2) - entrainrate*(thmix - td(k2))

                 Dd(k2+1)=state1(Sd(k2+1),Td(k2+1),k2+1)

                   write(0,*)"Dd,De,entr,k ",Dd(k2),De(k2),entrainrate,k2

     .             (dz1*(Dd(k2)-De(k2))/De(k2)

     .             +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1))

                     maxdepth = k2

                      write(0,*)"Plume velocity went to zero at ",k2

                 Wd(k2+1)=sqrt(wsqr)

                 radius=sqrt(newflux/(Wd(k2)*Dd(k2)))

                 maxdepth=k2

             DO k2=k+1,MaxDepth

               Md(k2)=Md(k)*PlumeEntrainment(k2)

                 write(0,*)"Md, Wd, and  k are ",Md(k2),Wd(k2),k2

             do k2=k,maxdepth-1

                dt = min(dt,dz(k2)/wd(k2))

            k2 = MIN(nlevels(ld),k1+19)

            WRITE(msgBuf,'(A,20F5.0)') suffix, (levs(k,ld),k=k1,k2)

      INTEGER i, j, k, k1, k2, mm, kLev

           k2 = MIN(nlevels(ld),k1+9)

           WRITE(dUnit,'(A,1P10E10.3)')' interp:', (levs(k,ld),k=k1,k2)

           k2 = MIN(nlevels(ld),k1+24)

           WRITE(dUnit,'(A,25I4)')' levels:',(NINT(levs(k,ld)),k=k1,k2)

      INTEGER j, k, k1, k2, l