source: LMDZ4/trunk/libf/phylmd/carbon_cycle_mod.F90 @ 3603

MODULE carbon_cycle_mod

! Author : Josefine GHATTAS, Patricia CADULE

  IMPLICIT NONE
  SAVE
  PRIVATE
  PUBLIC :: carbon_cycle_init, carbon_cycle

! Variables read from parmeter file physiq.def
  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_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES) 
!$OMP THREADPRIVATE(carbon_cycle_cpl)
  LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible

! Scalare values when no transport, from physiq.def
  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)

  INTEGER :: ntr_co2                ! Number of tracers concerning the carbon cycle
  INTEGER :: id_fco2_tot            ! Tracer index
  INTEGER :: id_fco2_ocn            !  - " -
  INTEGER :: id_fco2_land           !  - " -
  INTEGER :: id_fco2_land_use       !  - " -
  INTEGER :: id_fco2_fos_fuel       !  - " -
!$OMP THREADPRIVATE(ntr_co2, id_fco2_tot, id_fco2_ocn, id_fco2_land, id_fco2_land_use, id_fco2_fos_fuel)  

  REAL, DIMENSION(:), ALLOCATABLE :: fos_fuel        ! CO2 fossil fuel emission from file [gC/m2/d]
!$OMP THREADPRIVATE(fos_fuel)
  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day ! flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]
!$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)

! Following 2 fields will be initialized in surf_land_orchidee at each time step
  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 
!$OMP THREADPRIVATE(co2_send)

CONTAINS
  
  SUBROUTINE carbon_cycle_init(tr_seri, aerosol, radio)
    USE dimphy
    USE infotrac
    USE IOIPSL
    USE surface_data, ONLY : ok_veget, type_ocean

    IMPLICIT NONE
    INCLUDE "clesphys.h"
    INCLUDE "iniprint.h"
 
! Input argument
    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA]  

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

! Local variables
    INTEGER               :: ierr, it, iiq
    REAL, DIMENSION(klon) :: tr_seri_sum


! 0) Test for compatibility
    IF (carbon_cycle_cpl .AND. type_ocean/='couple') &
         CALL abort_gcm('carbon_cycle_init', 'Coupling with ocean model is needed for carbon_cycle_cpl',1)
    IF (carbon_cycle_cpl .AND..NOT. ok_veget) &
         CALL abort_gcm('carbon_cycle_init', 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl',1)


! 1) Check if transport of one tracer flux CO2 or 4 separated tracers
    IF (carbon_cycle_tr) THEN
       id_fco2_tot=0
       id_fco2_ocn=0
       id_fco2_land=0
       id_fco2_land_use=0
       id_fco2_fos_fuel=0
       
       ! Search in tracer list
       DO it=1,nbtr
          iiq=niadv(it+2)
          IF (tname(iiq) == "fCO2" ) THEN
             id_fco2_tot=it
          ELSE IF (tname(iiq) == "fCO2_ocn" ) THEN
             id_fco2_ocn=it
          ELSE IF (tname(iiq) == "fCO2_land" ) THEN
             id_fco2_land=it
          ELSE IF (tname(iiq) == "fCO2_land_use" ) THEN
             id_fco2_land_use=it
          ELSE IF (tname(iiq) == "fCO2_fos_fuel" ) THEN
             id_fco2_fos_fuel=it
          END IF
       END DO

       ! Count tracers found
       IF (id_fco2_tot /= 0 .AND. &
            id_fco2_ocn==0 .AND. id_fco2_land==0 .AND. id_fco2_land_use==0 .AND. id_fco2_fos_fuel==0) THEN
          
          ! transport  1 tracer flux CO2
          ntr_co2 = 1
          
       ELSE IF (id_fco2_tot==0 .AND. &
            id_fco2_ocn /=0 .AND. id_fco2_land/=0 .AND. id_fco2_land_use/=0 .AND. id_fco2_fos_fuel/=0) THEN
          
          ! transport 4 tracers seperatively
          ntr_co2 = 4
          
       ELSE
          CALL abort_gcm('carbon_cycle_init', 'error in coherence between traceur.def and gcm.def',1)
       END IF

       ! Definition of control varaiables for the tracers
       DO it=1,nbtr
          IF (it==id_fco2_tot .OR. it==id_fco2_ocn .OR. it==id_fco2_land .OR. &
               it==id_fco2_land_use .OR. it==id_fco2_fos_fuel) THEN
             aerosol(it) = .FALSE.
             radio(it)   = .FALSE.
          END IF
       END DO

    ELSE 
       ! No transport of CO2
       ntr_co2 = 0
    END IF ! carbon_cycle_tr


