Module arguments
   !
   Implicit None
   !
   Character(len=256),dimension(5) :: c_arg = ""
   !
   Character(len=5) :: EltIon = ""
   Real :: lambda = Huge(0E0)  ! Angstroem
   Real :: Teff   = Huge(0E0)  ! Kelvin    
   Real :: logg   = Huge(0E0)  ! log(cm/s2)
   Real :: Fe_H   = Huge(0E0)  ! logarithmic scale with respect to the solar value and where H is 12.
   !
   Integer :: lambda_ind_i, lambda_ind_j
   !
   Integer, parameter :: Nelt=3 ! Number of elements in given of ionization stage
   Character(len=5),dimension(Nelt),parameter :: Data_EltIon=(/"MgI  ","CaI  ","CaII "/)
   Real,dimension(Nelt,3),parameter :: Data_lambda=Reshape((/8473.688,8715.255,8736.012,&
                                                             8525.707,8583.318,8633.929,&
                                                             8498.018,8542.086,8662.135/),(/Nelt,3/),order=(/2,1/))
   Real,parameter :: Tol_lambda = 1.0
   Real,parameter :: Teffmin = 3500, Teffmax = 5500
   Real,parameter :: loggmin =  0.5, loggmax = 2.0
   Real,parameter :: Fe_Hmax = 0.5
   Real,dimension(Nelt,3),parameter :: Fe_Hmin = Reshape((/-1.0,-1.5,-2.0, -1.0,-1.0,-1.0, -4.0,-4.0,-4.0/),(/3,3/),order=(/2,1/)) 
   
End Module arguments
!
!***********************************************************************
!
Module coeff_fits
   !
   USE arguments, ONLY: Nelt
   !
   Implicit None
   !
   Real,dimension(Nelt,3,20), parameter :: mat_coeff=Reshape((/&
    1.03E+00,3.32E-02,-1.79E-02,-1.48E-01,1.50E-02,0.,0.,2.24E-02,4.39E-02,0.,-2.31E-02,-6.92E-03,0.,0.,4.43E-02,-1.23E-02,0., &   ! MgI  8473 
      -2.26E-02,1.90E-02,3.75E-02,                                                                                             &   ! MgI  8473 (END)
    1.22E+00,1.43E-01,-9.26E-02,-3.90E-02,-1.04E-01,6.79E-02,-2.76E-02,-1.86E-02,9.46E-02,-2.95E-01,0.,-2.49E-02,0.,0.,        &   ! MgI  8715 
       1.84E-01,-6.67E-02,2.82E-02,-7.15E-02,0.,1.72E-01,                                                                      &   ! MgI  8715 (END)
    1.31E+00,-6.84E-02,-2.36E-02,1.61E-01,-6.34E-02,0.,0.,8.05E-02,-2.16E-02,-2.42E-01,3.44E-02,-2.50E-02,0.,0.,1.20E-01,0.,   &   ! MgI  8736 
       0.,-1.18E-01,5.08E-02,-5.34E-02,                                                                                        &   ! MgI  8736 (END)
    5.95E-01,-2.81E-02,-1.67E-02,7.11E-02,4.74E-02,2.41E-02,1.22E-02,2.09E-01,5.01E-02,-2.49E-02,-9.98E-02,1.61E-02,0.,0.,0.,  &   ! CaI  8525
       0.,0.,1.04E-02,1.06E-02,6.04E-02,                                                                                       &   ! CaI  8525 (END)
    7.05E-01,0.,0.,2.76E-02,5.13E-02,2.12E-02,0.,1.84E-01,3.87E-02,-3.69E-02,-1.02E-01,0.,0.,0.,0.,0.,0.,1.14E-02,6.51E-03,    &   ! CaI  8583
       6.25E-02,                                                                                                               &   ! CaI  8583 (END)
    7.97E-01,6.00E-02,0.,0.,6.63E-02,2.02E-02,0.,1.32E-01,2.81E-02,-3.44E-02,-1.20E-01,0.,0.,0.,0.,0.,0.,1.53E-02,0.,4.85E-02, &   ! CaI  8633
    0.988,0.0368,0.00915,-0.0844,0.0422,0.,0.,-0.0387,-0.0381,0.112,-0.0254,0.0115,0.,0.,-0.0322,0.,0.,0.0417,0.0293,-0.0747,  &   ! CaII 8498 A
    0.972,0.0371,0.00305,-0.0274,0.0466,0.,0.00645,0.00533,-0.0293,0.0777,-0.0243,0.0142,-0.0108,0.,-0.0532,0.,-0.00581,0.,    &   ! CaII 8542 A
       0.0295,-0.0774,                                                                                                         &   ! CaII 8542 A (END)
    0.973,0.0382,0.00221,-0.0318,0.0460,0.,0.,-0.0167,-0.0291,0.0938,-0.0232,0.0134,-0.00903,0.,-0.0547,0.,0.,0.0207,0.0321,   &   ! CaII 8662 A
       -0.0814/),                                                                                                              &   ! CaII 8662 A (END)
 (/Nelt,3,20/),order=(/3,2,1/))
