!
      program mkfactors
!
! 2/10/06 btf: 
! Make glbmean runs varying 4 factors with 3 levels each, for total of 
! 3**4 = 81 runs. This is the first experiment with D. Nychka.
!
      implicit none
      integer,parameter ::
     |  nfac = 4,         ! number of factors
     |  nlev = 3,         ! number of valid levels for each factor
     |  nrun = nlev**nfac ! number of unique permutations (total runs) 
      real :: flevs(nfac,nlev)   ! valid levels (values) for factors
      real :: factors(nfac,nrun) ! output set of nfac levels for all runs
      character(len=16) :: facnames(nfac)
      character(len=24) :: filename
      integer :: n,i
      integer,parameter :: lu=7
!
! Output text file for use in multiinp:
      filename = 'glbmean_factors.dat     '
!
! Factor names:
      facnames = 
     |  (/'F107            ',   ! 10.7 cm solar flux
     |    'AJOUL_INP       ',   ! joule heating (denmod.F)
     |    'ALP_INP         ',   ! char energy (aurora.F, solheat.F)
     |    'EAUR_INP        '/)  ! goes into auroral prod (solheat.F)
!
! Assign valid values (levels) to each factor:
      flevs(1,:) = (/70. ,115. ,160./)
      flevs(2,:) = (/1.e10, 2.5e10, 5.e10/)
      flevs(3,:) = (/1., 2.5, 5./)
      flevs(4,:) = (/.01, .025, .05/)
!
      write(6,"(/,'mkfactors: nfac=',i4,' nlev=',i4,' nrun=',i4)")
     |  nfac,nlev,nrun
      write(6,"(/,'Valid levels for each factor:')")
      do n=1,nfac
        write(6,"('Factor ',i4,': ',(4e12.4))") n,flevs(n,:)
      enddo
!
! Make unique set of nfac levels for each run:
!
      call mkfac(nfac,nlev,nrun,flevs,factors) 
!
! Report to stdout:
      write(6,"(/,'Ordered set of factors:')")
      do n=1,nrun
        write(6,"('Run ',i3,':')",advance='no') n
        do i=1,nfac
          write(6,"(' ',a,' = ',e10.3)",advance='no') 
     |      trim(facnames(i)),factors(i,n)
        enddo
        write(6,"(' ')")
      enddo
!
! Scramble the set:
!
!     call scramble(nfac,nrun,factors)
!     write(6,"(/,'Scrambled factors:')")
!     do n=1,nrun
!       write(6,"('n=',i4,' factors=',6e12.4)") n,factors(:,n)
!     enddo
!
! Write text file suitable for multiinp:
!
      call mkinp(nfac,nrun,factors,facnames,filename,lu)

      end program mkfactors
!-----------------------------------------------------------------------
      subroutine mkfac(nfac,nlev,nrun,flevs,factors)
      implicit none
!
! Args:
      integer,intent(in) :: nfac,nlev,nrun
      real,intent(in) :: flevs(nfac,nlev)
      real,intent(out) :: factors(nfac,nrun)
!
! Local:
      integer :: k,i,n,ncon(nfac),ncyc(nfac)
!
      do i=1,nfac
        ncon(i) = nlev**(i-1) ! number of runs in which levs is constant
        ncyc(i) = nlev**i     ! number of runs to cycle through levs
      enddo
!
! Construct ordered set of all combinations of factors:
! 
      do k=1,nlev
        do i=1,nfac
          do n=(k-1)*ncon(i)+1,nrun,ncyc(i)
            factors(i,n:n+ncon(i)-1) = flevs(i,k)
          enddo
        enddo
      enddo
      end subroutine mkfac
!-----------------------------------------------------------------------
      subroutine scramble(nfac,nrun,factors)
      implicit none
!
! Args:
      integer,intent(in) :: nfac,nrun
      real,intent(inout) :: factors(nfac,nrun)
!
! Local:
      integer :: n,i,ii,iran,nscram,iscram(nrun)
      real :: facsave(nfac,nrun)
!
      facsave = factors
      iscram(:) = 0
      nscram = 5     ! scramble 5 times (iscram will contain dups)
      do n=1,nscram
        do i=1,nrun
          call random_int(iscram(i),1,nrun)
        enddo
      enddo
!
! Replace dups with numbers not appearing, so all numbers between
! 1 and nrun appear exactly once.
!
      do n=1,nrun
        do i=n+1,nrun
          if (iscram(i)==iscram(n)) then ! dup of iscram(n) at i
            do ii=1,nrun
              if (.not.any(iscram==ii)) then
                iscram(i) = ii ! replace dup w/ number not yet used
                exit
              endif
            enddo
          endif
        enddo  
      enddo
!
! Confirm that all numbers are present:
      do n=1,nrun
        if (.not.any(iscram==n)) then
          write(6,"('>>> The number ',i4,' does not appear in iscram')")
     |      n
          write(6,"('iscram=',/,(8i8))") iscram
          stop 'iscram'
        endif
      enddo
!
! Assign factor levels from scrambled indices:
      do i=1,nrun
        factors(:,i) = facsave(:,iscram(i))
      enddo
      end subroutine scramble
!-----------------------------------------------------------------------
      subroutine random_int(result,low,high)
      integer,intent(out) :: result
      integer,intent(in) :: low,high
      real :: harvest
!
      call random_number(harvest)
      result = int((high-low+1)*harvest+low)
      end subroutine random_int
!-----------------------------------------------------------------------
      subroutine mkinp(nfac,nrun,factors,facnames,filename,lu)
!
! Write factors to filename, suitable for use in perl script multiinp. 
!
      implicit none
!
! Args:
      integer,intent(in) :: nfac,nrun,lu
      real,intent(inout) :: factors(nfac,nrun)
      character(len=*) :: filename,facnames(nfac)
!
! Local:
      integer :: i,n,nn,nline=5,modnline
!
! $keyval{'F107'}   = [ ('150.','160.','170.') ];
!
      modnline = mod(nrun,nline)
      open(lu,file=filename)
      do i=1,nfac
        write(lu,"('$keyval{''',a,'''} = [(')") trim(facnames(i))
        do n=1,nrun-modnline,nline
          do nn=n,n+(nline-1)
            if (nn < nrun) then
              write(lu,"('''',e10.4,''',')",advance="no")
     |          factors(i,nn)
            else
              write(lu,"('''',e10.4,''')];')")
     |          factors(i,nn)
            endif
          enddo
          write(lu,"(' ')")
        enddo
        if (modnline > 0) then
          do n=1,modnline
            if (n < modnline) then
              write(lu,"('''',e10.4,''',')",advance="no")
     |          factors(i,nrun-modnline+n)
            else
              write(lu,"('''',e10.4,''')];')")
     |          factors(i,nrun-modnline+n)
            endif
          enddo
        endif
        write(lu,"(' ')")
      enddo  ! 1,nfac
      close(lu)
      end subroutine mkinp