program plot11

      parameter(nx=64,ny=64,nz=64)
      parameter(nxyz=nx*ny*nz)
c super droplets (SDs): only in cloudy volumes...
c nppg - the number of SDs per gridbox
      parameter(nppg=20)
      parameter(npp=(nx-1)*(ny-1)*(nz-1)*nppg)

      parameter(nt0=100)
      dimension time(nt0),rad(nt0),sup(nt0),std(nt0)

      real*8 u(nx,ny,nz),v(nx,ny,nz),w(nx,ny,nz)
      real*8 th(nx,ny,nz),qv(nx,ny,nz),ss(nx,ny,nz)
      real*8 xp(npp),yp(npp),zp(npp)
      real*8 qcp(npp),plic(npp),radp(npp)
      real*8 sum,sum1,sum2
      dimension x(2),y(2)

      call opngks
      call gsclip(0)
      call dfclrs

      n=nx
      m=ny
      l=nz

         do isss=1,1
         do icas=1,3

c        r_act=5.39 ! real
        r_act=1. ! fake

         if(isss.eq.1) then
       if(icas.eq.1) then
       open(10,file='../fort.9.lg.20.1mps',form='unformatted',
     1        status='old')
       endif
       if(icas.eq.2) then
       open(10,file='../fort.9.lg.20.3mps',form='unformatted',
     1        status='old')
       endif
       if(icas.eq.3) then
       open(10,file='../fort.9.lg.20.0.333mps',form='unformatted',
     1        status='old')
       endif
       endif

      ee0=1227.2
      t00=283.16
      hlat=2.47e6
      rv=461.
      rho0=1.
      rg=287.
      gg=9.81

        if(icas.eq.1) nt=100
        if(icas.eq.2) nt=33
        if(icas.eq.3) nt=100

        updr=1.  ! other w are included by time not really being height...

        do it=1,nt

          if(icas.eq.1) time(it)=it*2.
          if(icas.eq.2) time(it)=it*6.
          if(icas.eq.3) time(it)=it*2.

          read(10)    ! (((u(i,j,k),i=1,n),j=1,m),k=1,l)
          read(10)    ! (((v(i,j,k),i=1,n),j=1,m),k=1,l)
          read(10)    ! (((w(i,j,k),i=1,n),j=1,m),k=1,l)
          read(10)    (((th(i,j,k),i=1,n),j=1,m),k=1,l)
          read(10)    (((qv(i,j,k),i=1,n),j=1,m),k=1,l)
          read(10)   xp,yp,zp,qcp,plic,radp 

       print*,'--- read for time: ',time(it)
c         print*,(radp(ii),ii=1,5)
         sum1=0.
         sum2=0.
      do i=1,n-1
      do j=1,m-1
      do k=1,l-1
      thi=1./th(i,j,k)
      ess=ee0*exp(hlat/rv * (1./t00 - thi))
      pre=0.9e5 - gg*rho0*updr*time(it)
      qvs=rg/rv*ess/(pre-ess)
      ss(i,j,k)=qv(i,j,k)/qvs - 1.
      sum1=sum1+1
      sum2=sum2+ss(i,j,k)*100.
      enddo
      enddo
      enddo
        sup(it)=sum2/sum1

         sum1=0.
         sum2=0.
      do i=1,n-1
      do j=1,m-1
      do k=1,l-1
      sum1=sum1+1
      sum2=sum2+(ss(i,j,k)*100.-sup(it))**2
      enddo
      enddo
      enddo
        std(it)=sqrt(sum2/sum1)

         sum1=0.
         sum2=0.
         sum3=0.
         sum4=0.
          do i=1,npp
          sum3=sum3+1.
          sum4=sum4+radp(i)*1.e6
          if(radp(i).ge.r_act*1.e-6) then
          sum1=sum1+1.
          sum2=sum2+radp(i)*1.e6
          endif
          enddo
          sum4=sum4/sum3
          sum22=0.
          if(sum1.gt.0.) sum22=sum2/sum1
          rad(it)=sum4
          if(sum22.gt.r_act) rad(it)=sum22
