MODULE tracco2i_mod

!

! This module does the work for the interactive CO2 tracers

! Authors: Patricia Cadule and Olivier Boucher

!

! Purpose and description:

!  -----------------------

! Main routine for the interactive carbon cycle

! Gather all carbon fluxes and emissions from ORCHIDEE, PISCES and fossil fuel

! Compute the net flux in source field which is used in phytrac

! Compute global CO2 mixing ratio for radiation scheme if option is activated

! Redistribute CO2 evenly over the atmosphere if transport is desactivated

!

CONTAINS

    ! This subroutine calls carbon_cycle_init needed to be done before first call to phys_output_write in physiq.

    USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl

    ! Initialize carbon_cycle_mod

    ENDIF

       xlat, xlon, pphis, pphi, &

    USE dimphy

    USE infotrac_phy, ONLY: nbtr

    USE geometry_mod, ONLY: cell_area

    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: read_fco2_ocean_cor,var_fco2_ocean_cor,fco2_ocean_cor

    USE carbon_cycle_mod, ONLY: read_fco2_land_cor,var_fco2_land_cor,fco2_land_cor

    USE carbon_cycle_mod, ONLY: co2_send

    USE carbon_cycle_mod, ONLY: fco2_land_nbp, fco2_land_nep, fco2_land_fLuc

    USE carbon_cycle_mod, ONLY: fco2_land_fwoodharvest, fco2_land_fHarvest

    USE carbon_cycle_mod, ONLY: carbon_cycle_cpl, carbon_cycle_tr, carbon_cycle_rad, RCO2_glo, RCO2_tot

    USE carbon_cycle_mod, ONLY: ocean_area_tot

    USE carbon_cycle_mod, ONLY: land_area_tot

    USE mod_grid_phy_lmdz

    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root

    USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter

    USE mod_phys_lmdz_omp_data, ONLY: is_omp_root

    USE phys_cal_mod

    USE phys_state_var_mod, ONLY: pctsrf

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

    IMPLICIT NONE

    INCLUDE "clesphys.h"

    INCLUDE "YOMCST.h"

! Input argument

!---------------

    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)

    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique

    REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point

    REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point

    REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! geopotentiel du sol

    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel de chaque couche

    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature

    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)

    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)

    REAL,DIMENSION(klon,nbtr),INTENT(INOUT):: source  ! flux de traceur [U/m2/s]

! Output argument

!----------------

    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/kgA]

! Local variables

!----------------

    INTEGER                        :: it, k, i, nb

    LOGICAL, SAVE :: check_fCO2_nbp_in_cfname

!$OMP THREADPRIVATE(check_fCO2_nbp_in_cfname)

    INTEGER, SAVE :: day_pre=-1

!$OMP THREADPRIVATE(day_pre)

    REAL, PARAMETER :: secinday=86400.

    ENDIF

!--initialisation of CO2 field if not in restart file

!--dirty way of doing, do it better later

!--convert 280 ppm into kg CO2 / kg air

! Initialization of tr_seri(id_CO2) If it is not initialized

      ENDIF

!--check if fCO2_nbp is in

      ENDDO

    ENDIF

!--calculate mass of air in every grid box in kg air

    ENDDO

    ENDDO

!--call CO2 emission routine

!--co2bb is zero for now

!--unit kg CO2 m-2 s-1

!--retrieving land and ocean CO2 flux

!--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")

      CASE("fCO2_fLuc")

      CASE("fCO2_fwoodharvest")

      CASE("fCO2_fHarvest")

      CASE("fCO2_nbp")

!--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")

      END SELECT

    ENDDO

!$OMP MASTER

          ENDDO

      ENDIF !--is_mpi_root and is_omp_root

!$OMP END MASTER

    ENDIF

!$OMP MASTER

          ENDDO

      ENDIF !--is_mpi_root and is_omp_root

!$OMP END MASTER

ENDIF

    ENDIF !-- debutphy

! var_fco2_ocean_cor: correction of the surface downward CO2 flux into the ocean fgco2 (PgC/yr)

! This is the correction of the the net air to ocean carbon flux. Positive flux is into the ocean.

!    PRINT *, 'tracco2i_mod.F90 --- var_fco2_ocean_cor (PgC/yr) ',var_fco2_ocean_cor

!var_fco2_ocean_cor: correction of the net air to ocean carbon flux (input data is a scalar in PgC/yr and must be converted in kg CO2 m-2 s-1)

! Factors for carbon and carbon dioxide

! 1 mole CO2 = 44.009 g CO2 = 12.011 g C

! 1 ppm by volume of atmosphere CO2 = 2.13 Gt C

! 1 gC = 44.009/12.011 gCO2

! ocean_area_tot: ocean area (m2)

! year_len: year length (in days)

! conversion: PgC/yr --> kg CO2 m-2 s-1

! fco2_ocean_cor  / (86400.*year_len): PgC/yr to PgC/s

! fco2_ocean_cor  / (86400.*year_len)*(pctsrf(i,is_oce)+pctsrf(i,is_sic))/ocean_area_tot: PgC/s to PgC/s/m2

! (fco2_ocean_cor / (86400.*year_len)*(pctsrf(i,is_oce)+pctsrf(i,is_sic))/ocean_area_tot) *1e12: PgC/s/m2 to kgC/s/m2

