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)
c      parameter(nppg=4)
      parameter(npp=(nx-1)*(ny-1)*(nz-1)*nppg)

      parameter(nt0=100)
      dimension time(nt0),rad(nt0),std(nt0),afra(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
        r_act=1.

          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

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

        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)

         sum1=0.
         sum2=0.
          do i=1,npp
          if(radp(i).gt.r_act*1.e-6) then
          sum1=sum1+1.
          sum2=sum2+radp(i)
          endif
          enddo
          if(sum1.gt.0.) then
          rad(it)=sum2/sum1
          else
          rad(it)=0.
          endif


         sum1=0.
         sum2=0.
          do i=1,npp
          if(radp(i).gt.r_act*1.e-6) then
          sum1=sum1+1.
          sum2=sum2+(radp(i)-rad(it))**2
          endif
          enddo
          if(sum1.gt.0.) then
          std(it)=sqrt(sum2/sum1)
          else
          std(it)=0.
          endif

          rad(it)=rad(it)*1.e6
          std(it)=std(it)*1.e6

         sum1=0.
         sum2=0.
          do i=1,npp
          if(radp(i).gt.r_act*1.e-6) sum1=sum1+1.
          enddo
          afra(it)=sum1/npp

           print*,it,rad(it),std(it),afra(it)

        enddo

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

         do it=1,nt
      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
         call points(time(it),afra(it),1,-4,0)
        enddo



      call gstxci(15)
      call gsplci(15) ! black
      call gspmci(15) ! black
      call setusv('LW',2000)
      call set(.2,.8,.52,.92,0.,200.,0.,1.5,1)
      call labmod('(e5.1)','(e5.1)',5,5,2,2,28,20,0)
      call perim (2,10,3,5)
         do it=1,nt
      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
         call points(time(it),std(it),1,-2,0)
        enddo


      close(20)

      enddo
      enddo

      call setusv('LW',2000)
      call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
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, '0.0', 0.016,0.,0)
      CALL plchhq(.16,0.267, '0.5', 0.016,0.,0)
      CALL plchhq(.16,0.433, '1.0', 0.016,0.,0)
      CALL plchhq(.11,0.3,'activated fraction (1)', 0.016,90.,0)

      CALL plchhq(.16,0.52, '0.0', 0.016,0.,0)
      CALL plchhq(.16,0.65, '0.5', 0.016,0.,0)
      CALL plchhq(.16,0.79, '1.0', 0.016,0.,0)
      CALL plchhq(.16,0.92, '1.5', 0.016,0.,0)
      CALL plchhq(.11,0.72,'spectral width (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