MODULE carbon_cycle_mod
!=======================================================================
!   Authors: Patricia Cadule and Laurent Fairhead
!            base sur un travail anterieur mene par Patricia Cadule et Josefine Ghattas
!
! Purpose:
!   Control module for interactive CO2/carbon-cycle coupling in LMDZ
!   Supports both emission-driven and concentration-driven modes
!   Supports both monthly and daily emission file reading with SZ98 interpolation
!   Supports 2D (lonxlat) correction fields for land and ocean
!
! Authors:
!   Patricia Cadule (IPSL), Laurent Fairhead (IPSL)
!   Contributors: Olivier Boucher (initial version)
!
! History:
!   2020: Initial version (emission-driven only)
!   2026-01: Added concentration-driven mode (P. Cadule)
!   2026-02: Added daily emission files and 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=F):
!       Surface fluxes feed back to atmospheric CO2
!   (2) Concentration-driven (carbon_cycle_conc_driven=T):
!       CO2 is prescribed from co2_ppm: fluxes are diagnostic only
!
! Emission Reading Modes:
!   (A) read_daily_co2ff = .TRUE.
!       Files have time = year_len (360/365/366)
!       Direct indexing by day_cur
!   (B) read_daily_co2ff = .FALSE.
!       Files have time = 12 (monthly)
!       Apply SZ98 interpolation to get daily values
!
! level_coupling_esm: level of coupling of the biogeochemical fields
! Definitions:
! level_coupling_esm = 0 ! No field exchange
! level_coupling_esm = 1 ! Field exchange between LMDZ and ORCHIDEE models
! level_coupling_esm = 2 ! Field exchange between LMDZ and NEMO models
! level_coupling_esm = 3 ! Field exchange between LMDZ, ORCHIDEE, and NEMO
!=======================================================================

  IMPLICIT NONE
  SAVE
  PRIVATE

  ! Public interfaces
  PUBLIC :: carbon_cycle_init, infocfields_init
  PUBLIC :: ensure_real_1d

  !=========================================================================
  ! CONFIGURATION FLAGS (read from physiq.def via conf_phys_m)
  !=========================================================================

  ! If TRUE, prescribed CO2 (co2_ppm) is sent to land/ocean models
  ! and atmospheric CO2 tracer does not evolve from fluxes
  LOGICAL, PUBLIC :: carbon_cycle_conc_driven = .FALSE.
!$OMP THREADPRIVATE(carbon_cycle_conc_driven)

  ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/PISCES
  LOGICAL, PUBLIC :: carbon_cycle_cpl = .FALSE.
!$OMP THREADPRIVATE(carbon_cycle_cpl)

  ! 3D transport of CO2 in the atmosphere
  LOGICAL, PUBLIC :: carbon_cycle_tr = .FALSE.
!$OMP THREADPRIVATE(carbon_cycle_tr)

  ! CO2 interactive radiatively
  LOGICAL, PUBLIC :: carbon_cycle_rad = .FALSE.
!$OMP THREADPRIVATE(carbon_cycle_rad)

  ! Level of ESM coupling (0, 1, 2, or 3)
  INTEGER, PUBLIC :: level_coupling_esm = 0
!$OMP THREADPRIVATE(level_coupling_esm)

  ! Method selection for inputs:
  ! .TRUE.  = Read daily file (365/366 steps) directly
  ! .FALSE. = Read monthly file (12 steps) and interpolate via Taylor/SZ98
  LOGICAL, PUBLIC :: read_daily_co2ff = .FALSE.
!$OMP THREADPRIVATE(read_daily_co2ff)

  !=========================================================================
  ! FLUX CORRECTION FLAGS AND VALUES
  !=========================================================================

  LOGICAL, PUBLIC :: read_fco2_ocean_cor = .FALSE.
!$OMP THREADPRIVATE(read_fco2_ocean_cor)
  REAL, PUBLIC :: var_fco2_ocean_cor = 0.0
!$OMP THREADPRIVATE(var_fco2_ocean_cor)

  REAL, PUBLIC :: ocean_area_tot = 0.0
!$OMP THREADPRIVATE(ocean_area_tot)

  LOGICAL, PUBLIC :: read_fco2_land_cor = .FALSE.
!$OMP THREADPRIVATE(read_fco2_land_cor)
  REAL, PUBLIC :: var_fco2_land_cor = 0.0
