! 2019/01/28 - Padova - Michele Trabucchi - contact: michele.trabucchi(at)unipd.it
! giraffe: Grid Interpolation Routine And Flexible Frequency Estimator

program giraffe

   use grid_structure
   use constants

   implicit none

   integer, parameter :: narg = 8 ! the exact number of argument that have to be passed in the most general case (case B below)

   integer :: narg_in, &                     ! number of arguments detected
              lMI, lMO, &  ! length of main input and output file names (paths), respectively
              lTF, & ! length of a table-file name (path)
              ltmp, jj, &
              ll, l00, l0, l, j, i, &   ! used as running indices and to record lengths and positions of various strings in main input
              jja, jjr, jji, i1, i2, &
              ret, &
              ierr   ! error code

   character(len = 7) :: gret ! this is a flag that will be printed to the output file for each interpolation
                              ! it is a 7-digit number, each digit corresponding to one of the interpolation parameters
                              ! 1st digit, log(L/Lsun):
                              !    is 0 if no extrapolation occurred in this parameter
                              !    is 1 if a downward extrapolation occurred at least once
                              !    is 2 if an upward extrapolation occurred at least once
                              !    is 3 if both an upward and downward extrapolation occurred
                              ! 2nd digit, log(T_eff / Kelvin):
                              ! 3rd digit, M/Msun
                              ! 4th digit, Mc/Msun
                              ! 5th digit, X
                              ! 6th digit, Z
                              ! 7th digit, C/O
                              !    is 0 if no extrapolation occurred
                              !    is 1 if a downward extrapolation occurred at least once
                              !    is 1 if an upward extrapolation occurred at least once

   integer :: icL  = 1, & ! default number of column with log(L/Lsun) in main input file
              icT  = 2, & ! default number of column with log(T_eff / Kelvin) in main input file
              icM  = 3, & ! default number of column with M/Msun in main input file
              icMc = 4, & ! default number of column with M_c/Msun in main input file
              icX  = 5, & ! default number of column with X in main input file
              icZ  = 6, & ! default number of column with Z in main input file
              icCO = 7    ! default number of column with C/O in main input file

   integer :: iML, iZ, iX, iCO, iM

   character(len = 256) :: main_input, & ! name (path) of main input file
                           arg           ! string used for temporary storage of input arguments
   character(len = 269) :: main_output   ! name (path) of main output file
   character(len = *), parameter :: output_ext = '_interpolated' ! this string is appended to the name (path) of the main input file
                                                                 ! to create the name (path) of the main output file
   character(len = 1024) :: line ! string used for temporary storage of lines read from the main input file
   character(len = 2048) :: oline ! string used for temporary storage of lines read from the main input file
   character(len = 256) :: tstr
   character(len = 256) :: tablefile

   real(8) :: tmp ! used to store temporary data
   real(8), dimension(200) :: cols ! used to store data from each column of a line from the main input file

   real(8) :: lLi = UNDEF, & ! log(L/Lsun) read from a line in the main input file, initialized to UNDEF
              lTi = UNDEF, & ! log(T_eff / Kelvin) read from a line in the main input file, initialized to UNDEF
              Mi  = UNDEF, & ! M/Msun read from a line in the main input file, initialized to UNDEF
              Mci = UNDEF, & ! M_c/Msun read from a line in the main input file, initialized to UNDEF
              Xi  = UNDEF, & ! X read from a line in the main input file, initialized to UNDEF
              Zi  = UNDEF, & ! Z read from a line in the main input file, initialized to UNDEF
              Zri = UNDEF, & ! Zref computed from Zi, initialized to UNDEF
              COi = UNDEF, & ! C/O read from a line in the main input file, initialized to UNDEF
              fCOi = UNDEF    ! = fC = fCO (assuming fixed oxygen), computed from COi and COSUN

   real(8), dimension(ntot) :: TT, TT1, TT2 !,TT0
   real(8), dimension(2, ntot) :: tZ, tX, tfCO, tM
   real(8), dimension(3, ntot) :: tML
   real(8) :: dn, w1, w2
   integer :: t1, t2, clock_rate, clock_max

   logical :: input_exist, comment ! booleans used for check of existence of main input file and whether its lines are comments
   logical :: use_given_Z

!-----------------------------------------------------------------------------------------------------------------------------------

