source: LMDZ6/trunk/libf/phylmd/carbon_cycle_mod.F90 @ 3490

Last change on this file since 3490 was 3459, checked in by oboucher, 5 years ago

Changing hard-coded value of 28.97 g/mol by RMD physical constant defined in suphel.F90
As rmd = 28.9644, this involves a lack of numerical convergence with previous commit.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 33.8 KB
RevLine 
[1227]1MODULE carbon_cycle_mod
[3387]2!=======================================================================
3!   Authors: Patricia Cadule and Laurent Fairhead
4!            base sur un travail anterieur mene par Patricia Cadule et Josefine Ghattas
5!
6!  Purpose and description:
7!  -----------------------
8! Control module for the carbon CO2 tracers :
[1454]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!
[3387]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)
16!
17! Le cas online/offline est defini par le flag carbon_cycle_cpl (y/n)
18! Le transport du traceur CO2 est defini par le flag carbon_cycle_tr (y/n)
19! la provenance des champs (termes de puits) est denini par le flag level_coupling_esm
20!
21! level_coupling_esm : level of coupling of the biogeochemical fields between
22! LMDZ, ORCHIDEE and NEMO
23! Definitions of level_coupling_esm in physiq.def
24! level_coupling_esm = 0  ! No field exchange between LMDZ and ORCHIDEE models
25!                         ! No field exchange between LMDZ and NEMO
26! level_coupling_esm = 1  ! Field exchange between LMDZ and ORCHIDEE models
27!                         ! No field exchange between LMDZ and NEMO models
28! level_coupling_esm = 2  ! No field exchange between LMDZ and ORCHIDEE models
29!                         ! Field exchange between LMDZ and NEMO models
30! level_coupling_esm = 3  ! Field exchange between LMDZ and ORCHIDEE models
31!                         ! Field exchange between LMDZ and NEMO models
32!=======================================================================
[1227]33
34  IMPLICIT NONE
35  SAVE
36  PRIVATE
[3387]37  PUBLIC :: carbon_cycle_init, carbon_cycle, infocfields_init
[1227]38
39! Variables read from parmeter file physiq.def
[3447]40  LOGICAL, PUBLIC :: carbon_cycle_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES)
41!$OMP THREADPRIVATE(carbon_cycle_cpl)
[1250]42  LOGICAL, PUBLIC :: carbon_cycle_tr        ! 3D transport of CO2 in the atmosphere, parameter read in conf_phys
[1249]43!$OMP THREADPRIVATE(carbon_cycle_tr)
[3447]44  LOGICAL, PUBLIC :: carbon_cycle_rad       ! CO2 interactive radiatively
45!$OMP THREADPRIVATE(carbon_cycle_rad)
[3387]46  INTEGER, PUBLIC :: level_coupling_esm ! Level of coupling for the ESM - 0, 1, 2, 3
[3384]47!$OMP THREADPRIVATE(level_coupling_esm)
[3447]48  REAL, PUBLIC :: RCO2_glo, RCO2_tot
49!$OMP THREADPRIVATE(RCO2_glo, RCO2_tot)
[1454]50
[1759]51  LOGICAL :: carbon_cycle_emis_comp_omp=.FALSE.
[1227]52  LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible
[1454]53!$OMP THREADPRIVATE(carbon_cycle_emis_comp)
[1227]54
[1759]55  LOGICAL :: RCO2_inter_omp
[1454]56  LOGICAL :: RCO2_inter  ! RCO2 interactive : if true calculate new value RCO2 for the radiation scheme
57!$OMP THREADPRIVATE(RCO2_inter)
58
[3385]59! Scalare values when no transport, from physiq.def
60  REAL :: fos_fuel_s_omp
61  REAL :: fos_fuel_s  ! carbon_cycle_fos_fuel dans physiq.def
62!$OMP THREADPRIVATE(fos_fuel_s)
[1227]63  REAL :: emis_land_s ! not yet implemented
[1250]64!$OMP THREADPRIVATE(emis_land_s)
[1227]65
[1454]66  REAL :: airetot     ! Total area of the earth surface
67!$OMP THREADPRIVATE(airetot)
[1250]68
[3385]69  INTEGER :: ntr_co2  ! Number of tracers concerning the carbon cycle
[1454]70!$OMP THREADPRIVATE(ntr_co2)
71
[3385]72! fco2_ocn_day : flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]. Allocation and initalization done in cpl_mod
73  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day
74!$OMP THREADPRIVATE(fco2_ocn_day)
[1454]75
[3385]76  REAL, DIMENSION(:), ALLOCATABLE :: fco2_land_day   ! flux CO2 from land for 1 day (cumulated)  [gC/m2/d]
77!$OMP THREADPRIVATE(fco2_land_day)
78  REAL, DIMENSION(:), ALLOCATABLE :: fco2_lu_day     ! Emission from land use change for 1 day (cumulated) [gC/m2/d]
79!$OMP THREADPRIVATE(fco2_lu_day)
[3421]80  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ff ! Emission from fossil fuel [kgCO2/m2/s]
81!$OMP THREADPRIVATE(fco2_ff)
82  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_bb ! Emission from biomass burning [kgCO2/m2/s]
83!$OMP THREADPRIVATE(fco2_bb)
[3453]84  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land  ! Net flux from terrestrial ecocsystems [kgCO2/m2/s]
85!$OMP THREADPRIVATE(fco2_land)
86  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocean ! Net flux from ocean [kgCO2/m2/s]
87!$OMP THREADPRIVATE(fco2_ocean)
[1227]88
[1454]89  REAL, DIMENSION(:,:), ALLOCATABLE :: dtr_add       ! Tracer concentration to be injected
90!$OMP THREADPRIVATE(dtr_add)
91
[3385]92! Following 2 fields will be allocated and initialized in surf_land_orchidee
93  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_inst  ! flux CO2 from land at one time step
94!$OMP THREADPRIVATE(fco2_land_inst)
95  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_lu_inst    ! Emission from land use change at one time step
96!$OMP THREADPRIVATE(fco2_lu_inst)
[1227]97
[3385]98! Calculated co2 field to be send to the ocean via the coupler and to ORCHIDEE
[1454]99  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send ! Field allocated in phyetat0
[1250]100!$OMP THREADPRIVATE(co2_send)
[1227]101
[3447]102  INTEGER, PARAMETER, PUBLIC :: id_CO2=1              !--temporaire OB -- to be changed
103
[3387]104! nbfields : total number of fields
105  INTEGER, PUBLIC :: nbcf
106!$OMP THREADPRIVATE(nbcf)
[1454]107
[3387]108! nbcf_in : number of fields IN
109  INTEGER, PUBLIC  :: nbcf_in
110!$OMP THREADPRIVATE(nbcf_in)
111
112! nbcf_in_orc : number of fields IN
113  INTEGER, PUBLIC  :: nbcf_in_orc
114!$OMP THREADPRIVATE(nbcf_in_orc)
115
116! nbcf_in_inca : number of fields IN (from INCA)
117  INTEGER, PUBLIC  :: nbcf_in_inca
118!$OMP THREADPRIVATE(nbcf_in_inca)
119
120! nbcf_in_nemo : number of fields IN (from nemo)
121  INTEGER, PUBLIC  :: nbcf_in_nemo
122!$OMP THREADPRIVATE(nbcf_in_nemo)
123
124! nbcf_in_ant : number of fields IN (from anthropogenic sources)
125  INTEGER, PUBLIC  :: nbcf_in_ant
126!$OMP THREADPRIVATE(nbcf_in_ant)
127
128! nbcf_out : number of fields OUT
129  INTEGER, PUBLIC :: nbcf_out
130!$OMP THREADPRIVATE(nbcf_out)
131
132! Name of variables
133  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname     ! coupling field short name for restart (?) and diagnostics
134!$OMP THREADPRIVATE(cfname)
135
136  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname_in  ! coupling field short name for restart (?) and diagnostics
137!$OMP THREADPRIVATE(cfname_in)
138
139  CHARACTER(len=25), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfname_out ! coupling field short name for restart (?) and diagnostics
140!$OMP THREADPRIVATE(cfname_out)
141
142  CHARACTER(len=15), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfunits_in  !  coupling field units for diagnostics
143!$OMP THREADPRIVATE(cfunits_in)
144
145  CHARACTER(len=15), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfunits_out !  coupling field units for diagnostics
146!$OMP THREADPRIVATE(cfunits_out)
147
148  CHARACTER(len=120), ALLOCATABLE, DIMENSION(:), PUBLIC :: cftext_in  ! coupling field long name for diagnostics
149!$OMP THREADPRIVATE(cftext_in)
150
151  CHARACTER(len=120), ALLOCATABLE, DIMENSION(:), PUBLIC :: cftext_out ! coupling field long name for diagnostics
152!$OMP THREADPRIVATE(cftext_out)
153
154  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfmod1 ! model 1 (rreference) : LMDz
155!$OMP THREADPRIVATE(cfmod1)
156
157  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:), PUBLIC :: cfmod2 ! model 2
158!$OMP THREADPRIVATE(cfmod2)
159
[3391]160  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), PUBLIC :: field_out_names
161!$OMP THREADPRIVATE(field_out_names)
[3387]162
[3391]163  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), PUBLIC :: field_in_names
164!$OMP THREADPRIVATE(field_in_names)
[3387]165
[3391]166  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC :: fields_in   !  klon,nbcf_in
167!$OMP THREADPRIVATE(fields_in)
168
169  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC :: yfields_in  !  knon,nbcf_in
170!$OMP THREADPRIVATE(yfields_in)
171
172  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC :: fields_out  !  klon,nbcf_out
173!$OMP THREADPRIVATE(fields_out)
174
175  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC :: yfields_out !  knon,nbcf_out
176!$OMP THREADPRIVATE(yfields_out)
177
[3390]178  TYPE, PUBLIC :: co2_trac_type
[1454]179     CHARACTER(len = 8) :: name       ! Tracer name in tracer.def
180     INTEGER            :: id         ! Index in total tracer list, tr_seri
181     CHARACTER(len=30)  :: file       ! File name
[3385]182     LOGICAL            :: cpl        ! True if this tracers is coupled from ORCHIDEE or PISCES.
[1454]183                                      ! False if read from file.
184     INTEGER            :: updatefreq ! Frequence to inject in second
185     INTEGER            :: readstep   ! Actual time step to read in file
186     LOGICAL            :: updatenow  ! True if this tracer should be updated this time step
187  END TYPE co2_trac_type
188  INTEGER,PARAMETER :: maxco2trac=5  ! Maximum number of different CO2 fluxes
189  TYPE(co2_trac_type), DIMENSION(maxco2trac) :: co2trac
190
[1227]191CONTAINS
192 
[1454]193  SUBROUTINE carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
194! This subroutine is called from traclmdz_init, only at first timestep.
195! - Read controle parameters from .def input file
196! - Search for carbon tracers and set default values
197! - Allocate variables
198! - Test for compatibility
199
[1227]200    USE dimphy
[2351]201    USE geometry_mod, ONLY : cell_area
[1454]202    USE mod_phys_lmdz_transfert_para
[2320]203    USE infotrac_phy, ONLY: nbtr, nqo, niadv, tname
[1227]204    USE IOIPSL
205    USE surface_data, ONLY : ok_veget, type_ocean
[1454]206    USE phys_cal_mod, ONLY : mth_len
[2311]207    USE print_control_mod, ONLY: lunout
[1227]208
209    IMPLICIT NONE
210    INCLUDE "clesphys.h"
211 
212! Input argument
213    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA] 
[1454]214    REAL,INTENT(IN)                           :: pdtphys ! length of time step in physiq (sec)
[1227]215
216! InOutput arguments
217    LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: aerosol
218    LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: radio
219
220! Local variables
[1454]221    INTEGER               :: ierr, it, iiq, itc
222    INTEGER               :: teststop
[1227]223
[1454]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
[1759]228!$OMP MASTER
[3385]229       fos_fuel_s_omp = 0.
230       CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s_omp)
[1759]231!$OMP END MASTER
232!$OMP BARRIER
[3385]233       fos_fuel_s=fos_fuel_s_omp
234       WRITE(lunout,*) 'carbon_cycle_fos_fuel = ', fos_fuel_s
[1227]235    END IF
236
[1454]237    ! Read parmeter for calculation compatible emission
238    IF (.NOT. carbon_cycle_tr) THEN
[1759]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
[1454]245       WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
246       IF (carbon_cycle_emis_comp) THEN
[2311]247          CALL abort_physic('carbon_cycle_init', 'carbon_cycle_emis_comp option not yet implemented!!',1)
[1454]248       END IF
249    END IF
[1227]250
[1454]251    ! Read parameter for interactive calculation of the CO2 value for the radiation scheme
[1759]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
[1454]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
[1227]263
[1454]264
[3385]265! 2) Search for carbon tracers and set default values
[1454]266! ---------------------------------------------------
267    itc=0
[3385]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")
[1454]274          itc = itc + 1
[3385]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
[1227]284          END IF
[3385]285       CASE("fCO2_land")
[1454]286          itc = itc + 1
[3385]287          co2trac(itc)%name='fCO2_land'
288          co2trac(itc)%id=it
289          co2trac(itc)%file='fl_co2_land.nc'
[1454]290          IF (carbon_cycle_cpl .AND. ok_veget) THEN
[3385]291             co2trac(itc)%cpl=.TRUE.
292             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
[1454]293          ELSE
[3385]294             co2trac(itc)%cpl=.FALSE.
295!             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
296             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
[1454]297          END IF
[3385]298       CASE("fCO2_land_use")
[1454]299          itc = itc + 1
[3385]300          co2trac(itc)%name='fCO2_land_use'
301          co2trac(itc)%id=it
302          co2trac(itc)%file='fl_co2_land_use.nc'
[1454]303          IF (carbon_cycle_cpl .AND. ok_veget) THEN
[3385]304             co2trac(it)%cpl=.TRUE.
305             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
[1454]306          ELSE
[3385]307             co2trac(itc)%cpl=.FALSE.
308             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
[1454]309          END IF
[3385]310       CASE("fCO2_fos_fuel")
[1454]311          itc = itc + 1
[3385]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")
[1454]319          itc = itc + 1
[3385]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
[1454]325       CASE("fCO2")
326          ! fCO2 : One tracer transporting the total CO2 flux
327          itc = itc + 1
[3385]328          co2trac(itc)%name='fCO2'
329          co2trac(itc)%id=it
330          co2trac(itc)%file='fl_co2.nc'
[1454]331          IF (carbon_cycle_cpl) THEN
[3385]332             co2trac(itc)%cpl=.TRUE.
[1454]333          ELSE
[3385]334             co2trac(itc)%cpl=.FALSE.
[1454]335          END IF
[3385]336          co2trac(itc)%updatefreq = 86400
[1454]337          ! DOES THIS WORK ???? Problematic due to implementation of the coupled fluxes...
[2311]338          CALL abort_physic('carbon_cycle_init','transport of total CO2 has to be implemented and tested',1)
[1454]339       END SELECT
340    END DO
[1227]341
[3385]342    ! Total number of carbon CO2 tracers
343    ntr_co2 = itc
344   
[1454]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   
[1227]355
[1454]356! 3) Allocate variables
357! ---------------------
358    ! Allocate vector for storing fluxes to inject
[3385]359    ALLOCATE(dtr_add(klon,maxco2trac), stat=ierr)
[2311]360    IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 11',1)       
[3385]361   
[1454]362    ! Allocate variables for cumulating fluxes from ORCHIDEE
[3385]363    IF (RCO2_inter) THEN
364       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
365          ALLOCATE(fco2_land_day(klon), stat=ierr)
[2311]366          IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 2',1)
[3385]367          fco2_land_day(1:klon) = 0.
[1454]368         
[3385]369          ALLOCATE(fco2_lu_day(klon), stat=ierr)
[2311]370          IF (ierr /= 0) CALL abort_physic('carbon_cycle_init', 'pb in allocation 3',1)
[3385]371          fco2_lu_day(1:klon)   = 0.
[1227]372       END IF
[3385]373    END IF
[1227]374
[1454]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'
[2311]380!       CALL abort_physic('carbon_cycle_init', 'coupled ocean is needed for carbon_cycle_cpl',1)
[1454]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'
[2311]385!       CALL abort_physic('carbon_cycle_init', 'ok_veget is needed for carbon_cycle_cpl',1)
[1454]386!    END IF
387
388    ! Compiler test : following should never happen
389    teststop=0
390    DO it=1,teststop
[2311]391       CALL abort_physic('carbon_cycle_init', 'Entering loop from 1 to 0',1)
[1454]392    END DO
393
[3385]394    IF (ntr_co2==0) THEN
[1454]395       ! No carbon tracers found in tracer.def. It is not possible to do carbon cycle
[3385]396       WRITE(lunout,*) 'No carbon tracers found in tracer.def. Not ok with carbon_cycle_tr and/or carbon_cycle_cp'
[2311]397       CALL abort_physic('carbon_cycle_init', 'No carbon tracers found in tracer.def',1)
[1227]398    END IF
[1454]399   
400! 5) Calculate total area of the earth surface
401! --------------------------------------------
[2351]402    CALL reduce_sum(SUM(cell_area),airetot)
[1454]403    CALL bcast(airetot)
[1227]404
405  END SUBROUTINE carbon_cycle_init
406
[1454]407  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
408! Subroutine for injection of co2 in the tracers
[1227]409!
[1454]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)
[1227]414
[2320]415    USE infotrac_phy, ONLY: nbtr
[1227]416    USE dimphy
[1454]417    USE mod_phys_lmdz_transfert_para
[1227]418    USE phys_cal_mod, ONLY : mth_cur, mth_len
419    USE phys_cal_mod, ONLY : day_cur
[1785]420    USE indice_sol_mod
[2311]421    USE print_control_mod, ONLY: lunout
[2351]422    USE geometry_mod, ONLY : cell_area
[1227]423
424    IMPLICIT NONE
425
426    INCLUDE "clesphys.h"
[1454]427    INCLUDE "YOMCST.h"
[1227]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)
[1454]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
[1227]435
436! Local variables
[3385]437    INTEGER :: it
[1227]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
[3385]442    REAL, PARAMETER :: fact=1.E-15/2.12  ! transformation factor from gC/m2/day => ppm/m2/day
[1454]443    REAL, DIMENSION(klon) :: fco2_tmp
[1227]444    REAL :: sumtmp
445    REAL :: delta_co2_ppm
446   
447
[1454]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
[3385]451    newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.
[1227]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.
[1454]456
457! 2)  For each carbon tracer find out if it is time to inject (update)
458! --------------------------------------------------------------------
[3385]459    DO it = 1, ntr_co2
460       IF ( MOD(nstep,INT(co2trac(it)%updatefreq/pdtphys)) == 1 ) THEN
461          co2trac(it)%updatenow = .TRUE.
[1454]462       ELSE
[3385]463          co2trac(it)%updatenow = .FALSE.
[1227]464       END IF
[1454]465    END DO
[1227]466
[1454]467! 3) Get tracer update
468! --------------------------------------
[3385]469    DO it = 1, ntr_co2
470       IF ( co2trac(it)%updatenow ) THEN
471          IF ( co2trac(it)%cpl ) THEN
[1454]472             ! Get tracer from coupled model
[3385]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]
[1454]480             CASE DEFAULT
[3385]481                WRITE(lunout,*) 'Error with tracer ',co2trac(it)%name
[2311]482                CALL abort_physic('carbon_cycle', 'No coupling implemented for this tracer',1)
[1454]483             END SELECT
484          ELSE
485             ! Read tracer from file
[3385]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))
[1227]489
[3385]490             ! Converte from kgC/m2/h to kgC/m2/s
491             dtr_add(:,it) = dtr_add(:,it)/3600
[1454]492             ! Add individual treatment of values read from file
[3385]493             SELECT CASE(co2trac(it)%name)
494             CASE('fCO2_land')
[1454]495                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
[3385]496             CASE('fCO2_land_use')
[1454]497                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
[3385]498             CASE('fCO2_ocn')
[1454]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
[3385]503!                co2trac(it)%readstep = 0 ! Always read same value for fossil fuel(Cadule case)
[1454]504             END SELECT
[1227]505          END IF
[1454]506       END IF
507    END DO
[1227]508
[1454]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
[3385]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
[1454]520       END DO
[1227]521    END IF
522
523
[1454]524! 5) Calculations for new CO2 value for the radiation scheme(instead of reading value from .def)
525! ----------------------------------------------------------------------------------------------
[3385]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
[1227]536
[3385]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
[1227]540
[1454]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
[3385]547             
548          ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
[1454]549             ! Sum co2 fluxes comming from coupled models and parameter for fossil fuel
[3385]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
[1227]552          END IF
553
[1454]554          ! Calculate a global mean value of delta CO2 flux
[3385]555          fco2_tmp(1:klon) = fco2_tmp(1:klon) * cell_area(1:klon)
[1227]556          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
[1454]557          CALL bcast(sumtmp)
[1227]558          delta_co2_ppm = sumtmp/airetot
559         
[1454]560          ! Add initial value for co2_ppm and delta value
561          co2_ppm = co2_ppm0 + delta_co2_ppm
[3385]562         
[1454]563          ! Transformation of atmospheric CO2 concentration for the radiation code
[3459]564          RCO2 = co2_ppm * 1.0e-06 * RMCO2 / RMD
[1454]565         
[3385]566          WRITE(lunout,*) 'RCO2 is now updated! RCO2 = ', RCO2
567       END IF ! endday
[1227]568
[3385]569    END IF ! RCO2_inter
[1227]570
571
[1454]572! 6) Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
573! ----------------------------------------------------------------------------
[3385]574    IF (carbon_cycle_cpl) THEN
[1454]575
[3385]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.
[1454]579          DO it = 1, ntr_co2
[3385]580             fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
[1454]581          END DO
582          co2_send(1:klon) = fco2_tmp(1:klon) + co2_ppm0
[1227]583       ELSE
[3385]584          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
[1454]585          co2_send(1:klon) = co2_ppm
[1227]586       END IF
587
[3385]588    END IF
[1227]589
590  END SUBROUTINE carbon_cycle
591 
[3387]592  SUBROUTINE infocfields_init
593
[3435]594!    USE control_mod, ONLY: planet_type
[3387]595    USE phys_cal_mod, ONLY : mth_cur
596    USE mod_synchro_omp
597    USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_omp_root
598    USE mod_phys_lmdz_transfert_para
599    USE mod_phys_lmdz_omp_transfert
600    USE dimphy, ONLY: klon
601
602    IMPLICIT NONE
603
604!=======================================================================
605!
606!   Authors: Patricia Cadule and Laurent Fairhead 
607!   -------
608!
609!  Purpose and description:
610!  -----------------------
611!
612! Infofields
613! this routine enables to define the field exchanges in both directions between
614! the atmospheric circulation model (LMDZ) and ORCHIDEE. In the future this
615! routing might apply to other models (e.g., NEMO, INCA, ...).
616! Therefore, currently with this routine, it is possible to define the coupling
617! fields only between LMDZ and ORCHIDEE.
618! The coupling_fields.def file enables to define the name of the exchanged
619! fields at the coupling interface.
620! field_in_names : the set of names of the exchanged fields in input to ORCHIDEE
621! (LMDZ to ORCHIDEE)
622! field_out_names : the set of names of the exchanged fields in output of
623! ORCHIDEE (ORCHIDEE to LMDZ)
624! n : the number of exchanged fields at th coupling interface
625! nb_fields_in : number of inputs fields to ORCHIDEE (LMDZ to ORCHIDEE)
626! nb_fields_out : number of ouput fields of ORCHIDEE (ORCHIDEE to LMDZ)
627!
628! The syntax for coupling_fields.def is as follows:
629! IMPORTANT: each column entry must be separated from the previous one by 3
630! spaces and only that
631! field name         coupling          model 1         model 2         long_name
632!                    direction
633!   10char  -3spaces-  3char  -3spaces- 4char -3spaces- 4char -3spaces-  30char
634!
635! n
636! FIELD1 IN LMDZ ORC
637! ....
638! FIELD(j) IN LMDZ ORC
639! FIELD(j+1) OUT LMDZ ORC
640! ...
641! FIELDn OUT LMDZ ORC   
642!
643!=======================================================================
644!   ... 22/12/2017 ....
645!-----------------------------------------------------------------------
646! Declarations
647
648  INCLUDE "clesphys.h"
649  INCLUDE "dimensions.h"
650  INCLUDE "iniprint.h"
651
652! Local variables
653
654  INTEGER :: iq,  ierr, stat, error
655
656  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), SAVE  :: cfname_root
657  CHARACTER(LEN=120), ALLOCATABLE, DIMENSION(:), SAVE :: cftext_root
658  CHARACTER(LEN=15), ALLOCATABLE, DIMENSION(:), SAVE  :: cfunits_root
659
660  CHARACTER(len=3), ALLOCATABLE, DIMENSION(:) :: cfintent_root
661  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:) :: cfmod1_root
662  CHARACTER(len=5), ALLOCATABLE, DIMENSION(:) :: cfmod2_root
663
664  LOGICAL, ALLOCATABLE, DIMENSION(:), SAVE :: mask_in_root
665  LOGICAL, ALLOCATABLE, DIMENSION(:), SAVE :: mask_out_root
666
667  CHARACTER(len=*),parameter :: modname="infocfields"
668
[3435]669  CHARACTER(len=10),SAVE :: planet_type="earth"
670
[3387]671!-----------------------------------------------------------------------
672
673nbcf=0
674nbcf_in=0
675nbcf_out=0
676
677IF (planet_type=='earth') THEN
678
679    IF (is_mpi_root .AND. is_omp_root) THEN
680
681       IF (level_coupling_esm.GT.0) THEN
682
683          OPEN(200,file='coupling_fields.def',form='formatted',status='old', iostat=ierr)
684
685          IF (ierr.EQ.0) THEN
686
687             WRITE(lunout,*) trim(modname),': Open coupling_fields.def : ok'
688             READ(200,*) nbcf
689             WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf=',nbcf
[3390]690             ALLOCATE(cfname_root(nbcf))
691             ALLOCATE(cfintent_root(nbcf))
692             ALLOCATE(cfmod1_root(nbcf))
693             ALLOCATE(cfmod2_root(nbcf))
694             ALLOCATE(cftext_root(nbcf))
695             ALLOCATE(cfunits_root(nbcf))
696             ALLOCATE(mask_in_root(nbcf))
697             ALLOCATE(mask_out_root(nbcf))
[3387]698
699             nbcf_in=0
700             nbcf_out=0
[3390]701
[3387]702             DO iq=1,nbcf
703                WRITE(lunout,*) 'infofields : field=',iq
704                READ(200,'(A15,3X,A3,3X,A5,3X,A5,3X,A120,3X,A15)',IOSTAT=ierr) &
705                   cfname_root(iq),cfintent_root(iq),cfmod1_root(iq),cfmod2_root(iq),cftext_root(iq),cfunits_root(iq)
706                cfname_root(iq)=TRIM(cfname_root(iq))
707                cfintent_root(iq)=TRIM(cfintent_root(iq))
708                cfmod1_root(iq)=TRIM(cfmod1_root(iq))
709                cfmod2_root(iq)=TRIM(cfmod2_root(iq))
710                cftext_root(iq)=TRIM(cftext_root(iq))
711                cfunits_root(iq)=TRIM(cfunits_root(iq))
712                WRITE(lunout,*) 'coupling field: ',cfname_root(iq), & 
713                               ', number: ',iq,', INTENT: ',cfintent_root(iq)
714                WRITE(lunout,*) 'coupling field: ',cfname_root(iq), &
715                               ', number: ',iq,', model 1 (ref): ',cfmod1_root(iq),', model 2: ',cfmod2_root(iq)
716                WRITE(lunout,*) 'coupling field: ',cfname_root(iq), &
717                               ', number: ',iq,', long name: ',cftext_root(iq),', units ',cfunits_root(iq)
[3390]718                IF (nbcf_in+nbcf_out.LT.nbcf) THEN
719                  IF (cfintent_root(iq).NE.'OUT') THEN
720                    nbcf_in=nbcf_in+1
721                    mask_in_root(iq)=.TRUE.
722                    mask_out_root(iq)=.FALSE.
723                  ELSE IF (cfintent_root(iq).EQ.'OUT') THEN
724                    nbcf_out=nbcf_out+1
725                    mask_in_root(iq)=.FALSE.
726                    mask_out_root(iq)=.TRUE.
727                  ENDIF
[3387]728                ELSE
[3390]729                  WRITE(lunout,*) 'abort_gcm --- nbcf    : ',nbcf
730                  WRITE(lunout,*) 'abort_gcm --- nbcf_in : ',nbcf_in
[3387]731                  WRITE(lunout,*) 'abort_gcm --- nbcf_out: ',nbcf_out
[3435]732                  CALL abort_physic('infocfields_init','Problem in the definition of the coupling fields',1)
[3387]733               ENDIF
734             ENDDO !DO iq=1,nbcf
735          ELSE
736             WRITE(lunout,*) trim(modname),': infocfields_mod.F90 --- Problem in opening coupling_fields.def'
737             WRITE(lunout,*) trim(modname),': infocfields_mod.F90 --- WARNING using defaut values'
738          ENDIF ! ierr
739          CLOSE(200)
740
741       ENDIF ! level_coupling_esm
742
743    ENDIF !   (is_mpi_root .AND. is_omp_root)
744!$OMP BARRIER
745
746    CALL bcast(nbcf)
747    CALL bcast(nbcf_in)
748    CALL bcast(nbcf_out)
749
[3390]750    WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf    =',nbcf
751    WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf_in =',nbcf_in
[3387]752    WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf_out=',nbcf_out
753
[3390]754    ALLOCATE(cfname(nbcf))
755    ALLOCATE(cfname_in(nbcf_in))
756    ALLOCATE(cftext_in(nbcf_in))
757    ALLOCATE(cfname_out(nbcf_out))
758    ALLOCATE(cftext_out(nbcf_out))
759    ALLOCATE(cfmod1(nbcf))
760    ALLOCATE(cfmod2(nbcf))
761    ALLOCATE(cfunits_in(nbcf_in))
762    ALLOCATE(cfunits_out(nbcf_out))
[3387]763       
764    IF (is_mpi_root .AND. is_omp_root) THEN
765
[3390]766        IF (nbcf.GT.0)     cfname=cfname_root
767        IF (nbcf_in.GT.0)  cfname_in=PACK(cfname_root,mask_in_root)
[3387]768        IF (nbcf_out.GT.0) cfname_out=PACK(cfname_root,mask_out_root)
[3390]769        IF (nbcf_in.GT.0)  cftext_in=PACK(cftext_root,mask_in_root)
[3387]770        IF (nbcf_out.GT.0) cftext_out=PACK(cftext_root,mask_out_root)
[3390]771        IF (nbcf.GT.0)     cfmod1=cfmod1_root
772        IF (nbcf.GT.0)     cfmod2=cfmod2_root
773        IF (nbcf_in.GT.0)  cfunits_in=PACK(cfunits_root,mask_in_root)
[3387]774        IF (nbcf_out.GT.0) cfunits_out=PACK(cfunits_root,mask_out_root)
775
776        nbcf_in_orc=0
777        nbcf_in_nemo=0
778        nbcf_in_inca=0
779        nbcf_in_ant=0
780
781        DO iq=1,nbcf
[3390]782            IF (cfmod1(iq) == "ORC")  nbcf_in_orc  = nbcf_in_orc  + 1 
783            IF (cfmod1(iq) == "NEMO") nbcf_in_nemo = nbcf_in_nemo + 1 
784            IF (cfmod1(iq) == "INCA") nbcf_in_inca = nbcf_in_inca + 1
785            IF (cfmod1(iq) == "ALL")  nbcf_in_orc  = nbcf_in_orc  + 1  ! ALL = ORC/NEMO/INCA
786            IF (cfmod1(iq) == "ALL")  nbcf_in_nemo = nbcf_in_nemo + 1  ! ALL = ORC/NEMO/INCA
787            IF (cfmod1(iq) == "ALL")  nbcf_in_inca = nbcf_in_inca + 1  ! ALL = ORC/NEMO/INCA
788            IF (cfmod1(iq) == "ANT")  nbcf_in_ant  = nbcf_in_ant  + 1 
[3387]789        ENDDO
790
791    ENDIF !   (is_mpi_root .AND. is_omp_root)
792!$OMP BARRIER
793
794    CALL bcast(nbcf_in_orc)
795    CALL bcast(nbcf_in_nemo)
796    CALL bcast(nbcf_in_inca)
797    CALL bcast(nbcf_in_ant)
798
[3390]799    WRITE(lunout,*) 'nbcf_in_orc  =',nbcf_in_orc
800    WRITE(lunout,*) 'nbcf_in_nemo =',nbcf_in_nemo
801    WRITE(lunout,*) 'nbcf_in_inca =',nbcf_in_inca
802    WRITE(lunout,*) 'nbcf_in_ant  =',nbcf_in_ant
[3387]803
804    IF (nbcf_in.GT.0) THEN
805        DO iq=1,nbcf_in
806          CALL bcast(cfname_in(iq))
807          CALL bcast(cftext_in(iq))
808          CALL bcast(cfunits_in(iq))
809        ENDDO
810    ENDIF
811
812    IF (nbcf_out.GT.0) THEN
813        DO iq=1,nbcf_out
814          CALL bcast(cfname_out(iq))
815          CALL bcast(cftext_out(iq))
816          CALL bcast(cfunits_out(iq))
817        ENDDO
818    ENDIF
819
820    IF (nbcf.GT.0) THEN
821        DO iq=1,nbcf
822          CALL bcast(cfmod1(iq))
823          CALL bcast(cfmod2(iq))
824        ENDDO
825    ENDIF
826
[3390]827    IF (nbcf_in.GT.0)  WRITE(lunout,*)'infocfields_mod --- cfname_in: ',cfname_in
[3387]828    IF (nbcf_out.GT.0) WRITE(lunout,*)'infocfields_mod --- cfname_out: ',cfname_out
829
[3390]830    IF (nbcf_in.GT.0)  WRITE(lunout,*)'infocfields_mod --- cftext_in: ',cftext_in
[3387]831    IF (nbcf_out.GT.0) WRITE(lunout,*)'infocfields_mod --- cftext_out: ',cftext_out
832
833    IF (nbcf.GT.0) WRITE(lunout,*)'infocfields_mod --- cfmod1: ',cfmod1
834    IF (nbcf.GT.0) WRITE(lunout,*)'infocfields_mod --- cfmod2: ',cfmod2
835
[3390]836    IF (nbcf_in.GT.0)  WRITE(lunout,*)'infocfunits_mod --- cfunits_in: ',cfunits_in
[3387]837    IF (nbcf_out.GT.0) WRITE(lunout,*)'infocfunits_mod --- cfunits_out: ',cfunits_out
838
[3390]839    IF (nbcf_in.GT.0)  WRITE(*,*)'infocfields_init --- number of fields in to LMDZ: ',nbcf_in
[3387]840    IF (nbcf_out.GT.0) WRITE(*,*)'infocfields_init --- number of fields out of LMDZ: ',nbcf_out
841
842 ELSE
843 ! Default values for other planets
844    nbcf=0
845    nbcf_in=0
846    nbcf_out=0
847 ENDIF ! planet_type
848
[3391]849 ALLOCATE(fields_in(klon,nbcf_in),stat=error)
[3435]850 IF (error /= 0)  CALL abort_physic(modname,'Pb in allocation fields_in',1)
[3391]851 ALLOCATE(yfields_in(klon,nbcf_in),stat=error)
[3435]852 IF (error /= 0)  CALL abort_physic(modname,'Pb in allocation yfields_in',1)
[3391]853 ALLOCATE(fields_out(klon,nbcf_out),stat=error)
[3435]854 IF (error /= 0)  CALL abort_physic(modname,'Pb in allocation fields_out',1)
[3391]855 ALLOCATE(yfields_out(klon,nbcf_out),stat=error)
[3435]856 IF (error /= 0)  CALL abort_physic(modname,'Pb in allocation yfields_out',1)
[3387]857
858END SUBROUTINE infocfields_init
859
[1227]860END MODULE carbon_cycle_mod
Note: See TracBrowser for help on using the repository browser.