program radint

!     This program is a modified version of the GFDL supplied program for 
!     producing the data for the radiation code block data.
!     Martin Dix 1991-04-00

   use getopt_m
   implicit none

   integer, parameter :: r8 = selected_real_kind(15)
   integer, parameter :: nmax = 200
   real(kind=r8), dimension(nmax) :: pd, pd8
   real(kind=r8), dimension(4) :: cstd = (/ 0.5, 1.0, 2.0, 4.0 /)
   character(len=10) :: suffix, cval, clow, chi
   integer :: nlev, nlev2, ntype, nmodel, n, ntap, k, ir, npurp, nkkk, iq
   integer :: itin, itin1, itin2, itin3
   real(kind=r8), parameter :: co2std = 330.0
   real(kind=r8) :: ratio, co2mix, ccomp
   integer :: nopt, opt, ios
   character(len=80) :: optarg

   do 
      call getopt("c:",nopt,opt,optarg)
      if ( opt == -1 ) then
         exit ! End of options
      end if
      select case ( char(opt) )
      case ( "c" )
         read(optarg,*) co2mix
      case default
         print*, " Error: Unknown option "
         print*, "Usage: radint -c co2conc (ppmv)"
         stop
      end select
   end do
        
   open(41,file='pdmethod2')   
   open(42,file='pdmethod1')   
   open(43,file='pd8method2')   
   open(44,file='pd8method1')   
   open(10,file='co2_data') 
   open(16,file='profiles')

!  Find the number of levels from the pdmethod1 file. This was written with
!  format 5e16.9.
   pd = 0.0
   do k=1,nmax,5
      read(42,fmt='(5e16.9)',iostat=ios) pd(k:k+4)
      if ( ios /= 0 ) exit
   end do
   rewind(42)

   ! Last non-zero value is the last real one
   do n=k,1,-1
      if ( pd(n) /= 0. ) exit
   end do

   nlev = n-1
   print*, "Calculating CO2 for", nlev," level model"

   nlev2=nlev+2
   ntype=0
   read (41,*) (pd(k),k=1,nlev2)
   read (43,*) (pd8(k),k=1,nlev2)
!
   ratio = co2mix / co2std
   if ( ratio < 0.5 .or. ratio > 4.0 ) then
      print*, "Error: Datafiles for concentations < 165 or > 1320 are not available"
      stop
   endif

   write (*,*)  'The CO2 mixing ratio is',co2mix,' ppmv'
   write(10,'(i3)')  nlev
!  Write the sigma levels
   write(10,'(6f10.5)') (pd(k)/pd(nlev2),k=2,nlev2-1)
   write(10,'(e15.5)') co2mix*1e-6

   call ptz(nlev,ntype,nmodel,pd,pd8)  
!---use 1 tape (no interpolation) if the mixing ratio is sufficienlly
!   close to that standard value
!---also dont interpolate if the mixing ratio is sufficiently close
!   to multiples (0.5,2.0,4.0) of the standard value
   ntap = 2
   do n=1,4
      ccomp = abs(co2mix-cstd(n)*co2std)
      if ( ccomp <= 1.0e-4 ) ntap=1
   end do
   print*,  ' NTAP ', ntap

