c    @(#)trcphys.f	1.6
      subroutine trcphys(itcall)

      include'mladv.inc'
      include'trc.inc'

      common/parm/kz,m,mm1,mm2,mp1,msq,tim,tmax,tstart,hy,rea,dt,
     +   xbasin,yy(mx0),icalc,kzp1,kzm1,iftop,ifbot,iftvv,ifrel,theta,
     +   ifsbl,ifeva

      dimension ystart(maxtrace+1)

      external trcrhs
      external rkqc

      data tmlold/-9999./

      if ((itphys .eq. 2) .or. (itphys.eq. 3)) then
          icurtrace = 1
          do 25 ip=1,msq

             call findkml(hdpth(ip,mlcur),zml,l,kml)

             if (itcall .eq. 1) out(ip) = 0.

             if (itphys .eq. 2) then
                krem = kml
               else
                krem = keuph
              endif

             do 20 k=1,krem

                xn = trace(ip,k,mlcur,1) + trcqgofk(ip,k)
                if (xn .lt. 0.0) then
                   xn = 0.0
                  endif
                if (k .eq. kml) then
                   z1 = zml(kml) + dzml(1) / 2.0
                   dz1 = hdpth(ip,mlcur) + z1
                   xn = xn * dz1 / dzml(1)
                 endif
                if (k .eq. kml+1) then
                   z1 = -hdpth(ip,mlcur) - zml(kml+1) 
                   dz1 = z1 + dzml(1) / 2.0
                   xn = xn * dz1 / dzml(1)
                 endif
                trace(ip,l+1,mlcur,1) = trace(ip,l+1,mlcur,1) + xn
                if (itcall .eq. 1)  then
                   out(ip) = xn
                  else
                   out(ip) = 0.5 * (out(ip) + xn)
                 endif
                trace(ip,k,mlcur,1) = -trcqgofk(ip,k)

  20          continue

             if (ifrest .eq. 1) then
                do 22 k=krem+1,l
                   trace(ip,k,mlcur,1) = 0.
  22             continue
              endif

  25      continue

          return

       endif

c
c
c

      if (kc .eq. 0.0) then
         call light(1)
        endif

      t1 = 0.0
      t2 = qgts
      st1 = qgts / rkh1
      rkmin = qgts / rkmax

      if (mltime .ne. tmlold) then
         iflxst = 1
         tmlold = mltime
        else
         iflxst = 0
       endif

      do 5 ip=1,msq

         if (kc .ne. 0.0) then
            call light(ip)
           endif

         ystart(5) = trace(ip,l+1,mlcur,izoo)

         do 10 k=1,l

            do 500 icurtrace=1,ntrace
               ystart(icurtrace) = trace(ip,k,mlcur,icurtrace) +
     .                                              trcqgofk(ip,k)
 500         continue

            enow = epar(k)

            call odeint(ystart,ntrace+1,t1,t2,rkeps,st1,rkmin,
     .                                     nok,nbad,trcrhs,rkqc)

            if ((dtrflx .gt. 0.) .and. (iflxst .eq. 1)) then
               call savflx(ystart,ip,k)
              endif

            do 510 icurtrace=1,ntrace
               trace(ip,k,mlcur,icurtrace) = ystart(icurtrace) -
     .                                       trcqgofk(ip,k)
 510         continue

  10       continue

         trace(ip,l+1,mlcur,izoo) = ystart(5)

   5   continue

      return
      end
