MODULE tracco2i_mod
!
! This module does the work for the interactive CO2 tracers
! Authors: Patricia Cadule and Olivier Boucher
!
! Purpose and description:
!  -----------------------
! Main routine for the interactive carbon cycle
! Gather all carbon fluxes and emissions from ORCHIDEE, PISCES and fossil fuel
! Compute the net flux in source field which is used in phytrac
! Compute global CO2 mixing ratio for radiation scheme if option is activated
! Redistribute CO2 evenly over the atmosphere if transport is desactivated
!
! Purpose:
!   Interactive CO2 tracer module for LMDZ atmospheric model
!   Supports emission-driven AND concentration-driven modes
!   Supports daily emission files, AND monthly with SZ98 interpolation
!   Supports 2D (lonxlat) correction fields for land and ocean
!
! Authors:
!   Patricia Cadule (IPSL) and Olivier Boucher (IPSL/LMD)
!
! History:
!   2020: Initial version (P. Cadule and O. Boucher)
!   2026-01: Added concentration-driven mode + SZ98 interpolation (P. Cadule)
!   2026-02: Added 2D (lonxlat) correction fields for land and ocean (P. Cadule)
!
!
! Operating Modes:
!   1. EMISSION-DRIVEN (carbon_cycle_conc_driven = .FALSE.)
!      - CO2 evolves from surface fluxes and emissions
!      - source(:,id_CO2) contains the net flux
!      - Requires CO2 in tracer.def (co2_tracer_available = .TRUE.)
!
!   2. CONCENTRATION-DRIVEN (carbon_cycle_conc_driven = .TRUE.)
!      - Atmospheric CO2 is prescribed from co2_ppm
!      - source(:,id_CO2) = 0 (no flux feedback)
!      - co2_send transmits prescribed CO2 to land/ocean models
!      - Can operate WITHOUT CO2 in tracer.def (co2_tracer_available = .FALSE.)
!
!
! IMPORTANT: Transport Control
!   carbon_cycle_tr does NOT control tracer transport in LMDZ dynamics
!   Transport is controlled by the presence of CO2 in traceur.def and tracer.def
!
!   For true concentration-driven mode (no transport CPU cost):
!     - Remove CO2 from traceur.def and tracer.def 
!     - co2_tracer_available will be .FALSE.
!     - co2_send will use co2_ppm directly
!
! Emission Reading Modes:
!   (A) read_daily_co2ff = .TRUE.  ==> Files have time = year_len (360/365/366)
!   (B) read_daily_co2ff = .FALSE. ==> Files have time = 12, use SZ98 interpolation
!
!
!
!===============================================================================
!
  IMPLICIT NONE

  PRIVATE

  PUBLIC :: tracco2i_init, tracco2i

  !-----------------------------------------------------------------------------
  ! Module-level variables
  !-----------------------------------------------------------------------------

  INTEGER, SAVE :: day_pre = -1
  !$OMP THREADPRIVATE(day_pre)

  LOGICAL, SAVE :: first_call = .TRUE.
  !$OMP THREADPRIVATE(first_call)

  LOGICAL, SAVE :: check_fCO2_nbp_in_cfname = .FALSE.
  !$OMP THREADPRIVATE(check_fCO2_nbp_in_cfname)


  !---------------------------------------------------------------------------
  ! CO2 tracer availability flag
  !---------------------------------------------------------------------------
  ! This flag indicates whether the CO2 tracer exists in the tracer arrays.
  ! When carbon_cycle_tr = .FALSE. or CO2 is not in tracer.def, the tracer
  ! arrays may not include CO2. Accessing tr_seri(:,:,id_CO2) with invalid
  ! id_CO2 would cause segmentation faults.
  !
  ! The flag is set ONCE at first call based on:
  !   id_CO2 >= 1 .AND. id_CO2 <= nbtr
  !---------------------------------------------------------------------------
  LOGICAL, SAVE :: co2_tracer_available = .FALSE.
  !$OMP THREADPRIVATE(co2_tracer_available)


CONTAINS


!===============================================================================
! SUBROUTINE: tracco2i_init
!
! Purpose:
!   Initialise the carbon cycle module. Called from phytrac_init.
!   Must be called BEFORE first call to phys_output_write.
!===============================================================================
SUBROUTINE tracco2i_init()

  USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl, &
                              carbon_cycle_tr, carbon_cycle_rad,   &
                              carbon_cycle_conc_driven

  IMPLICIT NONE

  ! Call carbon_cycle_init if any carbon cycle pathway is active
  ! This allocates all necessary arrays
  IF (carbon_cycle_cpl .OR. carbon_cycle_tr .OR. &
      carbon_cycle_rad .OR. carbon_cycle_conc_driven) THEN
     CALL carbon_cycle_init()
  END IF

END SUBROUTINE tracco2i_init


