Ignore:
Timestamp:
Oct 10, 2019, 2:35:59 PM (5 years ago)
Author:
oboucher
Message:

Big update to the interactive carbon cycle
from Patricia's code

Location:
LMDZ6/trunk/libf/phylmd
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/carbon_cycle_mod.F90

    r3459 r3581  
    77!  -----------------------
    88! Control module for the carbon CO2 tracers :
    9 !   - Identification
    10 !   - Get concentrations comming from coupled model or read from file to tracers
    11 !   - Calculate new RCO2 for radiation scheme
    12 !   - Calculate new carbon flux for sending to coupled models (PISCES and ORCHIDEE)
    13 !
    14 ! Module permettant de mettre a jour les champs (puits et sources) pour le
    15 ! transport de CO2 en online (IPSL-CM et LMDZOR) et offline (lecture de carte)
     9!   - Initialisation of carbon cycle fields
     10!   - Definition of fluxes to be exchanged
     11!
     12! Rest of code is in tracco2i.F90
    1613!
    1714! Le cas online/offline est defini par le flag carbon_cycle_cpl (y/n)
    … …  
    3532  SAVE
    3633  PRIVATE
    37   PUBLIC :: carbon_cycle_init, carbon_cycle, infocfields_init
     34  PUBLIC :: carbon_cycle_init, infocfields_init
    3835
    3936! Variables read from parmeter file physiq.def
    … …  
    4643  INTEGER, PUBLIC :: level_coupling_esm ! Level of coupling for the ESM - 0, 1, 2, 3
    4744!$OMP THREADPRIVATE(level_coupling_esm)
    48   REAL, PUBLIC :: RCO2_glo, RCO2_tot
    49 !$OMP THREADPRIVATE(RCO2_glo, RCO2_tot)
     45  REAL, PUBLIC :: RCO2_glo
     46!$OMP THREADPRIVATE(RCO2_glo)
     47  REAL, PUBLIC :: RCO2_tot
     48!$OMP THREADPRIVATE(RCO2_tot)
    5049
    5150  LOGICAL :: carbon_cycle_emis_comp_omp=.FALSE.
    … …  
    8483  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
    8584!$OMP THREADPRIVATE(fco2_land)
     85  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_nbp  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
     86!$OMP THREADPRIVATE(fco2_land_nbp)
     87  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_nep  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
     88!$OMP THREADPRIVATE(fco2_land_nep)
     89  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fLuc  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
     90!$OMP THREADPRIVATE(fco2_land_fLuc)
     91  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fwoodharvest  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
     92!$OMP THREADPRIVATE(fco2_land_fwoodharvest)
     93  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_fHarvest  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
     94!$OMP THREADPRIVATE(fco2_land_fHarvest)
    8695  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocean ! Net flux from ocean [kgCO2/m2/s]
    8796!$OMP THREADPRIVATE(fco2_ocean)
    … …  
    191200CONTAINS
    192201 
    193   SUBROUTINE carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
     202  SUBROUTINE carbon_cycle_init()
    194203! This subroutine is called from traclmdz_init, only at first timestep.
    195204! - Read controle parameters from .def input file
    … …  
    199208
    200209    USE dimphy
    201     USE geometry_mod, ONLY : cell_area
    202     USE mod_phys_lmdz_transfert_para
    203     USE infotrac_phy, ONLY: nbtr, nqo, niadv, tname
    204210    USE IOIPSL
    205     USE surface_data, ONLY : ok_veget, type_ocean
    206     USE phys_cal_mod, ONLY : mth_len
    207211    USE print_control_mod, ONLY: lunout
    208212
    … …  
    210214    INCLUDE "clesphys.h"
    211215 
    212 ! Input argument
    213     REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA] 
    214     REAL,INTENT(IN)                           :: pdtphys ! length of time step in physiq (sec)
    215 
    216 ! InOutput arguments
    217     LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: aerosol
    218     LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: radio
    219 
    220216! Local variables
    221     INTEGER               :: ierr, it, iiq, itc
    222     INTEGER               :: teststop
    223 
    224 ! 1) Read controle parameters from .def input file
    225 ! ------------------------------------------------
    226     ! Read fosil fuel value if no transport
    227     IF (.NOT. carbon_cycle_tr) THEN
    228 !$OMP MASTER
    229        fos_fuel_s_omp = 0.
    230        CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s_omp)
    231 !$OMP END MASTER
    232 !$OMP BARRIER
    233        fos_fuel_s=fos_fuel_s_omp
    234        WRITE(lunout,*) 'carbon_cycle_fos_fuel = ', fos_fuel_s
    235     END IF
    236 
    237     ! Read parmeter for calculation compatible emission
    238     IF (.NOT. carbon_cycle_tr) THEN
    239 !$OMP MASTER
    240        carbon_cycle_emis_comp_omp=.FALSE.
    241        CALL getin('carbon_cycle_emis_comp',carbon_cycle_emis_comp_omp)
    242 !$OMP END MASTER
    243 !$OMP BARRIER
    244        carbon_cycle_emis_comp=carbon_cycle_emis_comp_omp
    245        WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
    246        IF (carbon_cycle_emis_comp) THEN
    247           CALL abort_physic('carbon_cycle_init', 'carbon_cycle_emis_comp option not yet implemented!!',1)
    248        END IF
    249     END IF
    250 
    251     ! Read parameter for interactive calculation of the CO2 value for the radiation scheme
    252 !$OMP MASTER
    253     RCO2_inter_omp=.FALSE.
    254     CALL getin('RCO2_inter',RCO2_inter_omp)
    255 !$OMP END MASTER
    256 !$OMP BARRIER
    257     RCO2_inter=RCO2_inter_omp
    258     WRITE(lunout,*) 'RCO2_inter = ', RCO2_inter
    259     IF (RCO2_inter) THEN
    260        WRITE(lunout,*) 'RCO2 will be recalculated once a day'
    261        WRITE(lunout,*) 'RCO2 initial = ', RCO2
    262     END IF
    263 
    264 
    265 ! 2) Search for carbon tracers and set default values
    266 ! ---------------------------------------------------
    267     itc=0
    268     DO it=1,nbtr
    269 !!       iiq=niadv(it+2)                                                            ! jyg
    270        iiq=niadv(it+nqo)                                                            ! jyg
    271        
    272        SELECT CASE(tname(iiq))
    273        CASE("fCO2_ocn")
    274           itc = itc + 1
    275           co2trac(itc)%name='fCO2_ocn'
    276           co2trac(itc)%id=it
    277           co2trac(itc)%file='fl_co2_ocean.nc'
    278           IF (carbon_cycle_cpl .AND. type_ocean=='couple') THEN
    279              co2trac(itc)%cpl=.TRUE.
    280              co2trac(itc)%updatefreq = 86400 ! Once a day as the coupling with OASIS/PISCES
    281           ELSE
    282              co2trac(itc)%cpl=.FALSE.
    283              co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
    284           END IF
    285        CASE("fCO2_land")
    286           itc = itc + 1
    287           co2trac(itc)%name='fCO2_land'
    288           co2trac(itc)%id=it
    289           co2trac(itc)%file='fl_co2_land.nc'
    290           IF (carbon_cycle_cpl .AND. ok_veget) THEN
    291              co2trac(itc)%cpl=.TRUE.
    292              co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
    293           ELSE
    294              co2trac(itc)%cpl=.FALSE.
    295 !             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
    296              co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
    297           END IF
    298        CASE("fCO2_land_use")
    299           itc = itc + 1
    300           co2trac(itc)%name='fCO2_land_use'
    301           co2trac(itc)%id=it
    302           co2trac(itc)%file='fl_co2_land_use.nc'
    303           IF (carbon_cycle_cpl .AND. ok_veget) THEN
    304              co2trac(it)%cpl=.TRUE.
    305              co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
    306           ELSE
    307              co2trac(itc)%cpl=.FALSE.
    308              co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
    309           END IF
    310        CASE("fCO2_fos_fuel")
    311           itc = itc + 1
    312           co2trac(itc)%name='fCO2_fos_fuel'
    313           co2trac(itc)%id=it
    314           co2trac(itc)%file='fossil_fuel.nc'
    315           co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
    316 !         co2trac(itc)%updatefreq = 86400  ! 86400sec = 24H Cadule case
    317           co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
    318        CASE("fCO2_bbg")
    319           itc = itc + 1
    320           co2trac(itc)%name='fCO2_bbg'
    321           co2trac(itc)%id=it
    322           co2trac(itc)%file='fl_co2_bbg.nc'
    323           co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
    324           co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
    325        CASE("fCO2")
    326           ! fCO2 : One tracer transporting the total CO2 flux
    327           itc = itc + 1
    328           co2trac(itc)%name='fCO2'
    329           co2trac(itc)%id=it
    330           co2trac(itc)%file='fl_co2.nc'
    331           IF (carbon_cycle_cpl) THEN
    332              co2trac(itc)%cpl=.TRUE.
    333           ELSE
    334              co2trac(itc)%cpl=.FALSE.
    335           END IF
    336           co2trac(itc)%updatefreq = 86400
    337           ! DOES THIS WORK ???? Problematic due to implementation of the coupled fluxes...
    338           CALL abort_physic('carbon_cycle_init','transport of total CO2 has to be implemented and tested',1)
    339        END SELECT
    340     END DO
    341 
    342     ! Total number of carbon CO2 tracers
    343     ntr_co2 = itc
    344    
    345     ! Definition of control varaiables for the tracers
    346     DO it=1,ntr_co2
    347        aerosol(co2trac(it)%id) = .FALSE.
    348        radio(co2trac(it)%id)   = .FALSE.
    349     END DO
    350    
    351     ! Vector indicating which timestep to read for each tracer
    352     ! Always start read in the beginning of the file
    353     co2trac(:)%readstep = 0
    354    
    355 
    356 ! 3) Allocate variables
    357 ! ---------------------
    358     ! Allocate vector for storing fluxes to inject
    359     ALLOCATE(dtr_add(klon,maxco2trac), stat=ierr)
    360     IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 11',1)       
    361    
    362     ! Allocate variables for cumulating fluxes from ORCHIDEE
    363     IF (RCO2_inter) THEN
    364        IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
    365           ALLOCATE(fco2_land_day(klon), stat=ierr)
    366           IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 2',1)
    367           fco2_land_day(1:klon) = 0.
    368          
    369           ALLOCATE(fco2_lu_day(klon), stat=ierr)
    370           IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 3',1)
    371           fco2_lu_day(1:klon)   = 0.
    372        END IF
    373     END IF
    374 
    375 
    376 ! 4) Test for compatibility
    377 ! -------------------------
    378 !    IF (carbon_cycle_cpl .AND. type_ocean/='couple') THEN
    379 !       WRITE(lunout,*) 'Coupling with ocean model is needed for carbon_cycle_cpl'
    380 !       CALL abort_physic('carbon_cycle_init', 'coupled ocean is needed for carbon_cycle_cpl',1)
    381 !    END IF
    382 !
    383 !    IF (carbon_cycle_cpl .AND..NOT. ok_veget) THEN
    384 !       WRITE(lunout,*) 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl'
    385 !       CALL abort_physic('carbon_cycle_init', 'ok_veget is needed for carbon_cycle_cpl',1)
    386 !    END IF
    387 
    388     ! Compiler test : following should never happen
    389     teststop=0
    390     DO it=1,teststop
    391        CALL abort_physic('carbon_cycle_init', 'Entering loop from 1 to 0',1)
    392     END DO
    393 
    394     IF (ntr_co2==0) THEN
    395        ! No carbon tracers found in tracer.def. It is not possible to do carbon cycle
    396        WRITE(lunout,*) 'No carbon tracers found in tracer.def. Not ok with carbon_cycle_tr and/or carbon_cycle_cp'
    397        CALL abort_physic('carbon_cycle_init', 'No carbon tracers found in tracer.def',1)
    398     END IF
    399    
    400 ! 5) Calculate total area of the earth surface
    401 ! --------------------------------------------
    402     CALL reduce_sum(SUM(cell_area),airetot)
    403     CALL bcast(airetot)
     217    INTEGER               :: ierr
     218
     219    IF (carbon_cycle_cpl) THEN
     220
     221       ierr=0
     222
     223       IF (.NOT.ALLOCATED(fco2_land)) ALLOCATE(fco2_land(klon), stat=ierr)
     224       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land',1)
     225       fco2_land(1:klon) = 0.
     226
     227       IF (.NOT.ALLOCATED(fco2_land_nbp)) ALLOCATE(fco2_land_nbp(klon), stat=ierr)
     228       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_nbp',1)
     229       fco2_land_nbp(1:klon) = 0.
     230
     231       IF (.NOT.ALLOCATED(fco2_land_nep)) ALLOCATE(fco2_land_nep(klon), stat=ierr)
     232       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_nep',1)
     233       fco2_land_nep(1:klon) = 0.
     234
     235       IF (.NOT.ALLOCATED(fco2_land_fLuc)) ALLOCATE(fco2_land_fLuc(klon), stat=ierr)
     236       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_fLuc',1)
     237       fco2_land_fLuc(1:klon) = 0.
     238
     239       IF (.NOT.ALLOCATED(fco2_land_fwoodharvest)) ALLOCATE(fco2_land_fwoodharvest(klon), stat=ierr)
     240       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_fwoodharvest',1)
     241       fco2_land_fwoodharvest(1:klon) = 0.
     242
     243       IF (.NOT.ALLOCATED(fco2_land_fHarvest)) ALLOCATE(fco2_land_fHarvest(klon), stat=ierr)
     244       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_land_fHarvest',1)
     245       fco2_land_fHarvest(1:klon) = 0.
     246
     247       IF (.NOT.ALLOCATED(fco2_ff)) ALLOCATE(fco2_ff(klon), stat=ierr)
     248       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_ff',1)
     249       fco2_ff(1:klon) = 0.
     250
     251       IF (.NOT.ALLOCATED(fco2_bb)) ALLOCATE(fco2_bb(klon), stat=ierr)
     252       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_bb',1)
     253       fco2_bb(1:klon) = 0.
     254
     255       IF (.NOT.ALLOCATED(fco2_ocean)) ALLOCATE(fco2_ocean(klon), stat=ierr)
     256       IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation fco2_ocean',1)
     257       fco2_bb(1:klon) = 0.
     258    ENDIF
    404259
    405260  END SUBROUTINE carbon_cycle_init
    406261
    407   SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
    408 ! Subroutine for injection of co2 in the tracers
    409 !
    410 ! - Find out if it is time to update
    411 ! - Get tracer from coupled model or from file
    412 ! - Calculate new RCO2 value for the radiation scheme
    413 ! - Calculate CO2 flux to send to ocean and land models (PISCES and ORCHIDEE)
    414 
    415     USE infotrac_phy, ONLY: nbtr
    416     USE dimphy
    417     USE mod_phys_lmdz_transfert_para
    418     USE phys_cal_mod, ONLY : mth_cur, mth_len
    419     USE phys_cal_mod, ONLY : day_cur
    420     USE indice_sol_mod
    421     USE print_control_mod, ONLY: lunout
    422     USE geometry_mod, ONLY : cell_area
    423 
    424     IMPLICIT NONE
    425 
    426     INCLUDE "clesphys.h"
    427     INCLUDE "YOMCST.h"
    428 
    429 ! In/Output arguments
    430     INTEGER,INTENT(IN) :: nstep      ! time step in physiq
    431     REAL,INTENT(IN)    :: pdtphys    ! length of time step in physiq (sec)
    432     REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf            ! Surface fraction
    433     REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri ! All tracers
    434     REAL, DIMENSION(klon,nbtr), INTENT(INOUT)       :: source  ! Source for all tracers
    435 
    436 ! Local variables
    437     INTEGER :: it
    438     LOGICAL :: newmonth ! indicates if a new month just started
    439     LOGICAL :: newday   ! indicates if a new day just started
    440     LOGICAL :: endday   ! indicated if last time step in a day
    441 
    442     REAL, PARAMETER :: fact=1.E-15/2.12  ! transformation factor from gC/m2/day => ppm/m2/day
    443     REAL, DIMENSION(klon) :: fco2_tmp
    444     REAL :: sumtmp
    445     REAL :: delta_co2_ppm
    446    
    447 
    448 ! 1) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
    449 ! -------------------------------------------------------------------------------------------------------
    450 
    451     newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.
    452 
    453     IF (MOD(nstep,INT(86400./pdtphys))==1) newday=.TRUE.
    454     IF (MOD(nstep,INT(86400./pdtphys))==0) endday=.TRUE.
    455     IF (newday .AND. day_cur==1) newmonth=.TRUE.
    456 
    457 ! 2)  For each carbon tracer find out if it is time to inject (update)
    458 ! --------------------------------------------------------------------
    459     DO it = 1, ntr_co2
    460        IF ( MOD(nstep,INT(co2trac(it)%updatefreq/pdtphys)) == 1 ) THEN
    461           co2trac(it)%updatenow = .TRUE.
    462        ELSE
    463           co2trac(it)%updatenow = .FALSE.
    464        END IF
    465     END DO
    466 
    467 ! 3) Get tracer update
    468 ! --------------------------------------
    469     DO it = 1, ntr_co2
    470        IF ( co2trac(it)%updatenow ) THEN
    471           IF ( co2trac(it)%cpl ) THEN
    472              ! Get tracer from coupled model
    473              SELECT CASE(co2trac(it)%name)
    474              CASE('fCO2_land')     ! from ORCHIDEE
    475                 dtr_add(:,it) = fco2_land_inst(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
    476              CASE('fCO2_land_use') ! from ORCHIDEE
    477                 dtr_add(:,it) = fco2_lu_inst(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
    478              CASE('fCO2_ocn')      ! from PISCES
    479                 dtr_add(:,it) = fco2_ocn_day(:)  *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
    480              CASE DEFAULT
    481                 WRITE(lunout,*) 'Error with tracer ',co2trac(it)%name
    482                 CALL abort_physic('carbon_cycle', 'No coupling implemented for this tracer',1)
    483              END SELECT
    484           ELSE
    485              ! Read tracer from file
    486              co2trac(it)%readstep = co2trac(it)%readstep + 1 ! increment time step in file
    487 ! Patricia   CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.FALSE.,dtr_add(:,it))
    488              CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.TRUE.,dtr_add(:,it))
    489 
    490              ! Converte from kgC/m2/h to kgC/m2/s
    491              dtr_add(:,it) = dtr_add(:,it)/3600
    492              ! Add individual treatment of values read from file
    493              SELECT CASE(co2trac(it)%name)
    494              CASE('fCO2_land')
    495                 dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
    496              CASE('fCO2_land_use')
    497                 dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
    498              CASE('fCO2_ocn')
    499                 dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_oce)
    500 ! Patricia :
    501 !             CASE('fCO2_fos_fuel')
    502 !                dtr_add(:,it) = dtr_add(:,it)/mth_len
    503 !                co2trac(it)%readstep = 0 ! Always read same value for fossil fuel(Cadule case)
    504              END SELECT
    505           END IF
    506        END IF
    507     END DO
    508 
    509 ! 4) Update co2 tracers :
    510 !    Loop over all carbon tracers and add source
    511 ! ------------------------------------------------------------------
    512     IF (carbon_cycle_tr) THEN
    513        DO it = 1, ntr_co2
    514           IF (.FALSE.) THEN
    515              tr_seri(1:klon,1,co2trac(it)%id) = tr_seri(1:klon,1,co2trac(it)%id) + dtr_add(1:klon,it)
    516              source(1:klon,co2trac(it)%id) = 0.
    517           ELSE
    518              source(1:klon,co2trac(it)%id) = dtr_add(1:klon,it)
    519           END IF
    520        END DO
    521     END IF
    522 
    523 
    524 ! 5) Calculations for new CO2 value for the radiation scheme(instead of reading value from .def)
    525 ! ----------------------------------------------------------------------------------------------
    526     IF (RCO2_inter) THEN
    527        ! Cumulate fluxes from ORCHIDEE at each timestep
    528        IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
    529           IF (newday) THEN ! Reset cumulative variables once a day
    530              fco2_land_day(1:klon) = 0.
    531              fco2_lu_day(1:klon)   = 0.
    532           END IF
    533           fco2_land_day(1:klon) = fco2_land_day(1:klon) + fco2_land_inst(1:klon) ![gC/m2/day]
    534           fco2_lu_day(1:klon)   = fco2_lu_day(1:klon)   + fco2_lu_inst(1:klon)   ![gC/m2/day]
    535        END IF
    536 
    537        ! At the end of a new day, calculate a mean scalare value of CO2
    538        ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
    539        IF (endday) THEN
    540 
    541           IF (carbon_cycle_tr) THEN
    542              ! Sum all co2 tracers to get the total delta CO2 flux
    543              fco2_tmp(:) = 0.
    544              DO it = 1, ntr_co2
    545                 fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
    546              END DO
    547              
    548           ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
    549              ! Sum co2 fluxes comming from coupled models and parameter for fossil fuel
    550              fco2_tmp(1:klon) = fos_fuel_s + ((fco2_lu_day(1:klon) + fco2_land_day(1:klon))*pctsrf(1:klon,is_ter) &
    551                   + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
    552           END IF
    553 
    554           ! Calculate a global mean value of delta CO2 flux
    555           fco2_tmp(1:klon) = fco2_tmp(1:klon) * cell_area(1:klon)
    556           CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    557           CALL bcast(sumtmp)
    558           delta_co2_ppm = sumtmp/airetot
    559          
    560           ! Add initial value for co2_ppm and delta value
    561           co2_ppm = co2_ppm0 + delta_co2_ppm
    562          
    563           ! Transformation of atmospheric CO2 concentration for the radiation code
    564           RCO2 = co2_ppm * 1.0e-06 * RMCO2 / RMD
    565          
    566           WRITE(lunout,*) 'RCO2 is now updated! RCO2 = ', RCO2
    567        END IF ! endday
    568 
    569     END IF ! RCO2_inter
    570 
    571 
    572 ! 6) Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
    573 ! ----------------------------------------------------------------------------
    574     IF (carbon_cycle_cpl) THEN
    575 
    576        IF (carbon_cycle_tr) THEN
    577           ! Sum all co2 tracers to get the total delta CO2 flux at first model layer
    578           fco2_tmp(:) = 0.
    579           DO it = 1, ntr_co2
    580              fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
    581           END DO
    582           co2_send(1:klon) = fco2_tmp(1:klon) + co2_ppm0
    583        ELSE
    584           ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
    585           co2_send(1:klon) = co2_ppm
    586        END IF
    587 
    588     END IF
    589 
    590   END SUBROUTINE carbon_cycle
    591  
    592262  SUBROUTINE infocfields_init
    593263

  • LMDZ6/trunk/libf/phylmd/phyetat0.F90

    r3505 r3581  
    114114  ENDIF
    115115
    116   ! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
    117   co2_ppm0   = tab_cntrl(16)
     116! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
     117! co2_ppm0   = tab_cntrl(16)
     118! initial value for interactive CO2 run when there is no tracer field for CO2 in restart
     119  co2_ppm0=284.32
    118120
    119121  solaire_etat0      = tab_cntrl(4)

  • LMDZ6/trunk/libf/phylmd/tracco2i_mod.F90

    r3549 r3581  
    22!
    33! This module does the work for the interactive CO2 tracers
     4! Authors: Patricia Cadule and Olivier Boucher
     5!
     6! Purpose and description:
     7!  -----------------------
     8! Main routine for the interactive carbon cycle
     9! Gather all carbon fluxes and emissions from ORCHIDEE, PISCES and fossil fuel
     10! Compute the net flux in source field which is used in phytrac
     11! Compute global CO2 mixing ratio for radiation scheme if option is activated
     12! Redistribute CO2 evenly over the atmosphere if transport is desactivated
    413!
    514CONTAINS
    … …  
    1221    USE infotrac_phy
    1322    USE geometry_mod, ONLY: cell_area
     23    USE carbon_cycle_mod, ONLY: carbon_cycle_init
    1424    USE carbon_cycle_mod, ONLY: id_CO2, nbcf_in, fields_in, cfname_in
    1525    USE carbon_cycle_mod, ONLY: fco2_ocn_day, fco2_ff, fco2_bb, fco2_land, fco2_ocean
    16     USE carbon_cycle_mod, ONLY: carbon_cycle_tr, carbon_cycle_rad, RCO2_glo, RCO2_tot
     26    USE carbon_cycle_mod, ONLY: fco2_land_nbp, fco2_land_nep, fco2_land_fLuc
     27    USE carbon_cycle_mod, ONLY: fco2_land_fwoodharvest, fco2_land_fHarvest
     28    USE carbon_cycle_mod, ONLY: carbon_cycle_cpl, carbon_cycle_tr, carbon_cycle_rad, RCO2_glo, RCO2_tot
    1729    USE mod_grid_phy_lmdz
    1830    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
    … …  
    5466    REAL, DIMENSION(klon_glo,klev) :: m_air_glo ! variable temporaire sur la grille global
    5567
    56     INTEGER, SAVE :: mth_pre=0, day_pre=0
    57 !$OMP THREADPRIVATE(mth_pre, day_pre)
     68    LOGICAL, SAVE :: check_fCO2_nbp_in_cfname
     69!$OMP THREADPRIVATE(check_fCO2_nbp_in_cfname)
     70    INTEGER, SAVE :: day_pre=-1
     71!$OMP THREADPRIVATE(day_pre)
    5872
    5973    IF (is_mpi_root) THEN
    … …  
    6680    IF (debutphy) THEN
    6781
     82! Initialisation de module carbon_cycle_mod
     83      IF (carbon_cycle_cpl) THEN
     84        CALL carbon_cycle_init()
     85      ENDIF
     86
     87! Initialisation de tr_seri(id_CO2) si pas initialise
    6888      IF (MAXVAL(tr_seri(:,:,id_CO2)).LT.1.e-15) THEN
    69         !!tr_seri(:,:,id_CO2)=280.e-6/RMD*RMCO2
    70         tr_seri(:,:,id_CO2)=400.e-6/RMD*RMCO2 !--initialised to 400 ppm for a test
     89        tr_seri(:,:,id_CO2)=co2_ppm0*1.e-6/RMD*RMCO2 !--initialised from co2_ppm0 in rdem
    7190      ENDIF
    7291
    73       ALLOCATE(fco2_ff(klon))
    74       ALLOCATE(fco2_bb(klon))
    75       ALLOCATE(fco2_land(klon))
    76       ALLOCATE(fco2_ocean(klon))
     92!--check if fCO2_nbp is in
     93      check_fCO2_nbp_in_cfname=.FALSE.
     94      DO nb=1, nbcf_in
     95        IF (cfname_in(nb)=="fCO2_nbp") check_fCO2_nbp_in_cfname=.TRUE.
     96      ENDDO
    7797
    7898    ENDIF
    … …  
    93113    fco2_land(:)=0.0
    94114    fco2_ocean(:)=0.0
     115    fco2_land_nbp(:)=0.
     116    fco2_land_nep(:)=0.
     117    fco2_land_fLuc(:)=0.
     118    fco2_land_fwoodharvest(:)=0.
     119    fco2_land_fHarvest(:)=0.
     120
    95121    DO nb=1, nbcf_in
    96       print *,'nb tracco2=', nb, cfname_in(nb)
    97 !--fCO2_nep comes in unit of kg C m-2 s-1
    98 !--converting to kg CO2 m-2 s-1
    99       IF (cfname_in(nb) == "fCO2_nbp" )   fco2_land(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
    100 !--fCO2_fgco2 comes in unit of mol C02 m-2 s-1
    101 !--converting to kg CO2 m-2 s-1 + change sign
    102       IF (cfname_in(nb) == "fCO2_fgco2" ) fco2_ocean(:)=-1.*fco2_ocn_day(:)*RMCO2/1.e3*(pctsrf(:,is_oce)+pctsrf(:,is_sic))
     122
     123      SELECT CASE(cfname_in(nb))
     124!--dealing with the different fluxes coming from ORCHIDEE
     125!--fluxes come in unit of kg C m-2 s-1 is converted into kg CO2 m-2 s-1
     126      CASE("fCO2_nep")
     127          fco2_land_nep(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
     128      CASE("fCO2_fLuc")
     129          fco2_land_fLuc(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
     130      CASE("fCO2_fwoodharvest")
     131          fco2_land_fwoodharvest(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
     132      CASE("fCO2_fHarvest")
     133          fco2_land_fHarvest(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
     134      CASE("fCO2_nbp")
     135          fco2_land_nbp(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
     136!--fCO2_fco2_ocn comes in unit of mol C02 m-2 s-1 is converted into kg CO2 m-2 s-1 + change sign
     137      CASE("fCO2_fgco2")
     138          fco2_ocean(:)=-1.*fco2_ocn_day(:)*RMCO2/1.e3*(pctsrf(:,is_oce)+pctsrf(:,is_sic))
     139      END SELECT
     140
    103141    ENDDO
    104142
    105 !--preparing the net anthropogenic flux at the surface for mixing layer
    106 !--unit kg CO2 / m2 / s
     143!--if fCO2_nbp is transferred we use it, otherwise we use the sum of what has been passed from ORCHIDEE
     144    IF (check_fCO2_nbp_in_cfname)  THEN
     145       fco2_land(:)=fco2_land_nbp(:)
     146    ELSE
     147       fco2_land(:)=fco2_land_nep(:)+fco2_land_fLuc(:)+fco2_land_fwoodharvest(:)+fco2_land_fHarvest(:)
     148    ENDIF
     149
     150!!--preparing the net anthropogenic flux at the surface for mixing layer
     151!!--unit kg CO2 / m2 / s
     152!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_ff) ',MAXVAL(fco2_ff)
     153!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_ff) ',MINVAL(fco2_ff)
     154!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_bb) ',MAXVAL(fco2_bb)
     155!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_bb) ',MINVAL(fco2_bb)
     156!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_land) ',MAXVAL(fco2_land)
     157!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_land) ',MINVAL(fco2_land)
     158!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_ocean) ',MAXVAL(fco2_ocean)
     159!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_ocean) ',MINVAL(fco2_ocean)
     160!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(source(:,id_CO2)) ',MAXVAL(source(:,id_CO2))
     161!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(source(:,id_CO2)) ',MINVAL(source(:,id_CO2))
     162!
     163!--build final source term for CO2
    107164    source(:,id_CO2)=fco2_ff(:)+fco2_bb(:)+fco2_land(:)+fco2_ocean(:)
    108165
    109166!--computing global mean CO2 for radiation
    110 !--every timestep for now but enough every day
    111 !    IF (debutphy.OR.mth_cur.NE.mth_pre) THEN
    112 !    IF (debutphy.OR.day_cur.NE.day_pre) THEN
     167!--for every timestep comment out the IF ENDIF statements
     168!--otherwise this is updated every day
     169    IF (debutphy.OR.day_cur.NE.day_pre) THEN
     170
    113171      CALL gather(tr_seri(:,:,id_CO2),co2_glo)
    114172      CALL gather(m_air,m_air_glo)
     173
    115174!$OMP MASTER
    116175
    … …  
    124183!$OMP END MASTER
    125184       CALL bcast(RCO2_glo)
    126        mth_pre=mth_cur
    127185       day_pre=day_cur
    128186!--if not carbon_cycle_tr, then we reinitialize the CO2 each day to its global mean value
    … …  
    130188         tr_seri(:,:,id_CO2)=RCO2_glo
    131189       ENDIF
    132 !    ENDIF
     190    ENDIF
    133191
    134192  END SUBROUTINE tracco2i
    … …  
    168226!! may be controlled via the .def later on
    169227!! also co2bb for now comes from ORCHIDEE
    170     LOGICAL, PARAMETER :: readco2ff=.TRUE., readco2bb=.FALSE.
     228    LOGICAL, PARAMETER :: readco2ff=.TRUE.
     229!! this should be left to FALSE for now
     230    LOGICAL, PARAMETER :: readco2bb=.FALSE.
    171231
    172232    CHARACTER (len = 20) :: modname = 'tracco2i.co2_emissions'
    … …  
    217277
    218278!--reading CO2 biomass burning emissions
     279!--using it will be inconsistent with treatment in ORCHIDEE
    219280      IF (readco2bb) THEN
    220281

  • LMDZ6/trunk/libf/phylmd/traclmdz_mod.F90

    r2320 r3581  
    9292    USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
    9393    USE press_coefoz_m, ONLY: press_coefoz
    94     USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl
    9594    USE mod_grid_phy_lmdz
    9695    USE mod_phys_lmdz_para
    … …  
    285284
    286285!
    287 ! Initialisation de module carbon_cycle_mod
    288 ! ----------------------------------------------
    289     IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
    290        CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
    291     END IF
    292 
    293286! Check if all tracers have restart values
    294287! ----------------------------------------------
    … …  
    346339    USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
    347340    USE o3_chem_m, ONLY: o3_chem
    348     USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
    349341    USE indice_sol_mod
    350342
    … …  
    612604    END IF
    613605
    614 !======================================================================
    615 !   Calcul de cycle de carbon
    616 !======================================================================
    617     IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
    618        CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
    619     END IF
    620 
    621606  END SUBROUTINE traclmdz
    622607
Note: See TracChangeset for help on using the changeset viewer.