C>>>> c:\christos\rdnec3c.for <<<<<<<<<<<<<
C>>> this should have the bug corrected  

C>>>>>>>>>>>>>c:\christos\rdnec3b.for <<<<<<<<<<<
C  try to fix bug in rdnec3a.for - some data taken from adjacent day -
c or days concatinated - pointed out by Nikos Christakis
C>> but I don't have rdnec3a-
c  one fix- initialize iymd





C>>>>>>>>>>>>>>>>>>>>rdnec3.for<<<<<<<<<<<<<<<<<<<<<<<<
C>>> and upgrade to rdnec2.for - to replace it when I get
C>>> auto file naming convention. ALso shear option


c   would option to add azim to bottom of plot; maybe hardware option to start



c>>>>>>>>>>>>>>>>>>>>>>>rdnec2.for<<<<<<<<<<<<<<<<<<<<
C>> daily plot, 3min avg; 4 days,per page
c>> 480 pixels across;275 down LT days, so 3hr tics not time labels                               

c>>>>>>>>>>>>>>>>>>>>rdnec1.for<<<<<<<<<<<<<<<<<<<<<<<<<<
c>>> an upgraded version of rdnec- to handle up to 1000 hts i/p, and > 1024 times
c here will treat as continuous - and assum could be than 1 run per file (with hdr line)

C>>>>>rdnec.for - the colour version -so far just 1 plot per o/p bmp

c>>>>>>>>>>.rdne.for<<<<<<<<<<<<<<
C>>> read arecibo Ne(h) (from Christos/Qihou Zou)
C file names NEyymm.asc
      SAVE
      PARAMETER (NHMAX=1000)
      character*128 prefix,fname,string,hdr,strng
      character*128 Lgnd
      dimension edens(NHMAX),month(12),oneday(700,480,2),avgaz(480,2)
      logical ex
      data month/0,31,59,90,120,151,181,212,243,273,304,334/
      data iskip/1/,iframe/0/,istop/0/,nframe/4/ !nframe per pg,chng to 2 if > 300 hts  
      character*3 ascale(2)
      
      data ascale/'Lin','Log'/,iazplt/1/ !iazplt =1 if want azplt tacked on
      common/diagno/iprt
      if(iazplt.ne.1)iazplt=0 !plot azim if =1 (and use to make space for !)
      print*,'diagnostic o/p(0=no,1=yes,2=more..)?'
      read*,iprt
      print*,'overlay time and ht. grid( 0=no,1=yes)?'
      read*,igrid
      if(igrid.ne.1)igrid=0
      if(iskip.eq.0)goto 30
      iskip=0
      print*,'what prefix for file(e.g. d:\christos\ne)?'
      read'(a128)',prefix
      Lnp=lentrim(prefix)
      print*,'NOTE!! THIS PROGRAM ASSUMES THAT data records in a file',
     1' *EQUALLY* spaced!!'
      print*,'what yymm (or yymmdd) to start(e.g. 9108 or 891128)?'
      read*,jdate1
      if(jdate1.lt.9999)then
        print*,'assume yymm type file names'
        ityp=1
        jyr1=jdate1/100
        jmon1=mod(jdate1,100)-1
      else
        print*,'assume yymmdd file names'
        ityp=2
        jyr1=jdate1/10000
        jmon1=mod(jdate1/100,100)
        jdm1=mod(jdate1,100)-1  
      endif  
      print*,'Noon centred plots (1) or Midnight(2) or ',
     1'Christos day(centre 18LT) (3)?'
      read*,jopt
      print*,'add azimuth plot at bottom (1=yes)?'
      read*,iazplt
      if(iazplt.ne.1)iazplt=0 
      nmiss=0
      print*,'What height range (km) to plot(e.g. 50,180)?'
      read*,ikm1,ikm2 

      print*,'what scale type (lin=1, log=2, lin grd(3),log +grd(4),',
     1' log-log grd (5))?'
      read*,iscl
      if(iscl.lt.1.or.iscl.gt.5)then
        print*,'no such scale type!'
        iscl=2
      endif  
      print*,'[NOTE NEXT: "good" scales are: 3,5 for option (2);',
     1' 0,4 for option (4); '
      print*,'and 0.,0.5 for option (5).]'
      print*
      print*,'Now enter log10 range(either scale, e.g. 3.1,5.1):'
      read*,vmin,vmax
      if(vmin.eq.vmax)then
        print*,'null range not allowed - auto scale not available'
        stop
      endif  
      if(vmin.gt.0.001.and.vmax.gt.0.001.and.
     1   (iscl.eq.1.or.iscl.eq.3))then  !linear  !BUT !!! if linear willnot go -ve!
        vmin=10.**vmin
        vmax=10.**vmax
      endif
      do it=1,480
        avgaz(it,1)=0.
        avgaz(it,2)=0
      do kh=1,700
        oneday(kh,it,1)=0.
        oneday(kh,it,2)=0.
      enddo
      enddo
      call plots(ia1,ia2,ia3,ia4,2)
   20 if(ityp.eq.1)then   !yymm names
        jmon1=jmon1+1
        if(jmon1.gt.12)then
          jmon1=1
          jyr1=jyr1+1
          if(jyr1.gt.99)then
            print*,'cant deal with yr rollover >99 yet'
            stop
          endif
        endif    
        write(string,'(2i2.2)')jyr1,jmon1
        fname=prefix(1:Lnp)//string(1:4)//'.asc'
        Lnf=lentrim(fname)
        print*,'Look for file=',fname(1:Lnf)
      else   !yymmdd names
        jdm1=jdm1+1
        if(jdm1.gt.month(jmon1))then  !non leap; so check below
          if(jmon1.eq.2.or.mod(jyr1,4).ne.0)then
            jmon1=jmon1+1
            jdm1=1
          elseif(jmon1.eq.2.and.mod(jyr1,4).eq.0.and.jdm1.gt.29)then
            jmon1=jmon1+1
            jdm1=1
          endif  
          if(jmon1.gt.12)then
            jmon1=1
            jyr1=jyr1+1
            if(jyr1.gt.99)then
              print*,'cant deal with yr rollover > 99 yet'
              stop
            endif
          endif
        endif
        write(string,'(3i2.2)')jyr1,jmon1,jdm1  
        fname=prefix(1:Lnp)//string(1:6)//'.asc'
        Lnf=lentrim(fname)
        print*,'Look for file=',fname(1:Lnf)
      endif  
      kskip=1
    5 inquire(file=fname,exist=ex)
      if(ex)then
        open(20,file=fname,form='formatted',status='old',recl=50000)
        nmiss=0
      else
        nmiss=nmiss+1
        Lnf=lentrim(fname)
        print*,'cant find file=',fname(1:Lnf)
        if(nmiss.gt.6)then
          print*,'.gt.6 missing files in a row:',
     1    ' get new file name(-1),stop(0),or look n further?'
          read,n
          if(n.eq.0)then
            if(igot.eq.1)then !something stored which need avg & plot,then stop
              print*,'dump last partial plot,and stop'
              istop=1
              iymd=9999999 !force day change
              jLian=9999999 ! and for midnight centred
              ixtos=9999999
              goto 40
            endif  