c           print*,it,pre,rad(it),sup(it)
           print*,it,rad(it),sup(it),std(it)

        enddo

      call setusv('LW',2000)
      call gstxci(15)
      call gsplci(15) ! black
      call set(.2,.8,.1,.3,0.,200.,-3.,3.,1)
      call labmod('(e5.1)','(e5.1)',5,5,2,2,28,20,0)
      call perim (2,10,2,3)
      call setusv('LW',1000)
        x(1)=0
        x(2)=200
        y(1)=0.15
        y(2)=0.15
        call dashdc('$''''$$''''$$''''$$''''$$''''$$''''$', 6,12)
        call curved(x,y,2)
        call dashdc('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$', 3,12)
        y(1)=0.
        y(2)=0.
        call curved(x,y,2)
      call setusv('LW',2000)

       ino=0.
      if(icas.eq.1) then
         call gsplci(12) ! red line
         call gspmci(12) ! red symbol
      endif
      if(icas.eq.2) then
         call gsplci(14) ! blue line
         call gspmci(14) ! blue symbol
      endif
      if(icas.eq.3) then
         call gsplci(15) ! black line
         call gspmci(15) ! black symbol
      endif

         do it=1,nt
        call points(time(it),sup(it),1,-4,0)
         enddo


      call setusv('LW',2000)
      call gstxci(15)
      call gsplci(15) ! black
      call set(.2,.8,.3,.5,0.,200.,0.,1.,1)
      call labmod('(e5.1)','(e5.1)',5,5,2,2,28,20,0)
      call perim (2,10,1,5)

      if(icas.eq.1) then
         call gsplci(12) ! red line
         call gspmci(12) ! red symbol
      endif
      if(icas.eq.2) then
         call gsplci(14) ! blue line
         call gspmci(14) ! blue symbol
      endif
      if(icas.eq.3) then
         call gsplci(15) ! black line
         call gspmci(15) ! black symbol
      endif

         do it=1,nt
        call points(time(it),std(it),1,-2,0)

        if(rad(it).ge.r_act .and.ino.eq.0)then
        ino=1
      call setusv('LW',1000)
        x(1)=time(it)
        x(2)=time(it)
        y(1)=-0.5
        y(2)=2.
c        call curved(x,y,2)
        endif
        enddo

      call setusv('LW',2000)
      call gstxci(15)
      call gsplci(15) ! black
      call set(.2,.8,.52,.92,0.,200.,1.e-2,1.e1,2)
      call labmod('(e5.1)','(e5.1)',5,5,2,2,28,20,0)
      call perim (2,10,1,9)
      if(icas.eq.1) then
         call gsplci(12) ! red line
         call gspmci(12) ! red symbol
      endif
      if(icas.eq.2) then
         call gsplci(14) ! blue line
         call gspmci(14) ! blue symbol
      endif
      if(icas.eq.3) then
         call gsplci(15) ! black line
         call gspmci(15) ! black symbol
      endif
         do it=1,nt
         call points(time(it),rad(it),1,-2,0)

        if(rad(it).ge.r_act .and.ino.eq.1)then
        ino=2
      call setusv('LW',1000)
        x(1)=time(it)
        x(2)=time(it)
        y(1)=0.1
        y(2)=7.
c        call curved(x,y,2)
        endif

      enddo

      close(10)

      enddo
      enddo

      call setusv('LW',2000)
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
      call gstxci(15)
      call gsplci(15) ! black
c
      call gstxci(15)
      call gsplci(15) ! black
      call gspmci(15) ! black
      CALL plchhq(.20,0.07, '0', 0.016,0.,0)
      CALL plchhq(.50,0.07, '100', 0.016,0.,0)
      CALL plchhq(.80,0.07, '200', 0.016,0.,0)
      CALL plchhq(.50,0.03,'height (m)', 0.016,0.,0)

      CALL plchhq(.16,0.1, '-3', 0.016,0.,0)
      CALL plchhq(.16,0.2, '0', 0.016,0.,0)
      CALL plchhq(.16,0.287, '3', 0.016,0.,0)
      CALL plchhq(.65,0.14,'mean', 0.014,0.,0)
      CALL plchhq(.16,0.313, '0', 0.016,0.,0)
      CALL plchhq(.16,0.49, '1', 0.016,0.,0)
      CALL plchhq(.65,0.46,'standard deviation', 0.014,0.,0)

      CALL plchhq(.11,0.3,'supersaturation (%)', 0.016,90.,0)

      CALL plchhq(.16,0.53, '10:S:-2', 0.016,0.,0)
      CALL plchhq(.16,0.65, '10:S:-1', 0.016,0.,0)
      CALL plchhq(.16,0.79, '10:S:0', 0.016,0.,0)
      CALL plchhq(.16,0.92, '10:S:1', 0.016,0.,0)
      CALL plchhq(.11,0.72,'radius (micron)', 0.016,90.,0)

      call frame
      call clsgks
      stop
      end

      subroutine dfclrs
cc as on page 4-31 in NCAR Graphics Guide to New Utilities
      dimension rgbv(3,15)
      data rgbv /
     +           1.00, 1.00, 1.00,
c dark blue into white
     +            .00,  .00,  1.0,
     +            1. ,   .5,   0.,
     +            1.0,  .50, 1.0,
     +            .92,  .92, .92,
     +            .88,  .88, .88,
     +            .84,  .84, .84,
     +            .80,  .80, .80,
     +            .76,  .76, .76,
     +            .72,  .72, .72,
     +            .68,  .68, .68,
cc
     +           1.00,  .00,  .00 ,
     +           0.00,  1.0,   .0 ,
     +           0.00,  0.0, 1.00 ,
     +            .00,  0.0, 0.00
     +                                /
c
cc background color
cc this is black      call gscr(1,0,0.,0.,0.)
cc this is Dons background, next is my modification
cc      call gscr(1,0,0.31,0.31,0.39)
c     call gscr(1,0,0.51,0.51,0.69)
      call gscr(1,0,1.00,1.00,1.00)
c     call gscr(1,0,0.00,0.00,0.25)
c
      do 101 i=1,15
       call gscr(1,i,rgbv(1,i),rgbv(2,i),rgbv(3,i))
 101  continue
      return
      end