!$OMP THREADPRIVATE(var_fco2_land_cor)

  REAL, PUBLIC :: land_area_tot = 0.0
!$OMP THREADPRIVATE(land_area_tot)

  !=========================================================================
  ! FLUX CORRECTION FLAGS AND VALUES (2D)
  !=========================================================================

  ! 2D FLUX CORRECTION: OCEAN

  ! Read a 2D correction field from NetCDF for the ocean
  LOGICAL, PUBLIC :: read_fco2_ocean_cor_2d = .FALSE.
!$OMP THREADPRIVATE(read_fco2_ocean_cor_2d)

  ! NetCDF file containing the ocean 2D correction
  CHARACTER(LEN=256), PUBLIC :: fco2_ocean_cor_2d_file = 'sflx_lmdz_co2_ocean_cor.nc'
!$OMP THREADPRIVATE(fco2_ocean_cor_2d_file)

  ! Variable name inside the NetCDF file
  CHARACTER(LEN=64), PUBLIC :: fco2_ocean_cor_2d_var = 'fco2_cor'
!$OMP THREADPRIVATE(fco2_ocean_cor_2d_var)

  ! 2D correction BASE field 
  ! This array stores the raw correction rate from the piControl NetCDF file
  ! At each time step, tracco2i computes:
  !   fco2_ocean_cor = fco2_ocean_cor_2d_base * pctsrf(is_oce + is_sic)
  ! Populated once at debutphy by co2_flux_corrections via scatter
  ! Dimension: (klon) on the local MPI/OMP domain
  REAL, ALLOCATABLE, PUBLIC :: fco2_ocean_cor_2d_base(:)
!$OMP THREADPRIVATE(fco2_ocean_cor_2d_base)

  ! 2D FLUX CORRECTION: LAND

  ! Read a 2D correction field from NetCDF for land
  LOGICAL, PUBLIC :: read_fco2_land_cor_2d = .FALSE.
!$OMP THREADPRIVATE(read_fco2_land_cor_2d)

  ! NetCDF file containing the land 2D correction
  CHARACTER(LEN=256), PUBLIC :: fco2_land_cor_2d_file = 'sflx_lmdz_co2_land_cor.nc'
!$OMP THREADPRIVATE(fco2_land_cor_2d_file)

  ! Variable name inside the NetCDF file
  CHARACTER(LEN=64), PUBLIC :: fco2_land_cor_2d_var = 'fco2_cor'
!$OMP THREADPRIVATE(fco2_land_cor_2d_var)

  ! 2D correction BASE field 
  ! This array stores the raw correction rate from the piControl NetCDF file
  ! At each time step, tracco2i computes:
  !   fco2_land_cor = fco2_land_cor_2d_base * pctsrf(is_ter)
  ! Populated once at debutphy by co2_flux_corrections via scatter
  ! Dimension: (klon) on the local MPI/OMP domain
  REAL, ALLOCATABLE, PUBLIC :: fco2_land_cor_2d_base(:)
!$OMP THREADPRIVATE(fco2_land_cor_2d_base)

  !=========================================================================
  ! ADDITIONAL VARIABLES 
  !=========================================================================

  ! For reduction operations
  REAL, PUBLIC :: ocean_area_tot_glo = 0.0
  REAL, PUBLIC :: land_area_tot_glo = 0.0
  REAL, PUBLIC :: RCO2_glo_shared = 0.0

  !=========================================================================
  ! ADDITIONAL FLAGS
  !=========================================================================

  LOGICAL :: carbon_cycle_emis_comp_omp = .FALSE.
  LOGICAL, PUBLIC :: carbon_cycle_emis_comp = .FALSE.
!$OMP THREADPRIVATE(carbon_cycle_emis_comp)

  LOGICAL, PUBLIC :: carbon_cycle_rco2_timestep = .FALSE.
!$OMP THREADPRIVATE(carbon_cycle_rco2_timestep)

  LOGICAL, PUBLIC :: carbon_cycle_diag_timestep = .FALSE.
!$OMP THREADPRIVATE(carbon_cycle_diag_timestep)

  !=========================================================================
  ! RADIATIVE CO2 VARIABLES
  !=========================================================================

  REAL, PUBLIC :: RCO2_glo = 0.0   ! Global mean mixing ratio [kg CO2 /kg air]
!$OMP THREADPRIVATE(RCO2_glo)
  REAL, PUBLIC :: RCO2_tot = 0.0   ! Total CO2 mass [kg]