c... if just 1 day in file, iframe will stay=0,because no day change; so count days too          
            if(iframe.ne.0)call plot(0,0,-1) !dump last partial;of 1 day(iframe -ve to start
            stop
          endif  
          if(n.lt.0)then
            print*,'what new file name e.g. f:\ne9011.asc (or q=quit)?'
            read'(a40)',fname
            if(fname.eq.'q'.or.fname.eq.'Q')stop
            goto 5
          endif
          nmiss=nmiss-n
        endif
        goto 20
      endif
        
   25  read(20,*,end=900)nrecs,nhts,hkm1,stepkm,zenang
        if(nhts.gt.NHMAX)then
          print*,'too many hts for dims!, max=',nhmax,
     1      ' but header says #=',nhts
          stop
        endif
        Lskip=1
        indh1=(float(ikm1)-hkm1)/stepkm
        if(indh1.lt.1)indh1=1
        indh2=(float(ikm2)-hkm1)/stepkm+1
        if(indh2.gt.nhts)indh2=nhts
        nh=indh2-indh1+1 ! # data hts
        if(nh.gt.300)nframe=2 !2 frames per page rom now on
        ny=(ikm2-ikm1)/stepkm !# pixels in height axis
        if(nh.gt.700)then
          print*,'Requested ht. range has requires > 700 pixels!'
          stop
        endif            

c... clear for optional filter  
        print 1311,nrecs,nhts,hkm1,stepkm,zenang
 1311 format('hdr:#recs=',i6,' #hts=',i3,' km1,step=',2f6.2,
     1  ' zen=',f6.1) 
c      kskip=1
      Lskip=1  !each hdr
      irec=0

c NOTE: for some reason azim =(~) 180 - 540, step 30 (at least for 891002)
c       and sweeps back and forth between these - i.e. doesn't go "around"
c   usual rec spacing about 1.2 min; some large gaps        

   30 read(20,*,end=100,err=800)irecno,idate,hrslt,azim,npulse,
     1 (edens(k),k=1,nhts)
     
      irec=irec+1
      idm=idate/10000
      imon=mod(idate/100,100)
      iyr=mod(idate,100)
      idy=month(imon)+idm
      iymd=iyr*10000+imon*100+idm  ! added April '08  fix?
      if(imon.gt.2.and.mod(iyr,4).eq.0)idy=idy+1 
      ihr=hrslt
      imin=mod(hrslt,1.)*60
      if(jopt.eq.1)it=(ihr*60+imin)/3+1
      if(jopt.eq.2)it=mod((ihr+12)*20+imin/3,480)+1 !midnight centred
      if(jopt.eq.3)it=mod((24+ihr-6)*20+imin/3,480)+1   !want:if ihr=6,it=1
      iymd=iyr*10000+imon*100+idm
      jLian=iyr*1000+idy+ihr/12 !not really Julian day,but it will do for midnight centred
      khr=ihr
      kdy=idy
      kyr=iyr
c...work out Christos day - starts at 6LT on named day
      khr=khr-6
      if(khr.lt.0)then
        khr=khr+24
        kdy=kdy-1
        if(kdy.lt.1)then
          Ly=365
          if(mod(kyr,4).eq.0)ly=366
          kyr=kyr-1
          kdy=Ly
        endif
      endif
      ixtos=kyr*1000+kdy          
      if(kskip.eq.1)then  !save 1st data and set up for plot
        iyrL=iyr
        imonL=imon
        idmL=idm
        idy1=idy
        ihr1=ihr
        imin1=imin
        iymdL=iymd
        kskip=0
        jLianL=jLian
        ixtosL=ixtos
      endif     
      if(lskip.eq.1)
     1  print 3344,iyrL,imonL,idmL,ihr1,imin1
 3344   format('Hdr d/t ',i2.2,2('/',i2.2),
     1    1x,2i2.2,'LT')   
      Lskip=0
      if(irec.ne.irecno)then
        print*,'wrong internal record# in rec#',irec,
     1   '; ascii rec# was ',irecno
        stop
      endif
   40 if((jopt.eq.1.and.iymdL.ne.iymd).or.
     1    (jopt.eq.2.and.jLianL.ne.jLian).or.
     1    (jopt.eq.3.and.ixtosL.ne.ixtos))then
        iframe=iframe+1
        ixf=15
        ix0=ixf+ny+2
        iy0=0
        iyf=480+1   !480 3min vals/day
        call setres(1,1) !make aLL CHARS 6x7 so I can easily squeeze more plots onto
                         !page later if I want
        ixoff=380*((iframe-1)/2)
        iyoff=36+mod(iframe-1,2)*490
        if(ix0+ixoff+30.gt.760)then !going to go off page,so dmp & restart frame count 
          call plot(0,0,-1)
          iframe=1
          ixoff=350*((iframe-1)/2)
          iyoff=36+mod(iframe-1,2)*490
        endif  
        print*,'test : print frame#,xoff,yoff',iframe,ixoff,iyoff
        call origin(ixoff,iyoff)
        call setcol(0)
        call rect(ix0,iy0,ixf,iyf)
        if(igrid.eq.1)then  !make it a manual ctrld option April 6 '08
c  try arbitrary grid(Jul'07)- later make 5 km step, 3hr step?
          do i=ix0,ixf,-50
            call line(i,iy0,i,iyf)
          enddo
          do i=iy0,iyf,100
            call line(ix0,i,ixf,i)
          enddo      
        endif
        if(jopt.eq.1)write(string,'(''One LT Day'')')  !convrtd to upper case for tiny chars
        if(jopt.eq.2)write(string,'(''One Midnight centred LT day'')')  !ditto
        if(jopt.eq.3)write(string,
     1         '(''One Christos day (06 LT - 06 LT)'')')
        nlen=128
        call strnln(string,nlen)
        call cymbol(ix0+10+15*iazplt,240,string,0,0,nlen)
        do iy=1,480,60 !every 3 hrs
         call line(ix0,iy,ix0+3,iy)        
        enddo
c..have to put ht tic here because at end too
        if(mod(iframe,2).eq.1)then !LHS plots
          write(string,'(''Height (Km)'')')
          call cymbol((ix0+ixf)/2,-30,string,0,3,11)
        endif  
          istep=(ikm2-ikm1)/8
          istep=(istep/10+1)*10
     
          do km=20,500,10
            ix=ix0-(km-ikm1)/stepkm
            if(ix.gt.ixf.and.ix.lt.ix0)then !on scale
               call line(ix,iy0,ix,iy0-2)
               call line(ix,iyf,ix,iyf+2)
               if(mod(iframe,2).eq.1.and.mod(km,20).eq.0)
     1             call rumber(ix+3,iy0-3,km,0,0,0)
            endif  
        enddo    
        wmin=vmin !for label
        wmax=vmax
        if(iscl.eq.1)then
          wmin=log10(wmin)
          wmax=log10(wmax)
        endif  
        if(jopt.eq.1)then
          if(iscl.lt.3)write(hdr,7719)iyrL,imonL,idmL,ascale(iscl),
     1      wmin,wmax
          if(iscl.eq.3)write(hdr,7719)iyrL,imonL,idmL,ascale(iscl-2),
     1      wmin,wmax,', Gradient    '
          if(iscl.eq.4)write(hdr,7719)iyrL,imonL,idmL,ascale(iscl-2),
     1      wmin,wmax,', Log_Grad    '
          if(iscl.eq.5)write(hdr,7719)iyrL,imonL,idmL,ascale(iscl-2),
     1      wmin,wmax,',Log_Log_Grad '   
 7719   format('Arecibo Ne(h):  day=',i2.2,'/',i2.2,'/',i2.2,
     1   1x,a3,' log10 rng:',f5.2,'-',f5.2,a14)
     
        elseif(jopt.eq.2.or.jopt.eq.3)then
              !midnight centred or "Christos day"; give 1st and 2nd
              !      LT day#s, 1st may be 0; but makes sense
              !   rather than getting y,m,d; e.g. 99/12/0-1 rather than 99/11/30-1
                
          if(iscl.lt.3)write(hdr,7720)iyr,imon,idm-1,idm,
     1      ascale(iscl),wmin,wmax
          if(iscl.eq.3)write(hdr,7720)iyr,imon,idm-1,idm,
     1      ascale(iscl-2),wmin,wmax,', Gradient   '
          if(iscl.eq.4)write(hdr,7720)iyr,imon,idm-1,idm,
     1      ascale(iscl-2),wmin,wmax,', Log_Grad   '
          if(iscl.eq.5)write(hdr,7720)iyr,imon,idm-1,idm,
     1      ascale(iscl-2),wmin,wmax,',Log_Log_Grad'
 7720   format('Arecibo Ne(h):  day=',i2.2,'/',i2.2,'/',i2.2,'-',i2.2,
     1   1x,a3,' log10 rng:',f5.2,'-',f5.2,a13)
        endif
        Lnh=lentrim(hdr)
        call setcol(0) !black
        call cymbol(ixf-6,(iy0+iyf)/2,hdr,0,0,Lnh)
        
      
        ioff=(hkm1-ikm1)/stepkm  !in case select 1st ht < than lowest data(=56Km)
        if(ioff.lt.0)ioff=0
        igot=1 !got something for this day
        if(iscl.lt.3)then
          do jt=1,480  !3 min index
            if(iazplt.eq.1)then
              v=avgaz(jt,2)
              if(v.gt.0.)then
                ind=12- mod(180.+avgaz(jt,1)/v,360.)/30.
                call setcol(0)
                call csym(ix0+5+ind,jt,1)
              endif  
            endif
          do kh=1,nh

            ix=ix0-ioff-kh
            iy=jt
            v=oneday(kh,jt,2)
            if(v.gt.0.)then
              v=oneday(kh,jt,1)/v !avg Ne
              if(iscl.eq.2)v=log10(v)
              icol=13-12*(v-vmin)/(vmax-vmin)+1
              if(icol.gt.12)icol=12
              if(icol.lt.1)icol=1    
              call setcol(icol)
              call csym(ix,iy,1)
            endif
            oneday(kh,jt,1)=0.
            oneday(kh,jt,2)=0. 
          enddo  
            avgaz(jt,1)=0.
            avgaz(jt,2)=0.
          enddo
        else
          do jt=1,480
            if(iazplt.eq.1)then
              v=avgaz(jt,2)
              if(v.gt.0.)then
                ind=12-mod(180.+avgaz(jt,1)/v,360.)/30.
                call setcol(0)
                call csym(ix0+5+ind,jt,1)
              endif  
            endif
          
          do kh=2,nh

            ix=ix0-ioff-kh
            iy=jt
            v1=oneday(kh-1,jt,2)
            v2=oneday(kh,jt,2)
            if(v1.gt.0..and.v2.gt.0.)then
              v=oneday(kh,jt,1)/v2 -oneday(kh-1,jt,1)/v1
              vv=-1
              if(v.gt.0.)vv=log10(v)
              if(iscl.eq.3)then
                icol=13-12*(v-vmin)/(vmax-vmin)+1
                if(icol.gt.12)icol=12
                if(icol.lt.1)icol=1    
                call setcol(icol)
                call csym(ix,iy,1)
              elseif(iscl.eq.4.and.v.gt.0.0)then
                v=log10(v)
                icol=13-12*(v-vmin)/(vmax-vmin)+1
                if(icol.gt.12)icol=12
                if(icol.lt.1)icol=1    
                call setcol(icol)
                call csym(ix,iy,1)
              elseif(iscl.eq.5.and.v.gt.0.0)then
                vv=log(v)
                if(vv.gt.0.0)then
                  v=log10(vv)
                  icol=13-12*(v-vmin)/(vmax-vmin)+1
                  if(icol.gt.12)icol=12
                  if(icol.lt.1)icol=1    
                  call setcol(icol)
                  call csym(ix,iy,1)
                endif  
              endif
            endif
            oneday(kh-1,jt,1)=0. !clear ready for next
            oneday(kh-1,jt,2)=0.
          enddo  
            avgaz(jt,1)=0.
            avgaz(jt,2)=0.
          enddo        
          oneday(nh,jt,1)=0.
          oneday(nh,jt,2)=0.  
        endif  
        if(iframe.eq.nframe.or.istop.eq.1)then
          iframe=0
          call plot(0,0,-1)
          if(istop.eq.1)stop
        endif  
        igot=0
        if(iymd.eq.9999999)goto 20   !APR'08 fix
      endif
      idmL=idm
      imonL=imon
      iyrL=iyr
      iymdL=iymd
      jLianL=jLian
      ixtosL=ixtos
      avgaz(it,2)=avgaz(it,2)+1
      avgaz(it,1)=avgaz(it,1)+azim
      do kh=1,nh
        v=edens(kh+indh1-1)
        if(v.gt.0.)then 
          oneday(kh,it,1)=oneday(kh,it,1)+v
          oneday(kh,it,2)=oneday(kh,it,2)+1
        endif
      enddo    
      igot=1  !show that there is some data to avg
c...skip legend for now ,.may put it at bottom
c        ix1= 20
c        iy1= iyf+20
c        if(iy1+50.gt.1020)then
c          print*,'legend off page!'
c        else
c          lgnd='Lg10Ne'
c        endif  
c        step=(wmax-wmin)/12.
c        start=wmin
c        L10=0 !dont do log10 in subr.
c        if(iscl.eq.1)then
c          start=vmin
c          step=(vmax-vmin)/12.
c          L10=1 !in legend get label value= log10( start + i-1 * step)
c        endif  
c        call legend(Lgnd,L10,start,step,12,ix1,iy1)
c        call plot(0,0,-1)
c        kskip=1
c      endif
      if(irec.eq.nrecs)goto 25 !try for next header
      goto 30 !continue reading this run
  100 print*,'EOF OK @ rec#',irec
      close(20)
      if(irec.ne.nrecs)then
        print*,'wrong# recs found! hdr=',nrecs,
     1  ' but read ',irec
      endif
      goto 20
      
  800 print*,'read err in data rec#',irec-1
      stop      
  900 print*,'EOF,generate next name'
      iymd=9999999  !APR '08 FIX
      close(20)    !APR '08 FIX
      goto 40      !APR '08 FIX
      end
      
c=============================================================      
      
      subroutine legend(lgnd,L10,start,step,ncol,ix1,iy1)
      SAVE
      
C.... iy1,ix1= LH-top cnr of legend; must have 5x6 + 15 + 5 pixels for legend at
c       RHS of plot
c>>>>>>>>> iyy=1,1024(across) ixx=1-768(down)
      
c.... colour scale.........
      character*131 lgnd,strng,hdr131
      call setres(1,1) !small char
      if(iy1.gt.960)then
        print*,'not enough room for legend,! max y>960,so make=960'
        iy1=960
      endif
      call setcol(0) !black
      ln=lentrim(lgnd)
      if(ln.gt.6)ln=6 !max 6 chr for legend title      
      call cymbol(ix1+10,iy1+7,lgnd(1:ln),0,0,ln)
      do k=1,ncol !going down page (max val to min)
        ixx=ix1+15*(k-1)+22
        call  setcol(k)
        call fill(ixx,iy1+5,ixx+15,iy1+20)
        call setcol(0)
c        call setres(1,1) !6 pixel wide char
        if(L10.eq.0)write(strng,'(f5.2)')(ncol-k)*step+start
        if(L10.eq.1)write(strng,'(f5.2)')log10((ncol-k)*step+start)
        call symbol(ixx+14,iy1+25,strng,0,0,5)
      enddo
      return
      end

      
      
      
      
      
      
      
      
      
      