!  Create componenents of file names
   if ( ntap == 1 ) then
      if ( abs(ratio-0.5) .lt. 1e-4 ) then
         write(cval,"(i4)") nint(0.5*co2std)
      else if ( abs(ratio-1.) .lt. 1e-4 ) then
         write(cval,"(i4)") nint(co2std)
      else if ( abs(ratio-2.) .lt. 1e-4 ) then
         write(cval,"(i4)") nint(2*co2std)
      else if ( abs(ratio-4.) .lt. 1e-4 ) then
         write(cval,"(i4)") nint(4*co2std)
      else
         print* , ' Error setting cval '
         stop
      end if
      cval = adjustl(cval)
      print*, " CVAL ", cval
   else if ( ntap == 2 ) then
      if ( ratio .gt. 0.5 .and. ratio .lt. 1. ) then
         write(clow,"(i4)") nint(0.5*co2std)
         write(chi,"(i4)")  nint(co2std)
      else if ( ratio .gt. 1. .and. ratio .lt. 2. ) then
         write(clow,"(i4)") nint(co2std)
         write(chi,"(i4)")  nint(2*co2std)
      else if ( ratio .gt. 2. .and. ratio .lt. 4. ) then
         write(clow,"(i4)") nint(2*co2std)
         write(chi,"(i4)")  nint(4*co2std)
      else
         print* , ' Error setting clow, chi '
         stop
      end if
      clow = adjustl(clow)
      chi = adjustl(chi)
      print*, " CLOW, HI ", clow, chi
   end if

!  Start of loop over ir
   do ir=1,4

      rewind 41
      rewind 42
      rewind 43
      rewind 44
      if ( ir == 1 ) then
         npurp=1
         nkkk=3
         suffix = '490850'
      else if ( ir == 2 ) then
         npurp=3 
         nkkk=3
         suffix = '490670'
      else if ( ir == 3 ) then
         npurp=3 
         nkkk=3
         suffix = '670850'
      else if ( ir == 4 ) then
         npurp=3
         nkkk=1
         suffix = '43um'
      end if
!

!***open lbl co2 transmission functions prior to executing the interpol-
!   ation pgm. this will be user-dependent. additionally, it depends
!   on ir (freq. range) and on ratio (co2 amt). the files below 
!   assume that ratio lies between 1 and 2, and that the file and
!   directory names are as specified below.

      if (ntap == 1) then
         open(20,file=trim(cval)//'ppmv'//trim(suffix),status='old')
      else if (ntap == 2) then
         open(20,file=trim(clow)//'ppmv'//trim(suffix),status='old')
         open(21,file=trim(chi)//'ppmv'//trim(suffix),status='old')
      endif

!***open output files for interpolated co2 transmission fctns
      if (npurp == 1) then
         open (22,file=trim(suffix)//'m2')
         open (23,file=trim(suffix)//'m1')
         open (24,file=trim(suffix)//'m2pd8')
         open (25,file=trim(suffix)//'m1pd8')
      endif
      if (npurp == 2) then
         open (22,file=trim(suffix)//'m1')
         open (23,file=trim(suffix)//'m1pd8')
      end if
      if (npurp == 3) then
         open (22,file=trim(suffix)//'m2')
         open (23,file=trim(suffix)//'m2pd8')
      endif
      call co2int(nlev,ir,npurp,nkkk,ntap,ratio)  
      close (20) 
      if (ntap == 2) then
         close (21)
      endif

      if (ir == 1) then
         rewind 22
         rewind 23
         rewind 24
         rewind 25
         itin=22
         itin1=24
         itin2=23
         itin3=25
         iq=ir
         open (12,file='bd3m2')
         open (11,file='bd3m1')
         call co2ins(nlev,itin,itin1,iq)
         call co2in1(nlev,itin2,itin3,iq)
      else
         rewind 22
         rewind 23
         itin=22
         itin1=23
         if (ir.ne.4) then
            iq=ir
         else
            iq=5 
         endif
         if (ir == 2) then
            open (12,file='bd2m2')
         endif
         if (ir == 3) then
            open (12,file='bd4m2')
         endif
         if (ir == 4) then
            open (12,file='bd5m2')
         endif
         call co2ins(nlev,itin,itin1,iq)
      endif
      ! Clean up the working files
      close(12,status='delete')
      close(22,status='delete')
      close(23,status='delete')
      if ( npurp == 1 ) then
         close(11,status='delete')
         close(24,status='delete')
         close(25,status='delete')
      end if
   end do  ! Loop over ir
   close(16,status='delete')

end program radint