MODULE carbon_cycle_mod
!=======================================================================
!   Authors: Patricia Cadule and Laurent Fairhead
!            base sur un travail anterieur mene par Patricia Cadule et Josefine Ghattas
!
!  Purpose and description:
!  -----------------------
! Control module for the carbon CO2 tracers :
!   - Identification
!   - Get concentrations comming from coupled model or read from file to tracers
!   - Calculate new RCO2 for radiation scheme
!   - Calculate new carbon flux for sending to coupled models (PISCES and ORCHIDEE)
!
! Module permettant de mettre a jour les champs (puits et sources) pour le
! transport de CO2 en online (IPSL-CM et LMDZOR) et offline (lecture de carte)
!
! Le cas online/offline est defini par le flag carbon_cycle_cpl (y/n)
! Le transport du traceur CO2 est defini par le flag carbon_cycle_tr (y/n)
! la provenance des champs (termes de puits) est denini par le flag level_coupling_esm
!
! level_coupling_esm : level of coupling of the biogeochemical fields between
! LMDZ, ORCHIDEE and NEMO
! Definitions of level_coupling_esm in physiq.def
! level_coupling_esm = 0  ! No field exchange between LMDZ and ORCHIDEE models
!                         ! No field exchange between LMDZ and NEMO
! level_coupling_esm = 1  ! Field exchange between LMDZ and ORCHIDEE models
!                         ! No field exchange between LMDZ and NEMO models
! level_coupling_esm = 2  ! No field exchange between LMDZ and ORCHIDEE models
!                         ! Field exchange between LMDZ and NEMO models
! level_coupling_esm = 3  ! Field exchange between LMDZ and ORCHIDEE models
!                         ! Field exchange between LMDZ and NEMO models
!=======================================================================

  IMPLICIT NONE
  SAVE
  PRIVATE
  PUBLIC :: carbon_cycle_init, carbon_cycle, infocfields_init

! Variables read from parmeter file physiq.def
  LOGICAL, PUBLIC :: carbon_cycle_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES) 
!$OMP THREADPRIVATE(carbon_cycle_cpl)
  LOGICAL, PUBLIC :: carbon_cycle_tr        ! 3D transport of CO2 in the atmosphere, parameter read in conf_phys
!$OMP THREADPRIVATE(carbon_cycle_tr)
  LOGICAL, PUBLIC :: carbon_cycle_rad       ! CO2 interactive radiatively 
!$OMP THREADPRIVATE(carbon_cycle_rad)
  INTEGER, PUBLIC :: level_coupling_esm ! Level of coupling for the ESM - 0, 1, 2, 3
!$OMP THREADPRIVATE(level_coupling_esm)
  REAL, PUBLIC :: RCO2_glo, RCO2_tot
!$OMP THREADPRIVATE(RCO2_glo, RCO2_tot)

  LOGICAL :: carbon_cycle_emis_comp_omp=.FALSE.
  LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible
!$OMP THREADPRIVATE(carbon_cycle_emis_comp)

  LOGICAL :: RCO2_inter_omp
  LOGICAL :: RCO2_inter  ! RCO2 interactive : if true calculate new value RCO2 for the radiation scheme
!$OMP THREADPRIVATE(RCO2_inter)

! Scalare values when no transport, from physiq.def
  REAL :: fos_fuel_s_omp
  REAL :: fos_fuel_s  ! carbon_cycle_fos_fuel dans physiq.def
!$OMP THREADPRIVATE(fos_fuel_s)
  REAL :: emis_land_s ! not yet implemented
!$OMP THREADPRIVATE(emis_land_s)

  REAL :: airetot     ! Total area of the earth surface
!$OMP THREADPRIVATE(airetot)

  INTEGER :: ntr_co2  ! Number of tracers concerning the carbon cycle
!$OMP THREADPRIVATE(ntr_co2)

! fco2_ocn_day : flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]. Allocation and initalization done in cpl_mod
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day 
!$OMP THREADPRIVATE(fco2_ocn_day)

  REAL, DIMENSION(:), ALLOCATABLE :: fco2_land_day   ! flux CO2 from land for 1 day (cumulated)  [gC/m2/d]