10 format('File not found: ' ,a)
11 format('Error while reading file: ',a)
12 format('GIRAFFE terminated normally - Elapsed CPU time: ')
20 format('BAD INPUT ARGUMENTS')
21 format('Use: giraffe [input file path] [col.logL] [col.logTe] [col.M] [col.Mc] [col.X] [col.Z] [col.C/O]')
22 format('where col.xyz is the number of the column, within the input file, where quantity xyz is listed')

77 format(a,'Zref',f6.4,'/X',f6.4,'/fCO',sp,f7.4,'/fC',f7.4,'/fN',f7.4,'/alpha',ss,f6.4,'/M',f6.4,'/Mc_AGB_',a5,'.tsq')

79 format(256x)
!80 format(3x,f9.2,1x,s1p,e10.3,ss,1x,0pf9.2)
81 format(1x,s1p,1(e14.6e3,1x,e14.6e3,1x,e14.6e3,1x),ss)

!-----------------------------------------------------------------------------------------------------------------------------------

   use_given_Z = .false.    !!! DON'T CHANGE
   call system_clock(t1, clock_rate, clock_max)

   snN(1) = 0.D0
   snN(2) = 0.D0
   snML(:) = gML(:)
   conversion = 2 ! DO NOT CHANGE!!!
   write(tstr, 79)

   ! Read input Arguments
   ! Two modes are allowed:
   ! A. pass the name (path) of the input file - it has to be structured in columns, where columns from 1 to 7 are as follows:
   !    col. 1 - contains luminosity: log(L/Lsun)
   !    col. 2 - contains effective temperature: log(T_eff / Kelvin)
   !    col. 3 - contains total mass: M/Msun
   !    col. 4 - contains core mass: M_c/Msun
   !    col. 5 - contains hydrogen mass fraction: X
   !    col. 6 - contains metallicity: Z
   !    col. 7 - contains carbon-to-oxygen ratio: C/O
   ! B. pass the name (path) of the input file AND the number of the column, within the input file,
   !    that contains values of the following parameters, in this order:
   !    log(L/Lsun), log(T_eff / Kelvin), M/Msun, M_c/Msun, X, Z, C/O
   ! Example:
   ! case A: ./giraffe.exe user/data/input_file.in
   ! case B: ./giraffe.exe user/data/mixed_input_file.in 1 5 3 7 10 11 14
   !   which means that log(L/Lsun) is in column 1, log(T_eff / Kelvin) is in column 5, M/Msun is in column 3, etc.
   narg_in = IArgC()
   if (narg_in == 1) then
      ! case A
      call GetArg(1, main_input)
      if (main_input(1:4) == 'help') then
         write(*,21)
         stop
      endif
   else if (narg_in == narg) then
      ! case B
      call GetArg(1, main_input)
      call GetArg(2, arg)
      ll = len_trim(arg)
      read(arg(1:ll),*) icL
      call GetArg(3, arg)
      ll = len_trim(arg)
      read(arg(1:ll),*) icT
      call GetArg(4, arg)
      ll = len_trim(arg)
      read(arg(1:ll),*) icM
      call GetArg(5, arg)
      ll = len_trim(arg)
      read(arg(1:ll),*) icMc
      call GetArg(6, arg)
      ll = len_trim(arg)
      read(arg(1:ll),*) icX
      call GetArg(7, arg)
      ll = len_trim(arg)
      read(arg(1:ll),*) icZ
      call GetArg(8, arg)
      ll = len_trim(arg)
      read(arg(1:ll),*) icCO
   else
      ! invalid arguments
      write(*,20)
      write(*,21)
      write(*,22)
      stop
   endif

   ! check existence of main input file
   lMI = len_trim(main_input)
   inquire(file = main_input(1:lMI), exist = input_exist)
   if (.not.input_exist) then
      write(*,10) main_input(1:lMI)
      stop
   endif

   ! make main output file name (path) from main input file name (path)
   main_output = main_input(1:lMI) // output_ext
   lMO = len_trim(main_output)

   ! open main input/output files for line-by-line reading/writing
   open(50, file = main_input(1:lMI), form = 'formatted', status = 'old')
   open(51, file = main_output(1:lMO), form = 'formatted', status = 'unknown')

   ! ****************** BEGIN MAIN LOOP ****************** !
   jj = 0
   do
      gret = '0000000'
      comment = .false.
      read(50,'(a1024)', iostat = ierr) line
      if (ierr < 0) then
         exit
      else if (ierr /= 0) then
         write(*,11)
         stop
      endif
      ! line parsing block - start
      ll = len_trim(line)
      l00 = 1
      j = 0
      do l = 1, ll
         if (line(l:l) == ' ') then
            l00 = l00 + 1
         else
            exit
         endif
      enddo
      l0 = l00
      do l=l00, ll - 1
         if (line(l:l) == '#') then
            comment = .true.
            exit
         endif
         if ((line(l:l) == ' ').and.(line(l+1:l+1) /= ' ')) then
            j = j + 1
            read(line(l0:l),*) tmp
            cols(j) = tmp
            l0 = l
         endif
      enddo
      if (comment) then
         do i = 1, nmodes
            write(tstr,"('  P_ad_',i1,'  GR_',i1,'  P_',i1)") i-1, i-1, i-1
            line = line(1:ll)//tstr
            ll = len_trim(line)
         enddo
         line = line(1:ll)//'  intpFLAG'
         ll = len_trim(line)
         write(51,'(a)') line(1:ll)
         cycle
      endif
      jj = jj + 1
      write(*,'(a,i10)') 'Processing line: ',jj
      j = j + 1
      read(line(l0:ll),*)  tmp
      cols(j) = tmp
      ! line parsing block - end
      ! set values to interpolate for
      lLi = cols(icL)
      lTi = cols(icT)
      Mi  = cols(icM)
      Mci = cols(icMc)
      Xi  = cols(icX)
      Zi  = cols(icZ)
      COi = cols(icCO)
      if (use_given_Z) then
         Zri = Zi
      else
         Zri = Zi / (1.D0 + CZSUN * (COi / COSUN - 1.D0))
      endif
      fCOi = LOG10(COi / COSUN) ! assuming fixed oxygen abundance !!!
      !Zri = (Z - XC - XN - XO) / (1 - CZSUN - NZSUN - OZSUN) --- this equation should be used when interpolating in XC, XN, XO
      ! select appropriate interpolation nodes in Zref, X, Mtot and fCO
      call set_interp_nodes(Zri,  gZ,  ngZ,  snZ,  ret)
      write(gret(6:6),'(i1)') ret
      call set_interp_nodes(Xi,   gX,  ngX,  snX,  ret)
      write(gret(5:5),'(i1)') ret
      call set_interp_nodes(Mi,   gM,  ngM,  snM,  ret)
      write(gret(3:3),'(i1)') ret
      call set_interp_nodes(fCOi, gCO, ngCO, snCO, ret)
      write(gret(7:7),'(i1)') ret
      do i = 1, ntot
         TT(i) = 0.D0
         do j = 1, 2
            tZ(i,j) = 0.D0
            tX(i,j) = 0.D0
            tfCO(i,j) = 0.D0
            tM(i,j) = 0.D0
         enddo
         do j = 1, 3
            tML(i,j) = 0.D0
         enddo
      enddo

      do iML = 1, 3
         do iZ = 1, 2
            do iX = 1, 2
               do iCO = 1, 2
                  do iM = 1, 2
                     write(tablefile,77) grid_dir, snZ(iZ), snX(iX), snCO(iCO), snCO(iCO), snN(1), snML(iML), snM(iM), 'lower'
                     lTF = len_trim(tablefile)
                     call interpolate_table_luminosity(lLi, tablefile(1:lTF), conversion, TT1, ret)
                     if (ret == -1) then
                        if (gret(1:1) == '0') then
                           write(gret(1:1),'(i1)') 1
                        else if (gret(1:1) == '2') then
                           write(gret(1:1),'(i1)') 3
                        endif
                     else if (ret == 1) then
                        if (gret(1:1) == '0') then
                           write(gret(1:1),'(i1)') 2
                        else if (gret(1:1) == '1') then
                           write(gret(1:1),'(i1)') 3
                        endif
                     else if (ret /= 0) then
                        stop
                     endif
                     write(tablefile,77) grid_dir, snZ(iZ), snX(iX), snCO(iCO), snCO(iCO), snN(1), snML(iML), snM(iM), 'upper'
                     lTF = len_trim(tablefile)
                     call interpolate_table_luminosity(lLi, tablefile(1:lTF), conversion, TT2, ret)
                     if (ret == -1) then
                        if (gret(1:1) == '0') then
                           write(gret(1:1),'(i1)') 1
                        else if (gret(1:1) == '2') then
                           write(gret(1:1),'(i1)') 3
                        endif
                     else if (ret == 1) then
                        if (gret(1:1) == '0') then
                           write(gret(1:1),'(i1)') 2
                        else if (gret(1:1) == '1') then
                           write(gret(1:1),'(i1)') 3
                        endif
                     else if (ret /= 0) then
                        stop
                     endif
                     if (Mci < TT1(jMc)) then
                        write(gret(4:4),'(i1)') 1
                     else if (Mci > TT2(jMc)) then
                        write(gret(4:4),'(i1)') 2
                     else
                        write(gret(4:4),'(i1)') 0
                     endif
                     dn = TT2(jMc) - TT1(jMc)
                     w1 = (TT2(jMc) - Mci) / dn
                     w2 = (Mci - TT1(jMc)) / dn
                     do i = 1, npar
                        tM(iM, i) = w1 * TT1(i) + w2 * TT2(i)
                     enddo
                     do i = 1, nmodes
                        jja = npar + npul * (i - 1) + joa
                        jjr = npar + npul * (i - 1) + jor
                        jji = npar + npul * (i - 1) + joi
                        if ((ABS(TT1(jjr) + 1.D0) <= 1.D-16).or.(ABS(TT2(jjr) + 1.D0) <= 1.D-16)) then
                           tM(iM, jja) = 1.D-10
                           tM(iM, jjr) = -1.D0
                           tM(iM, jji) = 1.D-10
                        else
                           tM(iM, jja) = w1 * TT1(jja) + w2 * TT2(jja)
                           tM(iM, jjr) = w1 * TT1(jjr) + w2 * TT2(jjr)
                           tM(iM, jji) = w1 * TT1(jji) + w2 * TT2(jji)
                        endif
                     enddo
                  enddo
                  dn = snM(2) - snM(1)
                  w1 = (snM(2) - Mi) / dn
                  w2 = (Mi - snM(1)) / dn
                  do i = 1, npar
                     tfCO(iCO, i) = w1 * tM(1, i) + w2 * tM(2, i)
                  enddo
                  do i = 1, nmodes
                     jja = npar + npul * (i - 1) + joa
                     jjr = npar + npul * (i - 1) + jor
                     jji = npar + npul * (i - 1) + joi
                     if ((ABS(tM(1, jjr) + 1.D0) <= 1.D-16).or.(ABS(tM(2, jjr) + 1.D0) <= 1.D-16)) then
                        tfCO(iCO, jja) = 1.D-10
                        tfCO(iCO, jjr) = -1.D0
                        tfCO(iCO, jji) = 1.D-10
                     else
                        tfCO(iCO, jja) = w1 * tM(1, jja) + w2 * tM(2, jja)
                        tfCO(iCO, jjr) = w1 * tM(1, jjr) + w2 * tM(2, jjr)
                        tfCO(iCO, jji) = w1 * tM(1, jji) + w2 * tM(2, jji)
                     endif
                  enddo
               enddo
               dn = snCO(2) - snCO(1)
               w1 = (snCO(2) - fCOi) / dn
               w2 = (fCOi - snCO(1)) / dn
               do i = 1, npar
                  tX(iX, i) = w1 * tfCO(1, i) + w2 * tfCO(2, i)
               enddo
               do i = 1, nmodes
                  jja = npar + npul * (i - 1) + joa
                  jjr = npar + npul * (i - 1) + jor
                  jji = npar + npul * (i - 1) + joi
                  if ((ABS(tfCO(1, jjr) + 1.D0) <= 1.D-16).or.(ABS(tfCO(2, jjr) + 1.D0) <= 1.D-16)) then
                     tX(iX, jja) = 1.D-10
                     tX(iX, jjr) = -1.D0
                     tX(iX, jji) = 1.D-10
                  else
                     tX(iX, jja) = w1 * tfCO(1, jja) + w2 * tfCO(2, jja)
                     tX(iX, jjr) = w1 * tfCO(1, jjr) + w2 * tfCO(2, jjr)
                     tX(iX, jji) = w1 * tfCO(1, jji) + w2 * tfCO(2, jji)
                  endif
               enddo
            enddo
            dn = snX(2) - snX(1)
            w1 = (snX(2) - Xi) / dn
            w2 = (Xi - snX(1)) / dn
            do i = 1, npar
               tZ(iZ, i) = w1 * tX(1, i) + w2 * tX(2, i)
            enddo
            do i = 1, nmodes
               jja = npar + npul * (i - 1) + joa
               jjr = npar + npul * (i - 1) + jor
               jji = npar + npul * (i - 1) + joi
               if ((ABS(tX(1, jjr) + 1.D0) <= 1.D-16).or.(ABS(tX(2, jjr) + 1.D0) <= 1.D-16)) then
                  tZ(iZ, jja) = 1.D-10
                  tZ(iZ, jjr) = -1.D0
                  tZ(iZ, jji) = 1.D-10
               else
                  tZ(iZ, jja) = w1 * tX(1, jja) + w2 * tX(2, jja)
                  tZ(iZ, jjr) = w1 * tX(1, jjr) + w2 * tX(2, jjr)
                  tZ(iZ, jji) = w1 * tX(1, jji) + w2 * tX(2, jji)
               endif
            enddo
         enddo
         dn = snZ(2) - snZ(1)
         w1 = (snZ(2) - Zri) / dn
         w2 = (Zri - snZ(1)) / dn
         do i = 1, npar
            tML(iML, i) = w1 * tZ(1, i) + w2 * tZ(2, i)
         enddo
         do i = 1, nmodes
            jja = npar + npul * (i - 1) + joa
            jjr = npar + npul * (i - 1) + jor
            jji = npar + npul * (i - 1) + joi
            if ((ABS(tZ(1, jjr) + 1.D0) <= 1.D-16).or.(ABS(tZ(2, jjr) + 1.D0) <= 1.D-16)) then
               tML(iML, jja) = 1.D-10
               tML(iML, jjr) = -1.D0
               tML(iML, jji) = 1.D-10
            else
               tML(iML, jja) = w1 * tZ(1, jja) + w2 * tZ(2, jja)
               tML(iML, jjr) = w1 * tZ(1, jjr) + w2 * tZ(2, jjr)
               tML(iML, jji) = w1 * tZ(1, jji) + w2 * tZ(2, jji)
            endif
         enddo
      enddo
      if (lTi < tML(1, jlogTe)) then
         write(gret(2:2),'(i1)') 1
      else if (lTi > tML(3, jlogTe)) then
         write(gret(2:2),'(i1)') 2
      else
         write(gret(2:2),'(i1)') 0
      endif
      write(*,*) gret
      if (lTi <= tML(1, jlogTe)) then
         i1 = 1
         i2 = 2
      else
         i1 = 2
         i2 = 3
      endif
      dn = tML(i2, jlogTe) - tML(i1, jlogTe)
      w1 = (tML(i2, jlogTe) - lTi) / dn
      w2 = (lTi - tML(i1, jlogTe)) / dn
      do i = 1, npar
         TT(i) = w1 * tML(i1, i) + w2 * tML(i2, i)
      enddo
      do i = 1, nmodes
         jja = npar + npul * (i - 1) + joa
         jjr = npar + npul * (i - 1) + jor
         jji = npar + npul * (i - 1) + joi
         if ((ABS(tML(i1, jjr) + 1.D0) <= 1.D-16).or.(ABS(tML(i2, jjr) + 1.D0) <= 1.D-16)) then
            TT(jja) = 1.D-10
            TT(jjr) = -1.D0
            TT(jji) = 1.D-10
         else
            TT(jja) = w1 * tML(i1, jja) + w2 * tML(i2, jja)
            TT(jjr) = w1 * tML(i1, jjr) + w2 * tML(i2, jjr)
            TT(jji) = w1 * tML(i1, jji) + w2 * tML(i2, jji)
         endif
      enddo

      ! write output
      oline(1:ll) = line(1:ll)
      do i = 1, nmodes
         jja = npar + npul * (i - 1) + joa
         jjr = npar + npul * (i - 1) + jor
         jji = npar + npul * (i - 1) + joi
         if (conversion == 2) then ! convert log(P_ad) and log(P) back to P_ad and P
            TT(jja) = 1.D1 ** TT(jja)
            TT(jji) = 1.D1 ** TT(jji)
         endif
         write(tstr, 81) TT(jja), TT(jjr), TT(jji)
         ltmp = len_trim(tstr)
         oline(ll + 1:ll + 1 + ltmp) = tstr(1:ltmp)
         ll = ll + 1 + ltmp
      enddo
      oline = oline(1:ll)//'  '//gret
      ll = len_trim(oline)
      write(51,'(a)') oline(1:ll)



   enddo
   close(50)
   close(51)



   call system_clock(t2, clock_rate, clock_max)
   write(*,12)
   write(*,*) 'Elapsed CPU time / seconds: ', real(t2 - t1) / real(clock_rate)
   write(*,*) 'Average per line / seconds: ', real(t2 - t1) / real(clock_rate) / real(jj)
   stop



!contains
!    subroutine interpolate_model()





end program giraffe