!$OMP THREADPRIVATE(RCO2_tot)

  !=========================================================================
  ! COMPATIBLE EMISSIONS 
  !=========================================================================

  REAL, PUBLIC :: emis_comp_kg_per_s_glo    = 0.0
  REAL, PUBLIC :: co2_mass_budget_resid_glo = 0.0
  REAL, PUBLIC :: co2_flux_surface_kg_s_glo = 0.0

  !=========================================================================
  ! INTERNAL CONFIGURATION
  !=========================================================================

  LOGICAL :: RCO2_inter_omp
  LOGICAL :: RCO2_inter
!$OMP THREADPRIVATE(RCO2_inter)

  REAL :: fos_fuel_s_omp
  REAL :: fos_fuel_s
!$OMP THREADPRIVATE(fos_fuel_s)

  REAL :: land_net_flux_s_omp
  REAL :: land_net_flux_s
!$OMP THREADPRIVATE(land_net_flux_s)

  REAL :: emis_land_s
!$OMP THREADPRIVATE(emis_land_s)

  REAL :: airetot
!$OMP THREADPRIVATE(airetot)

  INTEGER :: ntr_co2
!$OMP THREADPRIVATE(ntr_co2)

  !=========================================================================
  ! FLUX ARRAYS 
  ! Dimensions: (klon) where klon is local domain size
  !=========================================================================

  ! Daily ocean CO2 flux [gC m-2 d-1] - allocated/initialised in cpl_mod
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day
!$OMP THREADPRIVATE(fco2_ocn_day)

  ! Land CO2 flux [gC m-2 d-1]
  REAL, DIMENSION(:), ALLOCATABLE :: fco2_land_day
!$OMP THREADPRIVATE(fco2_land_day)

  ! Land-use change [gC m-2 d-1]
  REAL, DIMENSION(:), ALLOCATABLE :: fco2_lu_day
!$OMP THREADPRIVATE(fco2_lu_day)

  ! Fossil fuel [kg CO2 m-2 s-1]
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ff
!$OMP THREADPRIVATE(fco2_ff)

  ! Biomass burning [kg CO2 m-2 s-1]
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_bb
!$OMP THREADPRIVATE(fco2_bb)

  ! Net land flux [kg CO2 m-2 s-1]
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land
!$OMP THREADPRIVATE(fco2_land)

  ! nbp flux
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_nbp
!$OMP THREADPRIVATE(fco2_land_nbp)

  ! nep flux
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_nep
!$OMP THREADPRIVATE(fco2_land_nep)

  ! Land-use change
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fLuc
!$OMP THREADPRIVATE(fco2_land_fLuc)

  
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fwoodharvest
!$OMP THREADPRIVATE(fco2_land_fwoodharvest)


  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fHarvest
!$OMP THREADPRIVATE(fco2_land_fHarvest)

  ! Net ocean flux [kg CO2 m-2 s-1]
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocean
!$OMP THREADPRIVATE(fco2_ocean)

  ! Ocean flux correction
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocean_cor
!$OMP THREADPRIVATE(fco2_ocean_cor)

  ! Land flux correction
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_cor
!$OMP THREADPRIVATE(fco2_land_cor)


  !=========================================================================
  ! INSTANTANEOUS FIELDS 
  !=========================================================================

  ! Instantaneous fields (allocated in surf_land_orchidee)
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_inst
!$OMP THREADPRIVATE(fco2_land_inst)

  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_inst
!$OMP THREADPRIVATE(fco2_ocn_inst)

  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_lu_inst
!$OMP THREADPRIVATE(fco2_lu_inst)

  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: emis_comp_2d
!$OMP THREADPRIVATE(emis_comp_2d)

  REAL, DIMENSION(:,:), ALLOCATABLE, PUBLIC :: emis_comp_3d
!$OMP THREADPRIVATE(emis_comp_3d)

  !=========================================================================
  ! CO2 FIELD FOR COUPLER
  !=========================================================================

  !  [ppm] - allocated in phyetat0  
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send
!$OMP THREADPRIVATE(co2_send)

  REAL, ALLOCATABLE, PUBLIC :: source_co2_glo(:)
!$OMP THREADPRIVATE(source_co2_glo)
  REAL, ALLOCATABLE, PUBLIC :: area_glo(:)