!$OMP THREADPRIVATE(fco2_land_day)
  REAL, DIMENSION(:), ALLOCATABLE :: fco2_lu_day     ! Emission from land use change for 1 day (cumulated) [gC/m2/d]
!$OMP THREADPRIVATE(fco2_lu_day)
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ff ! Emission from fossil fuel [kgCO2/m2/s]
!$OMP THREADPRIVATE(fco2_ff)
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_bb ! Emission from biomass burning [kgCO2/m2/s]
!$OMP THREADPRIVATE(fco2_bb)
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
!$OMP THREADPRIVATE(fco2_land)
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocean ! Net flux from ocean [kgCO2/m2/s]
!$OMP THREADPRIVATE(fco2_ocean)

  REAL, DIMENSION(:,:), ALLOCATABLE :: dtr_add       ! Tracer concentration to be injected 
!$OMP THREADPRIVATE(dtr_add)

! Following 2 fields will be allocated and initialized in surf_land_orchidee
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_inst  ! flux CO2 from land at one time step
!$OMP THREADPRIVATE(fco2_land_inst)
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_lu_inst    ! Emission from land use change at one time step
!$OMP THREADPRIVATE(fco2_lu_inst)

! Calculated co2 field to be send to the ocean via the coupler and to ORCHIDEE 
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send ! Field allocated in phyetat0
!$OMP THREADPRIVATE(co2_send)

  INTEGER, PARAMETER, PUBLIC :: id_CO2=1              !--temporaire OB -- to be changed

! nbfields : total number of fields 
  INTEGER, PUBLIC :: nbcf
!$OMP THREADPRIVATE(nbcf)

! nbcf_in : number of fields IN
  INTEGER, PUBLIC  :: nbcf_in
!$OMP THREADPRIVATE(nbcf_in)

! nbcf_in_orc : number of fields IN 
  INTEGER, PUBLIC  :: nbcf_in_orc
!$OMP THREADPRIVATE(nbcf_in_orc)

! nbcf_in_inca : number of fields IN (from INCA)
  INTEGER, PUBLIC  :: nbcf_in_inca
!$OMP THREADPRIVATE(nbcf_in_inca)

! nbcf_in_nemo : number of fields IN (from nemo)
  INTEGER, PUBLIC  :: nbcf_in_nemo
!$OMP THREADPRIVATE(nbcf_in_nemo)

! nbcf_in_ant : number of fields IN (from anthropogenic sources)
  INTEGER, PUBLIC  :: nbcf_in_ant
!$OMP THREADPRIVATE(nbcf_in_ant)

! nbcf_out : number of fields OUT
  INTEGER, PUBLIC :: nbcf_out
!$OMP THREADPRIVATE(nbcf_out)

! Name of variables
  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname     ! coupling field short name for restart (?) and diagnostics
!$OMP THREADPRIVATE(cfname)

  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname_in  ! coupling field short name for restart (?) and diagnostics
!$OMP THREADPRIVATE(cfname_in)

  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname_out ! coupling field short name for restart (?) and diagnostics
!$OMP THREADPRIVATE(cfname_out)

  CHARACTER(len=15), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfunits_in  !  coupling field units for diagnostics
!$OMP THREADPRIVATE(cfunits_in)

  CHARACTER(len=15), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfunits_out !  coupling field units for diagnostics
!$OMP THREADPRIVATE(cfunits_out)

  CHARACTER(len=120), ALLOCATABLE, DIMENSION(:), PUBLIC :: cftext_in  ! coupling field long name for diagnostics
!$OMP THREADPRIVATE(cftext_in)

  CHARACTER(len=120), ALLOCATABLE, DIMENSION(:), PUBLIC :: cftext_out ! coupling field long name for diagnostics
!$OMP THREADPRIVATE(cftext_out)

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

  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfmod2 ! model 2
