c
c  @(#)sbldiag.f	1.1
c
      subroutine sbldiag(iopt)

      include 'mladv.inc'

      dimension prevtot(mn0,2),delqg(mn0,2),prevqg(mn0,2),
     .          deltot(mn0,2)

      dimension delsbl(mn0),deladv(mn0),delbal(mn0),delmix2(mn0),
     .          dtphys2(mn0),prevsbl(mn0),delphys(mn0)

      dimension delmix(mn0,2,2),dtphys(mn0,2,2)

      common/inform/ifdiff,ifpert,titl(20,2),runid(6),restid(9),
     +              idstring,ident

      character*80 ident,fname
      character*70 idstring
      character*4 restid,runid

      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 /extpar/ rlat0,rlng0,t0,v0,dhor,ht,time0,r0,f0star

      data icall/0/,mixcur/1/,mixpast/2/

      if (icall .eq. 0) then
         icall = 1
         iflen=lnblk(ident,80)
         fname=ident(:iflen)//'.sbl'
         open(unit=isbld,file=fname,form='formatted',status='new')
         call ndinfo(isbld)
        endif

      tmptim=(mltime-time0)*86400./t0
      zdl = zml(idl)
      icurtrace = 1


      if (iopt .eq. 1) then

         do 206 ip=1,msq
            deltot (ip,mlcur) = trace(ip,idl,mlcur,1) 
     .                        + trcqgofk(ip,idl) - prevtot (ip,mlcur)
            prevtot(ip,mlcur) = trace(ip,idl,mlcur ,1) 
     .                        + trcqgofk(ip,idl)
            delqg  (ip,mlcur) = trcqgofk(ip,idl) - prevqg(ip,mlcur)
            prevqg (ip,mlcur) = trcqgofk(ip,idl)
            delmix2(ip) = delmix(ip,mlcur ,mixcur ) 
     .                  + delmix(ip,mlpast,mixpast)
            dtphys2(ip) = dtphys(ip,mlcur ,mixcur ) 
     .                  + dtphys(ip,mlpast,mixpast)
            delsbl (ip) = deladv(ip) + delmix2(ip)
            delphys(ip) = delqg(ip,mlcur) + delsbl(ip)
            delbal (ip) = delqg(ip,mlcur) + delsbl(ip)
     .                  + dtphys2(ip) - deltot(ip,mlcur)
            prevsbl(ip) = trace(ip,idl,mlpast,1)
  206     continue

         itmp    = mixcur
         mixcur  = mixpast
         mixpast = itmp

         if (mladd .ne. mlpast) then
            write(*,*)'SBL diags performed on reverse time step'
            write(*,*)'Time :',tmptim
           endif

         if (ifsbld .ne. 0) then
            call xwrite(prevtot(1,mlcur),
     .                            msq,tmptim,20,idl,zdl,zdl,isbld)
            call xwrite(delqg  (1,mlcur),
     .                            msq,tmptim,21,idl,zdl,zdl,isbld)
            call xwrite(deltot  (1,mlcur),
     .                            msq,tmptim,25,idl,zdl,zdl,isbld)
            call xwrite(delsbl ,msq,tmptim,22,idl,zdl,zdl,isbld)
            call xwrite(delphys,msq,tmptim,23,idl,zdl,zdl,isbld)
            call xwrite(delbal ,msq,tmptim,26,idl,zdl,zdl,isbld)
            call xwrite(deladv ,msq,tmptim,27,idl,zdl,zdl,isbld)
            call xwrite(delmix2,msq,tmptim,31,idl,zdl,zdl,isbld)
            call xwrite(dtphys2,msq,tmptim,24,idl,zdl,zdl,isbld)
           endif

        else if (iopt .eq. 2) then
           do 205 ip=1,msq
              deladv  (ip) = trace(ip,idl,mlpast,1) - prevsbl(ip)
              prevsbl (ip) = trace(ip,idl,mlcur,1)
  205       continue

        else if (iopt .eq. 3) then
           do 208 ip=1,msq
              delmix(ip,mlcur,mixcur) = trace(ip,idl,mlcur,1) 
     .                                - prevsbl(ip)
              prevsbl(ip) = trace(ip,idl,mlcur,1)
  208       continue

        else if (iopt .eq. 4) then
           do 207 ip=1,msq
              dtphys(ip,mlcur,mixcur) = trace(ip,idl,mlcur,1) 
     .                                - prevsbl(ip)
  207       continue
 
        else if (iopt .eq. 5) then
           do 209 ip=1,msq
              prevsbl(ip) = trace(ip,idl,mlcur,1)
  209       continue

        else if (iopt .eq. 6) then
           do 202 ip=1,msq
              delmix(ip,mlcur,mixcur) = trace(ip,idl,mlcur,1) 
     .                                - prevsbl(ip)
              prevsbl(ip) = trace(ip,idl,mlcur,1)
  202       continue

        else if (iopt .eq. 7) then
           do 201 ip=1,msq
              dtphys(ip,mlcur,mixcur) = trace(ip,idl,mlcur,1) 
     .                                - prevsbl(ip)
  201       continue
 
       endif

       return
       end