!===============================================================================
! SUBROUTINE: tracco2i
!
! Purpose:
!   Main routine for interactive carbon cycle processing
!   Handles both emission-driven and concentration-driven modes
!
! Description:
!   1. Gathers surface CO2 fluxes from Land (ORCHIDEE), Ocean (PISCES),
!      Fossil Fuel (and Biomass Burning).
!   2. Builds the tracer source term [kg CO2 m-2 s-1].
!   3. Computes global mean CO2 mixing ratio and atmospheric mass.
!   4. Handles scenario-based target CO2 (diagnostic or emission-comp).
!
! Authors:
!   2026: Patricia Cadule
!   2020: Patricia Cadule, Olivier Boucher (initial)
! 
!
! Called from:
!   phytrac 
!
!===============================================================================
SUBROUTINE tracco2i(pdtphys, debutphy, &
                    xlat, xlon, pphis, pphi, &
                    t_seri, pplay, paprs, tr_seri, source)

  USE dimphy
  USE infotrac_phy, ONLY: nbtr
  USE geometry_mod, ONLY: cell_area

  ! Explicit import of needed variables from carbon_cycle_mod
  USE carbon_cycle_mod, ONLY: id_CO2
  USE carbon_cycle_mod, ONLY:  nbcf_in, fields_in, cfname_in, cfunits_in
  USE carbon_cycle_mod, ONLY : nbcf_out, fields_out, yfields_out, cfname_out
  USE carbon_cycle_mod, ONLY: fco2_ocn_day, fco2_ff, fco2_bb
  USE carbon_cycle_mod, ONLY: fco2_land, fco2_ocean
  USE carbon_cycle_mod, ONLY: read_fco2_ocean_cor, var_fco2_ocean_cor, fco2_ocean_cor
  USE carbon_cycle_mod, ONLY: read_fco2_land_cor, var_fco2_land_cor, fco2_land_cor
  USE carbon_cycle_mod, ONLY: read_fco2_ocean_cor_2d, read_fco2_land_cor_2d
  USE carbon_cycle_mod, ONLY: fco2_ocean_cor_2d_base, fco2_land_cor_2d_base
  USE carbon_cycle_mod, ONLY: co2_send
  USE carbon_cycle_mod, ONLY: fco2_land_nbp, fco2_land_nep, fco2_land_fLuc
  USE carbon_cycle_mod, ONLY: fco2_land_fwoodharvest, fco2_land_fHarvest
  USE carbon_cycle_mod, ONLY: carbon_cycle_cpl, carbon_cycle_tr
  USE carbon_cycle_mod, ONLY: carbon_cycle_rad, RCO2_glo, RCO2_tot
  
  USE carbon_cycle_mod, ONLY: carbon_cycle_conc_driven
  USE carbon_cycle_mod, ONLY: ocean_area_tot_glo, land_area_tot_glo
  USE carbon_cycle_mod, ONLY: ocean_area_tot, land_area_tot
  

  USE mod_grid_phy_lmdz
  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
  USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter, reduce_sum
  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root

  USE phys_cal_mod

  USE phys_state_var_mod, ONLY: pctsrf
  USE indice_sol_mod, ONLY: nbsrf, is_ter, is_lic, is_oce, is_sic

  USE yomcst_mod_h
  USE clesphys_mod_h
  USE print_control_mod, ONLY: lunout

  IMPLICIT NONE

  !-----------------------------------------------------------------------------
  ! Arguments
  !-----------------------------------------------------------------------------
  REAL, INTENT(IN)    :: pdtphys              ! Physics time step [s]
  LOGICAL, INTENT(IN) :: debutphy             ! TRUE at first physics call

  REAL, DIMENSION(klon), INTENT(IN)          :: xlat, xlon, pphis
  REAL, DIMENSION(klon,klev), INTENT(IN)     :: pphi, t_seri, pplay
  REAL, DIMENSION(klon,klev+1), INTENT(IN)   :: paprs

  REAL, DIMENSION(klon,nbtr), INTENT(INOUT)       :: source   ! Surface flux [U/m2/s]
  REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri  ! Tracer conc [U/kgA]

  !-----------------------------------------------------------------------------
  ! Local variables
  !-----------------------------------------------------------------------------
  CHARACTER(LEN=20), PARAMETER  :: modname = 'tracco2i'
  INTEGER :: it, k, i, nb
  REAL, DIMENSION(klon,klev) :: m_air     ! Mass of air in each grid box [kg]

  ! Global arrays for gather
  REAL, DIMENSION(klon_glo,klev)  :: co2_glo, m_air_glo
  REAL, DIMENSION(klon_glo,nbsrf) :: pctsrf_glo
  REAL, DIMENSION(klon_glo)       :: pctsrf_ter_glo, pctsrf_oce_glo, pctsrf_sic_glo
  REAL, DIMENSION(klon_glo)       :: cell_area_glo

  ! Local sums for MPI reduction
  REAL :: co2_mass_local, air_mass_local
  REAL :: co2_mass_global, air_mass_global

  REAL, PARAMETER :: secinday = 86400.0
  

  !===========================================================================
  ! 0. SAFETY CHECK: Validate tracer index bounds
  !
  ! When CO2 is not in tracer.def, id_CO2 may be invalid.
  ! We MUST check before any access to tr_seri or source with id_CO2.
  !===========================================================================

  IF (first_call) THEN
     ! Check if id_CO2 is a valid tracer index
     IF (id_CO2 >= 1 .AND. id_CO2 <= nbtr) THEN
        co2_tracer_available = .TRUE.
     ELSE
        co2_tracer_available = .FALSE.
     END IF

     IF (is_omp_root .AND. is_mpi_root) THEN
        WRITE(lunout,*) '=============================================='
        WRITE(lunout,*) modname, ': SAFETY CHECK'
        WRITE(lunout,*) '  nbtr                     = ', nbtr
        WRITE(lunout,*) '  id_CO2                   = ', id_CO2
        WRITE(lunout,*) '  carbon_cycle_tr          = ', carbon_cycle_tr
        WRITE(lunout,*) '  carbon_cycle_conc_driven = ', carbon_cycle_conc_driven
        WRITE(lunout,*) '  co2_tracer_available     = ', co2_tracer_available
        IF (.NOT. co2_tracer_available) THEN
           WRITE(lunout,*) '  WARNING: CO2 tracer NOT in tracer list!'
           WRITE(lunout,*) '           Tracer operations will be SKIPPED.'
           WRITE(lunout,*) '           co2_send will use co2_ppm directly.'
        END IF
        WRITE(lunout,*) '=============================================='
     END IF

     first_call = .FALSE.
  END IF

  !===========================================================================
  ! 1. INITIALISATION AT FIRST PHYSICS STEP
  !===========================================================================
  IF (debutphy) THEN

     ! Initialise tracer field ONLY if CO2 tracer is available
     IF (co2_tracer_available) THEN
        IF (carbon_cycle_conc_driven) THEN
           ! Concentration-driven: force tracer to prescribed value
           tr_seri(:,:,id_CO2) = co2_ppm * 1.0e-6 * (RMCO2 / RMD)
        ELSE
           ! Emission-driven: initialise only if essentially zero
           IF (MAXVAL(ABS(tr_seri(:,:,id_CO2))) < 1.0e-15) THEN
              tr_seri(:,:,id_CO2) = co2_ppm0 * 1.0e-6 * (RMCO2 / RMD)
           END IF
        END IF
     END IF

     ! Check if fCO2_nbp is in coupling fields
     check_fCO2_nbp_in_cfname = .FALSE.
     DO nb = 1, nbcf_in
        IF (cfname_in(nb) == "fCO2_nbp") check_fCO2_nbp_in_cfname = .TRUE.
     END DO

     ! Gather surface fractions for area calculations
     CALL gather(pctsrf, pctsrf_glo)
     CALL gather(pctsrf(:,is_ter), pctsrf_ter_glo)
     CALL gather(pctsrf(:,is_oce), pctsrf_oce_glo)
     CALL gather(pctsrf(:,is_sic), pctsrf_sic_glo)
     CALL gather(cell_area(:), cell_area_glo)

     ! Compute total areas for flux corrections
     ! IMPORTANT: areas are computed whenever ANY correction mode is
     ! active (scalar OR 2D), not only in scalar mode. This ensures:
     !   (a) Mixed modes (e.g. scalar ocean + 2D land) work correctly
     !   (b) Diagnostic output is meaningful in all modes
     !   (c) Global integral verification of 2D fields is possible

     IF (read_fco2_ocean_cor .OR. read_fco2_ocean_cor_2d) THEN
!$OMP MASTER
        IF (is_mpi_root .AND. is_omp_root) THEN
           ocean_area_tot = 0.0
           DO i = 1, klon_glo
              ocean_area_tot = ocean_area_tot + &
                   (pctsrf_oce_glo(i) + pctsrf_sic_glo(i)) * cell_area_glo(i)
           END DO
        END IF
!$OMP END MASTER
        CALL bcast(ocean_area_tot)
     END IF

     IF (read_fco2_land_cor .OR. read_fco2_land_cor_2d) THEN
!$OMP MASTER
        IF (is_mpi_root .AND. is_omp_root) THEN
           land_area_tot = 0.0
           DO i = 1, klon_glo
              land_area_tot = land_area_tot + pctsrf_ter_glo(i) * cell_area_glo(i)
           END DO
        END IF
