source: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/tracco2i_mod.F90 @ 3751

Last change on this file since 3751 was 3651, checked in by jghattas, 5 years ago

Correction pour tourner experience avec type_trac=co2i sur Jean-Zay : allocation necessaire sur tout les procs, sinon scatter se plain.

File size: 12.0 KB
Line 
1MODULE tracco2i_mod
2!
3! This module does the work for the interactive CO2 tracers
4! Authors: Patricia Cadule and Olivier Boucher
5!
6! Purpose and description:
7!  -----------------------
8! Main routine for the interactive carbon cycle
9! Gather all carbon fluxes and emissions from ORCHIDEE, PISCES and fossil fuel
10! Compute the net flux in source field which is used in phytrac
11! Compute global CO2 mixing ratio for radiation scheme if option is activated
12! Redistribute CO2 evenly over the atmosphere if transport is desactivated
13!
14CONTAINS
15
16  SUBROUTINE tracco2i_init()
17    ! This subroutine calls carbon_cycle_init needed to be done before first call to phys_output_write in physiq.
18    USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl
19
20    ! Initialize carbon_cycle_mod
21    IF (carbon_cycle_cpl) THEN
22       CALL carbon_cycle_init()
23    ENDIF
24
25  END SUBROUTINE tracco2i_init
26
27  SUBROUTINE tracco2i(pdtphys, debutphy, &
28       xlat, xlon, pphis, pphi, &
29       t_seri, pplay, paprs, tr_seri, source)
30
31    USE dimphy
32    USE infotrac_phy
33    USE geometry_mod, ONLY: cell_area
34    USE carbon_cycle_mod, ONLY: id_CO2, nbcf_in, fields_in, cfname_in
35    USE carbon_cycle_mod, ONLY: fco2_ocn_day, fco2_ff, fco2_bb, fco2_land, fco2_ocean
36    USE carbon_cycle_mod, ONLY: fco2_land_nbp, fco2_land_nep, fco2_land_fLuc
37    USE carbon_cycle_mod, ONLY: fco2_land_fwoodharvest, fco2_land_fHarvest
38    USE carbon_cycle_mod, ONLY: carbon_cycle_cpl, carbon_cycle_tr, carbon_cycle_rad, RCO2_glo, RCO2_tot
39    USE mod_grid_phy_lmdz
40    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
41    USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter
42    USE phys_cal_mod
43    USE phys_state_var_mod, ONLY: pctsrf
44    USE indice_sol_mod, ONLY: nbsrf, is_ter, is_lic, is_oce, is_sic
45
46    IMPLICIT NONE
47
48    INCLUDE "clesphys.h"
49    INCLUDE "YOMCST.h"
50
51! Input argument
52!---------------
53    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
54    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
55
56    REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point
57    REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point
58    REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! geopotentiel du sol
59    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel de chaque couche
60
61    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
62    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
63    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
64    REAL,DIMENSION(klon,nbtr),INTENT(INOUT):: source  ! flux de traceur [U/m2/s]
65
66! Output argument
67!----------------
68    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/kgA] 
69
70! Local variables
71!----------------
72
73    INTEGER                        :: it, k, i, nb
74    REAL, DIMENSION(klon,klev)     :: m_air     ! mass of air in every grid box [kg]
75    REAL, DIMENSION(klon_glo,klev) :: co2_glo   ! variable temporaire sur la grille global
76    REAL, DIMENSION(klon_glo,klev) :: m_air_glo ! variable temporaire sur la grille global
77
78    LOGICAL, SAVE :: check_fCO2_nbp_in_cfname
79!$OMP THREADPRIVATE(check_fCO2_nbp_in_cfname)
80    INTEGER, SAVE :: day_pre=-1
81!$OMP THREADPRIVATE(day_pre)
82
83    IF (is_mpi_root) THEN
84      PRINT *,'in tracco2i: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
85    ENDIF
86
87!--initialisation of CO2 field if not in restart file
88!--dirty way of doing, do it better later
89!--convert 280 ppm into kg CO2 / kg air
90    IF (debutphy) THEN
91
92! Initialisation de tr_seri(id_CO2) si pas initialise
93      IF (MAXVAL(tr_seri(:,:,id_CO2)).LT.1.e-15) THEN
94        tr_seri(:,:,id_CO2)=co2_ppm0*1.e-6/RMD*RMCO2 !--initialised from co2_ppm0 in rdem
95      ENDIF
96
97!--check if fCO2_nbp is in
98      check_fCO2_nbp_in_cfname=.FALSE.
99      DO nb=1, nbcf_in
100        IF (cfname_in(nb)=="fCO2_nbp") check_fCO2_nbp_in_cfname=.TRUE.
101      ENDDO
102
103    ENDIF
104
105!--calculate mass of air in every grid box in kg air
106    DO i=1, klon
107    DO k=1, klev
108      m_air(i,k)=(paprs(i,k)-paprs(i,k+1))/RG*cell_area(i)
109    ENDDO
110    ENDDO
111
112!--call CO2 emission routine
113!--co2bb is zero for now
114!--unit kg CO2 m-2 s-1
115    CALL co2_emissions(debutphy)
116
117!--retrieving land and ocean CO2 flux
118    fco2_land(:)=0.0
119    fco2_ocean(:)=0.0
120    fco2_land_nbp(:)=0.
121    fco2_land_nep(:)=0.
122    fco2_land_fLuc(:)=0.
123    fco2_land_fwoodharvest(:)=0.
124    fco2_land_fHarvest(:)=0.
125
126    DO nb=1, nbcf_in
127
128      SELECT CASE(cfname_in(nb))
129!--dealing with the different fluxes coming from ORCHIDEE
130!--fluxes come in unit of kg C m-2 s-1 is converted into kg CO2 m-2 s-1
131      CASE("fCO2_nep")
132          fco2_land_nep(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
133      CASE("fCO2_fLuc")
134          fco2_land_fLuc(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
135      CASE("fCO2_fwoodharvest")
136          fco2_land_fwoodharvest(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
137      CASE("fCO2_fHarvest")
138          fco2_land_fHarvest(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
139      CASE("fCO2_nbp")
140          fco2_land_nbp(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
141!--fCO2_fco2_ocn comes in unit of mol C02 m-2 s-1 is converted into kg CO2 m-2 s-1 + change sign
142      CASE("fCO2_fgco2")
143          fco2_ocean(:)=-1.*fco2_ocn_day(:)*RMCO2/1.e3*(pctsrf(:,is_oce)+pctsrf(:,is_sic))
144      END SELECT
145
146    ENDDO
147
148!--if fCO2_nbp is transferred we use it, otherwise we use the sum of what has been passed from ORCHIDEE
149    IF (check_fCO2_nbp_in_cfname)  THEN
150       fco2_land(:)=fco2_land_nbp(:)
151    ELSE
152       fco2_land(:)=fco2_land_nep(:)+fco2_land_fLuc(:)+fco2_land_fwoodharvest(:)+fco2_land_fHarvest(:)
153    ENDIF
154
155!!--preparing the net anthropogenic flux at the surface for mixing layer
156!!--unit kg CO2 / m2 / s
157!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_ff) ',MAXVAL(fco2_ff)
158!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_ff) ',MINVAL(fco2_ff)
159!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_bb) ',MAXVAL(fco2_bb)
160!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_bb) ',MINVAL(fco2_bb)
161!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_land) ',MAXVAL(fco2_land)
162!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_land) ',MINVAL(fco2_land)
163!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(fco2_ocean) ',MAXVAL(fco2_ocean)
164!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(fco2_ocean) ',MINVAL(fco2_ocean)
165!    PRINT *, 'tracco2i_mod.F90 --- MAXVAL(source(:,id_CO2)) ',MAXVAL(source(:,id_CO2))
166!    PRINT *, 'tracco2i_mod.F90 --- MINVAL(source(:,id_CO2)) ',MINVAL(source(:,id_CO2))
167!
168!--build final source term for CO2
169    source(:,id_CO2)=fco2_ff(:)+fco2_bb(:)+fco2_land(:)+fco2_ocean(:)
170
171!--computing global mean CO2 for radiation
172!--for every timestep comment out the IF ENDIF statements
173!--otherwise this is updated every day
174    IF (debutphy.OR.day_cur.NE.day_pre) THEN
175
176      CALL gather(tr_seri(:,:,id_CO2),co2_glo)
177      CALL gather(m_air,m_air_glo)
178
179!$OMP MASTER
180
181!--compute a global mean CO2 value and print its value in ppm
182       IF (is_mpi_root) THEN
183         RCO2_tot=SUM(co2_glo*m_air_glo)  !--unit kg CO2
184         RCO2_glo=RCO2_tot/SUM(m_air_glo) !--unit kg CO2 / kg air
185         PRINT *,'tracco2i: global CO2 in ppm =', RCO2_glo*1.e6*RMD/RMCO2
186         PRINT *,'tracco2i: total CO2 in kg =', RCO2_tot
187       ENDIF
188!$OMP END MASTER
189       CALL bcast(RCO2_glo)
190       day_pre=day_cur
191!--if not carbon_cycle_tr, then we reinitialize the CO2 each day to its global mean value
192       IF (.NOT.carbon_cycle_tr) THEN
193         tr_seri(:,:,id_CO2)=RCO2_glo
194       ENDIF
195    ENDIF
196
197  END SUBROUTINE tracco2i
198
199  SUBROUTINE co2_emissions(debutphy)
200
201    USE dimphy
202    USE infotrac_phy
203    USE geometry_mod, ONLY : cell_area
204    USE mod_grid_phy_lmdz
205    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
206    USE mod_phys_lmdz_para, ONLY: gather, scatter
207    USE phys_cal_mod
208
209    USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_varid, nf95_open
210    USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
211
212    USE carbon_cycle_mod, ONLY : fco2_ff, fco2_bb, fco2_land, fco2_ocean
213
214    IMPLICIT NONE
215
216    INCLUDE "YOMCST.h"
217    LOGICAL,INTENT(IN) :: debutphy
218
219! For NetCDF:
220    INTEGER ncid_in  ! IDs for input files
221    INTEGER varid, ncerr
222
223    INTEGER :: n_glo, n_month
224    REAL, POINTER:: vector(:), time(:)
225    REAL,ALLOCATABLE       :: flx_co2ff_glo(:,:) !  fossil-fuel CO2
226    REAL,ALLOCATABLE       :: flx_co2bb_glo(:,:) !  biomass-burning CO2
227    REAL,ALLOCATABLE, SAVE :: flx_co2ff(:,:)     !  fossil-fuel CO2
228    REAL,ALLOCATABLE, SAVE :: flx_co2bb(:,:)     !  biomass-burning CO2
229!$OMP THREADPRIVATE(flx_co2ff,flx_co2bb)
230
231!! may be controlled via the .def later on
232!! also co2bb for now comes from ORCHIDEE
233    LOGICAL, PARAMETER :: readco2ff=.TRUE.
234!! this should be left to FALSE for now
235    LOGICAL, PARAMETER :: readco2bb=.FALSE.
236
237    CHARACTER (len = 20) :: modname = 'tracco2i.co2_emissions'
238    CHARACTER (len = 80) :: abort_message
239
240    IF (debutphy) THEN
241
242    ALLOCATE(flx_co2ff(klon,12))
243    ALLOCATE(flx_co2bb(klon,12))
244
245!$OMP MASTER
246    IF (is_mpi_root) THEN
247   
248      IF (.NOT.ALLOCATED(flx_co2ff_glo)) ALLOCATE(flx_co2ff_glo(klon_glo,12))
249      IF (.NOT.ALLOCATED(flx_co2bb_glo)) ALLOCATE(flx_co2bb_glo(klon_glo,12))
250
251!--reading CO2 fossil fuel emissions
252      IF (readco2ff) THEN
253
254        ! ... Open the COZff file
255        CALL nf95_open("sflx_lmdz_co2_ff.nc", nf90_nowrite, ncid_in)
256
257        CALL nf95_inq_varid(ncid_in, "vector", varid)
258        CALL nf95_gw_var(ncid_in, varid, vector)
259        n_glo = size(vector)
260        IF (n_glo.NE.klon_glo) THEN
261           abort_message='sflx_lmdz_co2_ff: le nombre de points n est pas egal a klon_glo'
262           CALL abort_physic(modname,abort_message,1)
263        ENDIF
264
265        CALL nf95_inq_varid(ncid_in, "time", varid)
266        CALL nf95_gw_var(ncid_in, varid, time)
267        n_month = size(time)
268        IF (n_month.NE.12) THEN
269           abort_message='sflx_lmdz_co2_ff: le nombre de month n est pas egal a 12'
270           CALL abort_physic(modname,abort_message,1)
271        ENDIF
272
273!--reading flx_co2 for fossil fuel
274        CALL nf95_inq_varid(ncid_in, "flx_co2", varid)
275        ncerr = nf90_get_var(ncid_in, varid, flx_co2ff_glo)
276
277        CALL nf95_close(ncid_in)
278   
279      ELSE  !--co2ff not to be read
280        flx_co2ff_glo(:,:)=0.0
281      ENDIF
282
283!--reading CO2 biomass burning emissions
284!--using it will be inconsistent with treatment in ORCHIDEE
285      IF (readco2bb) THEN
286
287      ! ... Open the CO2bb file
288      CALL nf95_open("sflx_lmdz_co2_bb.nc", nf90_nowrite, ncid_in)
289
290      CALL nf95_inq_varid(ncid_in, "vector", varid)
291      CALL nf95_gw_var(ncid_in, varid, vector)
292      n_glo = size(vector)
293      IF (n_glo.NE.klon_glo) THEN
294         abort_message='sflx_lmdz_co2_bb: le nombre de points n est pas egal a klon_glo'
295         CALL abort_physic(modname,abort_message,1)
296      ENDIF
297
298      CALL nf95_inq_varid(ncid_in, "time", varid)
299      CALL nf95_gw_var(ncid_in, varid, time)
300      n_month = size(time)
301      IF (n_month.NE.12) THEN
302         abort_message='sflx_lmdz_co2_bb: le nombre de month n est pas egal a 12'
303         CALL abort_physic(modname,abort_message,1)
304      ENDIF
305
306!--reading flx_co2 for biomass burning
307      CALL nf95_inq_varid(ncid_in, "flx_co2", varid)
308      ncerr = nf90_get_var(ncid_in, varid, flx_co2bb_glo)
309
310      CALL nf95_close(ncid_in)
311   
312      ELSE  !--co2bb not to be read
313        flx_co2bb_glo(:,:)=0.0
314      ENDIF
315
316    ENDIF
317!$OMP END MASTER
318
319    ! Allocation needed for all proc otherwise scatter might complain
320    IF (.NOT.ALLOCATED(flx_co2ff_glo)) ALLOCATE(flx_co2ff_glo(0,0))
321    IF (.NOT.ALLOCATED(flx_co2bb_glo)) ALLOCATE(flx_co2bb_glo(0,0))
322
323!--scatter on all proc
324    CALL scatter(flx_co2ff_glo,flx_co2ff)
325    CALL scatter(flx_co2bb_glo,flx_co2bb)
326
327   IF (ALLOCATED(flx_co2ff_glo)) DEALLOCATE(flx_co2ff_glo)
328   IF (ALLOCATED(flx_co2bb_glo)) DEALLOCATE(flx_co2bb_glo)
329
330  ENDIF !--end debuthy
331
332!---select the correct month
333  IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
334    PRINT *,'probleme avec le mois dans co2_ini =', mth_cur
335  ENDIF
336
337  fco2_ff(:) = flx_co2ff(:,mth_cur)
338  fco2_bb(:) = flx_co2bb(:,mth_cur)
339
340  END SUBROUTINE co2_emissions
341
342END MODULE tracco2i_mod
Note: See TracBrowser for help on using the repository browser.