!$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)

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

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

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

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

  TYPE, PUBLIC :: co2_trac_type
     CHARACTER(len = 8) :: name       ! Tracer name in tracer.def
     INTEGER            :: id         ! Index in total tracer list, tr_seri
     CHARACTER(len=30)  :: file       ! File name
     LOGICAL            :: cpl        ! True if this tracers is coupled from ORCHIDEE or PISCES. 
                                      ! False if read from file.
     INTEGER            :: updatefreq ! Frequence to inject in second
     INTEGER            :: readstep   ! Actual time step to read in file
     LOGICAL            :: updatenow  ! True if this tracer should be updated this time step
  END TYPE co2_trac_type
  INTEGER,PARAMETER :: maxco2trac=5  ! Maximum number of different CO2 fluxes
  TYPE(co2_trac_type), DIMENSION(maxco2trac) :: co2trac

CONTAINS
  
  SUBROUTINE carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
! This subroutine is called from traclmdz_init, only at first timestep.
! - Read controle parameters from .def input file
! - Search for carbon tracers and set default values
! - Allocate variables
! - Test for compatibility

    USE dimphy
    USE geometry_mod, ONLY : cell_area 
    USE mod_phys_lmdz_transfert_para
    USE infotrac_phy, ONLY: nbtr, nqo, niadv, tname
    USE IOIPSL
    USE surface_data, ONLY : ok_veget, type_ocean
    USE phys_cal_mod, ONLY : mth_len
    USE print_control_mod, ONLY: lunout

    IMPLICIT NONE
    INCLUDE "clesphys.h"
 
! Input argument
    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA]  
    REAL,INTENT(IN)                           :: pdtphys ! length of time step in physiq (sec)

! InOutput arguments
    LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: aerosol
    LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: radio

! Local variables
    INTEGER               :: ierr, it, iiq, itc
    INTEGER               :: teststop

! 1) Read controle parameters from .def input file
! ------------------------------------------------
    ! Read fosil fuel value if no transport
    IF (.NOT. carbon_cycle_tr) THEN
!$OMP MASTER
       fos_fuel_s_omp = 0.
       CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s_omp)
!$OMP END MASTER
!$OMP BARRIER
       fos_fuel_s=fos_fuel_s_omp
       WRITE(lunout,*) 'carbon_cycle_fos_fuel = ', fos_fuel_s 
    END IF

    ! Read parmeter for calculation compatible emission
    IF (.NOT. carbon_cycle_tr) THEN
!$OMP MASTER
       carbon_cycle_emis_comp_omp=.FALSE.
       CALL getin('carbon_cycle_emis_comp',carbon_cycle_emis_comp_omp)
!$OMP END MASTER
!$OMP BARRIER
       carbon_cycle_emis_comp=carbon_cycle_emis_comp_omp
       WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
       IF (carbon_cycle_emis_comp) THEN
          CALL abort_physic('carbon_cycle_init', 'carbon_cycle_emis_comp option not yet implemented!!',1)
       END IF
    END IF

    ! Read parameter for interactive calculation of the CO2 value for the radiation scheme
!$OMP MASTER
    RCO2_inter_omp=.FALSE.
    CALL getin('RCO2_inter',RCO2_inter_omp)
!$OMP END MASTER
!$OMP BARRIER
    RCO2_inter=RCO2_inter_omp
    WRITE(lunout,*) 'RCO2_inter = ', RCO2_inter
    IF (RCO2_inter) THEN
       WRITE(lunout,*) 'RCO2 will be recalculated once a day'
       WRITE(lunout,*) 'RCO2 initial = ', RCO2
    END IF


! 2) Search for carbon tracers and set default values
! ---------------------------------------------------
    itc=0
    DO it=1,nbtr