!$OMP END MASTER
        CALL bcast(land_area_tot)
     END IF

     ! Diagnostic output: show correction mode and areas
     IF (is_omp_root .AND. is_mpi_root) THEN
        WRITE(lunout,*) modname, ': --- Flux correction configuration ---'
        WRITE(lunout,*) modname, ':   read_fco2_ocean_cor    = ', read_fco2_ocean_cor
        WRITE(lunout,*) modname, ':   read_fco2_ocean_cor_2d = ', read_fco2_ocean_cor_2d
        WRITE(lunout,*) modname, ':   read_fco2_land_cor     = ', read_fco2_land_cor
        WRITE(lunout,*) modname, ':   read_fco2_land_cor_2d  = ', read_fco2_land_cor_2d
        WRITE(lunout,*) modname, ':   ocean_area_tot =', ocean_area_tot, ' m2'
        WRITE(lunout,*) modname, ':   land_area_tot  =', land_area_tot, ' m2'
        IF (read_fco2_ocean_cor) THEN
           WRITE(lunout,*) modname, ':   var_fco2_ocean_cor =', &
                var_fco2_ocean_cor, ' PgC/yr (scalar mode)'
        END IF
        IF (read_fco2_land_cor) THEN
           WRITE(lunout,*) modname, ':   var_fco2_land_cor  =', &
                var_fco2_land_cor, ' PgC/yr (scalar mode)'
        END IF
     END IF

  END IF ! debutphy

  !===========================================================================
  ! 2. COMPUTE MASS OF AIR
  !===========================================================================
  DO k = 1, klev
     m_air(:,k) = (paprs(:,k) - paprs(:,k+1)) / RG * cell_area(:)
  END DO

  !===========================================================================
  ! 3. READ/UPDATE EMISSIONS
  !===========================================================================
  CALL co2_emissions(debutphy)

  !===========================================================================
  ! 3b. READ 2D FLUX CORRECTIONS (once at debutphy)
  !     Populates fco2_ocean_cor / fco2_land_cor 
  !     Only active when read_fco2_ocean_cor_2d or read_fco2_land_cor_2d
  !     is set to .TRUE. in physiq.def.
  !===========================================================================
  CALL co2_flux_corrections(debutphy)
  
  !===========================================================================
  ! 4. PROCESS COUPLING FIELDS
  !===========================================================================
  fco2_land(:)              = 0.0
  fco2_ocean(:)             = 0.0
  fco2_land_nbp(:)          = 0.0
  fco2_land_nep(:)          = 0.0
  fco2_land_fLuc(:)         = 0.0
  fco2_land_fwoodharvest(:) = 0.0
  fco2_land_fHarvest(:)     = 0.0

  ! Map coupling fields from ORCHIDEE/PISCES
  DO nb = 1, nbcf_in
     SELECT CASE (cfname_in(nb))
     CASE ("fCO2_nep")
        fco2_land_nep(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter)
     CASE ("fCO2_fLuc")
        fco2_land_fLuc(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter)
     CASE ("fCO2_fwoodharvest")
        fco2_land_fwoodharvest(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter)
     CASE ("fCO2_fHarvest")
        fco2_land_fHarvest(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter)
     CASE ("fCO2_nbp")
        fco2_land_nbp(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter)
     CASE ("fCO2_fgco2")
        fco2_ocean(:) = -1.0 * fco2_ocn_day(:) * RMCO2/1.0e3 * &
                        (pctsrf(:,is_oce) + pctsrf(:,is_sic))
     END SELECT
  END DO

  !===========================================================================
  ! 5. APPLY FLUX CORRECTIONS
  !
  !    The 2D base rate field is multiplied by the CURRENT
  !    pctsrf at EVERY time step, so the correction tracks sea-ice
  !    evolution dynamically. This makes the 2D path mathematically
  !    identical to the scalar path (both use runtime pctsrf).
  !
  !    Use cases:
  !      - 2D alone:    base_rate * pctsrf (from piControl correction file)
  !      - Scalar alone: scalar_rate * pctsrf (from physiq.def parameter)
  !      - Both:         additive combination (fine-tuning)
  !
  !    SIGN CONVENTION (source term, Section 7):
  !      source = fco2_ff + fco2_bb + fco2_land + fco2_ocean
  !             - fco2_ocean_cor - fco2_land_cor
  !
  !    Positive correction REDUCES atmospheric CO2 (net sink).
  !    Negative correction INCREASES atmospheric CO2 (net source).
  !===========================================================================

  ! --- Ocean correction ---
  ! Always start from zero, then add active components
  fco2_ocean_cor(:) = 0.0

  ! (a) 2D base rate * current pctsrf
  IF (read_fco2_ocean_cor_2d) THEN
     fco2_ocean_cor(:) = fco2_ocean_cor_2d_base(:) * &
                         (pctsrf(:,is_oce) + pctsrf(:,is_sic))
  END IF

  ! (b) Scalar uniform correction (additive if 2D also active)
  IF (read_fco2_ocean_cor .AND. ocean_area_tot > 0.0) THEN
     fco2_ocean_cor(:) = fco2_ocean_cor(:) + &
                         var_fco2_ocean_cor * 1.0e12 * RMCO2/RMC / &
                         ocean_area_tot / secinday / REAL(year_len) * &
                         (pctsrf(:,is_oce) + pctsrf(:,is_sic))
  END IF

  ! --- Land correction ---
  fco2_land_cor(:) = 0.0

  ! (a) 2D base rate * current pctsrf
  IF (read_fco2_land_cor_2d) THEN
     fco2_land_cor(:) = fco2_land_cor_2d_base(:) * pctsrf(:,is_ter)
  END IF

  ! (b) Scalar uniform correction (additive if 2D also active)
  IF (read_fco2_land_cor .AND. land_area_tot > 0.0) THEN
     fco2_land_cor(:) = fco2_land_cor(:) + &
                        var_fco2_land_cor * 1.0e12 * RMCO2/RMC / &
                        land_area_tot / secinday / REAL(year_len) * &
                        pctsrf(:,is_ter)
  END IF


  !===========================================================================
  ! 6. AGGREGATE LAND FLUXES
  !===========================================================================
  IF (check_fCO2_nbp_in_cfname) THEN
     fco2_land(:) = fco2_land_nbp(:)
  ELSE
     fco2_land(:) = fco2_land_nep(:) + fco2_land_fLuc(:) + &
                    fco2_land_fwoodharvest(:) + fco2_land_fHarvest(:)
  END IF

  !===========================================================================
  ! 7. BUILD FINAL SOURCE TERM
  !===========================================================================

  IF (co2_tracer_available) THEN
     ! Tracer exists - can write to source(:,id_CO2)
     IF (carbon_cycle_conc_driven) THEN
        !---------------------------------------------------------------------
        ! CONCENTRATION-DRIVEN MODE
        ! Source term is zero (no flux feedback to atmosphere)
        ! Fluxes computed above are for diagnostic outputs only
        !---------------------------------------------------------------------
        source(:,id_CO2) = 0.0
     ELSE
        !---------------------------------------------------------------------
        ! EMISSION-DRIVEN MODE
        ! Sum all flux components 
        !---------------------------------------------------------------------
        source(:,id_CO2) = fco2_ff(:) + fco2_bb(:) + fco2_land(:) + fco2_ocean(:) &
                         - fco2_ocean_cor(:) - fco2_land_cor(:)
     END IF
  END IF
  ! If co2_tracer_available = .FALSE., we do NOT touch source(:,id_CO2)

  !===========================================================================
  ! 8. COMPUTE GLOBAL MEAN CO2 (daily update)
  !===========================================================================
  IF (debutphy .OR. day_cur /= day_pre) THEN

     IF (carbon_cycle_conc_driven) THEN
        !---------------------------------------------------------------------
        ! Concentration-driven: prescribe from co2_ppm
        !---------------------------------------------------------------------
        RCO2_glo = co2_ppm * 1.0e-6 * (RMCO2/RMD)
        RCO2_glo = REAL(INT(RCO2_glo * 1.0e8)) / 1.0e8

        ! Force tracer to prescribed value (ONLY if available)
        IF (co2_tracer_available) THEN
           tr_seri(:,:,id_CO2) = RCO2_glo
        END IF

        IF (is_mpi_root .AND. is_omp_root) THEN
           WRITE(lunout,*) modname, ': Concentration-driven CO2 =', co2_ppm, ' ppm'
        END IF

     ELSE IF (co2_tracer_available) THEN
        !---------------------------------------------------------------------
        ! Emission-driven with tracer available: compute from tracer field
        !---------------------------------------------------------------------
        co2_mass_local = SUM(tr_seri(:,:,id_CO2) * m_air(:,:))
        air_mass_local = SUM(m_air(:,:))

        co2_mass_global = 0.0
        air_mass_global = 0.0
        CALL reduce_sum(co2_mass_local, co2_mass_global)
        CALL reduce_sum(air_mass_local, air_mass_global)

!$OMP BARRIER
        IF (air_mass_global > 0.0) THEN
           RCO2_glo = co2_mass_global / air_mass_global
        ELSE
           RCO2_glo = co2_ppm0 * 1.0e-6 * (RMCO2/RMD)
        END IF
        RCO2_glo = REAL(INT(RCO2_glo * 1.0e8)) / 1.0e8
!$OMP BARRIER

        IF (is_mpi_root .AND. is_omp_root) THEN
           WRITE(lunout,*) modname, ': global CO2 =', RCO2_glo*1.0e6*RMD/RMCO2, ' ppm'
        END IF

        ! If no transport flag, homogenise tracer field
        IF (.NOT. carbon_cycle_tr) THEN
           tr_seri(:,:,id_CO2) = RCO2_glo
        END IF

     ELSE
        !---------------------------------------------------------------------
        ! Tracer not available: use prescribed value
        !---------------------------------------------------------------------
        RCO2_glo = co2_ppm * 1.0e-6 * (RMCO2/RMD)
     END IF

     day_pre = day_cur
  END IF

  !===========================================================================
  ! 9. UPDATE CO2 TO SEND TO SURFACE MODELS
  !
  ! In concentration-driven mode OR when tracer not available:
  !   ==> send co2_ppm (prescribed)
  ! In emission-driven mode with tracer:
  !   ==> send computed atmospheric value
  !===========================================================================
  IF (ALLOCATED(co2_send)) THEN
     IF (carbon_cycle_conc_driven .OR. .NOT. co2_tracer_available) THEN
        co2_send(:) = co2_ppm
     ELSE
        co2_send(:) = tr_seri(:,1,id_CO2) * 1.0e6 * RMD / RMCO2
     END IF
  END IF

  IF (is_mpi_root .AND. is_omp_root) THEN
     IF (co2_tracer_available) THEN
        WRITE(lunout,*) modname, ': tr_seri L1 min/max (ppm) =', &
             MINVAL(tr_seri(:,1,id_CO2)*1.0e6*RMD/RMCO2), &
             MAXVAL(tr_seri(:,1,id_CO2)*1.0e6*RMD/RMCO2)
     END IF
     IF (ALLOCATED(co2_send)) THEN
        WRITE(lunout,*) modname, ': co2_send min/max =', &
             MINVAL(co2_send), MAXVAL(co2_send)
     END IF
  END IF

END SUBROUTINE tracco2i

!===============================================================================
! FUNCTION: to_lower
!===============================================================================
PURE FUNCTION to_lower(s) RESULT(t)
  IMPLICIT NONE
  CHARACTER(*), INTENT(IN) :: s
  CHARACTER(LEN(s))        :: t
  INTEGER :: i, ia

  DO i = 1, LEN(s)
     ia = IACHAR(s(i:i))
     IF (ia >= IACHAR('A') .AND. ia <= IACHAR('Z')) THEN
        t(i:i) = ACHAR(ia + 32)
     ELSE
        t(i:i) = s(i:i)
     END IF
  END DO
END FUNCTION to_lower


