source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/carbon_cycle_mod.F90 @ 1240

Last change on this file since 1240 was 1227, checked in by jghattas, 16 years ago
  • Inclusion d'un premier version du cycle de carbon dans LMDZ. Attention

!! Il s'agit d'un version ou les nouveaux cles cycle_carbon_tr et
cycle_carbon_cpl ne sont pas teste. Avec les ancinenes parametres le
modele donne les memes resultats qu'avant. L'interface avec ORCHIDEE n'a
pas encore etait modifie.

  • physiq.F, phys_cal_mod.F90 : ajout d'un nouveau module qui contient qq parametres pour le calendrier et le pas de temps acutelle de la physiq. Ce module pourrait etre elargie plus tard / LF + JG


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