!!       iiq=niadv(it+2)                                                            ! jyg
       iiq=niadv(it+nqo)                                                            ! jyg
       
       SELECT CASE(tname(iiq))
       CASE("fCO2_ocn")
          itc = itc + 1
          co2trac(itc)%name='fCO2_ocn'
          co2trac(itc)%id=it
          co2trac(itc)%file='fl_co2_ocean.nc'
          IF (carbon_cycle_cpl .AND. type_ocean=='couple') THEN 
             co2trac(itc)%cpl=.TRUE.
             co2trac(itc)%updatefreq = 86400 ! Once a day as the coupling with OASIS/PISCES
          ELSE
             co2trac(itc)%cpl=.FALSE.
             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
          END IF
       CASE("fCO2_land")
          itc = itc + 1
          co2trac(itc)%name='fCO2_land'
          co2trac(itc)%id=it
          co2trac(itc)%file='fl_co2_land.nc'
          IF (carbon_cycle_cpl .AND. ok_veget) THEN 
             co2trac(itc)%cpl=.TRUE.
             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
          ELSE
             co2trac(itc)%cpl=.FALSE.
!             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
          END IF
       CASE("fCO2_land_use")
          itc = itc + 1
          co2trac(itc)%name='fCO2_land_use'
          co2trac(itc)%id=it
          co2trac(itc)%file='fl_co2_land_use.nc'
          IF (carbon_cycle_cpl .AND. ok_veget) THEN 
             co2trac(it)%cpl=.TRUE.
             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
          ELSE
             co2trac(itc)%cpl=.FALSE.
             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
          END IF
       CASE("fCO2_fos_fuel")
          itc = itc + 1
          co2trac(itc)%name='fCO2_fos_fuel'
          co2trac(itc)%id=it
          co2trac(itc)%file='fossil_fuel.nc'
          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
!         co2trac(itc)%updatefreq = 86400  ! 86400sec = 24H Cadule case
          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
       CASE("fCO2_bbg")
          itc = itc + 1
          co2trac(itc)%name='fCO2_bbg'
          co2trac(itc)%id=it
          co2trac(itc)%file='fl_co2_bbg.nc'
          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
       CASE("fCO2")
          ! fCO2 : One tracer transporting the total CO2 flux
          itc = itc + 1
          co2trac(itc)%name='fCO2'
          co2trac(itc)%id=it
          co2trac(itc)%file='fl_co2.nc'
          IF (carbon_cycle_cpl) THEN 
             co2trac(itc)%cpl=.TRUE.
          ELSE
             co2trac(itc)%cpl=.FALSE.
          END IF
          co2trac(itc)%updatefreq = 86400
          ! DOES THIS WORK ???? Problematic due to implementation of the coupled fluxes...
          CALL abort_physic('carbon_cycle_init','transport of total CO2 has to be implemented and tested',1)
       END SELECT
    END DO

    ! Total number of carbon CO2 tracers
    ntr_co2 = itc 
    
    ! Definition of control varaiables for the tracers
    DO it=1,ntr_co2
       aerosol(co2trac(it)%id) = .FALSE.
       radio(co2trac(it)%id)   = .FALSE.
    END DO
    
    ! Vector indicating which timestep to read for each tracer
    ! Always start read in the beginning of the file
    co2trac(:)%readstep = 0
   

! 3) Allocate variables
! ---------------------
    ! Allocate vector for storing fluxes to inject
    ALLOCATE(dtr_add(klon,maxco2trac), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 11',1)       
    
    ! Allocate variables for cumulating fluxes from ORCHIDEE
    IF (RCO2_inter) THEN
       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
          ALLOCATE(fco2_land_day(klon), stat=ierr)
          IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 2',1)
          fco2_land_day(1:klon) = 0.
          
          ALLOCATE(fco2_lu_day(klon), stat=ierr)
          IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 3',1)
          fco2_lu_day(1:klon)   = 0.
       END IF
    END IF


! 4) Test for compatibility
! -------------------------
!    IF (carbon_cycle_cpl .AND. type_ocean/='couple') THEN
!       WRITE(lunout,*) 'Coupling with ocean model is needed for carbon_cycle_cpl'
!       CALL abort_physic('carbon_cycle_init', 'coupled ocean is needed for carbon_cycle_cpl',1)
!    END IF
!
!    IF (carbon_cycle_cpl .AND..NOT. ok_veget) THEN
!       WRITE(lunout,*) 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl'
!       CALL abort_physic('carbon_cycle_init', 'ok_veget is needed for carbon_cycle_cpl',1)
!    END IF

    ! Compiler test : following should never happen
    teststop=0
    DO it=1,teststop
       CALL abort_physic('carbon_cycle_init', 'Entering loop from 1 to 0',1)
    END DO

    IF (ntr_co2==0) THEN
       ! No carbon tracers found in tracer.def. It is not possible to do carbon cycle 
       WRITE(lunout,*) 'No carbon tracers found in tracer.def. Not ok with carbon_cycle_tr and/or carbon_cycle_cp'
       CALL abort_physic('carbon_cycle_init', 'No carbon tracers found in tracer.def',1)
    END IF
    