!$OMP THREADPRIVATE(area_glo)

  REAL, ALLOCATABLE, PUBLIC :: co2_glo(:,:)      ! (klon_glo, klev)
!$OMP THREADPRIVATE(co2_glo)

  REAL, ALLOCATABLE, PUBLIC :: m_air_glo(:,:)    ! (klon_glo, klev)
!$OMP THREADPRIVATE(m_air_glo)

  !=========================================================================
  ! TRACER INDEX
  !=========================================================================

   !--OB: CO2 tracer index (temporary hardcode)
  INTEGER, PARAMETER, PUBLIC :: id_CO2 = 1  ! CO2 tracer index (temporary)

  !=========================================================================


  INTEGER :: iq,  ierr, stat, error

  ! Total number of coupling fields
  INTEGER, PUBLIC :: nbcf = 0
!$OMP THREADPRIVATE(nbcf)

  ! Number of IN fields (to LMDZ)
  INTEGER, PUBLIC :: nbcf_in = 0
!$OMP THREADPRIVATE(nbcf_in)

  ! IN from ORCHIDEE
  INTEGER, PUBLIC :: nbcf_in_orc = 0
!$OMP THREADPRIVATE(nbcf_in_orc)

  ! IN from INCA
  INTEGER, PUBLIC :: nbcf_in_inca = 0
!$OMP THREADPRIVATE(nbcf_in_inca)

  ! IN from NEMO
  INTEGER, PUBLIC :: nbcf_in_nemo = 0
!$OMP THREADPRIVATE(nbcf_in_nemo)

  ! IN from anthropogenic sources
  INTEGER, PUBLIC :: nbcf_in_ant = 0
!$OMP THREADPRIVATE(nbcf_in_ant)

  ! Number of OUT fields (from LMDZ)
  INTEGER, PUBLIC :: nbcf_out = 0
!$OMP THREADPRIVATE(nbcf_out)

  !=========================================================================
  ! COUPLING FIELD METADATA
  !=========================================================================

  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname
!$OMP THREADPRIVATE(cfname)

  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname_in
!$OMP THREADPRIVATE(cfname_in)

  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname_out
!$OMP THREADPRIVATE(cfname_out)

  CHARACTER(len=15), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfunits_in
!$OMP THREADPRIVATE(cfunits_in)

  CHARACTER(len=15), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfunits_out
!$OMP THREADPRIVATE(cfunits_out)

  CHARACTER(len=120), ALLOCATABLE, DIMENSION(:), PUBLIC :: cftext_in
!$OMP THREADPRIVATE(cftext_in)

  CHARACTER(len=120), ALLOCATABLE, DIMENSION(:), PUBLIC :: cftext_out
!$OMP THREADPRIVATE(cftext_out)

  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfmod1
!$OMP THREADPRIVATE(cfmod1)

  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfmod2
!$OMP THREADPRIVATE(cfmod2)

  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), PUBLIC :: field_out_names
!$OMP THREADPRIVATE(field_out_names)

  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), PUBLIC :: field_in_names
!$OMP THREADPRIVATE(field_in_names)

  !===========================================================================
  ! COUPLING FIELD DATA ARRAYS
  !===========================================================================

  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC :: fields_in
!$OMP THREADPRIVATE(fields_in)

  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC :: yfields_in
!$OMP THREADPRIVATE(yfields_in)

  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC :: fields_out
!$OMP THREADPRIVATE(fields_out)

  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC :: yfields_out
!$OMP THREADPRIVATE(yfields_out)

  !=========================================================================
  ! CO2 TRACER TYPE DEFINITION
  !=========================================================================

  TYPE, PUBLIC :: co2_trac_type
     CHARACTER(len = 8) :: name
     INTEGER            :: id
     CHARACTER(len=30)  :: file
     LOGICAL            :: cpl
     INTEGER            :: updatefreq
     INTEGER            :: readstep
     LOGICAL            :: updatenow
  END TYPE co2_trac_type

  INTEGER, PARAMETER :: maxco2trac = 5
  TYPE(co2_trac_type), DIMENSION(maxco2trac) :: co2trac

CONTAINS

