source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phylmd/carbon_cycle_mod.F90 @ 5412

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

Transformation de l'include indicesol.h en un module indice_sol_mod et modification des appels dans tous les fichiers concernés.
Aucun changement des résultats ni des sorties du modèle vs 1784.
UG

...................................................

Replacement of the indicesol.h include by a module named indice_sol_mod. Modification of the calls in every affected files.
Results and outputs of simulations are unchanged in comparison with rev 1784.
UG

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.