source: LMDZ5/branches/testing/libf/phylmd/carbon_cycle_mod.F90 @ 1864

Last change on this file since 1864 was 1795, checked in by Ehouarn Millour, 11 years ago

Version testing basee sur la r1794


Testing release based on r1794

File size: 18.6 KB
Line 
1MODULE 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
78CONTAINS
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    USE indice_sol_mod
312
313    IMPLICIT NONE
314
315    INCLUDE "clesphys.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 
482END MODULE carbon_cycle_mod
Note: See TracBrowser for help on using the repository browser.