Context Navigation

carbon_cycle_mod.F90 @ 3271

Last change on this file since 3271 was 1444, checked in by jghattas, 14 years ago

Modifications for carbon tracers :

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

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

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

Download in other formats:

Original Format