!===============================================================================
! SUBROUTINE: ocean_flux_to_atm_source
!===============================================================================
SUBROUTINE ocean_flux_to_atm_source(flux_in, units_in, frac_ocean, flux_out)

  USE yomcst_mod_h, ONLY: RMCO2, RMC
  USE clesphys_mod_h

  IMPLICIT NONE

  REAL, INTENT(IN)          :: flux_in(:)
  CHARACTER(*), INTENT(IN)  :: units_in
  REAL, INTENT(IN)          :: frac_ocean(:)
  REAL, INTENT(OUT)         :: flux_out(:)

  REAL :: factor
  CHARACTER(LEN=64) :: u
  REAL, PARAMETER :: secinday = 86400.0
  CHARACTER(LEN=20), PARAMETER :: modname = 'ocean_flux_conv'

  u = to_lower(ADJUSTL(TRIM(units_in)))
  factor = -1.0

  SELECT CASE (u)
  CASE ('gc m-2 day-1', 'gc/m2/day', 'gc m-2 d-1', 'gc.m-2.day-1', 'gc.m-2.d-1')
     flux_out(:) = factor * flux_in(:) * 1.0e-3 / secinday * (RMCO2/RMC) * frac_ocean(:)
  CASE ('gc m-2 s-1', 'gc/m2/s', 'gc.m-2.s-1')
     flux_out(:) = factor * flux_in(:) * 1.0e-3 * (RMCO2/RMC) * frac_ocean(:)
  CASE ('kgc m-2 s-1', 'kgc/m2/s', 'kgc.m-2.s-1')
     flux_out(:) = factor * flux_in(:) * (RMCO2/RMC) * frac_ocean(:)
  CASE ('molc.m-2.s-1', 'molc/m2/s', 'molc m-2 s-1')
     flux_out(:) = factor * flux_in(:) * (RMCO2 / 1000.0) * frac_ocean(:)
  CASE DEFAULT
     CALL abort_physic(modname, 'Unsupported ocean flux units: '//TRIM(units_in), 1)
  END SELECT

END SUBROUTINE ocean_flux_to_atm_source

!===============================================================================
! "SZ98" ROUTINES 
!
! These routines implement the Sheng & Zwiers (1998) interpolation method
! for converting monthly mean emissions to smooth daily values whilst
! exactly preserving the monthly means.
!
! Reference:
!   Sheng, J., & Zwiers, F. (1998). An improved scheme for time-dependent
!   boundary conditions in atmospheric general circulation models.
!   Climate Dynamics, 14, 609-613.
!===============================================================================

!===============================================================================
! SUBROUTINE: co2_emissions
!
! Purpose:
!   Read CO2 fossil-fuel (and optionally biomass-burning) emissions from
!   NetCDF files and provide daily flux values to the atmospheric model.
!
! Authors:
!   2020: Olivier Boucher (IPSL/LMD) - initial
!   2026: Patricia Cadule (IPSL) - daily mode + SZ98 interpolation
!
!
! Operating Modes:
!   (A) read_daily_co2ff = .TRUE.
!       Input file has time = year_len (360/365/366) daily values.
!       Direct indexing by day-of-year.
!
!   (B) read_daily_co2ff = .FALSE.
!       Input file has time = 12 monthly mean values.
!       SZ98 interpolation provides smooth daily values that preserve
!       the monthly means exactly.
!!
! SZ98 Interpolation:
!   Sheng & Zwiers (1998) method ensures monthly means are preserved
!   when interpolating to daily values. The tridiagonal system is solved
!   ONCE per year.
!   Sheng, J. & Zwiers, F.W. (1998). An improved scheme for time-dependent
!   boundary conditions in atmospheric general circulation models.
!   Climate Dynamics, 14, 609-613.
!   DOI: 10.1007/s003820050244

!
! Units:
!   flx_co2ff/bb : kg CO2 m-2 s-1
!   fco2_ff/bb   : kg CO2 m-2 s-1 (output to source term)
!
!===============================================================================
SUBROUTINE co2_emissions(debutphy)

  USE dimphy
  USE geometry_mod, ONLY: cell_area
  USE mod_grid_phy_lmdz
  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root, omp_rank
  USE mod_phys_lmdz_para, ONLY: scatter
  USE print_control_mod, ONLY: lunout
  USE phys_cal_mod, ONLY: year_len, day_cur, mth_cur

  USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_varid, nf95_open
  USE netcdf,   ONLY: nf90_get_var, nf90_nowrite, nf90_noerr, nf90_strerror, &
                    nf90_inq_dimid, nf90_inquire_dimension

  ! Import variables from module
  USE carbon_cycle_mod, ONLY: fco2_ff, fco2_bb
  USE carbon_cycle_mod, ONLY: read_daily_co2ff


  USE clesphys_mod_h
  USE yomcst_mod_h

  IMPLICIT NONE


  !-----------------------------------------------------------------------------
  ! Arguments
  !-----------------------------------------------------------------------------
  LOGICAL, INTENT(IN) :: debutphy   ! .TRUE. at first physics call

  !-----------------------------------------------------------------------------
  !
  ! In monthly mode: (klon, 12) - 12 monthly means
  ! In daily mode:   (klon, year_len) - one value per day
  !
  ! These are allocated once at debutphy and persist across time steps.
  !-----------------------------------------------------------------------------
  REAL, ALLOCATABLE, SAVE :: flx_co2ff(:,:)     ! Fossil-fuel CO2 [kgCO2/m2/s]
  REAL, ALLOCATABLE, SAVE :: flx_co2bb(:,:)     ! Biomass-burning CO2 [kgCO2/m2/s]
!$OMP THREADPRIVATE(flx_co2ff, flx_co2bb)

  !-- SZ98 interpolation targets (monthly mode only, THREADPRIVATE)
  REAL, ALLOCATABLE, SAVE :: target_ff(:,:)     ! (klon, 12) mid-month targets
  REAL, ALLOCATABLE, SAVE :: target_bb(:,:)     ! (klon, 12) mid-month targets
!$OMP THREADPRIVATE(target_ff, target_bb)

  !-- Calendar information for SZ98 interpolation
  INTEGER, SAVE :: days_in_mth(12)              ! Days per month for current year
  REAL,    SAVE :: month_midday(12)             ! Day-of-year of month mid-points
!$OMP THREADPRIVATE(days_in_mth, month_midday)

  !-- Daily update guard: avoid recalculating every time step
  INTEGER, SAVE :: day_pre_emis = -1
!$OMP THREADPRIVATE(day_pre_emis)

  !-----------------------------------------------------------------------------
  ! Local variables for NetCDF reading
  !-----------------------------------------------------------------------------
  REAL, ALLOCATABLE :: flx_co2ff_glo(:,:)       ! Global field (MPI root only)
  REAL, ALLOCATABLE :: flx_co2bb_glo(:,:)       ! Global field (MPI root only)

  INTEGER :: ncid_in, varid, ncerr
  INTEGER :: n_glo, n_time_file, m, n_times
  INTEGER :: current_year_len, doy

  REAL, ALLOCATABLE :: vector_coord(:), time_coord(:)

  !-- Control flags 
  LOGICAL, PARAMETER :: readco2ff = .TRUE.
  LOGICAL, PARAMETER :: readco2bb = .FALSE.

  CHARACTER(LEN=30), PARAMETER :: modname = 'co2_emissions'
  CHARACTER(LEN=80) :: abort_message

  !=============================================================================
  ! 1. INITIALISATION 
  !=============================================================================
  IF (debutphy) THEN

     current_year_len = INT(year_len)

     !-- Determine number of time steps in the file
     IF (read_daily_co2ff) THEN
        n_times = current_year_len      ! 360, 365, or 366
     ELSE
        n_times = 12                     ! Monthly
     END IF

     !------------------------------------------------------------------
     ! A. Allocate arrays
     !------------------------------------------------------------------
     ALLOCATE(flx_co2ff(klon, n_times))
     ALLOCATE(flx_co2bb(klon, n_times))
     flx_co2ff(:,:) = 0.0
     flx_co2bb(:,:) = 0.0

     IF (.NOT. read_daily_co2ff) THEN
        ALLOCATE(target_ff(klon, 12))
        ALLOCATE(target_bb(klon, 12))
        target_ff(:,:) = 0.0
        target_bb(:,:) = 0.0
     END IF

     !------------------------------------------------------------------
     ! B. NetCDF reading
     !------------------------------------------------------------------
!$OMP MASTER
     IF (is_mpi_root) THEN

        !-- Allocate global buffer on MPI root
        IF (.NOT. ALLOCATED(flx_co2ff_glo)) ALLOCATE(flx_co2ff_glo(klon_glo, n_times))
        IF (.NOT. ALLOCATED(flx_co2bb_glo)) ALLOCATE(flx_co2bb_glo(klon_glo, n_times))
        flx_co2ff_glo(:,:) = 0.0
        flx_co2bb_glo(:,:) = 0.0

        !---------------------------------------------------------------
        ! B.1 Read fossil-fuel emissions
        !---------------------------------------------------------------
        IF (readco2ff) THEN

           CALL nf95_open('sflx_lmdz_co2_ff.nc', nf90_nowrite, ncid_in)

           !-- Validate spatial dimension
           CALL nf95_inq_varid(ncid_in, 'vector', varid)
           CALL nf95_gw_var(ncid_in, varid, vector_coord)
           n_glo = SIZE(vector_coord)
           DEALLOCATE(vector_coord)
           IF (n_glo /= klon_glo) THEN
              abort_message = 'sflx_lmdz_co2_ff: vector /= klon_glo'
              CALL abort_physic(modname, abort_message, 1)
           END IF

           !-- Validate time dimension
           CALL nf95_inq_varid(ncid_in, 'time', varid)
           CALL nf95_gw_var(ncid_in, varid, time_coord)
           n_time_file = SIZE(time_coord)
           DEALLOCATE(time_coord)
           IF (n_time_file /= n_times) THEN
              WRITE(abort_message, '(A,I4,A,I4)') &
                   'sflx_lmdz_co2_ff: time=', n_time_file, ' exp=', n_times
              CALL abort_physic(modname, abort_message, 1)
           END IF

           !-- Read emission field directly into (klon_glo, n_times) array
           CALL nf95_inq_varid(ncid_in, 'flx_co2', varid)
           ncerr = nf90_get_var(ncid_in, varid, flx_co2ff_glo)

           CALL nf95_close(ncid_in)

           !-- Reorder from file convention to LMDZ internal ordering
           !
           !   Longitude convention:
           !     Determined by LMDZ grid definition file used in preprocessing.
           !     If grid uses [0,360]: shift_lon = 0
           !     If grid uses [-180,+180]: shift_lon = nbp_lon / 2 (default)
           !
           !   If a different grid definition is used (e.g. xfirst=0), this value
           !   must be adjusted accordingly: shift_lon = 0 for xfirst=0.
           !   CSHIFT by nbp_lon/2 = 72 positions rotates each band by 180 degrees,
           !   aligning the Python [0,360] convention with the LMDZ [-180,+180] grid.
           !---------------------------------------------------------------
           DO m = 1, n_times
              CALL reorder_global_1d(flx_co2ff_glo(:,m), nbp_lon, nbp_lat, &
                                     flip_lat=.FALSE., shift_lon=nbp_lon/2)
           END DO

           WRITE(lunout,*) modname, ':   Reorder checked: flip_lat=F, shift_lon=', nbp_lon/2
           WRITE(*,*) modname, ': Read FF emissions OK, n_times=', n_times
           IF (read_daily_co2ff) THEN
              WRITE(*,*) '  Mode: DAILY (direct indexing by day-of-year)'
           ELSE
              WRITE(*,*) '  Mode: MONTHLY (SZ98 interpolation to daily)'
           END IF

        ELSE
           flx_co2ff_glo(:,:) = 0.0
        END IF

        !---------------------------------------------------------------
        ! B.2 Read biomass-burning emissions (disabled by default)
        !---------------------------------------------------------------
        IF (readco2bb) THEN

           CALL nf95_open('sflx_lmdz_co2_bb.nc', nf90_nowrite, ncid_in)

           CALL nf95_inq_varid(ncid_in, 'vector', varid)
           CALL nf95_gw_var(ncid_in, varid, vector_coord)
           n_glo = SIZE(vector_coord)
           DEALLOCATE(vector_coord)
           IF (n_glo /= klon_glo) THEN
              abort_message = 'sflx_lmdz_co2_bb: vector /= klon_glo'
              CALL abort_physic(modname, abort_message, 1)
           END IF

           CALL nf95_inq_varid(ncid_in, 'time', varid)
           CALL nf95_gw_var(ncid_in, varid, time_coord)
           n_time_file = SIZE(time_coord)
           DEALLOCATE(time_coord)
           IF (n_time_file /= n_times) THEN
              abort_message = 'sflx_lmdz_co2_bb: time mismatch'
              CALL abort_physic(modname, abort_message, 1)
           END IF

           CALL nf95_inq_varid(ncid_in, 'flx_co2', varid)
           ncerr = nf90_get_var(ncid_in, varid, flx_co2bb_glo)

           CALL nf95_close(ncid_in)
           
           DO m = 1, n_times
              CALL reorder_global_1d(flx_co2ff_glo(:,m), nbp_lon, nbp_lat, &
                                     flip_lat=.FALSE., shift_lon=nbp_lon/2)
           END DO
        ELSE
           flx_co2bb_glo(:,:) = 0.0
        END IF

     END IF ! is_mpi_root
!$OMP END MASTER

     !------------------------------------------------------------------
     ! C. Distribute global field 
     !------------------------------------------------------------------

     !-- Dummy allocation for non-MPI-root ranks
     IF (.NOT. ALLOCATED(flx_co2ff_glo)) ALLOCATE(flx_co2ff_glo(0, 0))
     IF (.NOT. ALLOCATED(flx_co2bb_glo)) ALLOCATE(flx_co2bb_glo(0, 0))

     !-- Scatter: distributes first dimension (klon_glo ==> klon)
     CALL scatter(flx_co2ff_glo, flx_co2ff)
     CALL scatter(flx_co2bb_glo, flx_co2bb)

     !-- Release global buffers 
     IF (ALLOCATED(flx_co2ff_glo)) DEALLOCATE(flx_co2ff_glo)
     IF (ALLOCATED(flx_co2bb_glo)) DEALLOCATE(flx_co2bb_glo)

     !------------------------------------------------------------------
     ! D. SZ98 TARGETS: compute mid-month interpolation targets
     !    (monthly mode only, each thread processes its own data)
     !
     !    The SZ98 tridiagonal solve finds target values T_j such that
     !    piecewise-linear interpolation between mid-month points
     !    preserves the original monthly means M_j exactly:
     !      M_j = a_j * T_{j-1} + b_j * T_j + c_j * T_{j+1}
     !
     !    This is done ONCE per year. Daily interpolation is O(1).
     !------------------------------------------------------------------
     IF (.NOT. read_daily_co2ff) THEN

        CALL get_days_in_month(current_year_len, days_in_mth)
        CALL compute_month_middays(days_in_mth, month_midday)

        CALL sz98_solve_tridiagonal(flx_co2ff, days_in_mth, 12, klon, target_ff)
        CALL sz98_solve_tridiagonal(flx_co2bb, days_in_mth, 12, klon, target_bb)

        IF (is_mpi_root) THEN
           WRITE(*,*) modname, ': SZ98 targets computed, year_len=', current_year_len
           WRITE(*,*) '  days_in_month =', days_in_mth
           WRITE(*,*) '  month_midday  =', month_midday
        END IF

     END IF

  END IF ! debutphy

  !=============================================================================
  ! 2. DAILY UPDATE: set fco2_ff and fco2_bb for the current day
  !
  !    Guard: only recalculate when the day changes (day_cur /= day_pre_emis).
  !    Within a single day, the emissions remain constant regardless of how
  !    many physics time steps occur. This saves unnecessary computation.
  !
  !    IMPORTANT: day_cur from phys_cal_mod is the day of the MONTH (1-31),
  !    NOT the day of year. For daily-file indexing and SZ98 interpolation,
  !    we need the day-of-year, computed from mth_cur and day_cur.
  !=============================================================================
  IF (day_cur /= day_pre_emis) THEN

     !-- Compute day-of-year from current month and day
     doy = compute_day_of_year(mth_cur, day_cur, INT(year_len))

     IF (read_daily_co2ff) THEN
        !---------------------------------------------------------------
        ! (A) DAILY MODE: direct indexing by day-of-year
        !---------------------------------------------------------------
        IF (doy >= 1 .AND. doy <= SIZE(flx_co2ff, 2)) THEN
           fco2_ff(:) = flx_co2ff(:, doy)
           fco2_bb(:) = flx_co2bb(:, doy)
        ELSE
           WRITE(*,*) modname, ': WARNING: doy out of range =', doy
           fco2_ff(:) = 0.0
           fco2_bb(:) = 0.0
        END IF

     ELSE
        !---------------------------------------------------------------
        ! (B) MONTHLY MODE: SZ98 piecewise-linear interpolation
        !     Interpolates between mid-month target values to give a
        !     smooth daily field that preserves monthly means.
        !---------------------------------------------------------------
        CALL sz98_interpolate_daily(target_ff, month_midday, doy, &
                                    klon, 12, fco2_ff)
        CALL sz98_interpolate_daily(target_bb, month_midday, doy, &
                                    klon, 12, fco2_bb)

     END IF

     day_pre_emis = day_cur

  END IF

END SUBROUTINE co2_emissions


!===============================================================================
! FUNCTION: compute_day_of_year
!
! Purpose:
!   Convert (month, day_of_month) to day-of-year (1 to year_len).
!   Uses the current calendar's month lengths.
!
! Arguments:
!   month       : Current month (1-12)
!   day         : Day of the month (1-31)
!   year_length : Total number of days in the year (360, 365, or 366)
!
! Returns:
!   Day-of-year (integer, 1 to year_length)
!===============================================================================
INTEGER FUNCTION compute_day_of_year(month, day, year_length)

  IMPLICIT NONE

  INTEGER, INTENT(IN) :: month, day, year_length
  INTEGER :: dim(12), im

  CALL get_days_in_month(year_length, dim)

  compute_day_of_year = day
  DO im = 1, month - 1
     compute_day_of_year = compute_day_of_year + dim(im)
  END DO

END FUNCTION compute_day_of_year


!===============================================================================
! SUBROUTINE: get_days_in_month
!
! Purpose:
!   Return the number of days in each month for a given year length.
!   Supports 360-day, 365-day (noleap), and 366-day (leap) calendars.
!
! Arguments:
!   year_length   : Total days in the year (360, 365, or 366)
!   days_in_month : Output array of 12 month lengths
!===============================================================================
SUBROUTINE get_days_in_month(year_length, days_in_month)

  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: year_length
  INTEGER, INTENT(OUT) :: days_in_month(12)

  SELECT CASE (year_length)

  CASE (360)
     !-- 360-day calendar: uniform 30-day months
     days_in_month(:) = 30

  CASE (365)
     !-- Standard non-leap year (noleap/365_day calendar)
     days_in_month = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)

  CASE (366)
     !-- Leap year (gregorian/standard calendar)
     days_in_month = (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)

  CASE DEFAULT
     !-- Fallback: assume non-leap year
     days_in_month = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)

  END SELECT

