source: LMDZ4/trunk/libf/phylmd/carbon_cycle_mod.F90 @ 5426

Last change on this file since 5426 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
File size: 14.5 KB
Line 
1MODULE carbon_cycle_mod
2
3! Author : Josefine GHATTAS, Patricia CADULE
4
5  IMPLICIT NONE
6  SAVE
7  PRIVATE
8  PUBLIC :: carbon_cycle_init, carbon_cycle
9
10! Variables read from parmeter file physiq.def
11  LOGICAL, PUBLIC :: carbon_cycle_tr        ! 3D transport of CO2 in the atmosphere, parameter read in conf_phys
12!$OMP THREADPRIVATE(carbon_cycle_tr)
13  LOGICAL, PUBLIC :: carbon_cycle_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES)
14!$OMP THREADPRIVATE(carbon_cycle_cpl)
15  LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible
16
17! Scalare values when no transport, from physiq.def
18  REAL :: fos_fuel_s  ! carbon_cycle_fos_fuel dans physiq.def
19!$OMP THREADPRIVATE(fos_fuel_s)
20  REAL :: emis_land_s ! not yet implemented
21!$OMP THREADPRIVATE(emis_land_s)
22
23  INTEGER :: ntr_co2                ! Number of tracers concerning the carbon cycle
24  INTEGER :: id_fco2_tot            ! Tracer index
25  INTEGER :: id_fco2_ocn            !  - " -
26  INTEGER :: id_fco2_land           !  - " -
27  INTEGER :: id_fco2_land_use       !  - " -
28  INTEGER :: id_fco2_fos_fuel       !  - " -
29!$OMP THREADPRIVATE(ntr_co2, id_fco2_tot, id_fco2_ocn, id_fco2_land, id_fco2_land_use, id_fco2_fos_fuel) 
30
31  REAL, DIMENSION(:), ALLOCATABLE :: fos_fuel        ! CO2 fossil fuel emission from file [gC/m2/d]
32!$OMP THREADPRIVATE(fos_fuel)
33  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day ! flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]
34!$OMP THREADPRIVATE(fco2_ocn_day)
35  REAL, DIMENSION(:), ALLOCATABLE :: fco2_land_day   ! flux CO2 from land for 1 day (cumulated)  [gC/m2/d]
36!$OMP THREADPRIVATE(fco2_land_day)
37  REAL, DIMENSION(:), ALLOCATABLE :: fco2_lu_day     ! Emission from land use change for 1 day (cumulated) [gC/m2/d]
38!$OMP THREADPRIVATE(fco2_lu_day)
39
40! Following 2 fields will be initialized in surf_land_orchidee at each time step
41  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_inst  ! flux CO2 from land at one time step
42!$OMP THREADPRIVATE(fco2_land_inst)
43  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_lu_inst    ! Emission from land use change at one time step
44!$OMP THREADPRIVATE(fco2_lu_inst)
45
46! Calculated co2 field to be send to the ocean via the coupler and to ORCHIDEE
47  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send
48!$OMP THREADPRIVATE(co2_send)
49
50CONTAINS
51 
52  SUBROUTINE carbon_cycle_init(tr_seri, aerosol, radio)
53    USE dimphy
54    USE infotrac
55    USE IOIPSL
56    USE surface_data, ONLY : ok_veget, type_ocean
57
58    IMPLICIT NONE
59    INCLUDE "clesphys.h"
60    INCLUDE "iniprint.h"
61 
62! Input argument
63    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA] 
64
65! InOutput arguments
66    LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: aerosol
67    LOGICAL,DIMENSION(nbtr), INTENT(INOUT) :: radio
68
69! Local variables
70    INTEGER               :: ierr, it, iiq
71    REAL, DIMENSION(klon) :: tr_seri_sum
72
73
74! 0) Test for compatibility
75    IF (carbon_cycle_cpl .AND. type_ocean/='couple') &
76         CALL abort_gcm('carbon_cycle_init', 'Coupling with ocean model is needed for carbon_cycle_cpl',1)
77    IF (carbon_cycle_cpl .AND..NOT. ok_veget) &
78         CALL abort_gcm('carbon_cycle_init', 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl',1)
79
80
81! 1) Check if transport of one tracer flux CO2 or 4 separated tracers
82    IF (carbon_cycle_tr) THEN
83       id_fco2_tot=0
84       id_fco2_ocn=0
85       id_fco2_land=0
86       id_fco2_land_use=0
87       id_fco2_fos_fuel=0
88       
89       ! Search in tracer list
90       DO it=1,nbtr
91          iiq=niadv(it+2)
92          IF (tname(iiq) == "fCO2" ) THEN
93             id_fco2_tot=it
94          ELSE IF (tname(iiq) == "fCO2_ocn" ) THEN
95             id_fco2_ocn=it
96          ELSE IF (tname(iiq) == "fCO2_land" ) THEN
97             id_fco2_land=it
98          ELSE IF (tname(iiq) == "fCO2_land_use" ) THEN
99             id_fco2_land_use=it
100          ELSE IF (tname(iiq) == "fCO2_fos_fuel" ) THEN
101             id_fco2_fos_fuel=it
102          END IF
103       END DO
104
105       ! Count tracers found
106       IF (id_fco2_tot /= 0 .AND. &
107            id_fco2_ocn==0 .AND. id_fco2_land==0 .AND. id_fco2_land_use==0 .AND. id_fco2_fos_fuel==0) THEN
108         
109          ! transport  1 tracer flux CO2
110          ntr_co2 = 1
111         
112       ELSE IF (id_fco2_tot==0 .AND. &
113            id_fco2_ocn /=0 .AND. id_fco2_land/=0 .AND. id_fco2_land_use/=0 .AND. id_fco2_fos_fuel/=0) THEN
114         
115          ! transport 4 tracers seperatively
116          ntr_co2 = 4
117         
118       ELSE
119          CALL abort_gcm('carbon_cycle_init', 'error in coherence between traceur.def and gcm.def',1)
120       END IF
121
122       ! Definition of control varaiables for the tracers
123       DO it=1,nbtr
124          IF (it==id_fco2_tot .OR. it==id_fco2_ocn .OR. it==id_fco2_land .OR. &
125               it==id_fco2_land_use .OR. it==id_fco2_fos_fuel) THEN
126             aerosol(it) = .FALSE.
127             radio(it)   = .FALSE.
128          END IF
129       END DO
130
131    ELSE
132       ! No transport of CO2
133       ntr_co2 = 0
134    END IF ! carbon_cycle_tr
135
136
137! 2) Allocate variable for CO2 fossil fuel emission
138    IF (carbon_cycle_tr) THEN
139       ! Allocate 2D variable
140       ALLOCATE(fos_fuel(klon), stat=ierr)
141       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 1',1)
142    ELSE
143       ! No transport : read value from .def
144       fos_fuel_s = 0.
145       CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s)
146       WRITE(lunout,*) 'carbon_cycle_fos_fuel = ', fos_fuel_s
147    END IF
148
149
150! 3) Allocate and initialize fluxes
151    IF (carbon_cycle_cpl) THEN
152       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
153       ALLOCATE(fco2_land_day(klon), stat=ierr)
154       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
155       ALLOCATE(fco2_lu_day(klon), stat=ierr)
156       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 4',1)
157
158       fco2_land_day(:) = 0.  ! JG : Doit prend valeur de restart
159       fco2_lu_day(:)   = 0.  ! JG : Doit prend valeur de restart
160
161       ! fco2_ocn_day is allocated in cpl_mod
162       ! fco2_land_inst and fco2_lu_inst are allocated in surf_land_orchidee
163       
164       ALLOCATE(co2_send(klon), stat=ierr)
165       IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 7',1)
166       
167       ! Calculate using restart tracer values
168       IF (carbon_cycle_tr) THEN
169          IF (ntr_co2==1) THEN
170             co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
171          ELSE ! ntr_co2==4
172             ! Calculate the delta CO2 flux
173             tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
174                  tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
175             co2_send(:) = tr_seri_sum(:) + co2_ppm0
176          END IF
177       ELSE
178          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
179          co2_send(:) = co2_ppm
180       END IF
181
182
183    ELSE
184       IF (carbon_cycle_tr) THEN
185          ! No coupling of CO2 fields :
186          ! corresponding fields will instead be read from files
187          ALLOCATE(fco2_ocn_day(klon), stat=ierr)
188          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 8',1)
189          ALLOCATE(fco2_land_day(klon), stat=ierr)
190          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 9',1)
191          ALLOCATE(fco2_lu_day(klon), stat=ierr)
192          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 10',1)       
193       END IF
194    END IF
195
196! 4) Read parmeter for calculation of emission compatible
197    IF (.NOT. carbon_cycle_tr) THEN
198       carbon_cycle_emis_comp=.FALSE.
199       CALL getin('carbon_cycle_emis_comp',carbon_cycle_emis_comp)
200       WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
201    END IF
202
203  END SUBROUTINE carbon_cycle_init
204
205!
206!
207!
208
209  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
210   
211    USE infotrac
212    USE dimphy
213    USE mod_phys_lmdz_transfert_para, ONLY : reduce_sum
214    USE phys_cal_mod, ONLY : mth_cur, mth_len
215    USE phys_cal_mod, ONLY : day_cur
216    USE comgeomphy
217
218    IMPLICIT NONE
219
220    INCLUDE "clesphys.h"
221    INCLUDE "indicesol.h"
222
223! In/Output arguments
224    INTEGER,INTENT(IN) :: nstep      ! time step in physiq
225    REAL,INTENT(IN)    :: pdtphys    ! length of time step in physiq (sec)
226    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
227    REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri
228
229! Local variables
230    LOGICAL :: newmonth ! indicates if a new month just started
231    LOGICAL :: newday   ! indicates if a new day just started
232    LOGICAL :: endday   ! indicated if last time step in a day
233
234    REAL, PARAMETER :: fact=1.E-15/2.12  ! transformation factor from gC/m2/day => ppm/m2/day
235    REAL, DIMENSION(klon) :: fco2_tmp, tr_seri_sum
236    REAL :: sumtmp
237    REAL :: airetot     ! Total area the earth
238    REAL :: delta_co2_ppm
239   
240! -) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
241
242    newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.
243
244    IF (MOD(nstep,INT(86400./pdtphys))==1) newday=.TRUE.
245    IF (MOD(nstep,INT(86400./pdtphys))==0) endday=.TRUE.
246    IF (newday .AND. day_cur==1) newmonth=.TRUE.
247   
248! -) Read new maps if new month started
249    IF (newmonth .AND. carbon_cycle_tr) THEN
250       CALL read_map2D('fossil_fuel.nc','fos_fuel', mth_cur, .FALSE., fos_fuel)
251       
252       ! division by month lenght to get dayly value
253       fos_fuel(:) = fos_fuel(:)/mth_len
254       
255       IF (.NOT. carbon_cycle_cpl) THEN
256          ! Get dayly values from monthly fluxes
257          CALL read_map2D('fl_co2_ocean.nc','CO2_OCN',mth_cur,.FALSE.,fco2_ocn_day)
258          CALL read_map2D('fl_co2_land.nc','CO2_LAND', mth_cur,.FALSE.,fco2_land_day)
259          CALL read_map2D('fl_co2_land_use.nc','CO2_LAND_USE',mth_cur,.FALSE.,fco2_lu_day)
260       END IF
261    END IF
262   
263
264
265! -) Update tracers at beginning of a new day. Beginning of a new day correspond to a new coupling period in cpl_mod.
266    IF (newday) THEN
267
268       IF (carbon_cycle_tr) THEN
269
270          ! Update tracers
271          IF (ntr_co2 == 1) THEN
272             ! Calculate the new flux CO2
273             tr_seri(:,1,id_fco2_tot) = tr_seri(:,1,id_fco2_tot) + &
274                  (fos_fuel(:) + &
275                  fco2_lu_day(:)  * pctsrf(:,is_ter) + &
276                  fco2_land_day(:)* pctsrf(:,is_ter) + &
277                  fco2_ocn_day(:) * pctsrf(:,is_oce)) * fact
278
279          ELSE ! ntr_co2 == 4
280             tr_seri(:,1,id_fco2_fos_fuel)  = tr_seri(:,1,id_fco2_fos_fuel) + fos_fuel(:) * fact ! [ppm/m2/day]
281
282             tr_seri(:,1,id_fco2_land_use)  = tr_seri(:,1,id_fco2_land_use) + &
283                  fco2_lu_day(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
284
285             tr_seri(:,1,id_fco2_land)      = tr_seri(:,1,id_fco2_land)     + &
286                  fco2_land_day(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
287
288             tr_seri(:,1,id_fco2_ocn)       = tr_seri(:,1,id_fco2_ocn)      + &
289                  fco2_ocn_day(:) *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
290
291          END IF
292
293       ELSE ! no transport
294          IF (carbon_cycle_cpl) THEN
295             IF (carbon_cycle_emis_comp) THEN
296                ! Calcul emission compatible a partir des champs 2D et co2_ppm
297                !!! TO DO!!
298                CALL abort_gcm('carbon_cycle', ' Option carbon_cycle_emis_comp not yet implemented',1)
299             END IF
300          END IF
301
302       END IF ! carbon_cycle_tr
303   
304       ! Reset cumluative variables
305       IF (carbon_cycle_cpl) THEN
306          fco2_land_day(:) = 0.
307          fco2_lu_day(:)   = 0.
308       END IF
309   
310    END IF ! newday
311       
312
313
314! -) Cumulate fluxes from ORCHIDEE at each timestep
315    IF (carbon_cycle_cpl) THEN
316       fco2_land_day(:) = fco2_land_day(:) + fco2_land_inst(:)
317       fco2_lu_day(:)   = fco2_lu_day(:)   + fco2_lu_inst(:)
318    END IF
319
320
321
322! -)  At the end of a new day, calculate a mean scalare value of CO2 to be used by
323!     the radiation scheme (instead of reading value from .def)
324
325! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
326    IF (endday) THEN
327       ! Calculte total area of the earth surface
328       CALL reduce_sum(SUM(airephy),airetot)
329       
330
331       IF (carbon_cycle_tr) THEN
332          IF (ntr_co2 == 1) THEN
333         
334             ! Calculate mean value of tracer CO2 to get an scalare value to be used in the
335             ! radiation scheme (instead of reading value from .def)
336             ! Mean value weighted with the grid cell area
337             
338             ! Calculate mean value
339             fco2_tmp(:) = tr_seri(:,1,id_fco2_tot) * airephy(:)
340             CALL reduce_sum(SUM(fco2_tmp),sumtmp)
341             co2_ppm = sumtmp/airetot + co2_ppm0
342             
343          ELSE ! ntr_co2 == 4
344             
345             ! Calculate the delta CO2 flux
346             tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
347                  tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
348             
349             ! Calculate mean value of delta CO2 flux
350             fco2_tmp(:) = tr_seri_sum(:) * airephy(:)
351             CALL reduce_sum(SUM(fco2_tmp),sumtmp)
352             delta_co2_ppm = sumtmp/airetot
353             
354             ! Add initial value for co2_ppm to delta value
355             co2_ppm = delta_co2_ppm + co2_ppm0
356          END IF
357
358       ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
359
360          ! Calculate the total CO2 flux and integrate to get scalare value for the radiation scheme
361          fco2_tmp(:) = (fos_fuel(:) + (fco2_lu_day(:) + fco2_land_day(:))*pctsrf(:,is_ter) &
362               + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
363         
364          ! Calculate mean value
365          fco2_tmp(:) = fco2_tmp(:) * airephy(:)
366          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
367          delta_co2_ppm = sumtmp/airetot
368
369          ! Update current value of the atmospheric co2_ppm
370          co2_ppm = co2_ppm + delta_co2_ppm
371         
372       END IF ! carbon_cycle_tr
373
374       ! transformation of the atmospheric CO2 concentration for the radiation code
375       RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
376
377    END IF
378
379    ! Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
380    IF (endday .AND. carbon_cycle_cpl) THEN
381       
382       IF (carbon_cycle_tr) THEN
383          IF (ntr_co2==1) THEN
384
385             co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
386
387          ELSE ! ntr_co2 == 4
388
389             co2_send(:) = tr_seri_sum(:) + co2_ppm0
390
391          END IF
392       ELSE
393          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
394          co2_send(:) = co2_ppm
395       END IF
396
397    END IF
398
399  END SUBROUTINE carbon_cycle
400 
401END MODULE carbon_cycle_mod
Note: See TracBrowser for help on using the repository browser.