c
c    @(#)mlfout.f	1.6
c
      subroutine mlfout
c
c --- routine to output the fields from the upper mixed layer model.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c variables for control of the output
c
c     nout      number of output quantities
c
c     dout      array which defines the output type, depth, and rate.
c               for some quantities the depth is not necessary.
c               if dout(1,1),(2,1),(3,1)=20,0.,1. then the mixed layer
c               depth will be plotted daily.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c include
c
      include'mladv.inc'
      include'trc.inc'

      dimension trcml(maxtrace)

c
c common
c
      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
c
      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)
c
      common/user/pram(16)
c
      common /extpar/ rlat0,rlng0,t0,v0,dhor,ht,time0,r0,f0star
c
      common /cnn/n,nm1,nm2,np1
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c common
c
      common /scrat/ r(mn0), r1(mn0), r2(mn0), r3(mn0)
c
      common /unit/iin,iout,idff,idif,iorr,iowr,iplt,ifndd
c
      common/inform/ifdiff,ifpert,titl(20,2),runid(6),restid(9),
     +              idstring,ident
c
      character*70 idstring
      character*80 ident,fname
      character*4 restid,runid
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      data icall/0/
      data icall2/0/
      data icall3/0/
      data icall4/0/
      data icall5/0/
      data icall6/0/
      data icall9/0/
      data icall10/0/
      data icalt1/0/
      data icalt2/0/
      data icalt3/0/
      data icalt4/0/

      if (nout.eq.0) return
c
      tmptim=(mltime-time0)*86400./t0
      do 1000 ipnt=1,nout
        it  =dout(1,ipnt)
        q   =dout(2,ipnt)
        rate=dout(3,ipnt)
        tstd = tstart * t0 / 86400.
        if (amod(mltime-tstd+qgts*.5/86400.,rate)
     .                           .ge.qgts/86400.) goto 1000
        if (it.eq.20) goto 200
        if (it.eq.21) goto 210
        if (it.eq.22) goto 220
        if (it.eq.23) goto 230
        if (it.eq.24) goto 240
        if (it.eq.25) goto 250

        if (it.eq.29) goto 290
        if (it.eq.30) goto 300

        if (it.eq.32) goto 320
        if (it.eq.34) goto 340
        if (it.eq.35) goto 350
c --- SB88 fields

        if (it.eq.37) goto 370
        if (it.eq.38) goto 380
c --- Tracer fields

        if (it.eq.51) goto 510
        if (it.eq.52) goto 520
        if (it.eq.53) goto 530
        if (it.eq.54) goto 540

        if (it.eq.72) goto 720
        if (it.eq.75) goto 750
        if (it.eq.79) goto 790
 
c --- Profile options

        if(it.eq.80) goto 800
        if(it.eq.81) goto 810

c --- Temperature and sigma t

        if(it.eq.82) goto 820
        if(it.eq.83) goto 830

c --- Tracers

        if(it.eq.84) goto 840
        if(it.eq.85) goto 850
        if(it.eq.91) goto 910
        if(it.eq.92) goto 920
        if(it.eq.93) goto 930
        if(it.eq.94) goto 940
        if(it.eq.98) goto 980

c
c --- mld time series
c
        if (it .eq. 99) goto 990
c
c
c      
        go to 1000
c
c
  200   continue
        call xwrite(hdpth(1,mlcur),msq,tmptim,20,0,0.,0.,iplt)
        go to 1000
c
c
  210   continue
        if(q.lt.0.0)go to 1000
        do 211 k=l,1,-1
          if(q.le.-zml(k))kml=k
  211   continue
        q=-zml(kml)
        do 215 ip=1,msq
          r(ip)=bqgofk(ip,kml)+buoy(ip,kml,mlcur)
  215   continue
        call xwrite(r,msq,tmptim,21,kml,q,q,iplt)
        go to 1000
c
c
  220   continue
        call xwrite(tbar,msq,tmptim,22,0,0.,0.,iplt)
        go to 1000

  230   continue
        kml = 1
        q = 0.0
        do 235 ip=1,msq
          r(ip)=tqgofk(ip,kml)+temp(ip,kml,mlcur)
  235   continue
        call xwrite(r,msq,tmptim,23,0,q,q,iplt)
        go to 1000
c
  240   continue
        if(q.lt.0.0)go to  1000
        kml=kmlofq(q,zml,l)
        q=-zml(kml)
        do 245 ip=1,msq
          r(ip)=bqgofk(ip,kml)+buoy(ip,kml,mlcur)
  245   continue
        call xwrite(r,msq,tmptim,24,kml,q,q,iplt)
        go to 1000