END SUBROUTINE get_days_in_month


!===============================================================================
! SUBROUTINE: compute_month_middays
!
! Purpose:
!   Compute the day-of-year corresponding to the mid-point of each month.
!   Used as interpolation nodes for the SZ98 piecewise-linear scheme.
!
! Example (non-leap year):
!   January:  mid-point = 0 + 31/2 = 15.5 (day 15-16 boundary)
!   February: mid-point = 31 + 28/2 = 45.0 (day 45)
!   ...
!
! Arguments:
!   days_in_month : Number of days in each month (12)
!   month_midday  : Output: day-of-year of each month mid-point (12)
!===============================================================================
SUBROUTINE compute_month_middays(days_in_month, month_midday)

  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: days_in_month(12)
  REAL,    INTENT(OUT)  :: month_midday(12)

  INTEGER :: cumul(13), im

  !-- Cumulative days at the START of each month
  cumul(1) = 0
  DO im = 1, 12
     cumul(im+1) = cumul(im) + days_in_month(im)
  END DO

  !-- Mid-point = start + half the month length
  DO im = 1, 12
     month_midday(im) = REAL(cumul(im)) + REAL(days_in_month(im)) / 2.0
  END DO

END SUBROUTINE compute_month_middays


!===============================================================================
! SUBROUTINE: sz98_solve_tridiagonal
!
! Purpose:
!   Solve the Sheng & Zwiers (1998) tridiagonal system to find mid-month
!   target values T that preserve monthly means M under piecewise-linear
!   interpolation:
!
!     M_j = a_j * T_{j-1} + b_j * T_j + c_j * T_{j+1}
!
! Mathematical Background:
!   The coefficients depend on month lengths d_j:
!     Interior months (j = 2, ..., 11):
!       a_j = d_j / (4 * (d_{j-1} + d_j))
!       c_j = d_j / (4 * (d_j + d_{j+1}))
!       b_j = 1 - a_j - c_j
!
!     Boundary months:
!       j=1:  b_1 = (3*d_1+4*d_2)/(4*(d_1+d_2)),  c_1 = d_1/(4*(d_1+d_2))
!       j=12: a_N = d_N/(4*(d_{N-1}+d_N)),  b_N = (4*d_{N-1}+3*d_N)/(4*(d_{N-1}+d_N))
!
!   When all d_j = 30 (360-day calendar): a=c=1/8, b=3/4 (interior)
!
! Algorithm:
!   Thomas algorithm (tridiagonal matrix solver) 
!
! Arguments:
!   monthly_means : Input monthly mean values (klon_loc, 12)
!   days_in_month : Number of days in each month (12)
!   n_months      : Number of months (always 12)
!   klon_loc      : Number of local grid points
!   targets       : Output mid-month target values (klon_loc, 12)
!
! Reference:
!   Sheng, J. & Zwiers, F.W. (1998), Climate Dynamics, 14, 609-613.
!===============================================================================
SUBROUTINE sz98_solve_tridiagonal(monthly_means, days_in_month, &
                                  n_months, klon_loc, targets)

  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: n_months
  INTEGER, INTENT(IN)  :: klon_loc
  REAL,    INTENT(IN)  :: monthly_means(klon_loc, n_months)
  INTEGER, INTENT(IN)  :: days_in_month(n_months)
  REAL,    INTENT(OUT) :: targets(klon_loc, n_months)

  !-- Tridiagonal coefficients
  REAL :: a(n_months)         ! Sub-diagonal
  REAL :: b(n_months)         ! Main diagonal
  REAL :: c(n_months)         ! Super-diagonal
  REAL :: d(n_months)         ! Month lengths as real

  !-- Thomas algorithm working arrays
  REAL :: cp(n_months)        ! Modified super-diagonal
  REAL :: dp(klon_loc, n_months) 

  INTEGER :: j, i
  REAL    :: denom

  !---------------------------------------------------------------------------
  ! Convert month lengths to real for coefficient computation
  !---------------------------------------------------------------------------
  DO j = 1, n_months
     d(j) = REAL(days_in_month(j))
  END DO

  !---------------------------------------------------------------------------
  ! Build tridiagonal matrix (SZ98 variable-length-month formulation)
  !---------------------------------------------------------------------------

  !-- First row (boundary condition: no month before January)
  a(1) = 0.0
  b(1) = (3.0*d(1) + 4.0*d(2)) / (4.0*(d(1) + d(2)))
  c(1) = d(1) / (4.0*(d(1) + d(2)))

  !-- Interior rows (February through November)
  DO j = 2, n_months - 1
     a(j) = d(j) / (4.0*(d(j-1) + d(j)))
     c(j) = d(j) / (4.0*(d(j)   + d(j+1)))
     b(j) = 1.0 - a(j) - c(j)
  END DO

  !-- Last row (boundary condition: no month after December)
  a(n_months) = d(n_months) / (4.0*(d(n_months-1) + d(n_months)))
  b(n_months) = (4.0*d(n_months-1) + 3.0*d(n_months)) / &
                (4.0*(d(n_months-1) + d(n_months)))
  c(n_months) = 0.0

  !---------------------------------------------------------------------------
  ! Thomas Algorithm: Forward elimination
  !---------------------------------------------------------------------------

  !-- First row
  cp(1) = c(1) / b(1)
  DO i = 1, klon_loc
     dp(i, 1) = monthly_means(i, 1) / b(1)
  END DO

  !-- Rows 2 to n_months
  DO j = 2, n_months
     denom = b(j) - a(j) * cp(j-1)
     IF (j < n_months) THEN
        cp(j) = c(j) / denom
     END IF
     DO i = 1, klon_loc
        dp(i, j) = (monthly_means(i, j) - a(j) * dp(i, j-1)) / denom
     END DO
  END DO

  !---------------------------------------------------------------------------
  ! Thomas Algorithm: Back substitution
  !---------------------------------------------------------------------------

  !-- Last row is already solved
  DO i = 1, klon_loc
     targets(i, n_months) = dp(i, n_months)
  END DO

  !-- Remaining rows (backwards)
  DO j = n_months - 1, 1, -1
     DO i = 1, klon_loc
        targets(i, j) = dp(i, j) - cp(j) * targets(i, j+1)
     END DO
  END DO

