Ignore:
Timestamp:
Oct 21, 2010, 4:25:20 PM (14 years ago)
Author:
jghattas
Message:

Modifications for carbon tracers :

  • add possibility to read source at different frequency
  • add dynamic varaiable to send fluxes in interface between ORCHIDEE and LMDZ : this is still under developpment under cpp key ORCH_NEW. No compatible official ORCHIDEE version does yet exist.
  • add variable co2_send in restart file.
  • clean up and bug fixes.


File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/branches/LMDZ5V1.0-dev/libf/phylmd/carbon_cycle_mod.F90

    r1279 r1444  
    11MODULE carbon_cycle_mod
    2 
     2! Controle module for the carbon CO2 tracers :
     3!   - Identification
     4!   - Get concentrations comming from coupled model or read from file to tracers
     5!   - Calculate new RCO2 for radiation scheme
     6!   - Calculate new carbon flux for sending to coupled models (PISCES and ORCHIDEE)
     7!
    38! Author : Josefine GHATTAS, Patricia CADULE
    49
     
    1318  LOGICAL, PUBLIC :: carbon_cycle_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES)
    1419!$OMP THREADPRIVATE(carbon_cycle_cpl)
     20
    1521  LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible
     22!$OMP THREADPRIVATE(carbon_cycle_emis_comp)
     23
     24  LOGICAL :: RCO2_inter  ! RCO2 interactive : if true calculate new value RCO2 for the radiation scheme
     25!$OMP THREADPRIVATE(RCO2_inter)
    1626
    1727! Scalare values when no transport, from physiq.def
     
    2131!$OMP THREADPRIVATE(emis_land_s)
    2232
    23   INTEGER :: ntr_co2                ! Number of tracers concerning the carbon cycle
    24   INTEGER :: id_fco2_tot            ! Tracer index
    25   INTEGER :: id_fco2_ocn            !  - " -
    26   INTEGER :: id_fco2_land           !  - " -
    27   INTEGER :: id_fco2_land_use       !  - " -
    28   INTEGER :: id_fco2_fos_fuel       !  - " -
    29 !$OMP THREADPRIVATE(ntr_co2, id_fco2_tot, id_fco2_ocn, id_fco2_land, id_fco2_land_use, id_fco2_fos_fuel) 
    30 
    31   REAL, DIMENSION(:), ALLOCATABLE :: fos_fuel        ! CO2 fossil fuel emission from file [gC/m2/d]
    32 !$OMP THREADPRIVATE(fos_fuel)
    33   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day ! flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]
     33  REAL :: airetot     ! Total area of the earth surface
     34!$OMP THREADPRIVATE(airetot)
     35
     36  INTEGER :: ntr_co2  ! Number of tracers concerning the carbon cycle
     37!$OMP THREADPRIVATE(ntr_co2)
     38
     39! fco2_ocn_day : flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]. Allocation and initalization done in cpl_mod
     40  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day
    3441!$OMP THREADPRIVATE(fco2_ocn_day)
     42
    3543  REAL, DIMENSION(:), ALLOCATABLE :: fco2_land_day   ! flux CO2 from land for 1 day (cumulated)  [gC/m2/d]
    3644!$OMP THREADPRIVATE(fco2_land_day)
     
    3846!$OMP THREADPRIVATE(fco2_lu_day)
    3947
    40 ! Following 2 fields will be initialized in surf_land_orchidee at each time step
     48  REAL, DIMENSION(:,:), ALLOCATABLE :: dtr_add       ! Tracer concentration to be injected
     49!$OMP THREADPRIVATE(dtr_add)
     50
     51! Following 2 fields will be allocated and initialized in surf_land_orchidee
    4152  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_inst  ! flux CO2 from land at one time step
    4253!$OMP THREADPRIVATE(fco2_land_inst)
     
    4556
    4657! Calculated co2 field to be send to the ocean via the coupler and to ORCHIDEE
    47   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send
     58  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send ! Field allocated in phyetat0
    4859!$OMP THREADPRIVATE(co2_send)
     60
     61
     62  TYPE, PUBLIC ::   co2_trac_type
     63     CHARACTER(len = 8) :: name       ! Tracer name in tracer.def
     64     INTEGER            :: id         ! Index in total tracer list, tr_seri
     65     CHARACTER(len=30)  :: file       ! File name
     66     LOGICAL            :: cpl        ! True if this tracers is coupled from ORCHIDEE or PISCES.
     67                                      ! False if read from file.
     68     INTEGER            :: updatefreq ! Frequence to inject in second
     69     INTEGER            :: readstep   ! Actual time step to read in file
     70     LOGICAL            :: updatenow  ! True if this tracer should be updated this time step
     71  END TYPE co2_trac_type
     72  INTEGER,PARAMETER :: maxco2trac=5  ! Maximum number of different CO2 fluxes
     73  TYPE(co2_trac_type), DIMENSION(maxco2trac) :: co2trac
    4974
    5075CONTAINS
    5176 
    52   SUBROUTINE carbon_cycle_init(tr_seri, aerosol, radio)
     77  SUBROUTINE carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
     78! This subroutine is called from traclmdz_init, only at first timestep.
     79! - Read controle parameters from .def input file
     80! - Search for carbon tracers and set default values
     81! - Allocate variables
     82! - Test for compatibility
     83
    5384    USE dimphy
     85    USE comgeomphy
     86    USE mod_phys_lmdz_transfert_para
    5487    USE infotrac
    5588    USE IOIPSL
    5689    USE surface_data, ONLY : ok_veget, type_ocean
     90    USE phys_cal_mod, ONLY : mth_len
    5791
    5892    IMPLICIT NONE
     
    6296! Input argument
    6397    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA] 
     98    REAL,INTENT(IN)                           :: pdtphys ! length of time step in physiq (sec)
    6499
    65100! InOutput arguments
     
    68103
    69104! Local variables
    70     INTEGER               :: ierr, it, iiq
    71     REAL, DIMENSION(klon) :: tr_seri_sum
    72 
    73 
    74 ! 0) Test for compatibility
    75     IF (carbon_cycle_cpl .AND. type_ocean/='couple') &
    76          CALL abort_gcm('carbon_cycle_init', 'Coupling with ocean model is needed for carbon_cycle_cpl',1)
    77     IF (carbon_cycle_cpl .AND..NOT. ok_veget) &
    78          CALL abort_gcm('carbon_cycle_init', 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl',1)
    79 
    80 
    81 ! 1) Check if transport of one tracer flux CO2 or 4 separated tracers
    82     IF (carbon_cycle_tr) THEN
    83        id_fco2_tot=0
    84        id_fco2_ocn=0
    85        id_fco2_land=0
    86        id_fco2_land_use=0
    87        id_fco2_fos_fuel=0
    88        
    89        ! Search in tracer list
    90        DO it=1,nbtr
    91           iiq=niadv(it+2)
    92           IF (tname(iiq) == "fCO2" ) THEN
    93              id_fco2_tot=it
    94           ELSE IF (tname(iiq) == "fCO2_ocn" ) THEN
    95              id_fco2_ocn=it
    96           ELSE IF (tname(iiq) == "fCO2_land" ) THEN
    97              id_fco2_land=it
    98           ELSE IF (tname(iiq) == "fCO2_land_use" ) THEN
    99              id_fco2_land_use=it
    100           ELSE IF (tname(iiq) == "fCO2_fos_fuel" ) THEN
    101              id_fco2_fos_fuel=it
    102           END IF
    103        END DO
    104 
    105        ! Count tracers found
    106        IF (id_fco2_tot /= 0 .AND. &
    107             id_fco2_ocn==0 .AND. id_fco2_land==0 .AND. id_fco2_land_use==0 .AND. id_fco2_fos_fuel==0) THEN
    108          
    109           ! transport  1 tracer flux CO2
    110           ntr_co2 = 1
    111          
    112        ELSE IF (id_fco2_tot==0 .AND. &
    113             id_fco2_ocn /=0 .AND. id_fco2_land/=0 .AND. id_fco2_land_use/=0 .AND. id_fco2_fos_fuel/=0) THEN
    114          
    115           ! transport 4 tracers seperatively
    116           ntr_co2 = 4
    117          
    118        ELSE
    119           CALL abort_gcm('carbon_cycle_init', 'error in coherence between traceur.def and gcm.def',1)
    120        END IF
    121 
    122        ! Definition of control varaiables for the tracers
    123        DO it=1,nbtr
    124           IF (it==id_fco2_tot .OR. it==id_fco2_ocn .OR. it==id_fco2_land .OR. &
    125                it==id_fco2_land_use .OR. it==id_fco2_fos_fuel) THEN
    126              aerosol(it) = .FALSE.
    127              radio(it)   = .FALSE.
    128           END IF
    129        END DO
    130 
    131     ELSE
    132        ! No transport of CO2
    133        ntr_co2 = 0
    134     END IF ! carbon_cycle_tr
    135 
    136 
    137 ! 2) Allocate variable for CO2 fossil fuel emission
    138     IF (carbon_cycle_tr) THEN
    139        ! Allocate 2D variable
    140        ALLOCATE(fos_fuel(klon), stat=ierr)
    141        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 1',1)
    142     ELSE
    143        ! No transport : read value from .def
     105    INTEGER               :: ierr, it, iiq, itc
     106    INTEGER               :: teststop
     107
     108
     109
     110! 1) Read controle parameters from .def input file
     111! ------------------------------------------------
     112    ! Read fosil fuel value if no transport
     113    IF (.NOT. carbon_cycle_tr) THEN
    144114       fos_fuel_s = 0.
    145115       CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s)
     
    148118
    149119
    150 ! 3) Allocate and initialize fluxes
    151     IF (carbon_cycle_cpl) THEN
    152        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
    153        ALLOCATE(fco2_land_day(klon), stat=ierr)
    154        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
    155        ALLOCATE(fco2_lu_day(klon), stat=ierr)
    156        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 4',1)
    157 
    158        fco2_land_day(:) = 0.  ! JG : Doit prend valeur de restart
    159        fco2_lu_day(:)   = 0.  ! JG : Doit prend valeur de restart
    160 
    161        ! fco2_ocn_day is allocated in cpl_mod
    162        ! fco2_land_inst and fco2_lu_inst are allocated in surf_land_orchidee
    163        
    164        ALLOCATE(co2_send(klon), stat=ierr)
    165        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 7',1)
    166        
    167        ! Calculate using restart tracer values
    168        IF (carbon_cycle_tr) THEN
    169           IF (ntr_co2==1) THEN
    170              co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
    171           ELSE ! ntr_co2==4
    172              ! Calculate the delta CO2 flux
    173              tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
    174                   tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
    175              co2_send(:) = tr_seri_sum(:) + co2_ppm0
    176           END IF
    177        ELSE
    178           ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
    179           co2_send(:) = co2_ppm
    180        END IF
    181 
    182 
    183     ELSE
    184        IF (carbon_cycle_tr) THEN
    185           ! No coupling of CO2 fields :
    186           ! corresponding fields will instead be read from files
    187           ALLOCATE(fco2_ocn_day(klon), stat=ierr)
    188           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 8',1)
    189           ALLOCATE(fco2_land_day(klon), stat=ierr)
    190           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 9',1)
    191           ALLOCATE(fco2_lu_day(klon), stat=ierr)
    192           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 10',1)       
    193        END IF
    194     END IF
    195 
    196 ! 4) Read parmeter for calculation of emission compatible
     120    ! Read parmeter for calculation compatible emission
    197121    IF (.NOT. carbon_cycle_tr) THEN
    198122       carbon_cycle_emis_comp=.FALSE.
    199123       CALL getin('carbon_cycle_emis_comp',carbon_cycle_emis_comp)
    200124       WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
    201     END IF
     125       IF (carbon_cycle_emis_comp) THEN
     126          CALL abort_gcm('carbon_cycle_init', 'carbon_cycle_emis_comp option not yet implemented!!',1)
     127       END IF
     128    END IF
     129
     130    ! Read parameter for interactive calculation of the CO2 value for the radiation scheme
     131    RCO2_inter=.FALSE.
     132    CALL getin('RCO2_inter',RCO2_inter)
     133    WRITE(lunout,*) 'RCO2_inter = ', RCO2_inter
     134    IF (RCO2_inter) THEN
     135       WRITE(lunout,*) 'RCO2 will be recalculated once a day'
     136       WRITE(lunout,*) 'RCO2 initial = ', RCO2
     137    END IF
     138
     139
     140! 2) Search for carbon tracers and set default values
     141! ---------------------------------------------------
     142    itc=0
     143    DO it=1,nbtr
     144       iiq=niadv(it+2)
     145       
     146       SELECT CASE(tname(iiq))
     147       CASE("fCO2_ocn")
     148          itc = itc + 1
     149          co2trac(itc)%name='fCO2_ocn'
     150          co2trac(itc)%id=it
     151          co2trac(itc)%file='fl_co2_ocean.nc'
     152          IF (carbon_cycle_cpl .AND. type_ocean=='couple') THEN
     153             co2trac(itc)%cpl=.TRUE.
     154             co2trac(itc)%updatefreq = 86400 ! Once a day as the coupling with OASIS/PISCES
     155          ELSE
     156             co2trac(itc)%cpl=.FALSE.
     157             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     158          END IF
     159       CASE("fCO2_land")
     160          itc = itc + 1
     161          co2trac(itc)%name='fCO2_land'
     162          co2trac(itc)%id=it
     163          co2trac(itc)%file='fl_co2_land.nc'
     164          IF (carbon_cycle_cpl .AND. ok_veget) THEN
     165             co2trac(itc)%cpl=.TRUE.
     166             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
     167          ELSE
     168             co2trac(itc)%cpl=.FALSE.
     169!             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
     170             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     171          END IF
     172       CASE("fCO2_land_use")
     173          itc = itc + 1
     174          co2trac(itc)%name='fCO2_land_use'
     175          co2trac(itc)%id=it
     176          co2trac(itc)%file='fl_co2_land_use.nc'
     177          IF (carbon_cycle_cpl .AND. ok_veget) THEN
     178             co2trac(it)%cpl=.TRUE.
     179             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
     180          ELSE
     181             co2trac(itc)%cpl=.FALSE.
     182             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
     183          END IF
     184       CASE("fCO2_fos_fuel")
     185          itc = itc + 1
     186          co2trac(itc)%name='fCO2_fos_fuel'
     187          co2trac(itc)%id=it
     188          co2trac(itc)%file='fossil_fuel.nc'
     189          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
     190!         co2trac(itc)%updatefreq = 86400  ! 86400sec = 24H Cadule case
     191          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     192       CASE("fCO2_bbg")
     193          itc = itc + 1
     194          co2trac(itc)%name='fCO2_bbg'
     195          co2trac(itc)%id=it
     196          co2trac(itc)%file='fl_co2_bbg.nc'
     197          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
     198          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     199       CASE("fCO2")
     200          ! fCO2 : One tracer transporting the total CO2 flux
     201          itc = itc + 1
     202          co2trac(itc)%name='fCO2'
     203          co2trac(itc)%id=it
     204          co2trac(itc)%file='fl_co2.nc'
     205          IF (carbon_cycle_cpl) THEN
     206             co2trac(itc)%cpl=.TRUE.
     207          ELSE
     208             co2trac(itc)%cpl=.FALSE.
     209          END IF
     210          co2trac(itc)%updatefreq = 86400
     211          ! DOES THIS WORK ???? Problematic due to implementation of the coupled fluxes...
     212          CALL abort_gcm('carbon_cycle_init','transport of total CO2 has to be implemented and tested',1)
     213       END SELECT
     214    END DO
     215
     216    ! Total number of carbon CO2 tracers
     217    ntr_co2 = itc
     218   
     219    ! Definition of control varaiables for the tracers
     220    DO it=1,ntr_co2
     221       aerosol(co2trac(it)%id) = .FALSE.
     222       radio(co2trac(it)%id)   = .FALSE.
     223    END DO
     224   
     225    ! Vector indicating which timestep to read for each tracer
     226    ! Always start read in the beginning of the file
     227    co2trac(:)%readstep = 0
     228   
     229
     230! 3) Allocate variables
     231! ---------------------
     232    ! Allocate vector for storing fluxes to inject
     233    ALLOCATE(dtr_add(klon,maxco2trac), stat=ierr)
     234    IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 11',1)       
     235   
     236    ! Allocate variables for cumulating fluxes from ORCHIDEE
     237    IF (RCO2_inter) THEN
     238       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
     239          ALLOCATE(fco2_land_day(klon), stat=ierr)
     240          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
     241          fco2_land_day(1:klon) = 0.
     242         
     243          ALLOCATE(fco2_lu_day(klon), stat=ierr)
     244          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
     245          fco2_lu_day(1:klon)   = 0.
     246       END IF
     247    END IF
     248
     249
     250! 4) Test for compatibility
     251! -------------------------
     252!    IF (carbon_cycle_cpl .AND. type_ocean/='couple') THEN
     253!       WRITE(lunout,*) 'Coupling with ocean model is needed for carbon_cycle_cpl'
     254!       CALL abort_gcm('carbon_cycle_init', 'coupled ocean is needed for carbon_cycle_cpl',1)
     255!    END IF
     256!
     257!    IF (carbon_cycle_cpl .AND..NOT. ok_veget) THEN
     258!       WRITE(lunout,*) 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl'
     259!       CALL abort_gcm('carbon_cycle_init', 'ok_veget is needed for carbon_cycle_cpl',1)
     260!    END IF
     261
     262    ! Compiler test : following should never happen
     263    teststop=0
     264    DO it=1,teststop
     265       CALL abort_gcm('carbon_cycle_init', 'Entering loop from 1 to 0',1)
     266    END DO
     267
     268    IF (ntr_co2==0) THEN
     269       ! No carbon tracers found in tracer.def. It is not possible to do carbon cycle
     270       WRITE(lunout,*) 'No carbon tracers found in tracer.def. Not ok with carbon_cycle_tr and/or carbon_cycle_cp'
     271       CALL abort_gcm('carbon_cycle_init', 'No carbon tracers found in tracer.def',1)
     272    END IF
     273   
     274! 5) Calculate total area of the earth surface
     275! --------------------------------------------
     276    CALL reduce_sum(SUM(airephy),airetot)
     277    CALL bcast(airetot)
    202278
    203279  END SUBROUTINE carbon_cycle_init
    204280
     281
     282  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
     283! Subroutine for injection of co2 in the tracers
    205284!
    206 !
    207 !
    208 
    209   SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
    210    
     285! - Find out if it is time to update
     286! - Get tracer from coupled model or from file
     287! - Calculate new RCO2 value for the radiation scheme
     288! - Calculate CO2 flux to send to ocean and land models (PISCES and ORCHIDEE)
     289
    211290    USE infotrac
    212291    USE dimphy
    213     USE mod_phys_lmdz_transfert_para, ONLY : reduce_sum
     292    USE mod_phys_lmdz_transfert_para
    214293    USE phys_cal_mod, ONLY : mth_cur, mth_len
    215294    USE phys_cal_mod, ONLY : day_cur
     
    220299    INCLUDE "clesphys.h"
    221300    INCLUDE "indicesol.h"
     301    INCLUDE "iniprint.h"
     302    INCLUDE "YOMCST.h"
    222303
    223304! In/Output arguments
    224305    INTEGER,INTENT(IN) :: nstep      ! time step in physiq
    225306    REAL,INTENT(IN)    :: pdtphys    ! length of time step in physiq (sec)
    226     REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
    227     REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri
     307    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf            ! Surface fraction
     308    REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri ! All tracers
     309    REAL, DIMENSION(klon,nbtr), INTENT(INOUT)       :: source  ! Source for all tracers
    228310
    229311! Local variables
     312    INTEGER :: it
    230313    LOGICAL :: newmonth ! indicates if a new month just started
    231314    LOGICAL :: newday   ! indicates if a new day just started
     
    233316
    234317    REAL, PARAMETER :: fact=1.E-15/2.12  ! transformation factor from gC/m2/day => ppm/m2/day
    235     REAL, DIMENSION(klon) :: fco2_tmp, tr_seri_sum
     318    REAL, DIMENSION(klon) :: fco2_tmp
    236319    REAL :: sumtmp
    237     REAL :: airetot     ! Total area the earth
    238320    REAL :: delta_co2_ppm
    239321   
    240 ! -) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
     322
     323! 1) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
     324! -------------------------------------------------------------------------------------------------------
    241325
    242326    newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.
     
    245329    IF (MOD(nstep,INT(86400./pdtphys))==0) endday=.TRUE.
    246330    IF (newday .AND. day_cur==1) newmonth=.TRUE.
    247    
    248 ! -) Read new maps if new month started
    249     IF (newmonth .AND. carbon_cycle_tr) THEN
    250        CALL read_map2D('fossil_fuel.nc','fos_fuel', mth_cur, .FALSE., fos_fuel)
    251        
    252        ! division by month lenght to get dayly value
    253        fos_fuel(:) = fos_fuel(:)/mth_len
    254        
    255        IF (.NOT. carbon_cycle_cpl) THEN
    256           ! Get dayly values from monthly fluxes
    257           CALL read_map2D('fl_co2_ocean.nc','CO2_OCN',mth_cur,.FALSE.,fco2_ocn_day)
    258           CALL read_map2D('fl_co2_land.nc','CO2_LAND', mth_cur,.FALSE.,fco2_land_day)
    259           CALL read_map2D('fl_co2_land_use.nc','CO2_LAND_USE',mth_cur,.FALSE.,fco2_lu_day)
    260        END IF
    261     END IF
    262    
    263 
    264 
    265 ! -) Update tracers at beginning of a new day. Beginning of a new day correspond to a new coupling period in cpl_mod.
    266     IF (newday) THEN
     331
     332! 2)  For each carbon tracer find out if it is time to inject (update)
     333! --------------------------------------------------------------------
     334    DO it = 1, ntr_co2
     335       IF ( MOD(nstep,INT(co2trac(it)%updatefreq/pdtphys)) == 1 ) THEN
     336          co2trac(it)%updatenow = .TRUE.
     337       ELSE
     338          co2trac(it)%updatenow = .FALSE.
     339       END IF
     340    END DO
     341
     342! 3) Get tracer update
     343! --------------------------------------
     344    DO it = 1, ntr_co2
     345       IF ( co2trac(it)%updatenow ) THEN
     346          IF ( co2trac(it)%cpl ) THEN
     347             ! Get tracer from coupled model
     348             SELECT CASE(co2trac(it)%name)
     349             CASE('fCO2_land')     ! from ORCHIDEE
     350                dtr_add(:,it) = fco2_land_inst(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
     351             CASE('fCO2_land_use') ! from ORCHIDEE
     352                dtr_add(:,it) = fco2_lu_inst(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
     353             CASE('fCO2_ocn')      ! from PISCES
     354                dtr_add(:,it) = fco2_ocn_day(:)  *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
     355             CASE DEFAULT
     356                WRITE(lunout,*) 'Error with tracer ',co2trac(it)%name
     357                CALL abort_gcm('carbon_cycle', 'No coupling implemented for this tracer',1)
     358             END SELECT
     359          ELSE
     360             ! Read tracer from file
     361             co2trac(it)%readstep = co2trac(it)%readstep + 1 ! increment time step in file
     362! Patricia   CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.FALSE.,dtr_add(:,it))
     363             CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.TRUE.,dtr_add(:,it))
     364
     365             ! Converte from kgC/m2/h to kgC/m2/s
     366             dtr_add(:,it) = dtr_add(:,it)/3600
     367             ! Add individual treatment of values read from file
     368             SELECT CASE(co2trac(it)%name)
     369             CASE('fCO2_land')
     370                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
     371             CASE('fCO2_land_use')
     372                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
     373             CASE('fCO2_ocn')
     374                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_oce)
     375! Patricia :
     376!             CASE('fCO2_fos_fuel')
     377!                dtr_add(:,it) = dtr_add(:,it)/mth_len
     378!                co2trac(it)%readstep = 0 ! Always read same value for fossil fuel(Cadule case)
     379             END SELECT
     380          END IF
     381       END IF
     382    END DO
     383
     384! 4) Update co2 tracers :
     385!    Loop over all carbon tracers and add source
     386! ------------------------------------------------------------------
     387    IF (carbon_cycle_tr) THEN
     388       DO it = 1, ntr_co2
     389          IF (.FALSE.) THEN
     390             tr_seri(1:klon,1,co2trac(it)%id) = tr_seri(1:klon,1,co2trac(it)%id) + dtr_add(1:klon,it)
     391             source(1:klon,co2trac(it)%id) = 0.
     392          ELSE
     393             source(1:klon,co2trac(it)%id) = dtr_add(1:klon,it)
     394          END IF
     395       END DO
     396    END IF
     397
     398
     399! 5) Calculations for new CO2 value for the radiation scheme(instead of reading value from .def)
     400! ----------------------------------------------------------------------------------------------
     401    IF (RCO2_inter) THEN
     402       ! Cumulate fluxes from ORCHIDEE at each timestep
     403       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
     404          IF (newday) THEN ! Reset cumulative variables once a day
     405             fco2_land_day(1:klon) = 0.
     406             fco2_lu_day(1:klon)   = 0.
     407          END IF
     408          fco2_land_day(1:klon) = fco2_land_day(1:klon) + fco2_land_inst(1:klon) ![gC/m2/day]
     409          fco2_lu_day(1:klon)   = fco2_lu_day(1:klon)   + fco2_lu_inst(1:klon)   ![gC/m2/day]
     410       END IF
     411
     412       ! At the end of a new day, calculate a mean scalare value of CO2
     413       ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
     414       IF (endday) THEN
     415
     416          IF (carbon_cycle_tr) THEN
     417             ! Sum all co2 tracers to get the total delta CO2 flux
     418             fco2_tmp(:) = 0.
     419             DO it = 1, ntr_co2
     420                fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
     421             END DO
     422             
     423          ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
     424             ! Sum co2 fluxes comming from coupled models and parameter for fossil fuel
     425             fco2_tmp(1:klon) = fos_fuel_s + ((fco2_lu_day(1:klon) + fco2_land_day(1:klon))*pctsrf(1:klon,is_ter) &
     426                  + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
     427          END IF
     428
     429          ! Calculate a global mean value of delta CO2 flux
     430          fco2_tmp(1:klon) = fco2_tmp(1:klon) * airephy(1:klon)
     431          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
     432          CALL bcast(sumtmp)
     433          delta_co2_ppm = sumtmp/airetot
     434         
     435          ! Add initial value for co2_ppm and delta value
     436          co2_ppm = co2_ppm0 + delta_co2_ppm
     437         
     438          ! Transformation of atmospheric CO2 concentration for the radiation code
     439          RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
     440         
     441          WRITE(lunout,*) 'RCO2 is now updated! RCO2 = ', RCO2
     442       END IF ! endday
     443
     444    END IF ! RCO2_inter
     445
     446
     447! 6) Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
     448! ----------------------------------------------------------------------------
     449    IF (carbon_cycle_cpl) THEN
    267450
    268451       IF (carbon_cycle_tr) THEN
    269 
    270           ! Update tracers
    271           IF (ntr_co2 == 1) THEN
    272              ! Calculate the new flux CO2
    273              tr_seri(:,1,id_fco2_tot) = tr_seri(:,1,id_fco2_tot) + &
    274                   (fos_fuel(:) + &
    275                   fco2_lu_day(:)  * pctsrf(:,is_ter) + &
    276                   fco2_land_day(:)* pctsrf(:,is_ter) + &
    277                   fco2_ocn_day(:) * pctsrf(:,is_oce)) * fact
    278 
    279           ELSE ! ntr_co2 == 4
    280              tr_seri(:,1,id_fco2_fos_fuel)  = tr_seri(:,1,id_fco2_fos_fuel) + fos_fuel(:) * fact ! [ppm/m2/day]
    281 
    282              tr_seri(:,1,id_fco2_land_use)  = tr_seri(:,1,id_fco2_land_use) + &
    283                   fco2_lu_day(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
    284 
    285              tr_seri(:,1,id_fco2_land)      = tr_seri(:,1,id_fco2_land)     + &
    286                   fco2_land_day(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
    287 
    288              tr_seri(:,1,id_fco2_ocn)       = tr_seri(:,1,id_fco2_ocn)      + &
    289                   fco2_ocn_day(:) *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
    290 
    291           END IF
    292 
    293        ELSE ! no transport
    294           IF (carbon_cycle_cpl) THEN
    295              IF (carbon_cycle_emis_comp) THEN
    296                 ! Calcul emission compatible a partir des champs 2D et co2_ppm
    297                 !!! TO DO!!
    298                 CALL abort_gcm('carbon_cycle', ' Option carbon_cycle_emis_comp not yet implemented',1)
    299              END IF
    300           END IF
    301 
    302        END IF ! carbon_cycle_tr
    303    
    304        ! Reset cumluative variables
    305        IF (carbon_cycle_cpl) THEN
    306           fco2_land_day(:) = 0.
    307           fco2_lu_day(:)   = 0.
    308        END IF
    309    
    310     END IF ! newday
    311        
    312 
    313 
    314 ! -) Cumulate fluxes from ORCHIDEE at each timestep
    315     IF (carbon_cycle_cpl) THEN
    316        fco2_land_day(:) = fco2_land_day(:) + fco2_land_inst(:)
    317        fco2_lu_day(:)   = fco2_lu_day(:)   + fco2_lu_inst(:)
    318     END IF
    319 
    320 
    321 
    322 ! -)  At the end of a new day, calculate a mean scalare value of CO2 to be used by
    323 !     the radiation scheme (instead of reading value from .def)
    324 
    325 ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
    326     IF (endday) THEN
    327        ! Calculte total area of the earth surface
    328        CALL reduce_sum(SUM(airephy),airetot)
    329        
    330 
    331        IF (carbon_cycle_tr) THEN
    332           IF (ntr_co2 == 1) THEN
    333          
    334              ! Calculate mean value of tracer CO2 to get an scalare value to be used in the
    335              ! radiation scheme (instead of reading value from .def)
    336              ! Mean value weighted with the grid cell area
    337              
    338              ! Calculate mean value
    339              fco2_tmp(:) = tr_seri(:,1,id_fco2_tot) * airephy(:)
    340              CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    341              co2_ppm = sumtmp/airetot + co2_ppm0
    342              
    343           ELSE ! ntr_co2 == 4
    344              
    345              ! Calculate the delta CO2 flux
    346              tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
    347                   tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
    348              
    349              ! Calculate mean value of delta CO2 flux
    350              fco2_tmp(:) = tr_seri_sum(:) * airephy(:)
    351              CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    352              delta_co2_ppm = sumtmp/airetot
    353              
    354              ! Add initial value for co2_ppm to delta value
    355              co2_ppm = delta_co2_ppm + co2_ppm0
    356           END IF
    357 
    358        ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
    359 
    360           ! Calculate the total CO2 flux and integrate to get scalare value for the radiation scheme
    361           fco2_tmp(:) = (fos_fuel(:) + (fco2_lu_day(:) + fco2_land_day(:))*pctsrf(:,is_ter) &
    362                + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
    363          
    364           ! Calculate mean value
    365           fco2_tmp(:) = fco2_tmp(:) * airephy(:)
    366           CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    367           delta_co2_ppm = sumtmp/airetot
    368 
    369           ! Update current value of the atmospheric co2_ppm
    370           co2_ppm = co2_ppm + delta_co2_ppm
    371          
    372        END IF ! carbon_cycle_tr
    373 
    374        ! transformation of the atmospheric CO2 concentration for the radiation code
    375        RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
    376 
    377     END IF
    378 
    379     ! Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
    380     IF (endday .AND. carbon_cycle_cpl) THEN
    381        
    382        IF (carbon_cycle_tr) THEN
    383           IF (ntr_co2==1) THEN
    384 
    385              co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
    386 
    387           ELSE ! ntr_co2 == 4
    388 
    389              co2_send(:) = tr_seri_sum(:) + co2_ppm0
    390 
    391           END IF
     452          ! Sum all co2 tracers to get the total delta CO2 flux at first model layer
     453          fco2_tmp(:) = 0.
     454          DO it = 1, ntr_co2
     455             fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
     456          END DO
     457          co2_send(1:klon) = fco2_tmp(1:klon) + co2_ppm0
    392458       ELSE
    393459          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
    394           co2_send(:) = co2_ppm
     460          co2_send(1:klon) = co2_ppm
    395461       END IF
    396462
Note: See TracChangeset for help on using the changeset viewer.