c    @(#)mldiag.f	1.3
c
c
c
      subroutine mldiag
c
c
c --- routine to calculate the diagnostics for the upper mixed layer
c     model.  
c
c
c --- flexibility
c
c     flexibility is obtained by defining the variable ndiag
c     which is is equal to the number of diagnostic quantities
c     and by writing the titles of the diagnostics.  this will
c     permit the addition, revision, and removal of diagnostics
c     without changing associated plotting software.
c
c
c
c --- variable definitions
c
c
c
c
c internal variables and constants:
c
c     ndiag     the number of diagnostic quantities which are calcu-
c               lated.  this should be the same as the number of
c               titles......
c
c     nchar     the number of characters in each title including
c               blanks.  it is generally not a good idea to change this
c               abbreviate the longes titles instead.
c
c     dtitl     an array of characters which must be at least
c               character*nchar dtitl(ndiag) 
c               used to store the titles of the diagnostics.
c
c
c     vdiag     an array at least (ndiag,max # of levels) which is used
c               to hold the calculated quantities.
c
c diagnostic quantities
c
c     bqsw      short wave radiation at the surface. (w/cm/cm)
c
c     bqsurf    long wave + sensible heat transfer at the
c               surface.                             (w/cm/cm)
c
c     bwflux    water flux at the surface.           (m/s)
c
c     bustar     surface friction velocity.          (m/s)
c
c     bhdpth     mean mixed layer depth.             (m)
c
c     bmdpth     mean mixing layer depth.            (m)
c
c     btbar      mean mixed layer temperature.       (c)
c
c     bbbar      mean mixing layer buoyancy.         (m/s/s)
c
c     ucfl       cfl for x direction.
c
c     vcfl       cfl for y direction.
c
c     wcfl       cfl for z direction. 
c
c     hrmsd      nrms difference in mixed layer depth  
c
c     trmsd      nrms difference in mixed layer temperature
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
c include
c
      include'mladv.inc'
c
c
c common
c
      character*70 idstring
      character*80 ident
      character*4 restid,runid
      common /inform/ifdiff,ifpert,titl(20,2),runid(6),restid(9),
     +               idstring,ident
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

      common /cnn/n,nm1,nm2,np1
c
      common/extpar/rlat0,rlng0,t0,v0,dhor,ht,time0,
     +                r0,f0star
c
      common /scrat/r(mn0),r1(mn0),r2(mn0),r3(mn0)
c
c
c dimension
c
      character*32 dtitl(17)
c
      dimension plim(2,17)
c
c_______________________________________________________________________________
c
c
c_______________________________________________________________________________
c
c
c data
c
      data jdiag/12/, nword/8/
c
      data dtitl(1)/'short wave radiation (langley/d)'/
      data dtitl(2)/'  surface heat flux  (langley/d)'/
      data dtitl(3)/'  surface water flux  (meters/d)'/
      data dtitl(4)/'  friction velocity   (meters/s)'/
      data dtitl(5)/'   mean mixed layer depth (m)   '/
      data dtitl(6)/'  mean mixing layer depth (m)   '/
      data dtitl(7)/'  mean mixing layer temp (c)    '/
      data dtitl(8)/'  mean mixing layer buoy (m/s/s)'/
      data dtitl(9)/'           u - cfl              '/
      data dtitl(10)/'                                '/
      data dtitl(11)/'                                '/
      data dtitl(12)/'                                '/
      data dtitl(13)/'                                '/
      data dtitl(14)/'                                '/
      data dtitl(15)/'                                '/
      data dtitl(16)/'                                '/
      data dtitl(17)/'                                '/
c    +            '  me','an m','ixin','g la','yer ',
c    +                   'temp',' (c)','    ',
c    +            'mean',' mix','ing ','laye','r bu',
c    +                   'oy (','m/s/','s)  ',
c    +            '    ','    ','   u',' - c','fl  ',
c    +                   '    ','    ','    ',
c    +            '    ','    ','   v',' - c','fl  ',
c    +                   '    ','    ','    ',
c    +            '    ','    ','   w',' - c','fl  ',
c    +                   '    ','    ','    ',
c    +            'miss','ed i','nter','ior ','entr',
c    +                   'ainm','ent ','(%) ',
c    +            'miss','ed b','ound','ary ','entr',
c    +                   'ainm','ent ','(%) ',
c    +            ' bou','yanc','y en','trai','nmen',
c    +                   't (m','/s/s','/s) ',
c    +            ' max','imum',' mix','ed l','ayer',
c    +                   ' dep','th (','m)  ',
c    +            '  nr','msd(','mixe','d la','yer ',
c    +                   'dept','h) (','%)  ',
c    +            '  nr','msd(','mixe','d la','yer ',
c    +                   'temp',')  (','%)  '/
c
      data plim/ 999., 999.,
     +           999., 999.,
     +           999., 999.,
     +           999., 999.,
     +            0.0, 999.,
     +            0.0, 999.,
     +           999., 999.,
     +           999., 999.,
     +            0.0,  2.0,
     +            0.0,  2.0,
     +            0.0,  2.0,
     +            0.0, 999.,
     +            0.0, 999.,
     +            0.0,  0.0,
     +            0.0, 999.,
     +            0.0, 999.,
     +            0.0, 999./
      data icall/0/
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
c if this is not the first call then skip the output.
c
      if(icall.ne.0)go to 99
      icall=1
c
c output titles
c
      nchar=32
      ndiag=17
      write(jdiag)runid,restid,titl
      write(jdiag)ndiag,nchar,nword
      write(jdiag)dtitl
c
c
c output combination instructions.
c
      write(jdiag)2,1,2,1,3,1,4,2,5,6,1,7,1,8,3,9,10,11,2,12,
     .            13,2,16,17,1,15
c
c
c output limits on the range of the plot
c
      write(jdiag)plim
c
c
c one time calculations.
c
      amp=0.0
      do 10 k=1,l
        tmp=-zml(l)/dzml(l)
        if(tmp.gt.amp)amp=tmp
   10 continue
c
c
c
   99 continue
c
c
c zero the diagnostic quantities.
c
      bqsw=0.0
      bqsurf=0.0
      bwflux=0.0
      bustar=0.0
      bhdpth=0.0
      bmdpth=0.0
      btbar=0.0
      bbbar=0.0
      ucfl=0.0
      vcfl=0.0
      wcfl=0.0
      bbe=0.0
      bhmax=0.0
c
c now loop through the field.
c
      do 100 ip=1,msq
        bqsw=bqsw+qsw(ip)
        bqsurf=bqsurf+qsurf(ip)
        bwflux=bwflux+wflux(ip)
        bustar=bustar+sqrt(sqrt(taux(ip)**2+tauy(ip)**2))
        bhdpth=bhdpth+out(ip)
        bmdpth=bmdpth+hdpth(ip,mlcur)
        btbar=btbar+tbar(ip)
        bbbar=bbbar+bbar(ip)
        bbe=bbe+be(ip)
        bhmax=bhmax+hmax(ip)
c
        tmp=abs(wzvel(ip)+wsurf(ip))
        if(tmp.gt.wcfl)wcfl=tmp
c
  100 continue
c
c
c
        tmptim=(mltime-time0)*86400./t0
c        call ref(23,0.0,tmptim,r)
c        hrmsd=rmsd(msq,hdpth(1,mlcur),r)
c
c
c        call ref(24,q,tmptim,r)
c        trmsd=rmsd(msq,tbar,r)
c
c
c
c
c  now normalize and change units.
c
      tmp=1./msq
      bqsw=bqsw*tmp*86400.
      bqsurf=bqsurf*tmp*86400.
      bwflux=bwflux*tmp*86400.
      bustar=bustar*tmp
      bhdpth=bhdpth*tmp
      bmdpth=bmdpth*tmp
      btbar=btbar*tmp
      bbbar=bbbar*tmp
      ucfl=ucfl*mlstep/dxyml
      vcfl=vcfl*mlstep/dxyml
      wcfl=wcfl*amp*mlstep
      bbe=bbe*tmp
      bhmax=tmp*bhmax
      hrmsd=hrmsd*100.
      trmsd=trmsd*100.
c
c now output
c
      write(jdiag)mltime,bqsw,bqsurf,bwflux,bustar,bhdpth,
     .            bmdpth,btbar,bbbar,ucfl,vcfl,wcfl,100*bimiss,
     .            100*bbmiss,bbe,bhmax,hrmsd,trmsd
c
c
c
      return
      end
