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

Last change on this file since 3513 was 3459, checked in by oboucher, 6 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
Line 
1MODULE carbon_cycle_mod
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 :
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)
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!=======================================================================
33
34  IMPLICIT NONE
35  SAVE
36  PRIVATE
37  PUBLIC :: carbon_cycle_init, carbon_cycle, infocfields_init
38
39! Variables read from parmeter file physiq.def
40  LOGICAL, PUBLIC :: carbon_cycle_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES)
41!$OMP THREADPRIVATE(carbon_cycle_cpl)
42  LOGICAL, PUBLIC :: carbon_cycle_tr        ! 3D transport of CO2 in the atmosphere, parameter read in conf_phys
43!$OMP THREADPRIVATE(carbon_cycle_tr)
44  LOGICAL, PUBLIC :: carbon_cycle_rad       ! CO2 interactive radiatively
45!$OMP THREADPRIVATE(carbon_cycle_rad)
46  INTEGER, PUBLIC :: level_coupling_esm ! Level of coupling for the ESM - 0, 1, 2, 3
47!$OMP THREADPRIVATE(level_coupling_esm)
48  REAL, PUBLIC :: RCO2_glo, RCO2_tot
49!$OMP THREADPRIVATE(RCO2_glo, RCO2_tot)
50
51  LOGICAL :: carbon_cycle_emis_comp_omp=.FALSE.
52  LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible
53!$OMP THREADPRIVATE(carbon_cycle_emis_comp)
54
55  LOGICAL :: RCO2_inter_omp
56  LOGICAL :: RCO2_inter  ! RCO2 interactive : if true calculate new value RCO2 for the radiation scheme
57!$OMP THREADPRIVATE(RCO2_inter)
58
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)
63  REAL :: emis_land_s ! not yet implemented
64!$OMP THREADPRIVATE(emis_land_s)
65
66  REAL :: airetot     ! Total area of the earth surface
67!$OMP THREADPRIVATE(airetot)
68
69  INTEGER :: ntr_co2  ! Number of tracers concerning the carbon cycle
70!$OMP THREADPRIVATE(ntr_co2)
71
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)
75
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)
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)
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)
88
89  REAL, DIMENSION(:,:), ALLOCATABLE :: dtr_add       ! Tracer concentration to be injected
90!$OMP THREADPRIVATE(dtr_add)
91
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)
97
98! Calculated co2 field to be send to the ocean via the coupler and to ORCHIDEE
99  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send ! Field allocated in phyetat0
100!$OMP THREADPRIVATE(co2_send)
101
102  INTEGER, PARAMETER, PUBLIC :: id_CO2=1              !--temporaire OB -- to be changed
103
104! nbfields : total number of fields
105  INTEGER, PUBLIC :: nbcf
106!$OMP THREADPRIVATE(nbcf)
107
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
160  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), PUBLIC :: field_out_names
161!$OMP THREADPRIVATE(field_out_names)
162
163  CHARACTER(LEN=20), ALLOCATABLE, DIMENSION(:), PUBLIC :: field_in_names
164!$OMP THREADPRIVATE(field_in_names)
165
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
178  TYPE, PUBLIC :: co2_trac_type
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
182     LOGICAL            :: cpl        ! True if this tracers is coupled from ORCHIDEE or PISCES.
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
191CONTAINS
192 
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
200    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
204    USE IOIPSL
205    USE surface_data, ONLY : ok_veget, type_ocean
206    USE phys_cal_mod, ONLY : mth_len
207    USE print_control_mod, ONLY: lunout
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] 
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
220! 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)
404
405  END SUBROUTINE carbon_cycle_init
406
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 
592  SUBROUTINE infocfields_init
593
594!    USE control_mod, ONLY: planet_type
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
669  CHARACTER(len=10),SAVE :: planet_type="earth"
670
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
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))
698
699             nbcf_in=0
700             nbcf_out=0
701
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)
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
728                ELSE
729                  WRITE(lunout,*) 'abort_gcm --- nbcf    : ',nbcf
730                  WRITE(lunout,*) 'abort_gcm --- nbcf_in : ',nbcf_in
731                  WRITE(lunout,*) 'abort_gcm --- nbcf_out: ',nbcf_out
732                  CALL abort_physic('infocfields_init','Problem in the definition of the coupling fields',1)
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
750    WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf    =',nbcf
751    WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf_in =',nbcf_in
752    WRITE(lunout,*) 'infocfields_mod.F90 --- nbcf_out=',nbcf_out
753
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))
763       
764    IF (is_mpi_root .AND. is_omp_root) THEN
765
766        IF (nbcf.GT.0)     cfname=cfname_root
767        IF (nbcf_in.GT.0)  cfname_in=PACK(cfname_root,mask_in_root)
768        IF (nbcf_out.GT.0) cfname_out=PACK(cfname_root,mask_out_root)
769        IF (nbcf_in.GT.0)  cftext_in=PACK(cftext_root,mask_in_root)
770        IF (nbcf_out.GT.0) cftext_out=PACK(cftext_root,mask_out_root)
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)
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
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 
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
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
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
827    IF (nbcf_in.GT.0)  WRITE(lunout,*)'infocfields_mod --- cfname_in: ',cfname_in
828    IF (nbcf_out.GT.0) WRITE(lunout,*)'infocfields_mod --- cfname_out: ',cfname_out
829
830    IF (nbcf_in.GT.0)  WRITE(lunout,*)'infocfields_mod --- cftext_in: ',cftext_in
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
836    IF (nbcf_in.GT.0)  WRITE(lunout,*)'infocfunits_mod --- cfunits_in: ',cfunits_in
837    IF (nbcf_out.GT.0) WRITE(lunout,*)'infocfunits_mod --- cfunits_out: ',cfunits_out
838
839    IF (nbcf_in.GT.0)  WRITE(*,*)'infocfields_init --- number of fields in to LMDZ: ',nbcf_in
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
849 ALLOCATE(fields_in(klon,nbcf_in),stat=error)
850 IF (error /= 0)  CALL abort_physic(modname,'Pb in allocation fields_in',1)
851 ALLOCATE(yfields_in(klon,nbcf_in),stat=error)
852 IF (error /= 0)  CALL abort_physic(modname,'Pb in allocation yfields_in',1)
853 ALLOCATE(fields_out(klon,nbcf_out),stat=error)
854 IF (error /= 0)  CALL abort_physic(modname,'Pb in allocation fields_out',1)
855 ALLOCATE(yfields_out(klon,nbcf_out),stat=error)
856 IF (error /= 0)  CALL abort_physic(modname,'Pb in allocation yfields_out',1)
857
858END SUBROUTINE infocfields_init
859
860END MODULE carbon_cycle_mod
Note: See TracBrowser for help on using the repository browser.