!=======================================================================
! ch4n2o_correction_mod.F90  -  Real-world ==> effective CH4 & N2O
!-----------------------------------------------------------------------
!
!   Author: Patricia Cadule
!
!   Concept:
!   Melnikova, I., Ciais, P., Tanaka, K., Shiogama, H., Tachiiri, K., Yokohata, T., and Boucher, O.: 
!   Carbon cycle and climate feedbacks under CO2 and non-CO2 overshoot pathways, 
!   Earth Syst. Dynam., 16, 257–273, 
!   https://doi.org/10.5194/esd-16-257-2025, 2025
!
!   Calibration inputs: 
!   Olivier Boucher, Jean-Louis Dufresne
!
!   Purpose
!   -------
!   Provide *effective* methane (CH4) and nitrous oxide (N2O)
!   concentrations to the radiation scheme such that the model's
!   radiative forcing is consistent with:
!
!   * Melnikova et al. (2025, ESD), Eq. (A4), Table A2  ("m2025"), and
!   * the Etminan (2016, GRL) based optimisation by O. Boucher &
!     J.-L. Dufresne ("mbd")
!
!   Both families are implemented here in the same direction:
!
!       C_real  (input4MIP / RFMIP)   [ppb]
!         | (transfer function)
!         v
!       C_eff   (model-effective)     [ppb]
!
!   Modes
!   -----
!   ghg_mode = "m2025"
!      Melnikova et al. (2025), inverse of Eq. (A4):
!
!        CH4: C_eff = CH4_PI + [ (C_real − CH4_PI)/A_CH4_M25 ]^(1/C_CH4_M25)
!        N2O: C_eff = N2O_PI + [ (C_real − N2O_PI)/B_N2O_M25 ]^(1/D_N2O_M25)
!
!      For C_real <= PI we use the identity mapping:
!        C_eff = C_real
!
!   ghg_mode = "mbd"
!      Melnikova-Boucher–Dufresne mapping:
!
!        CH4_mod2rw:
!          C_real = CH4_PI
!                   + A_CH4_MBD * (C_eff − CH4_MIN)^P_CH4_MBD
!                   − A_CH4_MBD * (CH4_PI − CH4_MIN)^P_CH4_MBD
!
!        Inverse (used here):
!          inner = (C_real − CH4_PI
!                   + A_CH4_MBD * (CH4_PI − CH4_MIN)^P_CH4_MBD) / A_CH4_MBD
!          if inner > 0:
!            C_eff = CH4_MIN + inner^(1/P_CH4_MBD)
!          else:
!            C_eff = CH4_MIN    ! saturate at lower effective bound
!
!        N2O is analogous, with N2O_MIN, A_N2O_MBD, P_N2O_MBD.
!
!      The valid real-world domain is C_real >= C_min_valid_rw, where
!      C_min_valid_rw = C_PI - A*(C_PI - C_min)^P. Below this the
!      analytic inverse is not defined; in that regime we clip to
!      C_eff = C_min and optionally print a single debug warning.
!
!   ghg_mode = "identity"
!      Pure pass-through for diagnostics:
!        C_eff = C_real
!
!   Notes
!   -------------------
!   * For very deep sub-PI values one may enforce a floor in effective
!     concentrations via:
!         ghg_min_eff_ppb_ch4, ghg_min_eff_ppb_n2o
!
!   Public API
!   ----------
!     CALL ch4n2o_init()
!     CALL ch4n2o_compute(year, ch4_in_ppb, n2o_in_ppb, ch4_eff_ppb, n2o_eff_ppb)
!     CALL ch4_eff_compute_only(year, ch4_in_ppb, ch4_eff_ppb)
!     CALL n2o_eff_compute_only(year, n2o_in_ppb, n2o_eff_ppb)
!
!=======================================================================
MODULE ch4n2o_correction_mod
  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
  USE mod_phys_lmdz_para,     ONLY: bcast
  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
  USE ioipsl_getin_p_mod,     ONLY: getin_p
  IMPLICIT NONE
  PRIVATE

  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 300)

  !---------------- Public routines ------------------------------------
  PUBLIC :: ch4n2o_init
  PUBLIC :: ch4n2o_compute
  PUBLIC :: ch4_eff_compute_only
  PUBLIC :: n2o_eff_compute_only

  ! Optional routing of effective ppb to outputs elsewhere
  LOGICAL, PUBLIC :: ok_CH4_eff_ppb = .FALSE.
!$OMP THREADPRIVATE(ok_CH4_eff_ppb)
  LOGICAL, PUBLIC :: ok_N2O_eff_ppb = .FALSE.