! 5) Calculate total area of the earth surface
! --------------------------------------------
    CALL reduce_sum(SUM(cell_area),airetot)
    CALL bcast(airetot)

  END SUBROUTINE carbon_cycle_init

  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
! Subroutine for injection of co2 in the tracers
!
! - Find out if it is time to update
! - Get tracer from coupled model or from file
! - Calculate new RCO2 value for the radiation scheme
! - Calculate CO2 flux to send to ocean and land models (PISCES and ORCHIDEE)

    USE infotrac_phy, ONLY: nbtr
    USE dimphy
    USE mod_phys_lmdz_transfert_para
    USE phys_cal_mod, ONLY : mth_cur, mth_len
    USE phys_cal_mod, ONLY : day_cur
    USE indice_sol_mod
    USE print_control_mod, ONLY: lunout
    USE geometry_mod, ONLY : cell_area

    IMPLICIT NONE

    INCLUDE "clesphys.h"
    INCLUDE "YOMCST.h"

! In/Output arguments
    INTEGER,INTENT(IN) :: nstep      ! time step in physiq
    REAL,INTENT(IN)    :: pdtphys    ! length of time step in physiq (sec)
    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf            ! Surface fraction
    REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri ! All tracers
    REAL, DIMENSION(klon,nbtr), INTENT(INOUT)       :: source  ! Source for all tracers

! Local variables
    INTEGER :: it
    LOGICAL :: newmonth ! indicates if a new month just started
    LOGICAL :: newday   ! indicates if a new day just started
    LOGICAL :: endday   ! indicated if last time step in a day

    REAL, PARAMETER :: fact=1.E-15/2.12  ! transformation factor from gC/m2/day => ppm/m2/day
    REAL, DIMENSION(klon) :: fco2_tmp
    REAL :: sumtmp
    REAL :: delta_co2_ppm
    

! 1) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
! -------------------------------------------------------------------------------------------------------

    newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.

    IF (MOD(nstep,INT(86400./pdtphys))==1) newday=.TRUE.
    IF (MOD(nstep,INT(86400./pdtphys))==0) endday=.TRUE.
    IF (newday .AND. day_cur==1) newmonth=.TRUE.

! 2)  For each carbon tracer find out if it is time to inject (update)
! --------------------------------------------------------------------
    DO it = 1, ntr_co2
       IF ( MOD(nstep,INT(co2trac(it)%updatefreq/pdtphys)) == 1 ) THEN
          co2trac(it)%updatenow = .TRUE.
       ELSE
          co2trac(it)%updatenow = .FALSE.
       END IF
    END DO

