c
c    @(#)trcinit.f	1.5
c
      subroutine trcinit

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

      parameter(npmax=100)
      dimension ztr(npmax),trcic(npmax,maxtrace),trcicz(lmax,maxtrace)

      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
      common /threed/hz(k0),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     .               amat(k0,k0)
      common/fields/vort(mn0,k0),stream(mn0,k0),prvjac(mn0,k0),
     +     topden(mn0),botden(mn0),relief(mn0),prvtop(mn0),
     +     prvbot(mn0),topfct(k0),botfct(k0),pvort(mn0)
      common /extpar/ rlat0, rlng0, t0, v0, dhor, ht, time0,
     +                r0,f0star
      common /vert/ rho0, tsurf, xc(k01), salsur, eta(k01)

      common/sbldef/qsw0,qsurf0,wflux0,u10m,angu10m,hdpth0

      real nitr

c ----------------------------------------------------------------------

      if (ifsbl .eq. 1) then

      do 1 k=1,l
         do 2 ip=1,msq

            if ((idinml.eq.0).and.(abs(zml(k)).le.hdpth(ip,mlcur))) then
               icurtrace = idin
               trace(ip,k,mlcur,idin) = -trcqgofk(ip,k)
              else
               trace(ip,k,mlcur,idin) = 0.
             endif

            trace(ip,k,mlcur,ipon) = 0.

            if (ibio .eq. ipnz1) then
               trace(ip,k,mlcur,icq) = 0.
               trace(ip,k,mlcur,izoo) = 0.
              endif

  2       continue
  1    continue

      else if (ifsbl .eq. 3) then

      dsigdb = -1000.0 * rho0 / 9.8
      dbdz = (f0star*dhor/ht)**2./(scsig*sigz(1))
      dsigdz = dsigdb * dbdz

      call findkml(hdpth0,zml,lmax,kml)
      icurtrace = idin

      do 3 k=1,l
         do 4 ip=1,msq

            if (k .le. kml) then
               trace(ip,k,mlcur,idin) = surno3 - trcqgofk(ip,k)
              else
               trace(ip,k,mlcur,idin) = 0.
             endif

            trace(ip,k,mlcur,ipon) = 0.

            trace(ip,k,mlcur,izoo) = 0.

            if (ibio .eq. ipnz1) then
               trace(ip,k,mlcur,icq) = 0.
              endif

  4       continue
  3    continue

      else if (ifsbl .eq. 4) then

      if (irdtrc .eq. 1) then
         open(unit=itrfile,file=trfile,form='formatted',status='old')
         i = 1
  8      read(itrfile,*,end=7)ztr(i),(trcic(i,j),j=1,ntrace)
         i = i + 1
         goto 8
  7      ntrd = i - 1
         do 21 j=1,ntrace
            do 20 i=1,l+1
               trcicz(i,j) = find(-zml(i),ztr,trcic(1,j),ntrd)
 20          continue
 21       continue
         pno31 = find(-zml(k),ztr,trcic(1,2),ntrd)
         zoo1  = find(-zml(k),ztr,trcic(1,3),ntrd)
         xnh41 = find(-zml(k),ztr,trcic(1,4),ntrd)
       endif

      bscale=f0star*v0*dhor/(ht*scsig*sigz(1))
      dbdz = (f0star*dhor/ht)**2./(scsig*sigz(1))

      call findkml(hdpth0,zml,lmax,kml)
      icurtrace = idin

      do 6 ip=1,msq

         delta = topden(ip) + artop
         bml = buoysur + bscale * delta - hdpth0 * dbdz
         xno3ml = nitr(stheta(bml))

         do 5 k=1,l+1

            if (irdtrc .eq. 1) then

               do 30 j=1,ntrace
                  icurtrace = j
                  trace(ip,k,mlcur,j) = trcicz(k,j) - trcqgofk(ip,k)
  30            continue

              else 

               if (k .le. kml) then
                  if (idinml .ne. 0) then
                     trace(ip,k,mlcur,idin) = xno3ml - trcqgofk(ip,k)
                    else
                     trace(ip,k,mlcur,idin) = - trcqgofk(ip,k)
                   endif
                 else
                  trace(ip,k,mlcur,idin) = 0.
                endif
               trace(ip,k,mlcur,ipon) = 0.
               trace(ip,k,mlcur,izoo) = 0.
               if (ibio .eq. ipnz1) then
                  trace(ip,k,mlcur,icq) = 0.
                 endif
               if (ibio .eq. ipnz3) then
                  trace(ip,k,mlcur,inh4) = 0.
                 endif

            endif

  5       continue
  6    continue

      else if (ifsbl .eq. 5) then

      call findkml(hdpth0,zml,lmax,kml)

      icurtrace = idin

      do 16 ip=1,msq

         do 15 k=1,l

            if (k .le. kml) then
               trace(ip,k,mlcur,idin) = xno3ml - trcqgofk(ip,k)
               trace(ip,k,mlcur,ipon) = float(ifsbl)
              else
               trace(ip,k,mlcur,idin) = 0.
               trace(ip,k,mlcur,ipon) = 0.
             endif

 15       continue

 16    continue

      endif

      return
      end