END SUBROUTINE sz98_solve_tridiagonal


!===============================================================================
! SUBROUTINE: sz98_interpolate_daily
!
! Purpose:
!   Interpolate mid-month target values to a daily value using
!   piecewise-linear interpolation between month mid-points.
!
! Method:
!   1. Find which two month mid-points bracket the current day-of-year
!   2. Linear interpolation between those targets
!   3. Constant extension before first / after last mid-point
!
!
! Arguments:
!   targets      : Mid-month target values from sz98_solve_tridiagonal
!   month_midday : Day-of-year of each month mid-point (12)
!   day_of_year  : Current day of year (1 to year_len)
!   klon_loc     : Number of local grid points
!   n_months     : Number of months (12)
!   daily_value  : Output: interpolated emission for current day (klon_loc)
!===============================================================================
SUBROUTINE sz98_interpolate_daily(targets, month_midday, day_of_year, &
                                  klon_loc, n_months, daily_value)

  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: klon_loc
  INTEGER, INTENT(IN)  :: n_months
  REAL,    INTENT(IN)  :: targets(klon_loc, n_months)
  REAL,    INTENT(IN)  :: month_midday(n_months)
  INTEGER, INTENT(IN)  :: day_of_year
  REAL,    INTENT(OUT) :: daily_value(klon_loc)

  REAL    :: day_real, weight
  INTEGER :: m_before, m_after, im, i

  day_real = REAL(day_of_year)

  !---------------------------------------------------------------------------
  ! Find bracketing month mid-points
  !---------------------------------------------------------------------------
  IF (day_real <= month_midday(1)) THEN
     !-- Before first mid-point: constant extension (use January target)
     DO i = 1, klon_loc
        daily_value(i) = targets(i, 1)
     END DO
     RETURN

  ELSE IF (day_real >= month_midday(n_months)) THEN
     !-- After last mid-point: constant extension (use December target)
     DO i = 1, klon_loc
        daily_value(i) = targets(i, n_months)
     END DO
     RETURN

  ELSE
     !-- Find the interval [month_midday(m_before), month_midday(m_after)]
     m_before = 1
     m_after  = 2
     DO im = 1, n_months - 1
        IF (day_real >= month_midday(im) .AND. &
            day_real <  month_midday(im+1)) THEN
           m_before = im
           m_after  = im + 1
           EXIT
        END IF
     END DO

     !-- Linear interpolation between the two bracketing targets
     weight = (day_real - month_midday(m_before)) / &
              (month_midday(m_after) - month_midday(m_before))

     DO i = 1, klon_loc
        daily_value(i) = (1.0 - weight) * targets(i, m_before) + &
                           weight        * targets(i, m_after)
     END DO

  END IF