!===============================================================================
! SUBROUTINE: carbon_cycle_init
!
! Purpose:
!   Initialise carbon cycle arrays and broadcast configuration flags
!   Called from tracco2i_init BEFORE the first physics time step
!
!===============================================================================
SUBROUTINE carbon_cycle_init()

  USE dimphy
  USE mod_phys_lmdz_para, ONLY: bcast
  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
  USE print_control_mod, ONLY: lunout

  IMPLICIT NONE

  CHARACTER(LEN=30), PARAMETER :: modname = 'carbon_cycle_init'
  INTEGER :: ierr
  INTEGER, PARAMETER :: MAX_YEAR_LEN = 366  ! Maximum days in a year


  !---------------------------------------------------------------------------
  ! 1. Broadcast configuration flags from OMP master to all threads
  !---------------------------------------------------------------------------
  CALL bcast(carbon_cycle_cpl)
  CALL bcast(carbon_cycle_tr)
  CALL bcast(carbon_cycle_rad)
  CALL bcast(carbon_cycle_conc_driven)
  CALL bcast(read_daily_co2ff)
  CALL bcast(level_coupling_esm)
  CALL bcast(read_fco2_ocean_cor)
  CALL bcast(read_fco2_land_cor)
  CALL bcast(var_fco2_ocean_cor)
  CALL bcast(var_fco2_land_cor)

  ! 2D flux correction flags
  CALL bcast(read_fco2_ocean_cor_2d)
  CALL bcast(fco2_ocean_cor_2d_file)
  CALL bcast(fco2_ocean_cor_2d_var)
  CALL bcast(read_fco2_land_cor_2d)
  CALL bcast(fco2_land_cor_2d_file)
  CALL bcast(fco2_land_cor_2d_var)

  CALL bcast(carbon_cycle_emis_comp)
  CALL bcast(carbon_cycle_rco2_timestep)
  CALL bcast(carbon_cycle_diag_timestep)

!$OMP BARRIER

  !---------------------------------------------------------------------------
  ! 2. Allocate THREADPRIVATE flux arrays 
  !---------------------------------------------------------------------------
  IF (carbon_cycle_cpl .OR. carbon_cycle_tr .OR. &
      carbon_cycle_rad .OR. carbon_cycle_conc_driven) THEN

     IF (klon <= 0) THEN
        CALL abort_physic(modname, 'Invalid klon value', 1)
     END IF

     ! Flux arrays (1D, size klon)
     IF (.NOT. ALLOCATED(fco2_land))              ALLOCATE(fco2_land(klon))
     IF (.NOT. ALLOCATED(fco2_land_nbp))          ALLOCATE(fco2_land_nbp(klon))
     IF (.NOT. ALLOCATED(fco2_land_nep))          ALLOCATE(fco2_land_nep(klon))
     IF (.NOT. ALLOCATED(fco2_land_fLuc))         ALLOCATE(fco2_land_fLuc(klon))
     IF (.NOT. ALLOCATED(fco2_land_fwoodharvest)) ALLOCATE(fco2_land_fwoodharvest(klon))
     IF (.NOT. ALLOCATED(fco2_land_fHarvest))     ALLOCATE(fco2_land_fHarvest(klon))
     IF (.NOT. ALLOCATED(fco2_ff))                ALLOCATE(fco2_ff(klon))
     IF (.NOT. ALLOCATED(fco2_bb))                ALLOCATE(fco2_bb(klon))
     IF (.NOT. ALLOCATED(fco2_ocean))             ALLOCATE(fco2_ocean(klon))
     IF (.NOT. ALLOCATED(fco2_ocean_cor))         ALLOCATE(fco2_ocean_cor(klon))
     IF (.NOT. ALLOCATED(fco2_land_cor))          ALLOCATE(fco2_land_cor(klon))
     IF (.NOT. ALLOCATED(fco2_ocn_day))           ALLOCATE(fco2_ocn_day(klon))
     IF (.NOT. ALLOCATED(co2_send))               ALLOCATE(co2_send(klon))

     ! 2D correction base arrays (only when 2D mode is active)
     ! These THREADPRIVATE arrays receive data via scatter in
     ! co2_flux_corrections. They store the raw correction rate
     ! (kgCO2/m2/s per unit total grid-cell area) from the piControl
     ! NetCDF file. At each time step, tracco2i multiplies by pctsrf
     IF (read_fco2_ocean_cor_2d) THEN
        IF (.NOT. ALLOCATED(fco2_ocean_cor_2d_base)) &
             ALLOCATE(fco2_ocean_cor_2d_base(klon))
        fco2_ocean_cor_2d_base(:) = 0.0
     END IF

     IF (read_fco2_land_cor_2d) THEN
        IF (.NOT. ALLOCATED(fco2_land_cor_2d_base)) &
             ALLOCATE(fco2_land_cor_2d_base(klon))
        fco2_land_cor_2d_base(:) = 0.0
     END IF

     ! Initialise to zero 
     fco2_land(:)              = 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
     fco2_ff(:)                = 0.0
     fco2_bb(:)                = 0.0
     fco2_ocean(:)             = 0.0
     fco2_ocean_cor(:)         = 0.0
     fco2_land_cor(:)          = 0.0
     fco2_ocn_day(:)           = 0.0
     co2_send(:)               = 0.0