c
c
  250   continue
        call xwrite(out,msq,tmptim,25,0,0.,0.,iplt)
        go to 1000
c
c
c
  290   continue
        if(q.gt.-zml(l))go to 1000

        do 295 ip=1,msq
          r(ip)= tauy(ip)/hdpth(ip,mlcur)/f0star
  295   continue
        call xwrite(r,msq,tmptim,29,0,q,q,iplt)
        go to 1000

  300   continue
        if(q.gt.-zml(l))go to 1000

        do 305 ip=1,msq
          r(ip)= -taux(ip)/hdpth(ip,mlcur)/f0star
  305   continue
        call xwrite(r,msq,tmptim,30,0,q,q,iplt)
        go to 1000
c
c

  320   continue
        if(q.gt.-zml(l))go to 1000

        do 325 ip=1,msq
          r(ip) = -wzvel(ip)*q + wsurf(ip)
  325   continue
        call xwrite(r,msq,tmptim,32,0,q,q,iplt)
        go to 1000
c
c
  340   continue
        q = 0.0
        fact = f0star * 1025.0
        do 345 ip=1,msq
          r(ip) = taux(ip)*fact
  345   continue
        call xwrite(r,msq,tmptim,34,0,q,q,iplt)
        go to 1000
c
c
  350   continue
        q = 0.0
        fact = f0star * 1025.0
        do 355 ip=1,msq
          r(ip) = tauy(ip)*fact
  355   continue
        call xwrite(r,msq,tmptim,35,0,q,q,iplt)
        go to 1000

c
c
  370   continue
        fact = 1.0
        do 375 ip=1,msq
          r(ip)=wflux(ip)*fact
  375   continue
        call xwrite(r,msq,tmptim,37,0,0.,0.,iplt)
        go to 1000
c
c
  380   continue
        do 385 ip=1,msq
          call mlbflux(hdpth(ip,mlcur),ip,r(ip))
  385   continue
        call xwrite(r,msq,tmptim,38,0,0.,0.,iplt)
        go to 1000