END SUBROUTINE sz98_interpolate_daily


!===============================================================================
! SUBROUTINE: co2_flux_corrections
!
! Purpose:
!   Read 2D spatially-varying flux correction fields from NetCDF files and
!   distribute them to all MPI ranks and OMP threads via scatter.
!
!   These corrections allow the user to apply gridpoint-level adjustments
!   to the ocean and/or land net CO2 flux, rather than a single global
!   scalar 
!
! NetCDF file format:
!   dimensions:
!     vector = klon_glo        (LMDZ 1D grid, same as emission files)
!   variables:
!     float vector(vector)     (coordinate variable)
!     float <varname>(vector)  (correction field)
!       units: "kgCO2/m2/s"   (same unit as fco2_ocean_cor / fco2_land_cor)
!
!   The variable name is configurable via fco2_ocean_cor_2d_var / fco2_land_cor_2d_var
!   in physiq.def (default: "fco2_cor").
!
! Unit convention:
!   The 2D field must be in kgCO2/m2/s per unit TOTAL grid-cell area
!   (i.e., already weighted by the appropriate surface-type fraction if needed).
!
! Sign convention:
!   Positive values REDUCE the source term (same as scalar mode):
!     source = fco2_ff + fco2_bb + fco2_land + fco2_ocean
!            - fco2_ocean_cor - fco2_land_cor
!
! Called from:
!   tracco2i, once at debutphy
!
! Authors:
!   2026-02: P. Cadule
!
!===============================================================================
SUBROUTINE co2_flux_corrections(debutphy)

  USE dimphy
  USE mod_grid_phy_lmdz, ONLY: klon_glo, nbp_lon, nbp_lat
  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
  USE mod_phys_lmdz_para, ONLY: scatter
  USE print_control_mod, ONLY: lunout

  USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_varid, nf95_open
  USE netcdf,   ONLY: nf90_get_var, nf90_nowrite, nf90_noerr, nf90_strerror

  ! Import 2D correction flags and filenames from carbon_cycle_mod
  USE carbon_cycle_mod, ONLY: read_fco2_ocean_cor_2d
  USE carbon_cycle_mod, ONLY: fco2_ocean_cor_2d_file
  USE carbon_cycle_mod, ONLY: fco2_ocean_cor_2d_var
  USE carbon_cycle_mod, ONLY: read_fco2_land_cor_2d
  USE carbon_cycle_mod, ONLY: fco2_land_cor_2d_file
  USE carbon_cycle_mod, ONLY: fco2_land_cor_2d_var

  ! Import target THREADPRIVATE arrays (already allocated in carbon_cycle_init)
  ! Section 5 computes: fco2_ocean_cor = base * pctsrf at every time step.
  USE carbon_cycle_mod, ONLY: fco2_ocean_cor_2d_base
  USE carbon_cycle_mod, ONLY: fco2_land_cor_2d_base

  IMPLICIT NONE

  !-----------------------------------------------------------------------------
  ! Arguments
  !-----------------------------------------------------------------------------
  LOGICAL, INTENT(IN) :: debutphy   ! .TRUE. at first physics call

  !-----------------------------------------------------------------------------
  ! Local variables for NetCDF reading 
  !-----------------------------------------------------------------------------
  REAL, ALLOCATABLE :: cor_glo(:)        ! Global correction field (klon_glo)
  REAL, ALLOCATABLE :: vector_coord(:)   ! Coordinate variable for validation

  INTEGER :: ncid_in, varid, ncerr
  INTEGER :: n_glo

  CHARACTER(LEN=30), PARAMETER :: modname = 'co2_flux_corrections'
  CHARACTER(LEN=80) :: abort_message

  !-----------------------------------------------------------------------------
  ! first call (debutphy)
  !-----------------------------------------------------------------------------
  IF (.NOT. debutphy) RETURN
  IF (.NOT. read_fco2_ocean_cor_2d .AND. .NOT. read_fco2_land_cor_2d) RETURN

  !=============================================================================
  ! A. OCEAN 2D CORRECTION
  !
  !    Read the ocean correction field
  !=============================================================================
  IF (read_fco2_ocean_cor_2d) THEN

     !------------------------------------------------------------------
     ! A.1 Read NetCDF 
     !------------------------------------------------------------------
!$OMP MASTER
     IF (is_mpi_root) THEN

        !-- Allocate global buffer on MPI root
        ALLOCATE(cor_glo(klon_glo))
        cor_glo(:) = 0.0

        !-- Open NetCDF file
        CALL nf95_open(TRIM(fco2_ocean_cor_2d_file), nf90_nowrite, ncid_in)

        !-- Validate spatial dimension
        CALL nf95_inq_varid(ncid_in, 'vector', varid)
        CALL nf95_gw_var(ncid_in, varid, vector_coord)
        n_glo = SIZE(vector_coord)
        DEALLOCATE(vector_coord)
        IF (n_glo /= klon_glo) THEN
           WRITE(abort_message, '(A,I8,A,I8)') &
                'ocean_cor_2d: vector=', n_glo, ' /= klon_glo=', klon_glo
           CALL abort_physic(modname, abort_message, 1)
        END IF

        !-- Read correction field
        CALL nf95_inq_varid(ncid_in, TRIM(fco2_ocean_cor_2d_var), varid)
        ncerr = nf90_get_var(ncid_in, varid, cor_glo)
        IF (ncerr /= nf90_noerr) THEN
           abort_message = 'ocean_cor_2d: '//TRIM(nf90_strerror(ncerr))
           CALL abort_physic(modname, abort_message, 1)
        END IF

        CALL nf95_close(ncid_in)

        !-- Diagnostic output
        WRITE(lunout,*) modname, ': --- Ocean 2D correction ---'
        WRITE(lunout,*) modname, ':   File = ', TRIM(fco2_ocean_cor_2d_file)
        WRITE(lunout,*) modname, ':   Var  = ', TRIM(fco2_ocean_cor_2d_var)
        WRITE(lunout,*) modname, ':   Min  = ', MINVAL(cor_glo), ' kgCO2/m2/s'
        WRITE(lunout,*) modname, ':   Max  = ', MAXVAL(cor_glo), ' kgCO2/m2/s'
        WRITE(lunout,*) modname, ':   Mean = ', SUM(cor_glo)/REAL(klon_glo), ' kgCO2/m2/s'

     END IF  ! is_mpi_root
!$OMP END MASTER

     !------------------------------------------------------------------
     ! A.2 Dummy allocation for non-root ranks
     !------------------------------------------------------------------
     IF (.NOT. ALLOCATED(cor_glo)) ALLOCATE(cor_glo(0))

     !------------------------------------------------------------------
     ! A.3 Scatter: distribute global field to all ranks/threads
     !     Section 5 computes fco2_ocean_cor = base * pctsrf every step.
     !------------------------------------------------------------------
     CALL scatter(cor_glo, fco2_ocean_cor_2d_base)

     !-- Release global buffer immediately
     IF (ALLOCATED(cor_glo)) DEALLOCATE(cor_glo)

  END IF  ! read_fco2_ocean_cor_2d

  !=============================================================================
  ! B. LAND 2D CORRECTION
  !
  !    Same pattern as ocean, targeting fco2_land_cor(:).
  !=============================================================================
  IF (read_fco2_land_cor_2d) THEN

     !------------------------------------------------------------------
     ! B.1 Read NetCDF 
     !------------------------------------------------------------------