! 2) Allocate variable for CO2 fossil fuel emission
    IF (carbon_cycle_tr) THEN
       ! Allocate 2D variable
       ALLOCATE(fos_fuel(klon), stat=ierr)
       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 1',1)
    ELSE
       ! No transport : read value from .def
       fos_fuel_s = 0.
       CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s)
       WRITE(lunout,*) 'carbon_cycle_fos_fuel = ', fos_fuel_s 
    END IF


! 3) Allocate and initialize fluxes
    IF (carbon_cycle_cpl) THEN
       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
       ALLOCATE(fco2_land_day(klon), stat=ierr)
       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
       ALLOCATE(fco2_lu_day(klon), stat=ierr)
       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 4',1)

       fco2_land_day(:) = 0.  ! JG : Doit prend valeur de restart
       fco2_lu_day(:)   = 0.  ! JG : Doit prend valeur de restart

       ! fco2_ocn_day is allocated in cpl_mod
       ! fco2_land_inst and fco2_lu_inst are allocated in surf_land_orchidee
       
       ALLOCATE(co2_send(klon), stat=ierr)
       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 7',1)
       
       ! Calculate using restart tracer values
       IF (carbon_cycle_tr) THEN
          IF (ntr_co2==1) THEN
             co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
          ELSE ! ntr_co2==4
             ! Calculate the delta CO2 flux
             tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
                  tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
             co2_send(:) = tr_seri_sum(:) + co2_ppm0
          END IF
       ELSE
          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
          co2_send(:) = co2_ppm
       END IF


    ELSE 
       IF (carbon_cycle_tr) THEN
          ! No coupling of CO2 fields : 
          ! corresponding fields will instead be read from files
          ALLOCATE(fco2_ocn_day(klon), stat=ierr)
          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 8',1)
          ALLOCATE(fco2_land_day(klon), stat=ierr)
          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 9',1)
          ALLOCATE(fco2_lu_day(klon), stat=ierr)
          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 10',1)       
       END IF
    END IF

! 4) Read parmeter for calculation of emission compatible
    IF (.NOT. carbon_cycle_tr) THEN
       carbon_cycle_emis_comp=.FALSE.
       CALL getin('carbon_cycle_emis_comp',carbon_cycle_emis_comp)
       WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
    END IF

  END SUBROUTINE carbon_cycle_init

!
!
!

  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
    
    USE infotrac
    USE dimphy
    USE mod_phys_lmdz_transfert_para, ONLY : reduce_sum
    USE phys_cal_mod, ONLY : mth_cur, mth_len
    USE phys_cal_mod, ONLY : day_cur
    USE comgeomphy

    IMPLICIT NONE

    INCLUDE "clesphys.h"
    INCLUDE "indicesol.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 ! Pourcentage de sol f(nature du sol)
    REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri

! Local variables
    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, tr_seri_sum
    REAL :: sumtmp
    REAL :: airetot     ! Total area the earth
    REAL :: delta_co2_ppm
    
! -) 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.
    
! -) Read new maps if new month started
    IF (newmonth .AND. carbon_cycle_tr) THEN
       CALL read_map2D('fossil_fuel.nc','fos_fuel', mth_cur, .FALSE., fos_fuel)
       
       ! division by month lenght to get dayly value
       fos_fuel(:) = fos_fuel(:)/mth_len
       
       IF (.NOT. carbon_cycle_cpl) THEN
          ! Get dayly values from monthly fluxes
          CALL read_map2D('fl_co2_ocean.nc','CO2_OCN',mth_cur,.FALSE.,fco2_ocn_day)
          CALL read_map2D('fl_co2_land.nc','CO2_LAND', mth_cur,.FALSE.,fco2_land_day)
          CALL read_map2D('fl_co2_land_use.nc','CO2_LAND_USE',mth_cur,.FALSE.,fco2_lu_day)
       END IF
    END IF
    


