source: LMDZ6/trunk/libf/phylmd/tracco2i_mod.F90 @ 3847

Last change on this file since 3847 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
RevLine 
[3361]1MODULE tracco2i_mod
2!
3! This module does the work for the interactive CO2 tracers
[3581]4! Authors: Patricia Cadule and Olivier Boucher
[3361]5!
[3581]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!
[3361]14CONTAINS
15
[3649]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
[3361]27  SUBROUTINE tracco2i(pdtphys, debutphy, &
28       xlat, xlon, pphis, pphi, &
29       t_seri, pplay, paprs, tr_seri, source)
30
31    USE dimphy
[3435]32    USE infotrac_phy
[3450]33    USE geometry_mod, ONLY: cell_area
[3453]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
[3581]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
[3361]39    USE mod_grid_phy_lmdz
[3450]40    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
[3361]41    USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter
42    USE phys_cal_mod
[3450]43    USE phys_state_var_mod, ONLY: pctsrf
44    USE indice_sol_mod, ONLY: nbsrf, is_ter, is_lic, is_oce, is_sic
[3361]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
[3421]73    INTEGER                        :: it, k, i, nb
74    REAL, DIMENSION(klon,klev)     :: m_air     ! mass of air in every grid box [kg]
[3361]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
[3581]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)
[3361]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
[3549]91
[3581]92! Initialisation de tr_seri(id_CO2) si pas initialise
[3361]93      IF (MAXVAL(tr_seri(:,:,id_CO2)).LT.1.e-15) THEN
[3581]94        tr_seri(:,:,id_CO2)=co2_ppm0*1.e-6/RMD*RMCO2 !--initialised from co2_ppm0 in rdem
[3361]95      ENDIF
[3549]96
[3581]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
[3549]102
[3361]103    ENDIF
104
[3421]105!--calculate mass of air in every grid box in kg air
[3361]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
[3421]114!--unit kg CO2 m-2 s-1
115    CALL co2_emissions(debutphy)
[3361]116
[3421]117!--retrieving land and ocean CO2 flux
[3453]118    fco2_land(:)=0.0
119    fco2_ocean(:)=0.0
[3581]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
[3421]126    DO nb=1, nbcf_in
[3581]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
[3421]146    ENDDO
147
[3581]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
[3453]169    source(:,id_CO2)=fco2_ff(:)+fco2_bb(:)+fco2_land(:)+fco2_ocean(:)
[3361]170
171!--computing global mean CO2 for radiation
[3581]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
[3361]176      CALL gather(tr_seri(:,:,id_CO2),co2_glo)
177      CALL gather(m_air,m_air_glo)
[3581]178
[3361]179!$OMP MASTER
[3549]180
[3450]181!--compute a global mean CO2 value and print its value in ppm
[3361]182       IF (is_mpi_root) THEN
[3450]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
[3361]187       ENDIF
188!$OMP END MASTER
189       CALL bcast(RCO2_glo)
[3450]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
[3581]195    ENDIF
[3361]196
197  END SUBROUTINE tracco2i
198
[3421]199  SUBROUTINE co2_emissions(debutphy)
[3361]200
201    USE dimphy
[3435]202    USE infotrac_phy
[3361]203    USE geometry_mod, ONLY : cell_area
204    USE mod_grid_phy_lmdz
[3366]205    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
[3361]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
[3549]212    USE carbon_cycle_mod, ONLY : fco2_ff, fco2_bb, fco2_land, fco2_ocean
[3421]213
[3361]214    IMPLICIT NONE
215
[3368]216    INCLUDE "YOMCST.h"
[3361]217    LOGICAL,INTENT(IN) :: debutphy
218
219! For NetCDF:
[3366]220    INTEGER ncid_in  ! IDs for input files
221    INTEGER varid, ncerr
[3361]222
223    INTEGER :: n_glo, n_month
224    REAL, POINTER:: vector(:), time(:)
[3366]225    REAL,ALLOCATABLE       :: flx_co2ff_glo(:,:) !  fossil-fuel CO2
[3421]226    REAL,ALLOCATABLE       :: flx_co2bb_glo(:,:) !  biomass-burning CO2
[3366]227    REAL,ALLOCATABLE, SAVE :: flx_co2ff(:,:)     !  fossil-fuel CO2
228    REAL,ALLOCATABLE, SAVE :: flx_co2bb(:,:)     !  biomass-burning CO2
[3421]229!$OMP THREADPRIVATE(flx_co2ff,flx_co2bb)
[3361]230
[3383]231!! may be controlled via the .def later on
232!! also co2bb for now comes from ORCHIDEE
[3581]233    LOGICAL, PARAMETER :: readco2ff=.TRUE.
234!! this should be left to FALSE for now
235    LOGICAL, PARAMETER :: readco2bb=.FALSE.
[3361]236
[3531]237    CHARACTER (len = 20) :: modname = 'tracco2i.co2_emissions'
238    CHARACTER (len = 80) :: abort_message
239
[3366]240    IF (debutphy) THEN
[3361]241
[3366]242    ALLOCATE(flx_co2ff(klon,12))
243    ALLOCATE(flx_co2bb(klon,12))
244
[3361]245!$OMP MASTER
[3366]246    IF (is_mpi_root) THEN
247   
[3383]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))
[3361]250
[3366]251!--reading CO2 fossil fuel emissions
252      IF (readco2ff) THEN
[3361]253
[3366]254        ! ... Open the COZff file
255        CALL nf95_open("sflx_lmdz_co2_ff.nc", nf90_nowrite, ncid_in)
[3361]256
[3366]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
[3531]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)
[3366]263        ENDIF
[3361]264
[3366]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
[3531]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)
[3366]271        ENDIF
[3361]272
[3366]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
[3581]284!--using it will be inconsistent with treatment in ORCHIDEE
[3366]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
[3531]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)
[3366]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
[3531]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)
[3366]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
[3361]317!$OMP END MASTER
318
[3651]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
[3361]323!--scatter on all proc
[3366]324    CALL scatter(flx_co2ff_glo,flx_co2ff)
325    CALL scatter(flx_co2bb_glo,flx_co2bb)
[3361]326
[3651]327   IF (ALLOCATED(flx_co2ff_glo)) DEALLOCATE(flx_co2ff_glo)
328   IF (ALLOCATED(flx_co2bb_glo)) DEALLOCATE(flx_co2bb_glo)
[3361]329
330  ENDIF !--end debuthy
331
332!---select the correct month
333  IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
[3366]334    PRINT *,'probleme avec le mois dans co2_ini =', mth_cur
[3361]335  ENDIF
[3549]336
[3421]337  fco2_ff(:) = flx_co2ff(:,mth_cur)
338  fco2_bb(:) = flx_co2bb(:,mth_cur)
[3361]339
340  END SUBROUTINE co2_emissions
341
342END MODULE tracco2i_mod
Note: See TracBrowser for help on using the repository browser.