Index: LMDZ6/trunk/libf/phylmd/carbon_cycle_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/carbon_cycle_mod.F90	(revision 3579)
+++ LMDZ6/trunk/libf/phylmd/carbon_cycle_mod.F90	(revision 3581)
@@ -7,11 +7,8 @@
 !  -----------------------
 ! 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)
+!   - Initialisation of carbon cycle fields
+!   - Definition of fluxes to be exchanged
+!
+! Rest of code is in tracco2i.F90
 !
 ! Le cas online/offline est defini par le flag carbon_cycle_cpl (y/n)
@@ -35,5 +32,5 @@
   SAVE
   PRIVATE
-  PUBLIC :: carbon_cycle_init, carbon_cycle, infocfields_init
+  PUBLIC :: carbon_cycle_init, infocfields_init
 
 ! Variables read from parmeter file physiq.def
@@ -46,6 +43,8 @@
   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)
+  REAL, PUBLIC :: RCO2_glo
+!$OMP THREADPRIVATE(RCO2_glo)
+  REAL, PUBLIC :: RCO2_tot
+!$OMP THREADPRIVATE(RCO2_tot)
 
   LOGICAL :: carbon_cycle_emis_comp_omp=.FALSE.
@@ -84,4 +83,14 @@
   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
 !$OMP THREADPRIVATE(fco2_land)
+  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_nbp  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
+!$OMP THREADPRIVATE(fco2_land_nbp)
+  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_nep  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
+!$OMP THREADPRIVATE(fco2_land_nep)
+  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fLuc  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
+!$OMP THREADPRIVATE(fco2_land_fLuc)
+  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fwoodharvest  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
+!$OMP THREADPRIVATE(fco2_land_fwoodharvest)
+  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fHarvest  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
+!$OMP THREADPRIVATE(fco2_land_fHarvest)
   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocean ! Net flux from ocean [kgCO2/m2/s]
 !$OMP THREADPRIVATE(fco2_ocean)
@@ -191,5 +200,5 @@
 CONTAINS
   
-  SUBROUTINE carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
+  SUBROUTINE carbon_cycle_init()
 ! This subroutine is called from traclmdz_init, only at first timestep.
 ! - Read controle parameters from .def input file
@@ -199,10 +208,5 @@
 
     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
 
@@ -210,384 +214,50 @@
     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)