!$OMP BARRIER
  END IF

  !---------------------------------------------------------------------------
  ! 3. Log initialisation status
  !---------------------------------------------------------------------------
  IF (is_omp_root .AND. is_mpi_root) THEN
     WRITE(lunout,*) '=============================================='
     WRITE(lunout,*) modname, ': Initialisation complete'
     WRITE(lunout,*) '  carbon_cycle_cpl         = ', carbon_cycle_cpl
     WRITE(lunout,*) '  carbon_cycle_tr          = ', carbon_cycle_tr
     WRITE(lunout,*) '  carbon_cycle_rad         = ', carbon_cycle_rad
     WRITE(lunout,*) '  carbon_cycle_conc_driven = ', carbon_cycle_conc_driven
     WRITE(lunout,*) '  read_daily_co2ff         = ', read_daily_co2ff
     WRITE(lunout,*) '  level_coupling_esm       = ', level_coupling_esm
     WRITE(lunout,*) '  klon (local domain)      = ', klon
     WRITE(lunout,*) '=============================================='
  END IF

END SUBROUTINE carbon_cycle_init
  
!===============================================================================
! SUBROUTINE: ensure_real_1d
!
! Purpose:
!   Allocate a 1D real array and initialise to zero
!   Handles both already-allocated and unallocated arrays
!
! Arguments:
!   arr  : Array to allocate (INTENT(INOUT))
!   n    : Size to allocate
!   name : Name for error messages
!
! Notes:
!   - If array is already allocated with correct size, does nothing
!   - If allocated with wrong size, reallocates
!   - Always initialises to zero
!===============================================================================
SUBROUTINE ensure_real_1d(arr, n, name)

  IMPLICIT NONE

  REAL, ALLOCATABLE, INTENT(INOUT) :: arr(:)
  INTEGER, INTENT(IN) :: n
  CHARACTER(LEN=*), INTENT(IN) :: name

  INTEGER :: ierr
  LOGICAL :: newly_allocated

  newly_allocated = .FALSE.

  ! Check if reallocation is needed
  IF (ALLOCATED(arr)) THEN
     IF (SIZE(arr) /= n) THEN
        DEALLOCATE(arr)
        newly_allocated = .TRUE.
     END IF
  ELSE
     newly_allocated = .TRUE.
  END IF

  ! Allocate if necessary
  IF (.NOT. ALLOCATED(arr)) THEN
     ALLOCATE(arr(n), STAT=ierr)
     IF (ierr /= 0) THEN
        CALL abort_physic('ensure_real_1d', 'Allocation failed: '//TRIM(name), 1)
     END IF
     newly_allocated = .TRUE.
  END IF

  ! Always initialise to zero 
  arr(:) = 0.0

END SUBROUTINE ensure_real_1d


!===============================================================================
! SUBROUTINE: infocfields_init
!
! Purpose:
!   Read coupling field definitions from coupling_fields.def
!   Allocate coupling field arrays (fields_in, fields_out)
!
! Called from:
!   phys_state_var_init (before carbon_cycle_init)
!===============================================================================
SUBROUTINE infocfields_init()

  USE dimphy
  USE mod_phys_lmdz_para, ONLY: bcast
  USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
  USE print_control_mod, ONLY: lunout

  IMPLICIT NONE

  INTEGER :: iq, ierr, error

  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), SAVE  :: cfname_root
!$OMP THREADPRIVATE(cfname_root)  
  CHARACTER(LEN=120), ALLOCATABLE, DIMENSION(:), SAVE :: cftext_root
!$OMP THREADPRIVATE(cftext_root)
  CHARACTER(LEN=15), ALLOCATABLE, DIMENSION(:), SAVE  :: cfunits_root