!$OMP THREADPRIVATE(ok_N2O_eff_ppb)

  INTEGER, PUBLIC :: ghg_year = 1850   ! for diagnostics only
  CHARACTER(LEN=16), PUBLIC :: ghg_mode = 'm2025'  ! 'm2025' | 'mbd' | 'identity'

  !---------------- Configuration (defaults) ---------------------------

  LOGICAL           :: ghg_debug = .FALSE.

  ! Pre-industrial pivots
  REAL(dp) :: ghg_ch4_zero_ppb = 808.25_dp
  REAL(dp) :: ghg_n2o_zero_ppb = 273.02_dp

  ! Melnikova et al. (2025) parameters (Table A2)
  REAL(dp) :: ghg_A_ch4_m25 = 19.44701978_dp   ! A_CH4_M25
  REAL(dp) :: ghg_C_ch4_m25 =  0.49593024_dp   ! C_CH4_M25
  REAL(dp) :: ghg_B_n2o_m25 =  0.84856644_dp   ! B_N2O_M25
  REAL(dp) :: ghg_D_n2o_m25 =  1.13865386_dp   ! D_N2O_M25

  ! Boucher–Dufresne parameters
  REAL(dp) :: ghg_ch4_min_ppb = 808.25_dp - 100.0_dp  ! CH4_MIN
  REAL(dp) :: ghg_n2o_min_ppb = 273.02_dp -  80.0_dp  ! N2O_MIN

  REAL(dp) :: ghg_A_ch4_mbd = 77.74824074_dp  ! A_CH4_MBD
  REAL(dp) :: ghg_P_ch4_mbd =  0.36559975_dp  ! P_CH4_MBD
  REAL(dp) :: ghg_A_n2o_mbd =  0.17792515_dp  ! A_N2O_MBD
  REAL(dp) :: ghg_P_n2o_mbd =  1.38444104_dp  ! P_N2O_MBD

  ! Real-world minima for which the MBD inverse is defined
  REAL(dp) :: ghg_ch4_min_valid_rw = 0.0_dp
  REAL(dp) :: ghg_n2o_min_valid_rw = 0.0_dp

  ! Optional floors for extreme sub-PI robustness
  REAL(dp) :: ghg_min_eff_ppb_ch4 = 0.0_dp
  REAL(dp) :: ghg_min_eff_ppb_n2o = 0.0_dp

  ! One-time flags to avoid flooding logs in debug mode
  LOGICAL :: warned_ch4_mbd_below = .FALSE.
  LOGICAL :: warned_n2o_mbd_below = .FALSE.

