Context Navigation

carbon_cycle_mod.F90 @ 1763

Last change on this file since 1763 was 1759, checked in by Ehouarn Millour, 12 years ago
IOIPSL routine getin is not threadsafe, so when running in OpenMP, it should be called by only one thread (and the result copied to other threads in the case of threadprivate variables). EM
File size: 18.6 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 comgeomphy
89	USE mod_phys_lmdz_transfert_para
90	USE infotrac
91	USE IOIPSL
92	USE surface_data, ONLY : ok_veget, type_ocean
93	USE phys_cal_mod, ONLY : mth_len
94
95	IMPLICIT NONE
96	INCLUDE "clesphys.h"
97	INCLUDE "iniprint.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_gcm('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)
160
161	SELECT CASE(tname(iiq))
162	CASE("fCO2_ocn")
163	itc = itc + 1
164	co2trac(itc)%name='fCO2_ocn'
165	co2trac(itc)%id=it
166	co2trac(itc)%file='fl_co2_ocean.nc'
167	IF (carbon_cycle_cpl .AND. type_ocean=='couple') THEN
168	co2trac(itc)%cpl=.TRUE.
169	co2trac(itc)%updatefreq = 86400 ! Once a day as the coupling with OASIS/PISCES
170	ELSE
171	co2trac(itc)%cpl=.FALSE.
172	co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
173	END IF
174	CASE("fCO2_land")
175	itc = itc + 1
176	co2trac(itc)%name='fCO2_land'
177	co2trac(itc)%id=it
178	co2trac(itc)%file='fl_co2_land.nc'
179	IF (carbon_cycle_cpl .AND. ok_veget) THEN
180	co2trac(itc)%cpl=.TRUE.
181	co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
182	ELSE
183	co2trac(itc)%cpl=.FALSE.
184	! co2trac(itc)%updatefreq = 10800 ! 10800sec = 3H
185	co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
186	END IF
187	CASE("fCO2_land_use")
188	itc = itc + 1
189	co2trac(itc)%name='fCO2_land_use'
190	co2trac(itc)%id=it
191	co2trac(itc)%file='fl_co2_land_use.nc'
192	IF (carbon_cycle_cpl .AND. ok_veget) THEN
193	co2trac(it)%cpl=.TRUE.
194	co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
195	ELSE
196	co2trac(itc)%cpl=.FALSE.
197	co2trac(itc)%updatefreq = 10800 ! 10800sec = 3H
198	END IF
199	CASE("fCO2_fos_fuel")
200	itc = itc + 1
201	co2trac(itc)%name='fCO2_fos_fuel'
202	co2trac(itc)%id=it
203	co2trac(itc)%file='fossil_fuel.nc'
204	co2trac(itc)%cpl=.FALSE. ! This tracer always read from file
205	! co2trac(itc)%updatefreq = 86400 ! 86400sec = 24H Cadule case
206	co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
207	CASE("fCO2_bbg")
208	itc = itc + 1
209	co2trac(itc)%name='fCO2_bbg'
210	co2trac(itc)%id=it
211	co2trac(itc)%file='fl_co2_bbg.nc'
212	co2trac(itc)%cpl=.FALSE. ! This tracer always read from file
213	co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
214	CASE("fCO2")
215	! fCO2 : One tracer transporting the total CO2 flux
216	itc = itc + 1
217	co2trac(itc)%name='fCO2'
218	co2trac(itc)%id=it
219	co2trac(itc)%file='fl_co2.nc'
220	IF (carbon_cycle_cpl) THEN
221	co2trac(itc)%cpl=.TRUE.
222	ELSE
223	co2trac(itc)%cpl=.FALSE.
224	END IF
225	co2trac(itc)%updatefreq = 86400
226	! DOES THIS WORK ???? Problematic due to implementation of the coupled fluxes...
227	CALL abort_gcm('carbon_cycle_init','transport of total CO2 has to be implemented and tested',1)
228	END SELECT
229	END DO
230
231	! Total number of carbon CO2 tracers
232	ntr_co2 = itc
233
234	! Definition of control varaiables for the tracers
235	DO it=1,ntr_co2
236	aerosol(co2trac(it)%id) = .FALSE.
237	radio(co2trac(it)%id) = .FALSE.
238	END DO
239
240	! Vector indicating which timestep to read for each tracer
241	! Always start read in the beginning of the file
242	co2trac(:)%readstep = 0
243
244
245	! 3) Allocate variables
246	! ---------------------
247	! Allocate vector for storing fluxes to inject
248	ALLOCATE(dtr_add(klon,maxco2trac), stat=ierr)
249	IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 11',1)
250
251	! Allocate variables for cumulating fluxes from ORCHIDEE
252	IF (RCO2_inter) THEN
253	IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
254	ALLOCATE(fco2_land_day(klon), stat=ierr)
255	IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
256	fco2_land_day(1:klon) = 0.
257
258	ALLOCATE(fco2_lu_day(klon), stat=ierr)
259	IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
260	fco2_lu_day(1:klon) = 0.
261	END IF
262	END IF
263
264
265	! 4) Test for compatibility
266	! -------------------------
267	! IF (carbon_cycle_cpl .AND. type_ocean/='couple') THEN
268	! WRITE(lunout,*) 'Coupling with ocean model is needed for carbon_cycle_cpl'
269	! CALL abort_gcm('carbon_cycle_init', 'coupled ocean is needed for carbon_cycle_cpl',1)
270	! END IF
271	!
272	! IF (carbon_cycle_cpl .AND..NOT. ok_veget) THEN
273	! WRITE(lunout,*) 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl'
274	! CALL abort_gcm('carbon_cycle_init', 'ok_veget is needed for carbon_cycle_cpl',1)
275	! END IF
276
277	! Compiler test : following should never happen
278	teststop=0
279	DO it=1,teststop
280	CALL abort_gcm('carbon_cycle_init', 'Entering loop from 1 to 0',1)
281	END DO
282
283	IF (ntr_co2==0) THEN
284	! No carbon tracers found in tracer.def. It is not possible to do carbon cycle
285	WRITE(lunout,*) 'No carbon tracers found in tracer.def. Not ok with carbon_cycle_tr and/or carbon_cycle_cp'
286	CALL abort_gcm('carbon_cycle_init', 'No carbon tracers found in tracer.def',1)
287	END IF
288
289	! 5) Calculate total area of the earth surface
290	! --------------------------------------------
291	CALL reduce_sum(SUM(airephy),airetot)
292	CALL bcast(airetot)
293
294	END SUBROUTINE carbon_cycle_init
295
296
297	SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
298	! Subroutine for injection of co2 in the tracers
299	!
300	! - Find out if it is time to update
301	! - Get tracer from coupled model or from file
302	! - Calculate new RCO2 value for the radiation scheme
303	! - Calculate CO2 flux to send to ocean and land models (PISCES and ORCHIDEE)
304
305	USE infotrac
306	USE dimphy
307	USE mod_phys_lmdz_transfert_para
308	USE phys_cal_mod, ONLY : mth_cur, mth_len
309	USE phys_cal_mod, ONLY : day_cur
310	USE comgeomphy
311
312	IMPLICIT NONE
313
314	INCLUDE "clesphys.h"
315	INCLUDE "indicesol.h"
316	INCLUDE "iniprint.h"
317	INCLUDE "YOMCST.h"
318
319	! In/Output arguments
320	INTEGER,INTENT(IN) :: nstep ! time step in physiq
321	REAL,INTENT(IN) :: pdtphys ! length of time step in physiq (sec)
322	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Surface fraction
323	REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT) :: tr_seri ! All tracers
324	REAL, DIMENSION(klon,nbtr), INTENT(INOUT) :: source ! Source for all tracers
325
326	! Local variables
327	INTEGER :: it
328	LOGICAL :: newmonth ! indicates if a new month just started
329	LOGICAL :: newday ! indicates if a new day just started
330	LOGICAL :: endday ! indicated if last time step in a day
331
332	REAL, PARAMETER :: fact=1.E-15/2.12 ! transformation factor from gC/m2/day => ppm/m2/day
333	REAL, DIMENSION(klon) :: fco2_tmp
334	REAL :: sumtmp
335	REAL :: delta_co2_ppm
336
337
338	! 1) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
339	! -------------------------------------------------------------------------------------------------------
340
341	newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.
342
343	IF (MOD(nstep,INT(86400./pdtphys))==1) newday=.TRUE.
344	IF (MOD(nstep,INT(86400./pdtphys))==0) endday=.TRUE.
345	IF (newday .AND. day_cur==1) newmonth=.TRUE.
346
347	! 2) For each carbon tracer find out if it is time to inject (update)
348	! --------------------------------------------------------------------
349	DO it = 1, ntr_co2
350	IF ( MOD(nstep,INT(co2trac(it)%updatefreq/pdtphys)) == 1 ) THEN
351	co2trac(it)%updatenow = .TRUE.
352	ELSE
353	co2trac(it)%updatenow = .FALSE.
354	END IF
355	END DO
356
357	! 3) Get tracer update
358	! --------------------------------------
359	DO it = 1, ntr_co2
360	IF ( co2trac(it)%updatenow ) THEN
361	IF ( co2trac(it)%cpl ) THEN
362	! Get tracer from coupled model
363	SELECT CASE(co2trac(it)%name)
364	CASE('fCO2_land') ! from ORCHIDEE
365	dtr_add(:,it) = fco2_land_inst(:)pctsrf(:,is_ter)fact ! [ppm/m2/day]
366	CASE('fCO2_land_use') ! from ORCHIDEE
367	dtr_add(:,it) = fco2_lu_inst(:) pctsrf(:,is_ter)fact ! [ppm/m2/day]
368	CASE('fCO2_ocn') ! from PISCES
369	dtr_add(:,it) = fco2_ocn_day(:) pctsrf(:,is_oce)fact ! [ppm/m2/day]
370	CASE DEFAULT
371	WRITE(lunout,*) 'Error with tracer ',co2trac(it)%name
372	CALL abort_gcm('carbon_cycle', 'No coupling implemented for this tracer',1)
373	END SELECT
374	ELSE
375	! Read tracer from file
376	co2trac(it)%readstep = co2trac(it)%readstep + 1 ! increment time step in file
377	! Patricia CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.FALSE.,dtr_add(:,it))
378	CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.TRUE.,dtr_add(:,it))
379
380	! Converte from kgC/m2/h to kgC/m2/s
381	dtr_add(:,it) = dtr_add(:,it)/3600
382	! Add individual treatment of values read from file
383	SELECT CASE(co2trac(it)%name)
384	CASE('fCO2_land')
385	dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
386	CASE('fCO2_land_use')
387	dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
388	CASE('fCO2_ocn')
389	dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_oce)
390	! Patricia :
391	! CASE('fCO2_fos_fuel')
392	! dtr_add(:,it) = dtr_add(:,it)/mth_len
393	! co2trac(it)%readstep = 0 ! Always read same value for fossil fuel(Cadule case)
394	END SELECT
395	END IF
396	END IF
397	END DO
398
399	! 4) Update co2 tracers :
400	! Loop over all carbon tracers and add source
401	! ------------------------------------------------------------------
402	IF (carbon_cycle_tr) THEN
403	DO it = 1, ntr_co2
404	IF (.FALSE.) THEN
405	tr_seri(1:klon,1,co2trac(it)%id) = tr_seri(1:klon,1,co2trac(it)%id) + dtr_add(1:klon,it)
406	source(1:klon,co2trac(it)%id) = 0.
407	ELSE
408	source(1:klon,co2trac(it)%id) = dtr_add(1:klon,it)
409	END IF
410	END DO
411	END IF
412
413
414	! 5) Calculations for new CO2 value for the radiation scheme(instead of reading value from .def)
415	! ----------------------------------------------------------------------------------------------
416	IF (RCO2_inter) THEN
417	! Cumulate fluxes from ORCHIDEE at each timestep
418	IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
419	IF (newday) THEN ! Reset cumulative variables once a day
420	fco2_land_day(1:klon) = 0.
421	fco2_lu_day(1:klon) = 0.
422	END IF
423	fco2_land_day(1:klon) = fco2_land_day(1:klon) + fco2_land_inst(1:klon) ![gC/m2/day]
424	fco2_lu_day(1:klon) = fco2_lu_day(1:klon) + fco2_lu_inst(1:klon) ![gC/m2/day]
425	END IF
426
427	! At the end of a new day, calculate a mean scalare value of CO2
428	! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
429	IF (endday) THEN
430
431	IF (carbon_cycle_tr) THEN
432	! Sum all co2 tracers to get the total delta CO2 flux
433	fco2_tmp(:) = 0.
434	DO it = 1, ntr_co2
435	fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
436	END DO
437
438	ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
439	! Sum co2 fluxes comming from coupled models and parameter for fossil fuel
440	fco2_tmp(1:klon) = fos_fuel_s + ((fco2_lu_day(1:klon) + fco2_land_day(1:klon))*pctsrf(1:klon,is_ter) &
441	+ fco2_ocn_day(:)pctsrf(:,is_oce)) fact
442	END IF
443
444	! Calculate a global mean value of delta CO2 flux
445	fco2_tmp(1:klon) = fco2_tmp(1:klon) * airephy(1:klon)
446	CALL reduce_sum(SUM(fco2_tmp),sumtmp)
447	CALL bcast(sumtmp)
448	delta_co2_ppm = sumtmp/airetot
449
450	! Add initial value for co2_ppm and delta value
451	co2_ppm = co2_ppm0 + delta_co2_ppm
452
453	! Transformation of atmospheric CO2 concentration for the radiation code
454	RCO2 = co2_ppm * 1.0e-06 * 44.011/28.97
455
456	WRITE(lunout,*) 'RCO2 is now updated! RCO2 = ', RCO2
457	END IF ! endday
458
459	END IF ! RCO2_inter
460
461
462	! 6) Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE
463	! ----------------------------------------------------------------------------
464	IF (carbon_cycle_cpl) THEN
465
466	IF (carbon_cycle_tr) THEN
467	! Sum all co2 tracers to get the total delta CO2 flux at first model layer
468	fco2_tmp(:) = 0.
469	DO it = 1, ntr_co2
470	fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
471	END DO
472	co2_send(1:klon) = fco2_tmp(1:klon) + co2_ppm0
473	ELSE
474	! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
475	co2_send(1:klon) = co2_ppm
476	END IF
477
478	END IF
479
480	END SUBROUTINE carbon_cycle
481
482	END MODULE carbon_cycle_mod

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

Download in other formats:

Original Format