c    @(#)ediag.f	1.1
c
c
c
      subroutine ediag
c
c --- this subroutine calculates the diagnostics for the code
c     and outputs the results. full field to idff, interior to idif.
c
c --- areaf is called for the full fields, areai is called
c --- for the interior fields.
c
      include 'qg_par.inc'
      common/scrat/r(mn0),r1(mn0),r2(mn0),r3(mn0)
      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/reffield/ref3(mn0),ref4(mn0,k0),ref5(mn0,k0)
      common/pref/irefs
      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/unit/iin,iout,idff,idif,iorr,iowr,iplt,ifndd
      common/cnn/n,nm1,nm2,np1
      common/extpar/lat0,lng0,t0,v0,dhor,ht,time0,
     +                r0gen,f0star
c
      dimension r0(mn0)
      dimension r4(mn0),r5(mn0),r6(mn0)
      dimension pers(mn0,k0),t1(mn0),t2(mn0),t3(mn0),t4(mn0)
      dimension t5(mn0),t6(mn0),t7(mn0),t8(mn0),t9(mn0)
      dimension dg(14,2)
c
      rndoff=.1e-7
c
c --- loop for all levels
c
        if(irefs.eq.1)call meshob(1)
c
      do 100 klevel=1,kz
c
c --- save the initial analysis field
c
        if(icalc.ne.0)goto 15
        do 10 ip=1,msq
          pers(ip,klevel)=stream(ip,klevel)
   10   continue
   15   continue
c
c --- get the 'real' data - r0 = psi, r1 = zeta
c
        if(irefs.eq.1)then
         do 18 ip=1,msq
          r0(ip)=ref4(ip,klevel)
          r1(ip)=ref5(ip,klevel)
   18   continue
        else
        call infld(r0,klevel,tim,3,0)
        call infld(r1,klevel,tim,4,0)
        endif
c
c --- get the statistics of the 'real' data
c
        call areaf(r0,1,1,1,s1af,s2af,s3af)
        call areai(r0,1,1,1,s1ai,s2ai,s3ai)
c
        call areaf(r1,0,1,0,dum,z2af,dum)
        call areai(r1,0,1,0,dum,z2ai,dum)
c
c --- get the statistics of the forecast data
c
        call areaf(stream(1,klevel),1,0,1,s1cf,dum,s3cf)
        call areai(stream(1,klevel),1,0,1,s1ci,dum,s3ci)
c
        call areaf(vort(1,klevel),1,1,0,z1cf,z2cf,dum)
        call areai(vort(1,klevel),1,1,0,z1ci,z2ci,dum)
c
c --- calculate the difference fields
c
        do 20 ip=1,msq
          r(ip)=stream(ip,klevel)-r0(ip)
          r1(ip)=vort(ip,klevel)-r1(ip)
          r2(ip)=stream(ip,klevel)-s1cf
          r3(ip)=r0(ip)-s1af
          r4(ip)=stream(ip,klevel)-pers(ip,klevel)
          r5(ip)=r0(ip)-pers(ip,klevel)
          r6(ip)=-r5(ip)
   20   continue
c
c --- get the statistics of the perturbation fields
c
        call areaf(r,1,1,0,s1df,s2df,dum)
        call areai(r,1,1,0,s1di,s2di,dum)
c
        call areaf(r1,0,1,0,dum,z2df,dum)
        call areai(r1,0,1,0,dum,z2di,dum)
c
        call areaf(r6,0,1,0,dum,s2ef,dum)
        call areai(r6,0,1,0,dum,s2ei,dum)
c
c --- calculate the standard deviation of the fcst error
c
        do 25 ip=1,msq
          t1(ip)=r(ip)-s1df
          t2(ip)=r(ip)-s1di
   25   continue
c
        call areaf(t1,0,1,0,dum,s2ff,dum)
        call areai(t2,0,1,0,dum,s2fi,dum)
c
c --- calculate absolute fcst correlation
c
        do 30 ip=1,msq
          t3(ip)=r2(ip)*r3(ip)
   30   continue
c
        call areaf(t3,1,0,0,s1gf,dum,dum)
        call areai(t3,1,0,0,s1gi,dum,dum)
c
        call areaf(r2,0,1,0,dum,s2hf,dum)
        call areai(r2,0,1,0,dum,s2hi,dum)
c
        call areaf(r3,0,1,0,dum,s2if,dum)
        call areai(r3,0,1,0,dum,s2ii,dum)
c 
c --- calculate fcst tendency correlation
c
        call areaf(r4,1,0,0,s1jf,dum,dum)
        call areai(r4,1,0,0,s1ji,dum,dum)
c
        call areaf(r5,1,0,0,s1kf,dum,dum)
        call areai(r5,1,0,0,s1ki,dum,dum)
c
        do 40 ip=1,msq
          t4(ip)=r4(ip)-s1jf
          t5(ip)=r4(ip)-s1ji
          t6(ip)=r5(ip)-s1kf
          t7(ip)=r5(ip)-s1ki
          t8(ip)=t4(ip)*t6(ip)
          t9(ip)=t5(ip)*t7(ip)
   40   continue
c
        call areaf(t4,0,1,0,dum,s2lf,dum)
        call areai(t5,0,1,0,dum,s2li,dum)
c
        call areaf(t6,0,1,0,dum,s2mf,dum)
        call areai(t7,0,1,0,dum,s2mi,dum)
c
        call areaf(t8,1,0,0,s1nf,dum,dum)
        call areai(t9,1,0,0,s1ni,dum,dum)
c
        if(z2af.eq.0.)z2af=rndoff
        if(z2ai.eq.0.)z2ai=rndoff
        if(s2af.eq.0.)s2af=rndoff
        if(s2ai.eq.0.)s2ai=rndoff
        if(s3af.eq.0.)s3af=rndoff
        if(s3ai.eq.0.)s3ai=rndoff
        if(s1ef.eq.0.)s1ef=rndoff
        if(s1ei.eq.0.)s1ei=rndoff
        if(s2ef.eq.0.)s2ef=rndoff
        if(s2ei.eq.0.)s2ei=rndoff
        if(s2hf.eq.0.)s2hf=rndoff
        if(s2hi.eq.0.)s2hi=rndoff
        if(s2if.eq.0.)s2if=rndoff
        if(s2ii.eq.0.)s2ii=rndoff
        if(s2lf.eq.0.)s2lf=rndoff
        if(s2li.eq.0.)s2li=rndoff
        if(s2mf.eq.0.)s2mf=rndoff
        if(s2mi.eq.0.)s2mi=rndoff
c
        e1=0.
        e2=0.
        e3=0.
        ei1=0.
        ei2=0.
        ei3=0.
        ml=mm1/4 + 1
        mu=m-ml + 1
	nl=nm1/4+1
	nu=n-nl+1
        do 60 j=2,nm1
          ip=(j-1)*m+1
          do 50 i=2,mm1
            ip=ip+1
            e4=(stream(ip+1,klevel)-stream(ip-1,klevel))**2+
     +          (stream(ip+m,klevel)-stream(ip-m,klevel))**2
            e1=amax1(e1,abs(r(ip)))
            e2=amax1(e2,abs(r1(ip)))
            e3=amax1(e3,e4)
            if(j.lt.nl .or. j.gt.nu .or. i.lt.ml .or. i.gt.mu)go to 50
            ei1=amax1(ei1,abs(r(ip)))
            ei2=amax1(ei2,abs(r1(ip)))
            ei3=amax1(ei3,e4)
   50     continue
   60   continue
c
c --- set the definitions of the diagnostics
c
        dg(1,1)=z1cf
        dg(1,2)=z1ci
        dg(2,1)=z2cf
        dg(2,2)=z2ci
        dg(3,1)=0.5*s3cf
        dg(3,2)=0.5*s3ci
        dg(4,1)=sqrt(e3)/(2.*hy)
        dg(4,2)=sqrt(ei3)/(2.*hy)
        dg(5,1)=s3cf/s3af-1.
        dg(5,2)=s3ci/s3ai-1.
        dg(6,1)=e1/sqrt(s2af)
        dg(6,2)=ei1/sqrt(s2ai)
        dg(7,1)=sqrt(s2df/s2af)
        dg(7,2)=sqrt(s2di/s2ai)
        dg(8,1)=sqrt(z2df/z2af)
        dg(8,2)=sqrt(z2di/z2ai)
        dg(9,1)=sqrt(s2ff/s2af)
        dg(9,2)=sqrt(s2fi/s2ai)
        dg(10,1)=dg(7,1)/sqrt(s2ef/s2af)
        dg(10,2)=dg(7,2)/sqrt(s2ei/s2ai)
        if(dg(10,1).gt.10.)dg(10,1)=10.
        if(dg(10,2).gt.10.)dg(10,2)=10.
        dg(11,1)=s1gf/sqrt(s2hf*s2if)
        dg(11,2)=s1gi/sqrt(s2hi*s2ii)
        dg(12,1)=s1nf/sqrt(s2lf*s2mf)
        dg(12,2)=s1ni/sqrt(s2li*s2mi)
c
c --- write the diagnostics to output files.
c
      if(ifndd.ne.0)go to 99
      ict=icalc
      if(icalc.gt.999)ict=icalc/10
      dtime=tim*t0/86400.+time0
      write(idff,801)klevel,ict,dtime,(dg(l,1),l=1,12)
      write(idif,801)klevel,ict,dtime,(dg(l,2),l=1,12)
      go to 100
c
   99 continue
      write(idff)klevel,icalc,tim,(dg(l,1),l=1,12)
      write(idif)klevel,icalc,tim,(dg(l,2),l=1,12)
  100 continue
      if(ifndd.ne.0)return
      write(idff,802)
      write(idif,802)
      return
c
c --- formats for numerical output.  headings are put on in the main.
c
  801 format(3x,i2,2x,i3,1x,0pf7.2,1p12e9.2)
  802 format(//)
      end