CONTAINS
!=======================================================================
  SUBROUTINE ch4n2o_init()
    CHARACTER(LEN=32) :: mode_in
    REAL(dp)          :: rbuf(18)
    INTEGER           :: ibuf(4), mode_code

    mode_in = ghg_mode
    CALL getin_p('ghg_mode',            mode_in)
    CALL getin_p('ghg_debug',           ghg_debug)
    CALL getin_p('ghg_year',            ghg_year)

    ! Allow overriding of parameters if needed
    CALL getin_p('ghg_ch4_zero_ppb',    ghg_ch4_zero_ppb)
    CALL getin_p('ghg_n2o_zero_ppb',    ghg_n2o_zero_ppb)

    CALL getin_p('ghg_A_ch4_m25',       ghg_A_ch4_m25)
    CALL getin_p('ghg_C_ch4_m25',       ghg_C_ch4_m25)
    CALL getin_p('ghg_B_n2o_m25',       ghg_B_n2o_m25)
    CALL getin_p('ghg_D_n2o_m25',       ghg_D_n2o_m25)

    CALL getin_p('ghg_ch4_min_ppb',     ghg_ch4_min_ppb)
    CALL getin_p('ghg_n2o_min_ppb',     ghg_n2o_min_ppb)
    CALL getin_p('ghg_A_ch4_mbd',       ghg_A_ch4_mbd)
    CALL getin_p('ghg_P_ch4_mbd',       ghg_P_ch4_mbd)
    CALL getin_p('ghg_A_n2o_mbd',       ghg_A_n2o_mbd)
    CALL getin_p('ghg_P_n2o_mbd',       ghg_P_n2o_mbd)

    CALL getin_p('ghg_min_eff_ppb_ch4', ghg_min_eff_ppb_ch4)
    CALL getin_p('ghg_min_eff_ppb_n2o', ghg_min_eff_ppb_n2o)

    CALL to_lower_inplace(mode_in)
    SELECT CASE (TRIM(mode_in))
    CASE ('m2025','mbd','identity')
      ghg_mode = TRIM(mode_in)
    CASE DEFAULT
      IF (is_mpi_root .AND. is_omp_root) THEN
        WRITE(*,*) '[ch4n2o] unknown ghg_mode="',TRIM(mode_in),'" => using m2025'
      END IF
      ghg_mode = 'm2025'
    END SELECT

    ! Broadcast scalars across MPI tasks
    rbuf = (/ ghg_ch4_zero_ppb, ghg_n2o_zero_ppb, &
              ghg_A_ch4_m25, ghg_C_ch4_m25, ghg_B_n2o_m25, ghg_D_n2o_m25, &
              ghg_ch4_min_ppb, ghg_n2o_min_ppb, &
              ghg_A_ch4_mbd, ghg_P_ch4_mbd, ghg_A_n2o_mbd, ghg_P_n2o_mbd, &
              ghg_min_eff_ppb_ch4, ghg_min_eff_ppb_n2o, &
              0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp /)
    CALL bcast(rbuf)

    ghg_ch4_zero_ppb      = rbuf(1)
    ghg_n2o_zero_ppb      = rbuf(2)
    ghg_A_ch4_m25         = rbuf(3)
    ghg_C_ch4_m25         = rbuf(4)
    ghg_B_n2o_m25         = rbuf(5)
    ghg_D_n2o_m25         = rbuf(6)
    ghg_ch4_min_ppb       = rbuf(7)
    ghg_n2o_min_ppb       = rbuf(8)
    ghg_A_ch4_mbd         = rbuf(9)
    ghg_P_ch4_mbd         = rbuf(10)
    ghg_A_n2o_mbd         = rbuf(11)
    ghg_P_n2o_mbd         = rbuf(12)
    ghg_min_eff_ppb_ch4   = rbuf(13)
    ghg_min_eff_ppb_n2o   = rbuf(14)

    ibuf(1) = MERGE(1,0,ghg_debug)
    ibuf(2) = ghg_year
    ibuf(3) = 0
    ibuf(4) = 0
    CALL bcast(ibuf)
    ghg_debug = (ibuf(1) == 1)
    ghg_year  = ibuf(2)

    ! Encode mode as integer for robustness
    SELECT CASE (TRIM(ghg_mode))
    CASE ('m2025');    mode_code = 1
    CASE ('mbd');      mode_code = 2
    CASE ('identity'); mode_code = 3
    END SELECT
    CALL bcast(mode_code)
    SELECT CASE (mode_code)
    CASE (1); ghg_mode = 'm2025'
    CASE (2); ghg_mode = 'mbd'
    CASE (3); ghg_mode = 'identity'
    CASE DEFAULT
      ghg_mode = 'm2025'
    END SELECT

    ! Compute minimum real-world values for which MBD inverse is defined
    ghg_ch4_min_valid_rw = ghg_ch4_zero_ppb - &
         ghg_A_ch4_mbd * (ghg_ch4_zero_ppb - ghg_ch4_min_ppb)**ghg_P_ch4_mbd
    ghg_n2o_min_valid_rw = ghg_n2o_zero_ppb - &
         ghg_A_n2o_mbd * (ghg_n2o_zero_ppb - ghg_n2o_min_ppb)**ghg_P_n2o_mbd

    IF (is_mpi_root .AND. is_omp_root) THEN
      WRITE(*,'(A,1X,A)') '[ch4n2o] mode =',TRIM(ghg_mode)
      WRITE(*,'(A,2(F9.2,1X))') '           pivots (ppb) CH4,N2O =', &
                                 ghg_ch4_zero_ppb, ghg_n2o_zero_ppb
      SELECT CASE (TRIM(ghg_mode))
      CASE ('m2025')
        WRITE(*,'(A,2(F10.5,1X))') '           CH4 A,C =',ghg_A_ch4_m25,ghg_C_ch4_m25
        WRITE(*,'(A,2(F10.5,1X))') '           N2O B,D =',ghg_B_n2o_m25,ghg_D_n2o_m25
      CASE ('mbd')
        WRITE(*,'(A,2(F10.6,1X))') '           CH4 A,P =',ghg_A_ch4_mbd,ghg_P_ch4_mbd
        WRITE(*,'(A,2(F10.6,1X))') '           N2O A,P =',ghg_A_n2o_mbd,ghg_P_n2o_mbd
        WRITE(*,'(A,2(F9.2,1X))')  '           mins (ppb) CH4,N2O =', &
                                   ghg_ch4_min_ppb, ghg_n2o_min_ppb
        WRITE(*,'(A,2(F9.2,1X))')  '           valid C_real >= (CH4,N2O) =', &
                                   ghg_ch4_min_valid_rw, ghg_n2o_min_valid_rw
      END SELECT
      WRITE(*,'(A,I6)') '           ghg_year =', ghg_year
    END IF
  END SUBROUTINE ch4n2o_init
