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

Last change on this file since 3549 was 3549, checked in by oboucher, 5 years ago

Allocate CO2 flux variables at first timestep
File size: 8.9 KB
Line 
1MODULE tracco2i_mod
2!
3! This module does the work for the interactive CO2 tracers
4!
5CONTAINS
6
7  SUBROUTINE tracco2i(pdtphys, debutphy, &
8       xlat, xlon, pphis, pphi, &
9       t_seri, pplay, paprs, tr_seri, source)
10
11    USE dimphy
12    USE infotrac_phy
13    USE geometry_mod, ONLY: cell_area
14    USE carbon_cycle_mod, ONLY: id_CO2, nbcf_in, fields_in, cfname_in
15    USE carbon_cycle_mod, ONLY: fco2_ocn_day, fco2_ff, fco2_bb, fco2_land, fco2_ocean
16    USE carbon_cycle_mod, ONLY: carbon_cycle_tr, carbon_cycle_rad, RCO2_glo, RCO2_tot
17    USE mod_grid_phy_lmdz
18    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
19    USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter
20    USE phys_cal_mod
21    USE phys_state_var_mod, ONLY: pctsrf
22    USE indice_sol_mod, ONLY: nbsrf, is_ter, is_lic, is_oce, is_sic
23
24    IMPLICIT NONE
25
26    INCLUDE "clesphys.h"
27    INCLUDE "YOMCST.h"
28
29! Input argument
30!---------------
31    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
32    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
33
34    REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point
35    REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point
36    REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! geopotentiel du sol
37    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel de chaque couche
38
39    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
40    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
41    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
42    REAL,DIMENSION(klon,nbtr),INTENT(INOUT):: source  ! flux de traceur [U/m2/s]
43
44! Output argument
45!----------------
46    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/kgA] 
47
48! Local variables
49!----------------
50
51    INTEGER                        :: it, k, i, nb
52    REAL, DIMENSION(klon,klev)     :: m_air     ! mass of air in every grid box [kg]
53    REAL, DIMENSION(klon_glo,klev) :: co2_glo   ! variable temporaire sur la grille global
54    REAL, DIMENSION(klon_glo,klev) :: m_air_glo ! variable temporaire sur la grille global
55
56    INTEGER, SAVE :: mth_pre=0, day_pre=0
57!$OMP THREADPRIVATE(mth_pre, day_pre)
58
59    IF (is_mpi_root) THEN
60      PRINT *,'in tracco2i: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
61    ENDIF
62
63!--initialisation of CO2 field if not in restart file
64!--dirty way of doing, do it better later
65!--convert 280 ppm into kg CO2 / kg air
66    IF (debutphy) THEN
67
68      IF (MAXVAL(tr_seri(:,:,id_CO2)).LT.1.e-15) THEN
69        !!tr_seri(:,:,id_CO2)=280.e-6/RMD*RMCO2
70        tr_seri(:,:,id_CO2)=400.e-6/RMD*RMCO2 !--initialised to 400 ppm for a test
71      ENDIF
72
73      ALLOCATE(fco2_ff(klon))
74      ALLOCATE(fco2_bb(klon))
75      ALLOCATE(fco2_land(klon))
76      ALLOCATE(fco2_ocean(klon))
77
78    ENDIF
79
80!--calculate mass of air in every grid box in kg air
81    DO i=1, klon
82    DO k=1, klev
83      m_air(i,k)=(paprs(i,k)-paprs(i,k+1))/RG*cell_area(i)
84    ENDDO
85    ENDDO
86
87!--call CO2 emission routine
88!--co2bb is zero for now
89!--unit kg CO2 m-2 s-1
90    CALL co2_emissions(debutphy)
91
92!--retrieving land and ocean CO2 flux
93    fco2_land(:)=0.0
94    fco2_ocean(:)=0.0
95    DO nb=1, nbcf_in
96      print *,'nb tracco2=', nb, cfname_in(nb)
97!--fCO2_nep comes in unit of kg C m-2 s-1
98!--converting to kg CO2 m-2 s-1
99      IF (cfname_in(nb) == "fCO2_nbp" )   fco2_land(:)=fields_in(:,nb)*RMCO2/RMC*pctsrf(:,is_ter)
100!--fCO2_fgco2 comes in unit of mol C02 m-2 s-1
101!--converting to kg CO2 m-2 s-1 + change sign
102      IF (cfname_in(nb) == "fCO2_fgco2" ) fco2_ocean(:)=-1.*fco2_ocn_day(:)*RMCO2/1.e3*(pctsrf(:,is_oce)+pctsrf(:,is_sic))
103    ENDDO
104
105!--preparing the net anthropogenic flux at the surface for mixing layer
106!--unit kg CO2 / m2 / s
107    source(:,id_CO2)=fco2_ff(:)+fco2_bb(:)+fco2_land(:)+fco2_ocean(:)
108
109!--computing global mean CO2 for radiation
110!--every timestep for now but enough every day
111!    IF (debutphy.OR.mth_cur.NE.mth_pre) THEN
112!    IF (debutphy.OR.day_cur.NE.day_pre) THEN
113      CALL gather(tr_seri(:,:,id_CO2),co2_glo)
114      CALL gather(m_air,m_air_glo)
115!$OMP MASTER
116
117!--compute a global mean CO2 value and print its value in ppm
118       IF (is_mpi_root) THEN
119         RCO2_tot=SUM(co2_glo*m_air_glo)  !--unit kg CO2
120         RCO2_glo=RCO2_tot/SUM(m_air_glo) !--unit kg CO2 / kg air
121         PRINT *,'tracco2i: global CO2 in ppm =', RCO2_glo*1.e6*RMD/RMCO2
122         PRINT *,'tracco2i: total CO2 in kg =', RCO2_tot
123       ENDIF
124!$OMP END MASTER
125       CALL bcast(RCO2_glo)
126       mth_pre=mth_cur
127       day_pre=day_cur
128!--if not carbon_cycle_tr, then we reinitialize the CO2 each day to its global mean value
129       IF (.NOT.carbon_cycle_tr) THEN
130         tr_seri(:,:,id_CO2)=RCO2_glo
131       ENDIF
132!    ENDIF
133
134  END SUBROUTINE tracco2i
135
136  SUBROUTINE co2_emissions(debutphy)
137
138    USE dimphy
139    USE infotrac_phy
140    USE geometry_mod, ONLY : cell_area
141    USE mod_grid_phy_lmdz
142    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
143    USE mod_phys_lmdz_para, ONLY: gather, scatter
144    USE phys_cal_mod
145
146    USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_varid, nf95_open
147    USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
148
149    USE carbon_cycle_mod, ONLY : fco2_ff, fco2_bb, fco2_land, fco2_ocean
150
151    IMPLICIT NONE
152
153    INCLUDE "YOMCST.h"
154    LOGICAL,INTENT(IN) :: debutphy
155
156! For NetCDF:
157    INTEGER ncid_in  ! IDs for input files
158    INTEGER varid, ncerr
159
160    INTEGER :: n_glo, n_month
161    REAL, POINTER:: vector(:), time(:)
162    REAL,ALLOCATABLE       :: flx_co2ff_glo(:,:) !  fossil-fuel CO2
163    REAL,ALLOCATABLE       :: flx_co2bb_glo(:,:) !  biomass-burning CO2
164    REAL,ALLOCATABLE, SAVE :: flx_co2ff(:,:)     !  fossil-fuel CO2
165    REAL,ALLOCATABLE, SAVE :: flx_co2bb(:,:)     !  biomass-burning CO2
166!$OMP THREADPRIVATE(flx_co2ff,flx_co2bb)
167
168!! may be controlled via the .def later on
169!! also co2bb for now comes from ORCHIDEE
170    LOGICAL, PARAMETER :: readco2ff=.TRUE., readco2bb=.FALSE.
171
172    CHARACTER (len = 20) :: modname = 'tracco2i.co2_emissions'
173    CHARACTER (len = 80) :: abort_message
174
175    IF (debutphy) THEN
176
177    ALLOCATE(flx_co2ff(klon,12))
178    ALLOCATE(flx_co2bb(klon,12))
179
180!$OMP MASTER
181    IF (is_mpi_root) THEN
182   
183      IF (.NOT.ALLOCATED(flx_co2ff_glo)) ALLOCATE(flx_co2ff_glo(klon_glo,12))
184      IF (.NOT.ALLOCATED(flx_co2bb_glo)) ALLOCATE(flx_co2bb_glo(klon_glo,12))
185
186!--reading CO2 fossil fuel emissions
187      IF (readco2ff) THEN
188
189        ! ... Open the COZff file
190        CALL nf95_open("sflx_lmdz_co2_ff.nc", nf90_nowrite, ncid_in)
191
192        CALL nf95_inq_varid(ncid_in, "vector", varid)
193        CALL nf95_gw_var(ncid_in, varid, vector)
194        n_glo = size(vector)
195        IF (n_glo.NE.klon_glo) THEN
196           abort_message='sflx_lmdz_co2_ff: le nombre de points n est pas egal a klon_glo'
197           CALL abort_physic(modname,abort_message,1)
198        ENDIF
199
200        CALL nf95_inq_varid(ncid_in, "time", varid)
201        CALL nf95_gw_var(ncid_in, varid, time)
202        n_month = size(time)
203        IF (n_month.NE.12) THEN
204           abort_message='sflx_lmdz_co2_ff: le nombre de month n est pas egal a 12'
205           CALL abort_physic(modname,abort_message,1)
206        ENDIF
207
208!--reading flx_co2 for fossil fuel
209        CALL nf95_inq_varid(ncid_in, "flx_co2", varid)
210        ncerr = nf90_get_var(ncid_in, varid, flx_co2ff_glo)
211
212        CALL nf95_close(ncid_in)
213   
214      ELSE  !--co2ff not to be read
215        flx_co2ff_glo(:,:)=0.0
216      ENDIF
217
218!--reading CO2 biomass burning emissions
219      IF (readco2bb) THEN
220
221      ! ... Open the CO2bb file
222      CALL nf95_open("sflx_lmdz_co2_bb.nc", nf90_nowrite, ncid_in)
223
224      CALL nf95_inq_varid(ncid_in, "vector", varid)
225      CALL nf95_gw_var(ncid_in, varid, vector)
226      n_glo = size(vector)
227      IF (n_glo.NE.klon_glo) THEN
228         abort_message='sflx_lmdz_co2_bb: le nombre de points n est pas egal a klon_glo'
229         CALL abort_physic(modname,abort_message,1)
230      ENDIF
231
232      CALL nf95_inq_varid(ncid_in, "time", varid)
233      CALL nf95_gw_var(ncid_in, varid, time)
234      n_month = size(time)
235      IF (n_month.NE.12) THEN
236         abort_message='sflx_lmdz_co2_bb: le nombre de month n est pas egal a 12'
237         CALL abort_physic(modname,abort_message,1)
238      ENDIF
239
240!--reading flx_co2 for biomass burning
241      CALL nf95_inq_varid(ncid_in, "flx_co2", varid)
242      ncerr = nf90_get_var(ncid_in, varid, flx_co2bb_glo)
243
244      CALL nf95_close(ncid_in)
245   
246      ELSE  !--co2bb not to be read
247        flx_co2bb_glo(:,:)=0.0
248      ENDIF
249
250    ENDIF
251!$OMP END MASTER
252
253!--scatter on all proc
254    CALL scatter(flx_co2ff_glo,flx_co2ff)
255    CALL scatter(flx_co2bb_glo,flx_co2bb)
256
257!$OMP MASTER
258    IF (is_mpi_root) THEN
259       DEALLOCATE(flx_co2ff_glo)
260       DEALLOCATE(flx_co2bb_glo)
261    ENDIF
262!$OMP END MASTER
263
264  ENDIF !--end debuthy
265
266!---select the correct month
267  IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
268    PRINT *,'probleme avec le mois dans co2_ini =', mth_cur
269  ENDIF
270
271  fco2_ff(:) = flx_co2ff(:,mth_cur)
272  fco2_bb(:) = flx_co2bb(:,mth_cur)
273
274  END SUBROUTINE co2_emissions
275
276END MODULE tracco2i_mod
Note: See TracBrowser for help on using the repository browser.