End Module coeff_fits
!
!***********************************************************************
!
Program NLTE_Gaia_lines
   !
   ! Give the NLTE/LTE Equivalent Width (EW), W/W*(line), interpolated at 
   ! Teff, logg and [Fe/H] for MgI, CaI and CaII lines in the RVS/Gaia,
   ! using the coefficient fits given in Merle et al. (2011) MNRAS
   ! The program needs 5 arguments:
   ! element with ionization stage, lambda (A), Teff (K), log g (log cgs) and [Fe/H] 
   ! If:     $ gfortran -o NLTE_Gaia_lines NLTE_Gaia_lines.f90
   ! Then:   $ ./NLTE_Gaia_lines CaII 8498 5200 1.35 -0.5
   !
   ! This means that we want to compute the NLTE/LTE EW for the CaII 8498.018 A line
   ! for the model atmosphere with Teff = 5200 K, log g = 1.35 and [Fe/H] = -0.5
   ! 
   !                             0.0       when        [Fe/H] >  0.0
   ! We remind that [alpha/Fe] = 0.4[Fe/H] when -1.0 < [Fe/H] <  0.0  
   !                             0.4       when        [Fe/H] < -1.0 
   !
   !History: 2011-07-28 (TM) Creation
   !
   USE arguments
   USE coeff_fits
   !
   Implicit None
   !
   Integer :: Narg,ios
   Real    :: WWS=Huge(0E0)
   !
   Narg = Command_argument_count()
   !
   Call CMD(Narg)
   Call VAL(EltIon,lambda,Teff,logg,Fe_H)
   Call FIT(EltIon,lambda,Teff,logg,Fe_H,WWS)
   !
   !Write(*,"(2A,F6.0,A,F5.2)")"WWS(",Trim(EltIon),lambda,") = ", WWS
   Write(*,"(F5.2)") WWS
   !
   CONTAINS
   !
   !--------------------------------------------------------------------
   !
   Real Function POLY_XYZ(x,y,z,coeff)
      Real,intent(in) :: x,y,z
      Real,dimension(:),intent(in) :: coeff
      !
      POLY_XYZ = coeff(1)         + coeff(2)*x       + coeff(3)*y       + coeff(4)*z       + &
                 coeff(5)*x**2    + coeff(6)*x*y     + coeff(7)*y**2    + coeff(8)*x*z     + &
                 coeff(9)*y*z     + coeff(10)*z**2   + coeff(11)*x**3   + coeff(12)*x**2*y + &
                 coeff(13)*x*y**2 + coeff(14)*y**3   + coeff(15)*x**2*z + coeff(16)*x*y*z  + &
                 coeff(17)*y**2*z + coeff(18)*x*z**2 + coeff(19)*y*z**2 + coeff(20)*z**3

   End Function POLY_XYZ
   !
   !--------------------------------------------------------------------
   !
   Subroutine FIT(EltIon,lambda,Teff,logg,Fe_H,WWS)
      Character(len=*),intent(in) :: EltIon
      Real,intent(in)             :: lambda,Teff,logg,Fe_H
      Real,intent(out)            :: WWS
      !
      Integer :: I
      Real :: x,y,z, xmid,ymid,zmid, xdev,ydev,zdev
      Real,dimension(20) :: coeff
      Logical :: blabla=.FALSE.
      !      
      xmid = 0.5*(Teffmax+Teffmin)
      xdev = 0.5*(Teffmax-Teffmin)
      ymid = 0.5*(loggmax+loggmin)
      ydev = 0.5*(loggmax-loggmin)
      zmid = 0.5*(Fe_Hmax+Fe_Hmin(lambda_ind_i,lambda_ind_j))
      zdev = 0.5*(Fe_Hmax-Fe_Hmin(lambda_ind_i,lambda_ind_j))
      !
      IF (blabla) Write(*,*)"xmid,xdev,ymid,ydev,zmid,zdev:",xmid,xdev,ymid,ydev,zmid,zdev
      !
      x = (Teff-xmid)/xdev
      y = (logg-ymid)/ydev
      z = (Fe_H-zmid)/zdev
      !
      IF (blabla) Write(*,*)"x,y,z:",x,y,z
      !
      coeff = mat_coeff(lambda_ind_i,lambda_ind_j,:)
      !
      IF (blabla) Write(*,'(ES10.2)') (coeff(I),I=1,20)
      !
      WWS = POLY_XYZ(x,y,z,coeff)
      !
   End Subroutine FIT
   !
   !--------------------------------------------------------------------
   !
   Subroutine CMD(Narg)
      Integer,intent(in) :: Narg
      Integer :: I
      Logical :: blabla = .FALSE.
      Character(Len=10) :: EltIon_line
      !
      IF (Narg/=4)  Then
         Call INFO()
         STOP "NLTE_GAIA_LINES PROGRAM: PB-00: MISSING ARGUMENTS."
      End If
      !
      Do I=1,Narg
         Call Get_command_argument(I,c_arg(I),status=ios)
         IF (ios/=0) STOP "NLTE_GAIA_LINES PROGRAM: PB-01: FAILED ARGUMENTS."
      End Do
      !
      IF(blabla) Write(*, *) "c_arg:", ( (Trim(c_arg(I))//" "),I=1,5 )
      ! 
      Read(c_arg(1),"(A)",iostat=ios) EltIon_line
      IF (ios/=0) STOP "NLTE_GAIA_LINES PROGRAM: PB-02: ELT WITH IONIZATION STAGE, ONLY MgI, CaI and CaII AVAILABLE."
      Select Case (EltIon_line)
         Case ("MgI_8473")
            EltIon  = "MgI  " ; lambda = 8473.0           
         Case ( "MgI_8715  " , "MgI_8710  " , "MgI_8712  " , "MgI_8717  " )   
            EltIon  = "MgI  " ; lambda = 8715.0
         Case ("MgI_8736")
            EltIon  = "MgI  " ; lambda = 8736.0
         Case ("CaI_8525")
            EltIon  = "CaI  " ; lambda = 8525.0
         Case ("CaI_8583")
            EltIon  = "CaI  " ; lambda = 8583.0
         Case ("CaI_8633")
            EltIon  = "CaI  " ; lambda = 8633.0               
         Case ("CaII_8498 ") 
            EltIon  = "CaII " ; lambda = 8498.0 
         Case ("CaII_8542 ")   
            EltIon  = "CaII " ; lambda = 8542.0 
         Case ("CaII_8662")
            EltIon  = "CaII " ; lambda = 8662.0
         Case Default 
            Write(*,*) " NLTE_GAIA_LINES PROGRAM: PB-03: WRONG ARGUMENT"
            Write(*,*) " EltIon_line AVAILABLE ARE : "
            Write(*,*) "    MgI_8473  , MgI_8715 (MgI_8710, MgI_8712, MgI_8717), MgI_8736 ,"
            Write(*,*) "    CaI_8525  , CaI_8583                               , CaI_8633 ,"
            Write(*,*) "    CaII_8498 , CaII_8542                              , CaII_8662."
            STOP
      End Select
      !
      Read(c_arg(2),"(F10.0)",iostat=ios) Teff
      IF (ios/=0) STOP "NLTE_GAIA_LINES PROGRAM: PB-04: EFFECTIVE TEMPERATURE IN KELVIN IS REQUIRED."
      Read(c_arg(3),"(F5.0)",iostat=ios) logg
      IF (ios/=0) STOP "NLTE_GAIA_LINES PROGRAM: PB-05: LOG OF SURFACE GRAVITY IN LOG(CGS) REQUIRED."
      Read(c_arg(4),"(F5.0)",iostat=ios) Fe_H
      IF (ios/=0) STOP "NLTE_GAIA_LINES PROGRAM: PB-11: GLOBAL METALLICITY WITH RESPECT TO THE SUN IS REQUIRED."
      !
      IF(blabla) Write(*,*) EltIon,lambda,Teff,logg,Fe_H  
   End Subroutine CMD
   !
   !--------------------------------------------------------------------
   !
   Subroutine VAL(EltIon,lambda,Teff,logg,Fe_H)
      Character(len=*),intent(in) :: EltIon
      Real,intent(in)             :: lambda,Teff,logg,Fe_H
      !
      Logical,dimension(5) :: larg=.FALSE.
      Logical,parameter  :: blabla=.FALSE.
      Integer :: I,J
      !
      larg(1)=Any(Data_EltIon==Trim(EltIon))
      !
      If (larg(1)) Then
         IF(blabla) Write(*,*)"Elt: ",EltIon," OK."
         !
         Do I=1,Size(Data_EltIon)
            IF (EltIon==Data_EltIon(I)) EXIT
         End Do
         !
         lambda_ind_i = I
         !
         larg(2) = Any(Abs(Data_lambda(I,:)-lambda)<=Tol_lambda)
         J=Minloc(Abs(Data_lambda(I,:)-lambda),dim=1)
         lambda_ind_j = J
         !
         If (EltIon=="MgI" .AND. (Abs(lambda-8710.169)<= Tol_lambda .OR. &
                                  Abs(lambda-8712.678)<= Tol_lambda .OR. &
                                  Abs(lambda-8717.812)<= Tol_lambda))         Then
            larg(2) = .TRUE.
            lambda_ind_j = 2
         End If      
         !
         If (.NOT.larg(2)) STOP " "
         IF (blabla) Write(*,*)"Line:    ",lambda," A, OK."  
         
         larg(3) = Teff>=Teffmin .AND. Teff<=Teffmax
         If (.NOT.larg(3)) STOP "NLTE_GAIA_LINES PROGRAM: PB-08: TEFF OUTSIDE OF THE VALID RANGE: 3500 < TEFF < 5500."
         IF (blabla) Write(*,*)"Teff   = ",Teff," K, OK." 

         larg(4) = logg>=loggmin .AND. logg<=loggmax
         If (.NOT.larg(4)) STOP "NLTE_GAIA_LINES PROGRAM: PB-09: LOGG OUTSIDE OF THE VALID RANGE: 0.5 < LOGG < 2.0."
         IF (blabla) Write(*,*)"log g  = ",logg,"   OK."
         !
         larg(5) = Fe_H>=Fe_Hmin(I,J) .AND. Fe_H<=Fe_Hmax
         If (.NOT.larg(5)) STOP "NLTE_GAIA_LINES PROGRAM: PB-10: [Fe/H] OUTSIDE OF THE VALID RANGE: [Fe/H]min < [Fe/H] < 0.5."   
         IF (blabla) Write(*,*)"[Fe/H] = ",Fe_H,"   OK."
         !
      Else
         STOP "NLTE_GAIA_LINES PROGRAM: PB-06: NO FITS FOR THIS ELEMENT."
      End If
      !
   End Subroutine VAL
   !
   !--------------------------------------------------------------------
   !
   Subroutine INFO()
      Write(*,*)""
      Write(*,*)" This  fortran  program  give  the  NLTE/LTE  Equivalent  Width  (EW),  W/W*(line),"
      Write(*,*)" interpolated at Teff, logg and [Fe/H] for MgI, CaI and CaII lines in the RVS/Gaia,"
      Write(*,*)" using the coefficient fits given in Merle et al. (2011) MNRAS"
      Write(*,*)" The program needs 4 arguments:"
      Write(*,*)" element  with  ionization stage and lambda (A), Teff (K), log g (log cgs) and [Fe/H]."
      Write(*,*)
      Write(*,*)" If:     $ gfortran -o NLTE_Gaia_lines NLTE_Gaia_lines.f90"
      Write(*,*)" Then:   $ ./NLTE_Gaia_lines CaII_8498 5200 1.35 -0.5"
      Write(*,*)
      Write(*,*)" This means that you want to compute the NLTE/LTE EW for the CaII 8498.018 A line"
      Write(*,*)" for the model atmosphere with Teff = 5200 K, log g = 1.35 and [Fe/H] = -0.5"
      Write(*,*)
      Write(*,*)" Elements and lines available:"
      Write(*,*)"   MgI:  8473, <8715> (8710, 8712, 8717), 8736"
      Write(*,*)"   CaI:  8525,          8583            , 8633"
      Write(*,*)"   CaII: 8498,          8542            , 8662"  
      Write(*,*)
      Write(*,*)" Thibault Merle -- 2011 07 28 --"
      Write(*,*) 
   End Subroutine INFO
   !
End Program NLTE_Gaia_lines