!$OMP MASTER
     IF (is_mpi_root) THEN

        !-- Allocate global buffer on MPI root
        ALLOCATE(cor_glo(klon_glo))
        cor_glo(:) = 0.0

        !-- Open NetCDF file
        CALL nf95_open(TRIM(fco2_land_cor_2d_file), nf90_nowrite, ncid_in)

        !-- Validate spatial dimension
        CALL nf95_inq_varid(ncid_in, 'vector', varid)
        CALL nf95_gw_var(ncid_in, varid, vector_coord)
        n_glo = SIZE(vector_coord)
        DEALLOCATE(vector_coord)
        IF (n_glo /= klon_glo) THEN
           WRITE(abort_message, '(A,I8,A,I8)') &
                'land_cor_2d: vector=', n_glo, ' /= klon_glo=', klon_glo
           CALL abort_physic(modname, abort_message, 1)
        END IF

        !-- Read correction field
        CALL nf95_inq_varid(ncid_in, TRIM(fco2_land_cor_2d_var), varid)
        ncerr = nf90_get_var(ncid_in, varid, cor_glo)
        IF (ncerr /= nf90_noerr) THEN
           abort_message = 'land_cor_2d: '//TRIM(nf90_strerror(ncerr))
           CALL abort_physic(modname, abort_message, 1)
        END IF

        CALL nf95_close(ncid_in)

        !-- Diagnostic output
        WRITE(lunout,*) modname, ': --- Land 2D correction ---'
        WRITE(lunout,*) modname, ':   File = ', TRIM(fco2_land_cor_2d_file)
        WRITE(lunout,*) modname, ':   Var  = ', TRIM(fco2_land_cor_2d_var)
        WRITE(lunout,*) modname, ':   Min  = ', MINVAL(cor_glo), ' kgCO2/m2/s'
        WRITE(lunout,*) modname, ':   Max  = ', MAXVAL(cor_glo), ' kgCO2/m2/s'
        WRITE(lunout,*) modname, ':   Mean = ', SUM(cor_glo)/REAL(klon_glo), ' kgCO2/m2/s'

     END IF  ! is_mpi_root
!$OMP END MASTER

     !------------------------------------------------------------------
     ! B.2 Dummy allocation
     !------------------------------------------------------------------
     IF (.NOT. ALLOCATED(cor_glo)) ALLOCATE(cor_glo(0))

     !------------------------------------------------------------------
     ! B.3 Scatter to fco2_land_cor_2d_base 
     !------------------------------------------------------------------
     CALL scatter(cor_glo, fco2_land_cor_2d_base)

     !-- Release global buffer immediately
     IF (ALLOCATED(cor_glo)) DEALLOCATE(cor_glo)

  END IF  ! read_fco2_land_cor_2d

END SUBROUTINE co2_flux_corrections

!===============================================================================
! SUBROUTINE: reorder_global_1d
!
! Purpose:
!   Reorder a 1D global physics field (klon_glo) to correct grid ordering
!   convention mismatches between external NetCDF files and the LMDZ internal
!   physics grid.
!
!   LMDZ Physics Grid Structure:
!   The LMDZ physics grid is NOT a regular (nlon x nlat) rectangle.
!   It has SINGULAR POLE POINTS:
!
!     field(1)                                     = North Pole (1 point)
!     field(2 : nlon+1)                            = 1st interior band (N ==> S)
!     field(nlon+2 : 2*nlon+1)                     = 2nd interior band
!     ...
!     field((nlat_int-1)*nlon+2 : nlat_int*nlon+1) = last interior band
!     field(klon_glo)                              = South Pole (1 point)
!
!   where:
!     nlat_int = nlat - 2         (number of interior latitude bands)
!     klon_glo = nlon * nlat_int + 2 = nlon * (nlat - 2) + 2
!
!   THIS IS NOT nlon * nlat.
!
!
! Arguments:
!   field(klon)  : INOUT - flat 1D field in LMDZ physics layout
!                  klon = nlon * (nlat - 2) + 2
!   nlon         : IN    - longitude points per interior band (nbp_lon)
!   nlat         : IN    - total latitude count including poles (nbp_lat)
!   flip_lat     : IN    - .TRUE. to reverse latitude band ordering
!                          (S=>N becomes N=>S, or vice versa)
!                          Also swaps pole values.
!   shift_lon    : IN    - circular shift in longitude (grid points, interior only)
!                          Typical: nlon/2 for [-180,+180] => [0,360] conversion
!                          Positive = left rotation (eastward shift of origin)
!                          Poles are scalar averages: unaffected by longitude shift.
!
! Algorithm:
!   Step 1: If flip_lat - swap pole values (indices 1 and klon_glo)
!   Step 2: Copy interior bands to 2D work array
!   Step 3: Write back with reversed band order (if flip_lat)
!           and/or CSHIFT (if shift_lon /= 0)
!   Step 4: the field is now in LMDZ internal ordering
!
!
! Authors:
!   2026-02: P. Cadule - IPSL
!===============================================================================
SUBROUTINE reorder_global_1d(field, nlon, nlat, flip_lat, shift_lon)

  IMPLICIT NONE

  !-- Arguments
  INTEGER, INTENT(IN)    :: nlon       ! Number of longitudes (nbp_lon)
  INTEGER, INTENT(IN)    :: nlat       ! Total latitudes incl. poles (nbp_lat)
  REAL,    INTENT(INOUT) :: field(nlon * (nlat - 2) + 2)  ! LMDZ physics 1D field
  LOGICAL, INTENT(IN)    :: flip_lat   ! Flip latitude band ordering
  INTEGER, INTENT(IN)    :: shift_lon  ! Longitude circular shift (grid points)

  !-- Local variables
  INTEGER :: nlat_int                     ! Number of interior latitude bands
  INTEGER :: klon_glo                     ! Total field size (for clarity)
  REAL    :: work(nlon, nlat - 2)         ! 2D work array for interior bands
  REAL    :: temp_pole                    ! Temporary for pole swap
  INTEGER :: j, idx_src, idx_dst

  !-- Quick exit if no reordering requested
  IF (.NOT. flip_lat .AND. shift_lon == 0) RETURN

  nlat_int = nlat - 2
  klon_glo = nlon * nlat_int + 2

  !---------------------------------------------------------------------------
  ! Step 1: If flipping latitudes, swap pole points
  !
  !   LMDZ convention:
  !     field(1)         = North Pole
  !     field(klon_glo)  = South Pole
  !
  !   File convention (CF, S=>N):
  !     field(1)         = South Pole (needs to become North Pole)
  !     field(klon_glo)  = North Pole (needs to become South Pole)
  !---------------------------------------------------------------------------
  IF (flip_lat) THEN
     temp_pole = field(1)
     field(1) = field(klon_glo)
     field(klon_glo) = temp_pole
  END IF

  !---------------------------------------------------------------------------
  ! Step 2: Copy interior bands to 2D work array
  !
  !   Interior points span indices 2 to klon_glo-1, arranged as:
  !     band j: field(2 + (j-1)*nlon) .. field(1 + j*nlon)
  !   for j = 1 .. nlat_int
  !---------------------------------------------------------------------------
  DO j = 1, nlat_int
     idx_src = 2 + (j - 1) * nlon
     work(1:nlon, j) = field(idx_src : idx_src + nlon - 1)
  END DO

  !---------------------------------------------------------------------------
  ! Step 3: Write back with reversed band order and/or longitude shift
  !
  !   If flip_lat:  output band j reads from work band (nlat_int - j + 1)
  !   If shift_lon: each band is circularly shifted by shift_lon positions
  !                 using Fortran intrinsic CSHIFT (left rotation)
  !---------------------------------------------------------------------------
  DO j = 1, nlat_int
     idx_dst = 2 + (j - 1) * nlon

     IF (flip_lat) THEN
        IF (shift_lon /= 0) THEN
           field(idx_dst : idx_dst + nlon - 1) = &
                CSHIFT(work(1:nlon, nlat_int - j + 1), shift_lon)
        ELSE
           field(idx_dst : idx_dst + nlon - 1) = work(1:nlon, nlat_int - j + 1)
        END IF
     ELSE
        ! No latitude flip, longitude shift only
        field(idx_dst : idx_dst + nlon - 1) = CSHIFT(work(1:nlon, j), shift_lon)
     END IF
  END DO

END SUBROUTINE reorder_global_1d

END MODULE tracco2i_mod
