program pset

!  This program creates the files pdmethod1, pdmethod2, pd8method1 and
!  pd8method2 that specify the model levels to the program RADINT
!
!  Command line options are the number of levels and the type of spacing,
!  -t n (default) for normal model formula
!  -t e for equally spaced
!  -t l for GFDL LBL set
   use getopt_m
   use vertutils_m
   implicit none
   integer :: n, np
   real, allocatable, dimension(:) :: sig, sigh
   integer :: nopt, opt, k
   character(len=80) :: optarg
   real, parameter :: p0 = 1013.25
   logical :: normal = .true., equal = .false., lbl=.false.

   n = 0

   do 
      call getopt("t:n:",nopt,opt,optarg)
      if ( opt == -1 ) then
         exit ! End of options
      end if
      select case ( char(opt) )
      case ( "n" )
         read(optarg,*) n
      case ( "t" )
         select case ( optarg )
         case ( "n" )
            normal = .true.
         case ( "e" )
            equal = .true.
            normal = .false.
         case ( "l" )
            lbl = .true.
            normal = .false.
         case default
            print*, " Unknown argument for -t "
            stop
         end select
      case default
         print*, " Unknown option "
         stop
      end select
   end do

   if ( n == 0 ) then
      print*, "Error, need to specify number of levels"
      stop
   end if

   np = n+1
   allocate ( sig(n), sigh(np) )

   if ( normal ) then
      call setsig ( sig, sigh, n )
   else if ( lbl ) then
      call setsig_lbl ( sig, sigh, n )
   else if ( equal ) then
      do k=1,np
         sigh(k) = 1.0 - (k-1)/real(n)
      end do
      sig = 0.5*(sigh(1:n)+sigh(2:np))
   end if

   print*, ' SIG ', sig
   print*, ' SIGH ', sigh
   open(unit=1,file='pdmethod1')
   open(unit=2,file='pd8method1')
   open(unit=3,file='pdmethod2')
   open(unit=4,file='pd8method2')
   write(1,1000) (sigh(k)*p0,k=np,1,-1)
   write(2,1000) (0.8*sigh(k)*p0,k=np,1,-1)
   write(3,1000) 0.,(sig(k)*p0,k=n,1,-1),p0
   write(3,1000) (sigh(k)*p0,k=np,1,-1)
   write(4,1000) 0.,(0.8*sig(k)*p0,k=n,1,-1),p0
   write(4,1000) (0.8*sigh(k)*p0,k=np,1,-1)
 1000 format(5e16.9)

end program pset