+    INTEGER               :: ierr
+
+    IF (carbon_cycle_cpl) THEN
+
+       ierr=0
+
+       IF (.NOT.ALLOCATED(fco2_land)) ALLOCATE(fco2_land(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land',1)
+       fco2_land(1:klon) = 0.
+
+       IF (.NOT.ALLOCATED(fco2_land_nbp)) ALLOCATE(fco2_land_nbp(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_nbp',1)
+       fco2_land_nbp(1:klon) = 0.
+
+       IF (.NOT.ALLOCATED(fco2_land_nep)) ALLOCATE(fco2_land_nep(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_nep',1)
+       fco2_land_nep(1:klon) = 0.
+
+       IF (.NOT.ALLOCATED(fco2_land_fLuc)) ALLOCATE(fco2_land_fLuc(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_fLuc',1)
+       fco2_land_fLuc(1:klon) = 0.
+
+       IF (.NOT.ALLOCATED(fco2_land_fwoodharvest)) ALLOCATE(fco2_land_fwoodharvest(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_fwoodharvest',1)
+       fco2_land_fwoodharvest(1:klon) = 0.
+
+       IF (.NOT.ALLOCATED(fco2_land_fHarvest)) ALLOCATE(fco2_land_fHarvest(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_fHarvest',1)
+       fco2_land_fHarvest(1:klon) = 0.
+
+       IF (.NOT.ALLOCATED(fco2_ff)) ALLOCATE(fco2_ff(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_ff',1)
+       fco2_ff(1:klon) = 0.
+
+       IF (.NOT.ALLOCATED(fco2_bb)) ALLOCATE(fco2_bb(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_bb',1)
+       fco2_bb(1:klon) = 0.
+
+       IF (.NOT.ALLOCATED(fco2_ocean)) ALLOCATE(fco2_ocean(klon), stat=ierr)
+       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_ocean',1)
+       fco2_bb(1:klon) = 0.
+    ENDIF
 
   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
 
Index: LMDZ6/trunk/libf/phylmd/phyetat0.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 3579)
+++ LMDZ6/trunk/libf/phylmd/phyetat0.F90	(revision 3581)
@@ -114,6 +114,8 @@
   ENDIF
 
-  ! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
-  co2_ppm0   = tab_cntrl(16)
+! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
+! co2_ppm0   = tab_cntrl(16)
+! initial value for interactive CO2 run when there is no tracer field for CO2 in restart
+  co2_ppm0=284.32
 
   solaire_etat0      = tab_cntrl(4)
Index: LMDZ6/trunk/libf/phylmd/tracco2i_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/tracco2i_mod.F90	(revision 3579)
+++ LMDZ6/trunk/libf/phylmd/tracco2i_mod.F90	(revision 3581)
@@ -2,4 +2,13 @@
 !
 ! 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
 !
 CONTAINS
@@ -12,7 +21,10 @@
     USE infotrac_phy
     USE geometry_mod, ONLY: cell_area
+    USE carbon_cycle_mod, ONLY: carbon_cycle_init
     USE carbon_cycle_mod, ONLY: id_CO2, nbcf_in, fields_in, cfname_in
     USE carbon_cycle_mod, ONLY: fco2_ocn_day, fco2_ff, fco2_bb, fco2_land, fco2_ocean
-    USE carbon_cycle_mod, ONLY: carbon_cycle_tr, carbon_cycle_rad, RCO2_glo, RCO2_tot
+    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, carbon_cycle_rad, RCO2_glo, RCO2_tot
     USE mod_grid_phy_lmdz
     USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
@@ -54,6 +66,8 @@
     REAL, DIMENSION(klon_glo,klev) :: m_air_glo ! variable temporaire sur la grille global
 
-    INTEGER, SAVE :: mth_pre=0, day_pre=0
-!$OMP THREADPRIVATE(mth_pre, day_pre)
+    LOGICAL, SAVE :: check_fCO2_nbp_in_cfname
+!$OMP THREADPRIVATE(check_fCO2_nbp_in_cfname)
+    INTEGER, SAVE :: day_pre=-1
+!$OMP THREADPRIVATE(day_pre)
 
     IF (is_mpi_root) THEN
@@ -66,13 +80,19 @@
     IF (debutphy) THEN
 
+! Initialisation de module carbon_cycle_mod
+      IF (carbon_cycle_cpl) THEN
+        CALL carbon_cycle_init()
+      ENDIF
+
+! Initialisation de tr_seri(id_CO2) si pas initialise
       IF (MAXVAL(tr_seri(:,:,id_CO2)).LT.1.e-15) THEN 
-        !!tr_seri(:,:,id_CO2)=280.e-6/RMD*RMCO2
-        tr_seri(:,:,id_CO2)=400.e-6/RMD*RMCO2 !--initialised to 400 ppm for a test
+        tr_seri(:,:,id_CO2)=co2_ppm0*1.e-6/RMD*RMCO2 !--initialised from co2_ppm0 in rdem
       ENDIF 
 
-      ALLOCATE(fco2_ff(klon))
-      ALLOCATE(fco2_bb(klon))
-      ALLOCATE(fco2_land(klon))
-      ALLOCATE(fco2_ocean(klon))
+!--check if fCO2_nbp is in 
+      check_fCO2_nbp_in_cfname=.FALSE.
+      DO nb=1, nbcf_in
+        IF (cfname_in(nb)=="fCO2_nbp") check_fCO2_nbp_in_cfname=.TRUE.
+      ENDDO
 
     ENDIF
@@ -93,24 +113,63 @@
     fco2_land(:)=0.0
     fco2_ocean(:)=0.0
+    fco2_land_nbp(:)=0.
+    fco2_land_nep(:)=0.
+    fco2_land_fLuc(:)=0.
+    fco2_land_fwoodharvest(:)=0.
+    fco2_land_fHarvest(:)=0.
+
     DO nb=1, nbcf_in
-      print *,'nb tracco2=', nb, cfname_in(nb)
-!--fCO2_nep comes in unit of kg C m-2 s-1
-!--converting to kg CO2 m-2 s-1
-      IF (cfname_in(nb) == "fCO2_nbp" )   fco2_land(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
-!--fCO2_fgco2 comes in unit of mol C02 m-2 s-1
-!--converting to kg CO2 m-2 s-1 + change sign
-      IF (cfname_in(nb) == "fCO2_fgco2" ) fco2_ocean(:)=-1.*fco2_ocn_day(:)*RMCO2/1.e3*(pctsrf(:,is_oce)+pctsrf(:,is_sic))
+
+      SELECT CASE(cfname_in(nb))
+!--dealing with the different fluxes coming from ORCHIDEE
+!--fluxes come in unit of kg C m-2 s-1 is converted into kg CO2 m-2 s-1
+      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)
+!--fCO2_fco2_ocn comes in unit of mol C02 m-2 s-1 is converted into kg CO2 m-2 s-1 + change sign
+      CASE("fCO2_fgco2") 
+          fco2_ocean(:)=-1.*fco2_ocn_day(:)*RMCO2/1.e3*(pctsrf(:,is_oce)+pctsrf(:,is_sic))
+      END SELECT
+
     ENDDO
 
-!--preparing the net anthropogenic flux at the surface for mixing layer
-!--unit kg CO2 / m2 / s
+!--if fCO2_nbp is transferred we use it, otherwise we use the sum of what has been passed from ORCHIDEE
+    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(:)
+    ENDIF
+
+!!--preparing the net anthropogenic flux at the surface for mixing layer
+!!--unit kg CO2 / m2 / s
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_ff) ',MAXVAL(fco2_ff)
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_ff) ',MINVAL(fco2_ff)
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_bb) ',MAXVAL(fco2_bb)
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_bb) ',MINVAL(fco2_bb)
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_land) ',MAXVAL(fco2_land)
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_land) ',MINVAL(fco2_land)
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_ocean) ',MAXVAL(fco2_ocean)
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_ocean) ',MINVAL(fco2_ocean)
+!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(source(:,id_CO2)) ',MAXVAL(source(:,id_CO2))
+!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(source(:,id_CO2)) ',MINVAL(source(:,id_CO2))
+!
+!--build final source term for CO2
     source(:,id_CO2)=fco2_ff(:)+fco2_bb(:)+fco2_land(:)+fco2_ocean(:)
 
 !--computing global mean CO2 for radiation
-!--every timestep for now but enough every day
-!    IF (debutphy.OR.mth_cur.NE.mth_pre) THEN
-!    IF (debutphy.OR.day_cur.NE.day_pre) THEN
+!--for every timestep comment out the IF ENDIF statements 
+!--otherwise this is updated every day
+    IF (debutphy.OR.day_cur.NE.day_pre) THEN
+
       CALL gather(tr_seri(:,:,id_CO2),co2_glo)
       CALL gather(m_air,m_air_glo)
+
 !$OMP MASTER
 
@@ -124,5 +183,4 @@
 !$OMP END MASTER
        CALL bcast(RCO2_glo)
-       mth_pre=mth_cur
        day_pre=day_cur
 !--if not carbon_cycle_tr, then we reinitialize the CO2 each day to its global mean value
@@ -130,5 +188,5 @@
          tr_seri(:,:,id_CO2)=RCO2_glo
        ENDIF
-!    ENDIF
+    ENDIF
 
   END SUBROUTINE tracco2i
@@ -168,5 +226,7 @@
 !! may be controlled via the .def later on
 !! also co2bb for now comes from ORCHIDEE
-    LOGICAL, PARAMETER :: readco2ff=.TRUE., readco2bb=.FALSE. 
+    LOGICAL, PARAMETER :: readco2ff=.TRUE.
+!! this should be left to FALSE for now
+    LOGICAL, PARAMETER :: readco2bb=.FALSE. 
 
     CHARACTER (len = 20) :: modname = 'tracco2i.co2_emissions'
@@ -217,4 +277,5 @@
 
 !--reading CO2 biomass burning emissions
+!--using it will be inconsistent with treatment in ORCHIDEE
       IF (readco2bb) THEN 
 
Index: LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90	(revision 3579)
+++ LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90	(revision 3581)
@@ -92,5 +92,4 @@
     USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
     USE press_coefoz_m, ONLY: press_coefoz
-    USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl
     USE mod_grid_phy_lmdz
     USE mod_phys_lmdz_para
@@ -285,10 +284,4 @@
 
 !
-! Initialisation de module carbon_cycle_mod
-! ----------------------------------------------
-    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
-       CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
-    END IF
-
 ! Check if all tracers have restart values
 ! ----------------------------------------------
@@ -346,5 +339,4 @@
     USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
     USE o3_chem_m, ONLY: o3_chem
-    USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
     USE indice_sol_mod
 
@@ -612,11 +604,4 @@
     END IF
 
-!======================================================================
-!   Calcul de cycle de carbon
-!======================================================================
-    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
-       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
-    END IF
-
   END SUBROUTINE traclmdz
 