! 3) Get tracer update
! --------------------------------------
    DO it = 1, ntr_co2
       IF ( co2trac(it)%updatenow ) THEN
          IF ( co2trac(it)%cpl ) THEN
             ! Get tracer from coupled model
             SELECT CASE(co2trac(it)%name)
             CASE('fCO2_land')     ! from ORCHIDEE
                dtr_add(:,it) = fco2_land_inst(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
             CASE('fCO2_land_use') ! from ORCHIDEE
                dtr_add(:,it) = fco2_lu_inst(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
             CASE('fCO2_ocn')      ! from PISCES
                dtr_add(:,it) = fco2_ocn_day(:)  *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
             CASE DEFAULT
                WRITE(lunout,*) 'Error with tracer ',co2trac(it)%name
                CALL abort_physic('carbon_cycle', 'No coupling implemented for this tracer',1)
             END SELECT
          ELSE
             ! Read tracer from file
             co2trac(it)%readstep = co2trac(it)%readstep + 1 ! increment time step in file
! Patricia   CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.FALSE.,dtr_add(:,it))
             CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.TRUE.,dtr_add(:,it))

             ! Converte from kgC/m2/h to kgC/m2/s
             dtr_add(:,it) = dtr_add(:,it)/3600
             ! Add individual treatment of values read from file
             SELECT CASE(co2trac(it)%name)
             CASE('fCO2_land')
                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
             CASE('fCO2_land_use')
                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
             CASE('fCO2_ocn')
                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_oce)
! Patricia :
!             CASE('fCO2_fos_fuel')
!                dtr_add(:,it) = dtr_add(:,it)/mth_len
!                co2trac(it)%readstep = 0 ! Always read same value for fossil fuel(Cadule case)
             END SELECT
          END IF
       END IF
    END DO

! 4) Update co2 tracers : 
!    Loop over all carbon tracers and add source
! ------------------------------------------------------------------
    IF (carbon_cycle_tr) THEN
       DO it = 1, ntr_co2
          IF (.FALSE.) THEN
             tr_seri(1:klon,1,co2trac(it)%id) = tr_seri(1:klon,1,co2trac(it)%id) + dtr_add(1:klon,it)
             source(1:klon,co2trac(it)%id) = 0.
          ELSE
             source(1:klon,co2trac(it)%id) = dtr_add(1:klon,it)
          END IF
       END DO
    END IF


! 5) Calculations for new CO2 value for the radiation scheme(instead of reading value from .def)
! ----------------------------------------------------------------------------------------------
    IF (RCO2_inter) THEN
       ! Cumulate fluxes from ORCHIDEE at each timestep
       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
          IF (newday) THEN ! Reset cumulative variables once a day 
             fco2_land_day(1:klon) = 0.
             fco2_lu_day(1:klon)   = 0.
          END IF
          fco2_land_day(1:klon) = fco2_land_day(1:klon) + fco2_land_inst(1:klon) ![gC/m2/day]
          fco2_lu_day(1:klon)   = fco2_lu_day(1:klon)   + fco2_lu_inst(1:klon)   ![gC/m2/day]
       END IF

       ! At the end of a new day, calculate a mean scalare value of CO2
       ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ? 
       IF (endday) THEN

          IF (carbon_cycle_tr) THEN
             ! Sum all co2 tracers to get the total delta CO2 flux
             fco2_tmp(:) = 0.
             DO it = 1, ntr_co2
                fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
             END DO
             
          ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
             ! Sum co2 fluxes comming from coupled models and parameter for fossil fuel
             fco2_tmp(1:klon) = fos_fuel_s + ((fco2_lu_day(1:klon) + fco2_land_day(1:klon))*pctsrf(1:klon,is_ter) &
                  + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
          END IF

          ! Calculate a global mean value of delta CO2 flux
          fco2_tmp(1:klon) = fco2_tmp(1:klon) * cell_area(1:klon)
          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
          CALL bcast(sumtmp)
          delta_co2_ppm = sumtmp/airetot
          
          ! Add initial value for co2_ppm and delta value
          co2_ppm = co2_ppm0 + delta_co2_ppm
          
          ! Transformation of atmospheric CO2 concentration for the radiation code
          RCO2 = co2_ppm * 1.0e-06 * RMCO2 / RMD
          
          WRITE(lunout,*) 'RCO2 is now updated! RCO2 = ', RCO2
       END IF ! endday

    END IF ! RCO2_inter


! 6) Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
! ----------------------------------------------------------------------------
    IF (carbon_cycle_cpl) THEN

       IF (carbon_cycle_tr) THEN
          ! Sum all co2 tracers to get the total delta CO2 flux at first model layer
          fco2_tmp(:) = 0.
          DO it = 1, ntr_co2
             fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
          END DO
          co2_send(1:klon) = fco2_tmp(1:klon) + co2_ppm0
       ELSE
          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
          co2_send(1:klon) = co2_ppm
       END IF

    END IF

  END SUBROUTINE carbon_cycle
  
  SUBROUTINE infocfields_init