! -) Update tracers at beginning of a new day. Beginning of a new day correspond to a new coupling period in cpl_mod.
    IF (newday) THEN

       IF (carbon_cycle_tr) THEN

          ! Update tracers
          IF (ntr_co2 == 1) THEN
             ! Calculate the new flux CO2
             tr_seri(:,1,id_fco2_tot) = tr_seri(:,1,id_fco2_tot) + &
                  (fos_fuel(:) + &
                  fco2_lu_day(:)  * pctsrf(:,is_ter) + &
                  fco2_land_day(:)* pctsrf(:,is_ter) + &
                  fco2_ocn_day(:) * pctsrf(:,is_oce)) * fact

          ELSE ! ntr_co2 == 4
             tr_seri(:,1,id_fco2_fos_fuel)  = tr_seri(:,1,id_fco2_fos_fuel) + fos_fuel(:) * fact ! [ppm/m2/day]

             tr_seri(:,1,id_fco2_land_use)  = tr_seri(:,1,id_fco2_land_use) + &
                  fco2_lu_day(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]

             tr_seri(:,1,id_fco2_land)      = tr_seri(:,1,id_fco2_land)     + &
                  fco2_land_day(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]

             tr_seri(:,1,id_fco2_ocn)       = tr_seri(:,1,id_fco2_ocn)      + &
                  fco2_ocn_day(:) *pctsrf(:,is_oce)*fact ! [ppm/m2/day]

          END IF

       ELSE ! no transport
          IF (carbon_cycle_cpl) THEN
             IF (carbon_cycle_emis_comp) THEN
                ! Calcul emission compatible a partir des champs 2D et co2_ppm
                !!! TO DO!!
                CALL abort_gcm('carbon_cycle', ' Option carbon_cycle_emis_comp not yet implemented',1)
             END IF
          END IF

       END IF ! carbon_cycle_tr
    
       ! Reset cumluative variables
       IF (carbon_cycle_cpl) THEN
          fco2_land_day(:) = 0.
          fco2_lu_day(:)   = 0.
       END IF
    
    END IF ! newday
       


! -) Cumulate fluxes from ORCHIDEE at each timestep
    IF (carbon_cycle_cpl) THEN
       fco2_land_day(:) = fco2_land_day(:) + fco2_land_inst(:)
       fco2_lu_day(:)   = fco2_lu_day(:)   + fco2_lu_inst(:)
    END IF



! -)  At the end of a new day, calculate a mean scalare value of CO2 to be used by 
!     the radiation scheme (instead of reading value from .def)

! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ? 
    IF (endday) THEN 
       ! Calculte total area of the earth surface
       CALL reduce_sum(SUM(airephy),airetot)
       

       IF (carbon_cycle_tr) THEN
          IF (ntr_co2 == 1) THEN
          
             ! Calculate mean value of tracer CO2 to get an scalare value to be used in the 
             ! radiation scheme (instead of reading value from .def)
             ! Mean value weighted with the grid cell area
             
             ! Calculate mean value
             fco2_tmp(:) = tr_seri(:,1,id_fco2_tot) * airephy(:) 
             CALL reduce_sum(SUM(fco2_tmp),sumtmp)
             co2_ppm = sumtmp/airetot + co2_ppm0
             
          ELSE ! ntr_co2 == 4
             
             ! Calculate the delta CO2 flux
             tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
                  tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
             
             ! Calculate mean value of delta CO2 flux
             fco2_tmp(:) = tr_seri_sum(:) * airephy(:)
             CALL reduce_sum(SUM(fco2_tmp),sumtmp)
             delta_co2_ppm = sumtmp/airetot
             
             ! Add initial value for co2_ppm to delta value
             co2_ppm = delta_co2_ppm + co2_ppm0
          END IF

       ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr

          ! Calculate the total CO2 flux and integrate to get scalare value for the radiation scheme
          fco2_tmp(:) = (fos_fuel(:) + (fco2_lu_day(:) + fco2_land_day(:))*pctsrf(:,is_ter) &
               + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
          
          ! Calculate mean value
          fco2_tmp(:) = fco2_tmp(:) * airephy(:)
          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
          delta_co2_ppm = sumtmp/airetot

          ! Update current value of the atmospheric co2_ppm
          co2_ppm = co2_ppm + delta_co2_ppm
          
       END IF ! carbon_cycle_tr

       ! transformation of the atmospheric CO2 concentration for the radiation code
       RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97 

    END IF

    ! Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
    IF (endday .AND. carbon_cycle_cpl) THEN
       
       IF (carbon_cycle_tr) THEN
          IF (ntr_co2==1) THEN

             co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0

          ELSE ! ntr_co2 == 4

             co2_send(:) = tr_seri_sum(:) + co2_ppm0

          END IF
       ELSE
          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
          co2_send(:) = co2_ppm
       END IF

    END IF

  END SUBROUTINE carbon_cycle
  
END MODULE carbon_cycle_mod