Context Navigation

carbon_cycle_mod.F90

Last change on this file was 3331, checked in by acozic, 7 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 18.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format