!    USE control_mod, ONLY: planet_type
    USE phys_cal_mod, ONLY : mth_cur
    USE mod_synchro_omp
    USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_omp_root
    USE mod_phys_lmdz_transfert_para
    USE mod_phys_lmdz_omp_transfert
    USE dimphy, ONLY: klon

    IMPLICIT NONE

!=======================================================================
!
!   Authors: Patricia Cadule and Laurent Fairhead  
!   -------
!
!  Purpose and description:
!  -----------------------
!
! Infofields
! this routine enables to define the field exchanges in both directions between
! the atmospheric circulation model (LMDZ) and ORCHIDEE. In the future this
! routing might apply to other models (e.g., NEMO, INCA, ...).
! Therefore, currently with this routine, it is possible to define the coupling
! fields only between LMDZ and ORCHIDEE.
! The coupling_fields.def file enables to define the name of the exchanged
! fields at the coupling interface.
! field_in_names : the set of names of the exchanged fields in input to ORCHIDEE
! (LMDZ to ORCHIDEE)
! field_out_names : the set of names of the exchanged fields in output of
! ORCHIDEE (ORCHIDEE to LMDZ)
! n : the number of exchanged fields at th coupling interface
! nb_fields_in : number of inputs fields to ORCHIDEE (LMDZ to ORCHIDEE)
! nb_fields_out : number of ouput fields of ORCHIDEE (ORCHIDEE to LMDZ)
!
! The syntax for coupling_fields.def is as follows:
! IMPORTANT: each column entry must be separated from the previous one by 3
! spaces and only that
! field name         coupling          model 1         model 2         long_name
!                    direction
!   10char  -3spaces-  3char  -3spaces- 4char -3spaces- 4char -3spaces-  30char 
!
! n
! FIELD1 IN LMDZ ORC
! ....
! FIELD(j) IN LMDZ ORC
! FIELD(j+1) OUT LMDZ ORC
! ...
! FIELDn OUT LMDZ ORC   
! 
!=======================================================================
!   ... 22/12/2017 ....
!-----------------------------------------------------------------------
! Declarations

  INCLUDE "clesphys.h"
  INCLUDE "dimensions.h"
  INCLUDE "iniprint.h"