! (fco2_ocean_cor / (86400.*year_len)*(pctsrf(i,is_oce)+pctsrf(i,is_sic))/ocean_area_tot) * 1e12 * (RMCO2/RMC): kgC/s/m2 to kgCO2/s/m2

         fco2_ocean_cor(i)=(var_fco2_ocean_cor*(RMCO2/RMC) &

              *(pctsrf(i,is_oce)+pctsrf(i,is_sic))/ocean_area_tot &

      ENDDO

    ELSE

    ENDIF

! var_fco2_land_cor: correction of the carbon Mass Flux out of Atmosphere Due to Net Biospheric Production on Land  (PgC/yr)

! This is the correction of the net mass flux of carbon between land and atmosphere calculated as

! photosynthesis MINUS the sum of plant and soil respiration, carbon fluxes from

! fire, harvest, grazing and land use change. Positive flux is into the land.

!    PRINT *, 'tracco2i_mod.F90 --- var_fco2_land_cor (m2) ',var_fco2_land_cor

!var_fco2_land_cor: correction of the et air to land carbon flux (input data is a scalar in PgC/yr and must be converted in kg CO2 m-2 s-1)

! Factors for carbon and carbon dioxide

! 1 mole CO2 = 44.009 g CO2 = 12.011 g C

! 1 ppm by volume of atmosphere CO2 = 2.13 Gt C

! 1 gC = 44.009/12.011 gCO2

! land_area_tot: land area (m2)

! year_len: year length (in days)

! conversion: PgC/yr --> kg CO2 m-2 s-1

! fco2_land_cor  / (86400.*year_len): PgC/yr to PgC/s

! fco2_land_cor  / (86400.*year_len)*pctsrf(i,is_ter)/land_area_tot: PgC/s to PgC/s/m2

! (fco2_land_cor / (86400.*year_len)*pctsrf(i,is_ter)/land_area_tot) *1e12: PgC/s/m2 to kgC/s/m2

! (fco2_land_cor / (86400.*year_len)*pctsrf(i,is_ter)/land_area_tot) * 1e12 * (RMCO2/RMC): kgC/s/m2 to kgCO2/s/m2

      ENDDO

    ELSE

    ENDIF

!--if fCO2_nbp is transferred we use it, otherwise we use the sum of what has been passed from ORCHIDEE

    ELSE

    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

!--computing global mean CO2 for radiation

!--for every timestep comment out the IF ENDIF statements

!--otherwise this is updated every day

!$OMP MASTER

!--compute a global mean CO2 value and print its value in ppm

         ! the following operation is only to maintain precision consistency

         ! of RCO2_glo which differs whether it is directly computed or read from

         ! a restart file (after having been computed)

       ENDIF

!$OMP END MASTER

!--if not carbon_cycle_tr, then we reinitialize the CO2 each day to its global mean value

       ENDIF

    ENDIF

!    USE infotrac_phy

    USE geometry_mod, ONLY : cell_area

    USE mod_grid_phy_lmdz

    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root

    USE mod_phys_lmdz_para, ONLY: gather, scatter

    USE phys_cal_mod

    USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_varid, nf95_open

    USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite

    USE carbon_cycle_mod, ONLY : fco2_ff, fco2_bb, fco2_land, fco2_ocean

    IMPLICIT NONE

    INCLUDE "YOMCST.h"

    LOGICAL,INTENT(IN) :: debutphy

! For NetCDF:

    INTEGER ncid_in  ! IDs for input files

    INTEGER varid, ncerr

    INTEGER :: n_glo, n_month

    REAL,ALLOCATABLE, SAVE :: flx_co2ff(:,:)     !  fossil-fuel CO2

    REAL,ALLOCATABLE, SAVE :: flx_co2bb(:,:)     !  biomass-burning CO2

!$OMP THREADPRIVATE(flx_co2ff,flx_co2bb)

!! may be controlled via the .def later on

!! also co2bb for now comes from ORCHIDEE

    LOGICAL, PARAMETER :: readco2ff=.TRUE.

!! this should be left to FALSE for now

    LOGICAL, PARAMETER :: readco2bb=.FALSE.

    CHARACTER (len = 20) :: modname = 'tracco2i.co2_emissions'

    CHARACTER (len = 80) :: abort_message

!$OMP MASTER

!--reading CO2 fossil fuel emissions

        ! ... Open the CO2ff file

        CALL nf95_gw_var(ncid_in, varid, vector)

        ENDIF

        CALL nf95_gw_var(ncid_in, varid, time)

        ENDIF

!--reading flx_co2 for fossil fuel

      ELSE  !--co2ff not to be read

        flx_co2ff_glo(:,:)=0.0

      ENDIF

!--reading CO2 biomass burning emissions

!--using it will be inconsistent with treatment in ORCHIDEE

      IF (readco2bb) THEN

      ! ... Open the CO2bb file

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

      CALL nf95_inq_varid(ncid_in, "vector", varid)

      CALL nf95_gw_var(ncid_in, varid, vector)

      n_glo = size(vector)

      IF (n_glo.NE.klon_glo) THEN

         abort_message='sflx_lmdz_co2_bb: le nombre de points n est pas egal a klon_glo'

         CALL abort_physic(modname,abort_message,1)

      ENDIF

      CALL nf95_inq_varid(ncid_in, "time", varid)

      CALL nf95_gw_var(ncid_in, varid, time)

      n_month = size(time)

      IF (n_month.NE.12) THEN

         abort_message='sflx_lmdz_co2_bb: le nombre de month n est pas egal a 12'

         CALL abort_physic(modname,abort_message,1)

      ENDIF

!--reading flx_co2 for biomass burning

      CALL nf95_inq_varid(ncid_in, "flx_co2", varid)

      ncerr = nf90_get_var(ncid_in, varid, flx_co2bb_glo)

      CALL nf95_close(ncid_in)

      ELSE  !--co2bb not to be read

      ENDIF

    ENDIF

!$OMP END MASTER

    ! Allocation needed for all proc otherwise scatter might complain

!--scatter on all proc

  ENDIF !--end debuthy

!---select the correct month

  ENDIF

END MODULE tracco2i_mod