!-----------------------------------------------------------------------

  SUBROUTINE ch4n2o_compute(year, ch4_in_ppb, n2o_in_ppb, ch4_eff_ppb, n2o_eff_ppb)
    INTEGER,  INTENT(IN)  :: year
    REAL(dp), INTENT(IN)  :: ch4_in_ppb, n2o_in_ppb
    REAL(dp), INTENT(OUT) :: ch4_eff_ppb, n2o_eff_ppb

    CALL ch4_eff_compute_only(year, ch4_in_ppb, ch4_eff_ppb)
    CALL n2o_eff_compute_only(year, n2o_in_ppb, n2o_eff_ppb)

    IF (ghg_debug .AND. is_mpi_root .AND. is_omp_root) THEN
      WRITE(*,'(A,I6,2(A,F12.3))') 'CH4: y=',year,' real=',ch4_in_ppb,' eff=',ch4_eff_ppb
      WRITE(*,'(A,I6,2(A,F12.3))') 'N2O: y=',year,' real=',n2o_in_ppb,' eff=',n2o_eff_ppb
    END IF
  END SUBROUTINE ch4n2o_compute
!-----------------------------------------------------------------------

  SUBROUTINE ch4_eff_compute_only(year, ch4_in_ppb, ch4_eff_ppb)
    INTEGER,  INTENT(IN)  :: year
    REAL(dp), INTENT(IN)  :: ch4_in_ppb
    REAL(dp), INTENT(OUT) :: ch4_eff_ppb

    SELECT CASE (TRIM(ghg_mode))
    CASE ('identity')
      ch4_eff_ppb = ch4_in_ppb

    CASE ('m2025')
      ch4_eff_ppb = ceff_ch4_m2025(ch4_in_ppb, ghg_ch4_zero_ppb, &
                                   ghg_A_ch4_m25, ghg_C_ch4_m25)

    CASE ('mbd')
      ch4_eff_ppb = ceff_ch4_mbd(ch4_in_ppb, ghg_ch4_zero_ppb, ghg_ch4_min_ppb, &
                                 ghg_A_ch4_mbd, ghg_P_ch4_mbd)

      IF (ghg_debug .AND. .NOT. warned_ch4_mbd_below) THEN
        IF (ch4_in_ppb < ghg_ch4_min_valid_rw .AND. is_mpi_root .AND. is_omp_root) THEN
          WRITE(*,'(A,F10.3,A,F10.3,A)') &
            '[ch4n2o] CH4 MBD: C_real below valid domain (', ch4_in_ppb, &
            ' < ', ghg_ch4_min_valid_rw, ' ppb), C_eff clipped to CH4_MIN.'
          warned_ch4_mbd_below = .TRUE.
        END IF
      END IF
    END SELECT

    IF (ghg_min_eff_ppb_ch4 > 0.0_dp) THEN
      ch4_eff_ppb = MAX(ch4_eff_ppb, ghg_min_eff_ppb_ch4)
    END IF
  END SUBROUTINE ch4_eff_compute_only
!-----------------------------------------------------------------------

  SUBROUTINE n2o_eff_compute_only(year, n2o_in_ppb, n2o_eff_ppb)
    INTEGER,  INTENT(IN)  :: year
    REAL(dp), INTENT(IN)  :: n2o_in_ppb
    REAL(dp), INTENT(OUT) :: n2o_eff_ppb

    SELECT CASE (TRIM(ghg_mode))
    CASE ('identity')
      n2o_eff_ppb = n2o_in_ppb

    CASE ('m2025')
      n2o_eff_ppb = ceff_n2o_m2025(n2o_in_ppb, ghg_n2o_zero_ppb, &
                                   ghg_B_n2o_m25, ghg_D_n2o_m25)

    CASE ('mbd')
      n2o_eff_ppb = ceff_n2o_mbd(n2o_in_ppb, ghg_n2o_zero_ppb, ghg_n2o_min_ppb, &
                                 ghg_A_n2o_mbd, ghg_P_n2o_mbd)

      IF (ghg_debug .AND. .NOT. warned_n2o_mbd_below) THEN
        IF (n2o_in_ppb < ghg_n2o_min_valid_rw .AND. is_mpi_root .AND. is_omp_root) THEN
          WRITE(*,'(A,F10.3,A,F10.3,A)') &
            '[ch4n2o] N2O MBD: C_real below valid domain (', n2o_in_ppb, &
            ' < ', ghg_n2o_min_valid_rw, ' ppb), C_eff clipped to N2O_MIN.'
          warned_n2o_mbd_below = .TRUE.
        END IF
      END IF
    END SELECT

    IF (ghg_min_eff_ppb_n2o > 0.0_dp) THEN
      n2o_eff_ppb = MAX(n2o_eff_ppb, ghg_min_eff_ppb_n2o)
    END IF
  END SUBROUTINE n2o_eff_compute_only