!$OMP THREADPRIVATE(cfunits_root)
  CHARACTER(len=3), ALLOCATABLE, DIMENSION(:)         :: cfintent_root
  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:)         :: cfmod1_root
  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:)         :: cfmod2_root
  LOGICAL, ALLOCATABLE, DIMENSION(:), SAVE            :: mask_in_root
!$OMP THREADPRIVATE(mask_in_root)  
  LOGICAL, ALLOCATABLE, DIMENSION(:), SAVE            :: mask_out_root
!$OMP THREADPRIVATE(mask_out_root)  

  CHARACTER(len=*), PARAMETER :: modname = "infocfields_init"
  CHARACTER(len=10), SAVE     :: planet_type = "earth"


  nbcf = 0
  nbcf_in = 0
  nbcf_out = 0

  IF (planet_type == 'earth') THEN

     IF (is_mpi_root .AND. is_omp_root) THEN

        IF (level_coupling_esm > 0) THEN

           OPEN(200, file='coupling_fields.def', form='formatted', status='old', iostat=ierr)

           IF (ierr == 0) THEN
              WRITE(lunout,*) TRIM(modname), ': Open coupling_fields.def : ok'
              READ(200,*) nbcf
              WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf=', nbcf

              ALLOCATE(cfname_root(nbcf))
              ALLOCATE(cfintent_root(nbcf))
              ALLOCATE(cfmod1_root(nbcf))
              ALLOCATE(cfmod2_root(nbcf))
              ALLOCATE(cftext_root(nbcf))
              ALLOCATE(cfunits_root(nbcf))
              ALLOCATE(mask_in_root(nbcf))
              ALLOCATE(mask_out_root(nbcf))

              nbcf_in = 0
              nbcf_out = 0

              DO iq = 1, nbcf
                 WRITE(lunout,*) 'infofields : field=', iq
                 READ(200,'(A15,3X,A3,3X,A5,3X,A5,3X,A120,3X,A15)', IOSTAT=ierr) &
                    cfname_root(iq), cfintent_root(iq), cfmod1_root(iq), &
                    cfmod2_root(iq), cftext_root(iq), cfunits_root(iq)
                 cfname_root(iq)   = TRIM(cfname_root(iq))
                 cfintent_root(iq) = TRIM(cfintent_root(iq))
                 cfmod1_root(iq)   = TRIM(cfmod1_root(iq))
                 cfmod2_root(iq)   = TRIM(cfmod2_root(iq))
                 cftext_root(iq)   = TRIM(cftext_root(iq))
                 cfunits_root(iq)  = TRIM(cfunits_root(iq))

                 WRITE(lunout,*) 'coupling field: ', cfname_root(iq), &
                                ', number: ', iq, ', INTENT: ', cfintent_root(iq)

                 IF (nbcf_in + nbcf_out < nbcf) THEN
                    IF (cfintent_root(iq) /= 'OUT') THEN
                       nbcf_in = nbcf_in + 1
                       mask_in_root(iq) = .TRUE.
                       mask_out_root(iq) = .FALSE.
                    ELSE IF (cfintent_root(iq) == 'OUT') THEN
                       nbcf_out = nbcf_out + 1
                       mask_in_root(iq) = .FALSE.
                       mask_out_root(iq) = .TRUE.
                    END IF
                 ELSE
                    WRITE(lunout,*) 'abort_gcm --- nbcf    : ', nbcf
                    WRITE(lunout,*) 'abort_gcm --- nbcf_in : ', nbcf_in
                    WRITE(lunout,*) 'abort_gcm --- nbcf_out: ', nbcf_out
                    CALL abort_physic(modname, 'Problem in coupling fields definition', 1)
                 END IF
              END DO
           ELSE
              WRITE(lunout,*) TRIM(modname), ': Problem opening coupling_fields.def'
              WRITE(lunout,*) TRIM(modname), ': WARNING using default values'
           END IF
           CLOSE(200)

        END IF ! level_coupling_esm

     END IF ! is_mpi_root .AND. is_omp_root