! Local variables

  INTEGER :: iq,  ierr, stat, error

  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), SAVE  :: cfname_root
  CHARACTER(LEN=120), ALLOCATABLE, DIMENSION(:), SAVE :: cftext_root
  CHARACTER(LEN=15), ALLOCATABLE, DIMENSION(:), SAVE  :: 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
  LOGICAL, ALLOCATABLE, DIMENSION(:), SAVE :: mask_out_root

  CHARACTER(len=*),parameter :: modname="infocfields"

  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.GT.0) THEN

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

          IF (ierr.EQ.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)
                WRITE(lunout,*) 'coupling field: ',cfname_root(iq), & 
                               ', number: ',iq,', model 1 (ref): ',cfmod1_root(iq),', model 2: ',cfmod2_root(iq)
                WRITE(lunout,*) 'coupling field: ',cfname_root(iq), & 
                               ', number: ',iq,', long name: ',cftext_root(iq),', units ',cfunits_root(iq)
                IF (nbcf_in+nbcf_out.LT.nbcf) THEN
                  IF (cfintent_root(iq).NE.'OUT') THEN
                    nbcf_in=nbcf_in+1
                    mask_in_root(iq)=.TRUE.
                    mask_out_root(iq)=.FALSE.
                  ELSE IF (cfintent_root(iq).EQ.'OUT') THEN
                    nbcf_out=nbcf_out+1
                    mask_in_root(iq)=.FALSE.
                    mask_out_root(iq)=.TRUE.
                  ENDIF
                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('infocfields_init','Problem in the definition of the coupling fields',1)
               ENDIF
             ENDDO !DO iq=1,nbcf
          ELSE
             WRITE(lunout,*) trim(modname),': infocfields_mod.F90 --- Problem in opening coupling_fields.def'
             WRITE(lunout,*) trim(modname),': infocfields_mod.F90 --- WARNING using defaut values'
          ENDIF ! ierr
          CLOSE(200)

       ENDIF ! level_coupling_esm

    ENDIF !   (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.GT.0)     cfname=cfname_root
        IF (nbcf_in.GT.0)  cfname_in=PACK(cfname_root,mask_in_root)
        IF (nbcf_out.GT.0) cfname_out=PACK(cfname_root,mask_out_root)
        IF (nbcf_in.GT.0)  cftext_in=PACK(cftext_root,mask_in_root)
        IF (nbcf_out.GT.0) cftext_out=PACK(cftext_root,mask_out_root)
        IF (nbcf.GT.0)     cfmod1=cfmod1_root
        IF (nbcf.GT.0)     cfmod2=cfmod2_root
        IF (nbcf_in.GT.0)  cfunits_in=PACK(cfunits_root,mask_in_root)
        IF (nbcf_out.GT.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  ! ALL = ORC/NEMO/INCA
            IF (cfmod1(iq) == "ALL")  nbcf_in_nemo = nbcf_in_nemo + 1  ! ALL = ORC/NEMO/INCA 
            IF (cfmod1(iq) == "ALL")  nbcf_in_inca = nbcf_in_inca + 1  ! ALL = ORC/NEMO/INCA 
            IF (cfmod1(iq) == "ANT")  nbcf_in_ant  = nbcf_in_ant  + 1  
        ENDDO

    ENDIF !   (is_mpi_root .AND. is_omp_root)
!$OMP BARRIER

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

    WRITE(lunout,*) 'nbcf_in_orc  =',nbcf_in_orc
    WRITE(lunout,*) 'nbcf_in_nemo =',nbcf_in_nemo
    WRITE(lunout,*) 'nbcf_in_inca =',nbcf_in_inca
    WRITE(lunout,*) 'nbcf_in_ant  =',nbcf_in_ant

    IF (nbcf_in.GT.0) THEN 
        DO iq=1,nbcf_in
          CALL bcast(cfname_in(iq))
          CALL bcast(cftext_in(iq))
          CALL bcast(cfunits_in(iq))
        ENDDO
    ENDIF

    IF (nbcf_out.GT.0) THEN
        DO iq=1,nbcf_out
          CALL bcast(cfname_out(iq))
          CALL bcast(cftext_out(iq))
          CALL bcast(cfunits_out(iq))
        ENDDO
    ENDIF

    IF (nbcf.GT.0) THEN
        DO iq=1,nbcf
          CALL bcast(cfmod1(iq))
          CALL bcast(cfmod2(iq))
        ENDDO
    ENDIF

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

    IF (nbcf_in.GT.0)  WRITE(lunout,*)'infocfields_mod --- cftext_in: ',cftext_in
    IF (nbcf_out.GT.0) WRITE(lunout,*)'infocfields_mod --- cftext_out: ',cftext_out

    IF (nbcf.GT.0) WRITE(lunout,*)'infocfields_mod --- cfmod1: ',cfmod1 
    IF (nbcf.GT.0) WRITE(lunout,*)'infocfields_mod --- cfmod2: ',cfmod2 

    IF (nbcf_in.GT.0)  WRITE(lunout,*)'infocfunits_mod --- cfunits_in: ',cfunits_in
    IF (nbcf_out.GT.0) WRITE(lunout,*)'infocfunits_mod --- cfunits_out: ',cfunits_out

    IF (nbcf_in.GT.0)  WRITE(*,*)'infocfields_init --- number of fields in to LMDZ: ',nbcf_in
    IF (nbcf_out.GT.0) WRITE(*,*)'infocfields_init --- number of fields out of LMDZ: ',nbcf_out

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

 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