      program formlex
c
c Read input data from ft05, calculate some
c stuff like ctpoten and hp from kp, and
c format in lexical reads for tgcm runs
c
c (use this when geophysical data obtained
c  from getndcs is not available, e.g. for
c  dates later than about 6/89)
c
      parameter(mxdays=15,msgun=6,rmiss=-32767.)
      dimension iyd(mxdays),f107(mxdays),f107a(mxdays),
     +  ap(mxdays),rkp(8,mxdays),hpfkp(8*mxdays),
     +  cpfkp(8*mxdays),ibegt(4),iendt(4),itm03g(4,8*mxdays),
     +  itmp1(4),i3ht(8),itm24g(4,mxdays),day3hr(mxdays*8)
      character*32 botlab
      character*3 months(12)
      data months /'Jan','Feb','Mar','Apr','May','Jun',
     +             'Jul','Aug','Sep','Oct','Nov','Dec'/
      data i3ht/0130 , 0430 , 0730 , 1030 , 1330 , 1630 , 1930 , 2230/
c
c Read cards from ft05:
c
      read(5,*) (ibegt(i),i=1,4),(iendt(j),j=1,4),lexbday  
      write(6,"(' ')")
      write(6,"(1X,'lexical day (LEXBDAY) =',I3,'  is begin day',2I5
     +          /,' Begin time (year moday hrmin centisec) =',4I5,/,
     +         '   End time (year moday hrmin centisec) =',4I5)")
     +        lexbday,(ibegt(k),k=1,2),(ibegt(i),i=1,4),(iendt(j),j=1,4)
      write(botlab,"(a3,i3,', ',i4,' to ',a3,i3,', ',i4)")
     +  months(ibegt(2)/100),ibegt(2)-ibegt(2)/100*100,ibegt(1),
     +  months(iendt(2)/100),iendt(2)-iendt(2)/100*100,iendt(1)
      write(6,"('botlab=',/a)") botlab
      itmp1(1) = ibegt(1)
      itmp1(2) = ibegt(2)
      itmp1(3) = 100*(ibegt(3)/100) + 30
      itmp1(4) = 0
      read(5,*) ndays
      if (ndays.gt.mxdays) then
        write(6,"('>>> formlex: ndays > mxdays: ', 
     +    ' ndays=',i4,' mxdays=',i4)") ndays,mxdays
        stop 'ndays'
      endif
      write(6,"(' ')")
      do 100 i=1,ndays
        read(5,*) iyd(i),f107(i),f107a(i),ap(i),
     +            (rkp(ii,i),ii=1,8)
        write(6,"(' day ',i2,': ',i5,2f7.2,f5.2,8f5.2)")
     +    i,iyd(i),f107(i),f107a(i),ap(i),(rkp(ii,i),ii=1,8)
 100  continue
c
c Power and ctpoten from KP:
c
      ih = 0
      do 200 i=1,ndays
        itm24g(1,i) = itmp1(1)
        itm24g(2,i) = itmp1(2)+(i-1)
        itm24g(3,i) = 1200
        itm24g(4,i) = 0000
        do 205 ii=1,8
          ih = ih+1
          hpfkp(ih) = amax1(3.0, -2.78+9.33*rkp(ii,i))
          cpfkp(ih) = 29. + 11. * rkp(ii,i)          
          itm03g(1,ih) = itmp1(1)
          itm03g(2,ih) = itmp1(2)+(i-1)
          itm03g(3,ih) = i3ht(ii)
          itm03g(4,ih) = 0000
 205    continue
 200  continue
c
c Write lex file, formatted for tgcm inputs:
c 
      open(61,file='formlex.lex')
      WRITE (61,"('C  lexical day ',I3,' = year',I5,' month',I3,
     | ' day',I3)") LEXBDAY,IBEGT(1),IBEGT(2)/100,MOD(IBEGT(2),100)
      IMODA1 = MOD(IBEGT(2),100)
      LEXADD = LEXBDAY - IMODA1
      n3hrg = 8*ndays
C
       CALL LEXWRT ('C  P from Kp, P = max(0.1,-2.78+9.33*Kp)',
     |  'POWER =', 1, HPFKP, ITM03G, N3HRG, 4, LEXADD)
       CALL LEXWRT ('C  CP from Kp, CP = 29 + 11 * Kp',
     |  'CTPOTEN =', 1, CPFKP, ITM03G, N3HRG, 4, LEXADD)
       CALL LEXWRT (' ','F107A =', 1, F107A, ITM24G, ndays, 4, LEXADD)
       CALL LEXWRT (' ','F107 =', 1, F107, ITM24G, ndays, 4, LEXADD)
       CALL LEXWRT (' ','APD =', 0, AP, ITM24G, ndays, 4, LEXADD)
       CALL LEXWRT (' ','KP =', 1, rkp, ITM03G, N3HRG, 4, LEXADD)
c
c Make line plots:
c
      call opngks
      write(6,"(' ')")
      call agsetc('LABEL/NAME.','B')
      call agseti('LINE/NUMBER.',-100)
      call agsetc('LINE/TEXT.',botlab)
c
c Daily Ap:
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','AP')
      call ezy(ap,ndays,'Daily Ap Index')
      write(6,"('Plotted daily Ap')")
c
c F107:
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','F10.7')
      call ezy(f107,ndays,'Daily 10.7 cm Solar Flux')
      write(6,"('Plotted f107')")
c
c F107a:
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','F10.7a')
      call ezy(f107a,ndays,'81-day average 10.7 cm Solar Flux')
      write(6,"('Plotted f107a')")
c
      iii = 0
      do 300 i=1,ndays
      do 300 ii=1,8
        iii = iii+1
        day3hr(iii) = float(i)+float(ii-1)*3./24.
 300  continue
c     write(6,"(' ')")
c     write(6,"(' day3hr=',/(6f12.4))") (day3hr(i),i=1,ndays*8)
c     write(6,"(' ')")
c
c Kp:
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','Kp')
      call ezxy(day3hr,rkp,8*ndays,
     +  '3-hr Kp')
      write(6,"('Plotted 3-hr Kp')")
c
c Power from kp:
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','power')
      call ezxy(day3hr,hpfkp,8*ndays,
     +  'Hem. Power = max(0.1,-2.78+9.33Kp)')
      write(6,"('Plotted hem power')")
c
c Potential from kp:
c
      call agsetc('LABEL/NAME.','L')
      call agseti('LINE/NUMBER.',100)
      call agsetc('LINE/TEXT.','ctpoten')
      call ezxy(day3hr,cpfkp,8*ndays,'Crosstail Potential (29.+11.*Kp)')
      write(6,"('Plotted crosstail potential')")
c
      stop
      end