!$OMP BARRIER

     CALL bcast(nbcf)
     CALL bcast(nbcf_in)
     CALL bcast(nbcf_out)

     WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf    =', nbcf
     WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf_in =', nbcf_in
     WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf_out=', nbcf_out

     ALLOCATE(cfname(nbcf))
     ALLOCATE(cfname_in(nbcf_in))
     ALLOCATE(cftext_in(nbcf_in))
     ALLOCATE(cfname_out(nbcf_out))
     ALLOCATE(cftext_out(nbcf_out))
     ALLOCATE(cfmod1(nbcf))
     ALLOCATE(cfmod2(nbcf))
     ALLOCATE(cfunits_in(nbcf_in))
     ALLOCATE(cfunits_out(nbcf_out))

     IF (is_mpi_root .AND. is_omp_root) THEN
        IF (nbcf > 0)     cfname = cfname_root
        IF (nbcf_in > 0)  cfname_in = PACK(cfname_root, mask_in_root)
        IF (nbcf_out > 0) cfname_out = PACK(cfname_root, mask_out_root)
        IF (nbcf_in > 0)  cftext_in = PACK(cftext_root, mask_in_root)
        IF (nbcf_out > 0) cftext_out = PACK(cftext_root, mask_out_root)
        IF (nbcf > 0)     cfmod1 = cfmod1_root
        IF (nbcf > 0)     cfmod2 = cfmod2_root
        IF (nbcf_in > 0)  cfunits_in = PACK(cfunits_root, mask_in_root)
        IF (nbcf_out > 0) cfunits_out = PACK(cfunits_root, mask_out_root)

        nbcf_in_orc = 0
        nbcf_in_nemo = 0
        nbcf_in_inca = 0
        nbcf_in_ant = 0

        DO iq = 1, nbcf
           IF (cfmod1(iq) == "ORC")  nbcf_in_orc  = nbcf_in_orc  + 1
           IF (cfmod1(iq) == "NEMO") nbcf_in_nemo = nbcf_in_nemo + 1
           IF (cfmod1(iq) == "INCA") nbcf_in_inca = nbcf_in_inca + 1
           IF (cfmod1(iq) == "ALL")  nbcf_in_orc  = nbcf_in_orc  + 1
           IF (cfmod1(iq) == "ALL")  nbcf_in_nemo = nbcf_in_nemo + 1
           IF (cfmod1(iq) == "ALL")  nbcf_in_inca = nbcf_in_inca + 1
           IF (cfmod1(iq) == "ANT")  nbcf_in_ant  = nbcf_in_ant  + 1
        END DO
     END IF

!$OMP BARRIER

     CALL bcast(nbcf_in_orc)
     CALL bcast(nbcf_in_nemo)
     CALL bcast(nbcf_in_inca)
     CALL bcast(nbcf_in_ant)

     IF (nbcf_in > 0) THEN
        DO iq = 1, nbcf_in
           CALL bcast(cfname_in(iq))
           CALL bcast(cftext_in(iq))
           CALL bcast(cfunits_in(iq))
        END DO
     END IF

     IF (nbcf_out > 0) THEN
        DO iq = 1, nbcf_out
           CALL bcast(cfname_out(iq))
           CALL bcast(cftext_out(iq))
           CALL bcast(cfunits_out(iq))
        END DO
     END IF

     IF (nbcf > 0) THEN
        DO iq = 1, nbcf
           CALL bcast(cfmod1(iq))
           CALL bcast(cfmod2(iq))
        END DO
     END IF

     IF (nbcf_in > 0)  WRITE(lunout,*) 'infocfields_mod --- cfname_in: ', cfname_in
     IF (nbcf_out > 0) WRITE(lunout,*) 'infocfields_mod --- cfname_out: ', cfname_out

  ELSE
     ! Default values for other planets
     nbcf = 0
     nbcf_in = 0
     nbcf_out = 0
  END IF ! planet_type

  ! Allocate coupling field data arrays
  ALLOCATE(fields_in(klon, nbcf_in), stat=error)
  IF (error /= 0) CALL abort_physic(modname, 'Pb in allocation fields_in', 1)
  ALLOCATE(yfields_in(klon, nbcf_in), stat=error)
  IF (error /= 0) CALL abort_physic(modname, 'Pb in allocation yfields_in', 1)
  ALLOCATE(fields_out(klon, nbcf_out), stat=error)
  IF (error /= 0) CALL abort_physic(modname, 'Pb in allocation fields_out', 1)
  ALLOCATE(yfields_out(klon, nbcf_out), stat=error)
  IF (error /= 0) CALL abort_physic(modname, 'Pb in allocation yfields_out', 1)

END SUBROUTINE infocfields_init

END MODULE carbon_cycle_mod