!-----------------------------------------------------------------------
!  Maths kernels (PURE ELEMENTAL)
!-----------------------------------------------------------------------

  PURE ELEMENTAL FUNCTION ceff_ch4_m2025(c_real, c_pi, A, Cexp) RESULT(c_eff)
    REAL(dp), INTENT(IN) :: c_real, c_pi, A, Cexp
    REAL(dp)             :: c_eff, delta
#ifdef USE_OMP_SIMD
!$OMP DECLARE SIMD(ceff_ch4_m2025)
#endif
    IF (c_real <= c_pi) THEN
      c_eff = c_real
    ELSE
      delta = (c_real - c_pi) / A
      c_eff = c_pi + delta**(1.0_dp/Cexp)
    END IF
  END FUNCTION ceff_ch4_m2025
!-----------------------------------------------------------------------

  PURE ELEMENTAL FUNCTION ceff_n2o_m2025(c_real, c_pi, B, Dexp) RESULT(c_eff)
    REAL(dp), INTENT(IN) :: c_real, c_pi, B, Dexp
    REAL(dp)             :: c_eff, delta
#ifdef USE_OMP_SIMD
!$OMP DECLARE SIMD(ceff_n2o_m2025)
#endif
    IF (c_real <= c_pi) THEN
      c_eff = c_real
    ELSE
      delta = (c_real - c_pi) / B
      c_eff = c_pi + delta**(1.0_dp/Dexp)
    END IF
  END FUNCTION ceff_n2o_m2025
!-----------------------------------------------------------------------

  PURE ELEMENTAL FUNCTION ceff_ch4_mbd(c_real, c_pi, c_min, A, Pexp) RESULT(c_eff)
    ! Boucher–Dufresne CH4 inverse mapping
    REAL(dp), INTENT(IN) :: c_real, c_pi, c_min, A, Pexp
    REAL(dp)             :: c_eff, inner
#ifdef USE_OMP_SIMD
!$OMP DECLARE SIMD(ceff_ch4_mbd)
#endif
    inner = (c_real - c_pi + A*(c_pi - c_min)**Pexp) / A
    IF (inner > 0.0_dp) THEN
      c_eff = c_min + inner**(1.0_dp/Pexp)
    ELSE
      c_eff = c_min
    END IF
  END FUNCTION ceff_ch4_mbd
!-----------------------------------------------------------------------

  PURE ELEMENTAL FUNCTION ceff_n2o_mbd(c_real, c_pi, c_min, A, Pexp) RESULT(c_eff)
    ! Boucher–Dufresne N2O inverse mapping
    REAL(dp), INTENT(IN) :: c_real, c_pi, c_min, A, Pexp
    REAL(dp)             :: c_eff, inner
#ifdef USE_OMP_SIMD
!$OMP DECLARE SIMD(ceff_n2o_mbd)
#endif
    inner = (c_real - c_pi + A*(c_pi - c_min)**Pexp) / A
    IF (inner > 0.0_dp) THEN
      c_eff = c_min + inner**(1.0_dp/Pexp)
    ELSE
      c_eff = c_min
    END IF
  END FUNCTION ceff_n2o_mbd

!-----------------------------------------------------------------------
! Utility: in-place lower-casing for mode strings
!-----------------------------------------------------------------------
  SUBROUTINE to_lower_inplace(s)
    CHARACTER(*), INTENT(INOUT) :: s
    INTEGER :: i, ia
    DO i = 1, LEN_TRIM(s)
      ia = IACHAR(s(i:i))
      IF (ia >= IACHAR('A') .AND. ia <= IACHAR('Z')) s(i:i) = ACHAR(ia+32)
    END DO
  END SUBROUTINE to_lower_inplace

END MODULE ch4n2o_correction_mod