MODULE MMP_GCM
  !! Interface to YAMMS for the LMDZ GCM.
  USE MM_LIB
  USE CFGPARSE
  USE DATASETS
  IMPLICIT NONE

  PUBLIC

  !> Alpha function parameters. 
  !!
  !! It stores the parameters of the inter-moments relation functions.
  !! 
  !! The inter-moments relation function is represented by the sum of exponential
  !! quadratic expressions:
  !!
  !! $$ 
  !! \displaystyle \alpha(k) = \sum_{i=1}^{n} \exp\left( a_{i}\times k^{2} + 
  !! b_{i}\times k^{2} +c_{i}\right) 
  !! $$
  TYPE, PUBLIC :: aprm
    !> Quadratic coefficients of the quadratic expressions.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: a
    !> Linear coefficients of the quadratic expressions.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: b
    !> Free term of the quadratic expressions.
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: c
  END TYPE
  !> Inter-moment relation set of parameters for the spherical mode.
  TYPE(aprm), PUBLIC, SAVE :: mmp_asp
  !> Inter-moment relation set of parameters for the fractal mode.
  TYPE(aprm), PUBLIC, SAVE :: mmp_afp

  !> Data set for @f$<Q>_{SF}^{M0}@f$.
  TYPE(dset2d), PUBLIC, SAVE, TARGET             :: mmp_qbsf0
  !> Extended values of [[mmp_globals(module):mmp_qbsf0(variable)]] dataset.
  REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mmp_qbsf0_e
  !> Data set for @f$<Q>_{SF}^{M3}@f$.
  TYPE(dset2d), PUBLIC, SAVE, TARGET             :: mmp_qbsf3
  !> Extended values of [[mmp_globals(module):mmp_qbsf3(variable)]] dataset.
  REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mmp_qbsf3_e
  !> Data set for @f$<Q>_{FF}^{M0}@f$.
  TYPE(dset2d), PUBLIC, SAVE, TARGET             :: mmp_qbff0
  !> Extended values of [[mmp_globals(module):mmp_qbff0(variable)]] dataset.
  REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mmp_qbff0_e
  
  !> Data set for linear interpolation of transfert probability (M0/CO).
  TYPE(dset1d), PUBLIC, SAVE, TARGET :: mmp_pco0p
  !> Data set for linear interpolation of transfert probability (M3/CO).
  TYPE(dset1d), PUBLIC, SAVE, TARGET :: mmp_pco3p
  !> Data set for linear interpolation of transfert probability (M0/FM).
  TYPE(dset1d), PUBLIC, SAVE, TARGET :: mmp_pfm0p
  !> Data set for linear interpolation of transfert probability (M3/FM). 
  TYPE(dset1d), PUBLIC, SAVE, TARGET :: mmp_pfm3p

  !> \(b_{0}^{t}\) coefficients for Free-molecular regime kernel approximation.
  REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(5) :: mmp_bt0 = (/1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp/)
  !> \(b_{3}^{t}\) coefficients for Free-molecular regime kernel approximation.
  REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(5) :: mmp_bt3 = (/1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp/)

  !> Spherical probability transfert control flag.
  LOGICAL, SAVE :: mmp_w_ps2s = .true.
  !> Aerosol electric charge correction control flag.
  LOGICAL, SAVE :: mmp_w_qe   = .true.


  CONTAINS 
    
  SUBROUTINE abort_program(err)
    !! Dump error message and abort the program.
    TYPE(error), INTENT(in) :: err !! Error object.
    WRITE(stderr,'(a)') "ERROR: "//TRIM(err%msg) 
    CALL EXIT(err%id)
  END SUBROUTINE abort_program


  SUBROUTINE mmp_initialize(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, air_rad,air_mmol,clouds,cfgpath)
    !! Initialize global parameters of the model.
    !! 
    !! The function initializes all the global parameters of the model from direct input.
    !! Boolean (and Fiadero) parameters are optional as they are rather testing parameters. Their 
    !! default values are suitable for production runs.  
    !! @note
    !! If the method fails to initialize parameters (i.e. returned error is not 0). Then the model
    !! should probably be aborted as the global variables of the model will not be correctly setup.
    !!
    !! @warning
    !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it 
    !! initializes global variable that are not thread private.
    !!
    !! '''
    !! !$OMP SINGLE
    !! call mmp_initialize(...)
    !! !$OMP END SINGLE
    !!
    REAL(kind=mm_wp), INTENT(in)           :: dt
      !! Microphysics timestep in seconds.
    REAL(kind=mm_wp), INTENT(in)           :: df
      !! Fractal dimension of fractal aerosol.
    REAL(kind=mm_wp), INTENT(in)           :: rm
      !! Monomer radius in meter.
    REAL(kind=mm_wp), INTENT(in)           :: rho_aer
      !! Aerosol density in \(kg.m^{-3}\).
    REAL(kind=mm_wp), INTENT(in)           :: p_prod
      !!  Aerosol production pressure level in Pa.
    REAL(kind=mm_wp), INTENT(in)           :: tx_prod
      !! Spherical aerosol mode production rate in \(kg.m^{-2}.s^{-1}\).
    REAL(kind=mm_wp), INTENT(in)           :: rc_prod
      !! Spherical mode characteristic radius for production in meter.
    REAL(kind=mm_wp), INTENT(in)           :: rplanet
      !! Planet radius in meter
    REAL(kind=mm_wp), INTENT(in)           :: g0
      !! Planet gravity acceleration at ground level in \(m.s^{-2}\).
    REAL(kind=mm_wp), INTENT(in)           :: air_rad
      !! Air molecules mean radius in meter.
    REAL(kind=mm_wp), INTENT(in)           :: air_mmol
      !! Air molecules mean molar mass in \(kg.mol^{-1}\).
    LOGICAL, INTENT(in)                    :: clouds
      !! Clouds microphysics control flag.
    CHARACTER(len=*), INTENT(in), OPTIONAL :: cfgpath
      !! Internal microphysic configuration file.

    INTEGER          :: coag_choice
    REAL(kind=mm_wp) :: fiad_max, fiad_min
    LOGICAL          :: w_h_prod, w_h_sed, w_h_coag, w_c_sed, w_c_nucond, &
                        no_fiadero, fwsed_m0, fwsed_m3
    TYPE(error)      :: err
    INTEGER                                           :: i
    TYPE(cfgparser)                                   :: cparser
    CHARACTER(len=st_slen)                            :: spcpath,pssfile,mqfile
    CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE       :: tmp

    w_h_prod    = .true.
    w_h_sed     = .true.
    w_h_coag    = .true.
    w_c_sed     = clouds
    w_c_nucond  = clouds
    fwsed_m0    = .true.
    fwsed_m3    = .false.
    no_fiadero  = .false.
    fiad_min    = 0.1_mm_wp
    fiad_max    = 10._mm_wp
    coag_choice = 7

    WRITE(*,'(a)') "##### MMP_GCM SPEAKING #####"
    WRITE(*,'(a)') "I will initialize ze microphysics model in moments YAMMS"
    WRITE(*,'(a)') "On error I will simply abort the program. Stay near your computer !"
    WRITE(*,*)
    WRITE(*,'(a)') "Reading muphys configuration file ("//trim(cfgpath)//")..."
    err = cfg_read_config(cparser,TRIM(cfgpath),.true.) 
    IF (err /= 0) THEN
      ! RETURN AN ERROR !!
      call abort_program(err)
    ENDIF
   
    ! YAMMS internal parameters:
    ! the following parameters are primarily used to test and debug YAMMS.
    ! They are set in an optional configuration file and default to suitable values for production runs.
    err = mm_check_opt(cfg_get_value(cparser,"haze_production",w_h_prod)    ,w_h_prod   ,.true.   ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"haze_sedimentation",w_h_sed)  ,w_h_sed    ,.true.   ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"haze_coagulation",w_h_coag)   ,w_h_coag   ,.true.   ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"clouds_sedimentation",w_c_sed),w_c_sed    ,clouds   ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"clouds_nucl_cond",w_c_nucond) ,w_c_nucond ,clouds   ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"wsed_m0",fwsed_m0)            ,fwsed_m0   ,.true.   ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"wsed_m3",fwsed_m3)            ,fwsed_m3   ,.false.  ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"no_fiadero",no_fiadero)       ,no_fiadero ,.false.  ,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"fiadero_min_ratio",fiad_min)  ,fiad_min   ,0.1_mm_wp,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"fiadero_max_ratio",fiad_max)  ,fiad_max   ,10._mm_wp,mm_log)
    err = mm_check_opt(cfg_get_value(cparser,"haze_coag_interactions",coag_choice),coag_choice,7,mm_log)

    ! Retrieve clouds species configuration file
    spcpath = ''
    IF (clouds) THEN
      err = mm_check_opt(cfg_get_value(cparser,"specie_cfg",spcpath), spcpath, wlog=mm_log)
      IF (err/=0) call abort_program(err)
    ENDIF

    ! YAMMS initialization.
    err = mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, &
                           air_rad,air_mmol,coag_choice,clouds,spcpath,  &
                           w_h_prod,w_h_sed,w_h_coag,w_c_nucond,  &
                           w_c_sed,fwsed_m0,fwsed_m3, &
                           no_fiadero,fiad_min,fiad_max)
    IF (err /= 0) call abort_program(err)

    ! Extra initialization (needed for YAMMS method interfaces)
    err = mm_check_opt(cfg_get_value(cparser, "transfert_probability", mmp_w_ps2s), mmp_w_ps2s, wlog=mm_log)
    IF (err/=0) call abort_program(err)
    err = mm_check_opt(cfg_get_value(cparser, "electric_charging"    , mmp_w_qe  ), mmp_w_qe, wlog=mm_log)
    IF (err/=0) call abort_program(err)

    ! initialize transfert probabilities look-up tables
    IF (mm_w_haze_coag .AND. mmp_w_ps2s) THEN
      err = mm_check_opt(cfg_get_value(cparser, "ps2s_file", pssfile), pssfile)
      IF (err /= 0) call abort_program(err)

      IF (.NOT.read_dset(pssfile,'p_m0_co',mmp_pco0p)) THEN
        call abort_program(error("Cannot get 'p_m0_co' from "//pssfile,-1)) 
      ENDIF
      IF (.NOT.read_dset(pssfile,'p_m3_co',mmp_pco3p)) THEN
        call abort_program(error("Cannot get 'p_m3_co' from "//pssfile,-1)) 
      ENDIF
      IF (.NOT.read_dset(pssfile,'p_m0_fm',mmp_pfm0p)) THEN
        call abort_program(error("Cannot get 'p_m0_fm' from "//pssfile,-1))
      ENDIF
      IF (.NOT.read_dset(pssfile,'p_m3_fm',mmp_pfm3p)) THEN
        call abort_program(error("Cannot get 'p_m3_fm' from "//pssfile,-1))
      ENDIF
    ENDIF
    ! initialize mean electric correction look-up tables
    IF (mm_w_haze_coag .AND. mmp_w_qe) THEN
      err = mm_check_opt(cfg_get_value(cparser, "mq_file", mqfile), mqfile)
      IF (err /= 0) call abort_program(err)

      IF (.NOT.read_dset(mqfile,'qbsf0',mmp_qbsf0)) THEN
        call abort_program(error("Cannot get 'qbsf0' from "//mqfile,-1))
      ELSE
        mmp_qbsf0_e(1,1) = MINVAL(mmp_qbsf0%x) 
        mmp_qbsf0_e(1,2) = MAXVAL(mmp_qbsf0%x)
        mmp_qbsf0_e(2,1) = MINVAL(mmp_qbsf0%y)
        mmp_qbsf0_e(2,2) = MAXVAL(mmp_qbsf0%y)
      ENDIF
      IF (.NOT.read_dset(mqfile,'qbsf3',mmp_qbsf3)) THEN
        call abort_program(error("Cannot get 'qbsf3' from "//mqfile,-1))
      ELSE
        mmp_qbsf3_e(1,1) = MINVAL(mmp_qbsf3%x) 
        mmp_qbsf3_e(1,2) = MAXVAL(mmp_qbsf3%x)
        mmp_qbsf3_e(2,1) = MINVAL(mmp_qbsf3%y)
        mmp_qbsf3_e(2,2) = MAXVAL(mmp_qbsf3%y)
      ENDIF
      IF (.NOT.read_dset(mqfile,'qbff0',mmp_qbff0)) THEN
        call abort_program(error("Cannot get 'qbff0' from "//mqfile,-1))
      ELSE
        mmp_qbff0_e(1,1) = MINVAL(mmp_qbff0%x) 
        mmp_qbff0_e(1,2) = MAXVAL(mmp_qbff0%x)
        mmp_qbff0_e(2,1) = MINVAL(mmp_qbff0%y)
        mmp_qbff0_e(2,2) = MAXVAL(mmp_qbff0%y)
      ENDIF
    ENDIF
    ! spherical mode inter-moments function parameters
    IF (.NOT.cfg_has_section(cparser,'alpha_s')) call abort_program(error("Cannot find [alpha_s] section",-1))
    err = read_aprm(cparser,'alpha_s',mmp_asp)
    IF (err /= 0) call abort_program(error("alpha_s: "//TRIM(err%msg),-1))
    ! fractal mode inter-moments function parameters
    IF (.NOT.cfg_has_section(cparser,'alpha_f')) call abort_program(error("Cannot find [alpha_f] section",-1))
    err = read_aprm(cparser,'alpha_f',mmp_afp)
    IF (err /= 0) call abort_program(error("alpha_s: "//TRIM(err%msg),-1))
    ! btk coefficients
    IF (.NOT.cfg_has_section(cparser,'btks')) call abort_program(error("Cannot find [btks] section",-1))
    err = cfg_get_value(cparser,"btks/bt0",tmp) ; IF (err/=0) call abort_program(error("bt0: "//TRIM(err%msg),-1))
    IF (SIZE(tmp) /= 5)  call abort_program(error("bt0: Inconsistent number of coefficients",-1))
    mmp_bt0 = tmp
    err = cfg_get_value(cparser,"btks/bt3",tmp) ; IF (err/=0) call abort_program(error("bt3: "//TRIM(err%msg),-1))
    IF (SIZE(tmp) /= 5)  call abort_program(error("bt3: Inconsistent number of coefficients",-1))
    mmp_bt3 = tmp

    ! dump parameters ...
    WRITE(*,'(a)')        "========= MUPHYS PARAMETERS ==========="
    WRITE(*,'(a,L2)')     "transfert_probability: ", mmp_w_ps2s
    WRITE(*,'(a,L2)')     "electric_charging    : ", mmp_w_qe
    call mm_dump_parameters()
      
  END SUBROUTINE mmp_initialize

  FUNCTION read_aprm(parser,sec,pp) RESULT(err)
    !! Read and store [[mmp_globals(module):aprm(type)]] parameters. 
    TYPE(cfgparser), INTENT(in)  :: parser !! Configuration parser 
    CHARACTER(len=*), INTENT(in) :: sec    !! Name of the section that contains the parameters.
    TYPE(aprm), INTENT(out)      :: pp     !! [[mmp_globals(module):aprm(type)]] object that stores the parameters values.
    TYPE(error) :: err                     !! Error status of the function.
    err = cfg_get_value(parser,TRIM(sec)//'/a',pp%a) ; IF (err /= 0) RETURN
    err = cfg_get_value(parser,TRIM(sec)//'/b',pp%b) ; IF (err /= 0) RETURN
    err = cfg_get_value(parser,TRIM(sec)//'/c',pp%c) ; IF (err /= 0) RETURN
    IF (SIZE(pp%a) /= SIZE(pp%b) .OR. SIZE(pp%a) /= SIZE(pp%c)) &
    err = error("Inconsistent number of coefficients (a,b, and c must have the same size)",-1)
    RETURN
  END FUNCTION read_aprm

END MODULE MMP_GCM