c
c
c
  510   continue
        kml=kmlofq(q,zml,l+1)
        q=-zml(kml)
        icurtrace = 1
        do 511 ip=1,msq
           if (kml .le. l) then
              r(ip) = trace(ip,kml,mlcur,1) + trcqgofk(ip,kml)
             else
              r(ip) = trace(ip,kml,mlcur,1)
            endif
  511   continue
        call xwrite(r,msq,tmptim,it,kml,q,q,iplt)
        go to 1000

  520   continue
        kml=kmlofq(q,zml,l)
        q=-zml(kml)
        icurtrace = 2 
        do 521 ip=1,msq
           r(ip) = trace(ip,kml,mlcur,2) + trcqgofk(ip,kml) 
  521   continue
        call xwrite(r,msq,tmptim,it,kml,q,q,iplt)
        go to 1000

  530   continue
        kml=kmlofq(q,zml,l)
        q=-zml(kml)
        icurtrace = 3
        do 531 ip=1,msq 
           r(ip) = trace(ip,kml,mlcur,3) + trcqgofk(ip,kml)  
  531   continue 
        call xwrite(r,msq,tmptim,it,kml,q,q,iplt)
        go to 1000

  540   continue
        kml=kmlofq(q,zml,l)
        q=-zml(kml)
        icurtrace = 4
        do 541 ip=1,msq 
           r(ip) = trace(ip,kml,mlcur,4) + trcqgofk(ip,kml)  
  541   continue 
        call xwrite(r,msq,tmptim,it,kml,q,q,iplt)
        go to 1000

  720   continue
        do 725 ip=1,msq
          r(ip)= -qsurf(ip)
  725   continue
        call xwrite(r,msq,tmptim,72,0,0.,0.,iplt)
        go to 1000
  750   continue
        do 755 ip=1,msq
          r(ip)= qsw(ip)
  755   continue
        call xwrite(r,msq,tmptim,75,0,0.,0.,iplt)
        go to 1000
  790   continue
        do 795 ip=1,msq
          r(ip)=qsw(ip) -qsurf(ip)
  795   continue
        call xwrite(r,msq,tmptim,79,0,0.,0.,iplt)
        go to 1000

  800   continue
        if (icall.eq.0) then
           icall = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.prf'
           open(unit=iuprf,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        do 805 k=1,l
           write(iuprf,806)zml(k),buoy(ip,k,mlcur)+bqgofk(ip,k)
  805   continue
  806   format(99e12.4)
        goto 1000

  810   continue
        if (icall2.eq.0) then
           icall2 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.prf2'
           open(unit=iuprf2,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        do 815 k=1,l
           write(iuprf2,806)zml(k),trace(ip,k,mlcur,2)
  815   continue
        goto 1000

  820   continue
        if (icall3.eq.0) then
           icall3 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.tem'
           open(unit=iutem,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        do 825 k=1,l
           write(iutem,806)zml(k),temp(ip,k,mlcur)+tqgofk(ip,k)
  825   continue
        goto 1000

  830   continue
        if (icall4.eq.0) then
           icall4 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.sgt'
           open(unit=iusgt,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        do 835 k=1,l
           write(iusgt,806)zml(k),stheta(bqgofk(ip,k)+buoy(ip,k,mlcur))
  835   continue
        goto 1000

  840   continue
        if (icall5.eq.0) then
           icall5 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.prd'
           open(unit=iuprd,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        call light(ip)
        yktot = 0.0
        do 845 k=1,l
 
           icurtrace = idin
           s = trace(ip,k,mlcur,idin) + trcqgofk(ip,k)
 
           icurtrace = ipon
           phyto = trace(ip,k,mlcur,ipon) + trcqgofk(ip,k)
 
           icurtrace = izoo
           zoo = trace(ip,k,mlcur,izoo) + trcqgofk(ip,k)
 
           icurtrace = inh4
           s2 = trace(ip,k,mlcur,inh4) + trcqgofk(ip,k)
 
           q1 = (s * exp(-kn * s2)) / (ks + s)
           q2 = s2 / (ks2 + s2)
           q  = q1 + q2

           pofi = pmax * (1.0 - exp(-p1*epar(k)/pmax) ) *
     .                            exp(-p2*epar(k)/pmax) * phyto
 

           xktot = kw + (kc * phyto)
           yktot = yktot + xktot

           write(iuprd,806)mltime,zml(k),epar(k),xktot,yktot/float(k),
     .                  q,q1,q2,pofi,phyto

  845   continue
        goto 1000

  850   continue
        if (icall6.eq.0) then
           icall6 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.trf'
           open(unit=iutrf,file=fname,form='formatted',status='new')
         endif
        ip = int(q)

        do 855 k=1,l
           write(iutrf,806)mltime,zml(k),(tfluxd(ip,k,i),i=1,ntflx)
  855   continue
        goto 1000

  910   continue
        if (icalt1.eq.0) then
           icalt1 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.tr1'
           open(unit=iutr1,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        icurtrace = 1
        do 915 k=1,l
           write(iutr1,806)zml(k),trace(ip,k,mlcur,1) + trcqgofk(ip,k)
  915   continue
        goto 1000

  920   continue
        if (icalt2.eq.0) then
           icalt2 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.tr2'
           open(unit=iutr2,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        icurtrace = 2
        do 925 k=1,l
           write(iutr2,806)zml(k),trace(ip,k,mlcur,2) + trcqgofk(ip,k)
  925   continue
        goto 1000

  930   continue
        if (icalt3.eq.0) then
           icalt3 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.tr3'
           open(unit=iutr3,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        icurtrace = 3
        do 935 k=1,l
           write(iutr3,806)zml(k),trace(ip,k,mlcur,3) + trcqgofk(ip,k)
  935   continue
        goto 1000

  940   continue
        if (icalt4.eq.0) then
           icalt4 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.tr4'
           open(unit=iutr4,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        icurtrace = 4
        do 945 k=1,l
           write(iutr4,806)zml(k),trace(ip,k,mlcur,4) + trcqgofk(ip,k)
  945   continue
        goto 1000

  980   continue
        if (icall10.eq.0) then
           icall10 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.mlt'
           open(unit=iumlt,file=fname,form='formatted',status='new')
         endif
        ip = int(q)
        do 981 icurtrace=1,ntrace
           trcml(icurtrace) = trace(ip,1,mlcur,icurtrace) + 
     .                                         trcqgofk(ip,1)
  981    continue
        write(iumlt,806)mltime,(trcml(i),i=1,ntrace)
        goto 1000

  990   continue
        if (icall9.eq.0) then
           icall9 = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.mld'
           open(unit=iumld,file=fname,form='formatted',status='new')
         endif
        call diag(hdpth(1,mlcur),msq,ave,sdv,rms,xmin,xmax)
        write(iumld,806)mltime,ave,xmin,xmax,sdv,rms
        goto 1000

c
c
c
 1000 continue
c
      return
      end
