1 | ! Routine to read the emissions of the different species |
---|
2 | ! |
---|
3 | SUBROUTINE read_newemissions(julien, jH_emi, edgar, flag_dms, & |
---|
4 | debutphy, & |
---|
5 | pdtphys, lafinphy, nbjour, pctsrf, & |
---|
6 | t_seri, xlat, xlon, & |
---|
7 | pmflxr, pmflxs, prfl, psfl, & |
---|
8 | u10m_ec, v10m_ec, dust_ec, & |
---|
9 | lmt_sea_salt, lmt_so2ff_l, & |
---|
10 | lmt_so2ff_h, lmt_so2nff, lmt_so2ba, & |
---|
11 | lmt_so2bb_l, lmt_so2bb_h, & |
---|
12 | lmt_so2volc_cont, lmt_altvolc_cont, & |
---|
13 | lmt_so2volc_expl, lmt_altvolc_expl, & |
---|
14 | lmt_dmsbio, lmt_h2sbio, lmt_dmsconc, & |
---|
15 | lmt_bcff, lmt_bcnff, lmt_bcbb_l, & |
---|
16 | lmt_bcbb_h, lmt_bcba, lmt_omff, & |
---|
17 | lmt_omnff, lmt_ombb_l, lmt_ombb_h, & |
---|
18 | lmt_omnat, lmt_omba) |
---|
19 | |
---|
20 | USE dimphy |
---|
21 | USE indice_sol_mod |
---|
22 | USE lmdz_grid_phy |
---|
23 | USE lmdz_phys_para |
---|
24 | USE lmdz_ssum_scopy, ONLY: scopy |
---|
25 | |
---|
26 | IMPLICIT NONE |
---|
27 | |
---|
28 | INCLUDE "dimensions.h" |
---|
29 | INCLUDE 'paramet.h' |
---|
30 | INCLUDE 'chem.h' |
---|
31 | INCLUDE 'chem_spla.h' |
---|
32 | |
---|
33 | LOGICAL :: debutphy, lafinphy, edgar |
---|
34 | INTEGER :: test_vent, test_day, step_vent, flag_dms, nbjour |
---|
35 | INTEGER :: julien, i, iday |
---|
36 | SAVE step_vent, test_vent, test_day, iday |
---|
37 | !$OMP THREADPRIVATE(step_vent, test_vent, test_day, iday) |
---|
38 | REAL :: pct_ocean(klon), pctsrf(klon, nbsrf) |
---|
39 | REAL :: pdtphys ! pas d'integration pour la physique (seconde) |
---|
40 | REAL :: t_seri(klon, klev) ! temperature |
---|
41 | |
---|
42 | REAL :: xlat(klon) ! latitudes pour chaque point |
---|
43 | REAL :: xlon(klon) ! longitudes pour chaque point |
---|
44 | |
---|
45 | ! |
---|
46 | ! Emissions: |
---|
47 | ! --------- |
---|
48 | ! |
---|
49 | !---------------------------- SEA SALT & DUST emissions ------------------------ |
---|
50 | REAL :: lmt_sea_salt(klon, ss_bins) !Sea salt 0.03-8.0 um !NOT SAVED OK |
---|
51 | REAL :: clyfac, avgdryrate, drying |
---|
52 | ! je REAL u10m_ec1(klon), v10m_ec1(klon), dust_ec1(klon) |
---|
53 | ! je REAL u10m_ec2(klon), v10m_ec2(klon), dust_ec2(klon) |
---|
54 | |
---|
55 | REAL, SAVE, ALLOCATABLE :: u10m_ec1(:), v10m_ec1(:), dust_ec1(:) |
---|
56 | REAL, SAVE, ALLOCATABLE :: u10m_ec2(:), v10m_ec2(:), dust_ec2(:) |
---|
57 | !$OMP THREADPRIVATE(u10m_ec1, v10m_ec1, dust_ec1) |
---|
58 | !$OMP THREADPRIVATE(u10m_ec2, v10m_ec2, dust_ec2) |
---|
59 | ! as REAL u10m_nc(iip1,jjp1), v10m_nc(iip1,jjp1) |
---|
60 | REAL :: u10m_ec(klon), v10m_ec(klon), dust_ec(klon) |
---|
61 | ! REAL cly(klon), wth(klon), zprecipinsoil(klon) |
---|
62 | REAL, SAVE, ALLOCATABLE :: cly(:), wth(:), zprecipinsoil(:) |
---|
63 | REAL :: cly_glo(klon_glo), wth_glo(klon_glo) |
---|
64 | REAL :: zprecipinsoil_glo(klon_glo) |
---|
65 | !$OMP THREADPRIVATE(cly,wth,zprecipinsoil) |
---|
66 | |
---|
67 | |
---|
68 | ! je SAVE u10m_ec2, v10m_ec2, dust_ec2 |
---|
69 | ! je SAVE u10m_ec1, v10m_ec1, dust_ec1 ! Added on titane |
---|
70 | ! je SAVE cly, wth, zprecipinsoil ! Added on titane |
---|
71 | ! SAVE cly, wth, zprecipinsoil, u10m_ec2, v10m_ec2, dust_ec2 |
---|
72 | !------------------------- BLACK CARBON emissions ---------------------- |
---|
73 | REAL :: lmt_bcff(klon) ! emissions de BC fossil fuels |
---|
74 | REAL :: lmt_bcnff(klon) ! emissions de BC non-fossil fuels |
---|
75 | REAL :: lmt_bcbb_l(klon) ! emissions de BC biomass basses |
---|
76 | REAL :: lmt_bcbb_h(klon) ! emissions de BC biomass hautes |
---|
77 | REAL :: lmt_bcba(klon) ! emissions de BC bateau |
---|
78 | !------------------------ ORGANIC MATTER emissions --------------------- |
---|
79 | REAL :: lmt_omff(klon) ! emissions de OM fossil fuels |
---|
80 | REAL :: lmt_omnff(klon) ! emissions de OM non-fossil fuels |
---|
81 | REAL :: lmt_ombb_l(klon) ! emissions de OM biomass basses |
---|
82 | REAL :: lmt_ombb_h(klon) ! emissions de OM biomass hautes |
---|
83 | REAL :: lmt_omnat(klon) ! emissions de OM Natural |
---|
84 | REAL :: lmt_omba(klon) ! emissions de OM bateau |
---|
85 | !------------------------- SULFUR emissions ---------------------------- |
---|
86 | REAL :: lmt_so2ff_l(klon) ! emissions so2 fossil fuels (low) |
---|
87 | REAL :: lmt_so2ff_h(klon) ! emissions so2 fossil fuels (high) |
---|
88 | REAL :: lmt_so2nff(klon) ! emissions so2 non-fossil fuels |
---|
89 | REAL :: lmt_so2bb_l(klon) ! emissions de so2 biomass burning basse |
---|
90 | REAL :: lmt_so2bb_h(klon) ! emissions de so2 biomass burning hautes |
---|
91 | REAL :: lmt_so2ba(klon) ! emissions de so2 bateau |
---|
92 | REAL :: lmt_so2volc_cont(klon) ! emissions so2 volcan continuous |
---|
93 | REAL :: lmt_altvolc_cont(klon) ! altitude so2 volcan continuous |
---|
94 | REAL :: lmt_so2volc_expl(klon) ! emissions so2 volcan explosive |
---|
95 | REAL :: lmt_altvolc_expl(klon) ! altitude so2 volcan explosive |
---|
96 | REAL :: lmt_dmsconc(klon) ! concentration de dms oceanique |
---|
97 | REAL :: lmt_dmsbio(klon) ! emissions de dms bio |
---|
98 | REAL :: lmt_h2sbio(klon) ! emissions de h2s bio |
---|
99 | |
---|
100 | REAL, SAVE, ALLOCATABLE :: lmt_dms(:) ! emissions de dms |
---|
101 | !$OMP THREADPRIVATE(lmt_dms) |
---|
102 | ! |
---|
103 | ! Lessivage |
---|
104 | ! --------- |
---|
105 | ! |
---|
106 | REAL :: pmflxr(klon, klev + 1), pmflxs(klon, klev + 1) !--convection |
---|
107 | REAL :: prfl(klon, klev + 1), psfl(klon, klev + 1) !--large-scale |
---|
108 | ! REAL pmflxr(klon,klev), pmflxs(klon,klev) !--convection |
---|
109 | ! REAL prfl(klon,klev), psfl(klon,klev) !--large-scale |
---|
110 | ! |
---|
111 | ! Variable interne |
---|
112 | ! ---------------- |
---|
113 | ! |
---|
114 | INTEGER :: icount |
---|
115 | REAL :: tau_1, tau_2 |
---|
116 | REAL :: max_flux, min_flux |
---|
117 | INTRINSIC MIN, MAX |
---|
118 | ! |
---|
119 | ! JE: Changes due to new pdtphys in new physics. |
---|
120 | ! REAL windintime ! time in hours of the wind input files resolution |
---|
121 | ! REAL dayemintime ! time in hours of the other emissions input files resolution |
---|
122 | REAL :: jH_init ! shift in the hour (count as days) respecto to |
---|
123 | ! realhour = (pdtphys*i)/3600/24 -days_elapsed |
---|
124 | REAL :: jH_emi, jH_vent, jH_day |
---|
125 | SAVE jH_init, jH_vent, jH_day |
---|
126 | !$OMP THREADPRIVATE(jH_init,jH_vent,jH_day) |
---|
127 | REAL, PARAMETER :: vent_resol = 6. ! resolution of winds in hours |
---|
128 | REAL, PARAMETER :: day_resol = 24. ! resolution of daily emmis. in hours |
---|
129 | ! INTEGER test_day1 |
---|
130 | ! SAVE test_day1 |
---|
131 | ! REAL tau_1j,tau_2j |
---|
132 | ! je |
---|
133 | ! allocate if necessary |
---|
134 | ! |
---|
135 | |
---|
136 | IF (.NOT. ALLOCATED(u10m_ec1)) ALLOCATE(u10m_ec1(klon)) |
---|
137 | IF (.NOT. ALLOCATED(v10m_ec1)) ALLOCATE(v10m_ec1(klon)) |
---|
138 | IF (.NOT. ALLOCATED(dust_ec1)) ALLOCATE(dust_ec1(klon)) |
---|
139 | IF (.NOT. ALLOCATED(u10m_ec2)) ALLOCATE(u10m_ec2(klon)) |
---|
140 | IF (.NOT. ALLOCATED(v10m_ec2)) ALLOCATE(v10m_ec2(klon)) |
---|
141 | IF (.NOT. ALLOCATED(dust_ec2)) ALLOCATE(dust_ec2(klon)) |
---|
142 | IF (.NOT. ALLOCATED(cly)) ALLOCATE(cly(klon)) |
---|
143 | IF (.NOT. ALLOCATED(wth)) ALLOCATE(wth(klon)) |
---|
144 | IF (.NOT. ALLOCATED(zprecipinsoil)) ALLOCATE(zprecipinsoil(klon)) |
---|
145 | IF (.NOT. ALLOCATED(lmt_dms)) ALLOCATE(lmt_dms(klon)) |
---|
146 | ! end je nov2013 |
---|
147 | ! |
---|
148 | !*********************************************************************** |
---|
149 | ! DUST EMISSIONS |
---|
150 | !*********************************************************************** |
---|
151 | ! |
---|
152 | IF (debutphy) THEN |
---|
153 | !---Fields are read only at the beginning of the period |
---|
154 | !--reading wind and dust |
---|
155 | iday = julien |
---|
156 | step_vent = 1 |
---|
157 | test_vent = 0 |
---|
158 | test_day = 0 |
---|
159 | CALL read_vent(.TRUE., step_vent, nbjour, u10m_ec2, v10m_ec2) |
---|
160 | print *, 'Read (debut) dust emissions: step_vent,julien,nbjour', & |
---|
161 | step_vent, julien, nbjour |
---|
162 | CALL read_dust(.TRUE., step_vent, nbjour, dust_ec2) |
---|
163 | ! Threshold velocity map |
---|
164 | !$OMP MASTER |
---|
165 | IF (is_mpi_root .AND. is_omp_root) THEN |
---|
166 | zprecipinsoil_glo(:) = 0.0 |
---|
167 | OPEN(51, file = 'wth.dat', status = 'unknown', form = 'formatted') |
---|
168 | READ(51, '(G18.10)') (wth_glo(i), i = 1, klon_glo) |
---|
169 | CLOSE(51) |
---|
170 | ! Clay content |
---|
171 | OPEN(52, file = 'cly.dat', status = 'unknown', form = 'formatted') |
---|
172 | READ(52, '(G18.10)') (cly_glo(i), i = 1, klon_glo) |
---|
173 | CLOSE(52) |
---|
174 | OPEN(53, file = 'precipinsoil.dat', & |
---|
175 | status = 'old', form = 'formatted', err = 999) |
---|
176 | READ(53, '(G18.10)') (zprecipinsoil_glo(i), i = 1, klon_glo) |
---|
177 | PRINT *, 'lecture precipinsoil.dat' |
---|
178 | 999 CONTINUE |
---|
179 | CLOSE(53) |
---|
180 | ENDIF |
---|
181 | !$OMP END MASTER |
---|
182 | !$OMP BARRIER |
---|
183 | CALL scatter(wth_glo, wth) |
---|
184 | CALL scatter(cly_glo, cly) |
---|
185 | CALL scatter(zprecipinsoil_glo, zprecipinsoil) |
---|
186 | |
---|
187 | !JE20140908<< GOTO 1000 |
---|
188 | ! DO i=1, klon |
---|
189 | ! zprecipinsoil(i)=0.0 |
---|
190 | ! ENDDO |
---|
191 | ! 1000 CLOSE(53) |
---|
192 | !JE20140908>> |
---|
193 | jH_init = jH_emi |
---|
194 | jH_vent = jH_emi |
---|
195 | jH_day = jH_emi |
---|
196 | ! test_day1=0 |
---|
197 | !JE end |
---|
198 | ! |
---|
199 | |
---|
200 | ENDIF !--- debutphy |
---|
201 | |
---|
202 | print *, 'READ_EMISSION: test_vent & test_day = ', test_vent, & |
---|
203 | test_day |
---|
204 | IF (test_vent==0) THEN !--on lit toutes les 6 h |
---|
205 | CALL SCOPY(klon, u10m_ec2, 1, u10m_ec1, 1) |
---|
206 | CALL SCOPY(klon, v10m_ec2, 1, v10m_ec1, 1) |
---|
207 | CALL SCOPY(klon, dust_ec2, 1, dust_ec1, 1) |
---|
208 | step_vent = step_vent + 1 |
---|
209 | ! !PRINT *,'step_vent=', step_vent |
---|
210 | CALL read_vent(.FALSE., step_vent, nbjour, u10m_ec2, v10m_ec2) |
---|
211 | print *, 'Reading dust emissions: step_vent, julien, nbjour ', & |
---|
212 | step_vent, julien, nbjour |
---|
213 | ! !print *,'test_vent, julien = ',test_vent, julien |
---|
214 | CALL read_dust(.FALSE., step_vent, nbjour, dust_ec2) |
---|
215 | |
---|
216 | ENDIF !--test_vent |
---|
217 | |
---|
218 | ! ubicacion original |
---|
219 | ! test_vent=test_vent+1 |
---|
220 | ! IF (test_vent.EQ.(6*2)) test_vent=0 !on remet a zero ttes les 6 h |
---|
221 | |
---|
222 | !JE tau_2=FLOAT(test_vent)/12. |
---|
223 | !JE tau_1=1.-tau_2 |
---|
224 | tau_2 = (jH_vent - jH_init) * 24. / (vent_resol) |
---|
225 | tau_1 = 1. - tau_2 |
---|
226 | ! PRINT*,'JEdec jHv,JHi,ventres',jH_vent,jH_init,vent_resol |
---|
227 | ! PRINT*,'JEdec tau2,tau1',tau_2,tau_1 |
---|
228 | ! PRINT*,'JEdec step_vent',step_vent |
---|
229 | DO i = 1, klon |
---|
230 | ! PRINT*,'JE tau_2,tau_2j',tau_2,tau_2j |
---|
231 | u10m_ec(i) = tau_1 * u10m_ec1(i) + tau_2 * u10m_ec2(i) |
---|
232 | v10m_ec(i) = tau_1 * v10m_ec1(i) + tau_2 * v10m_ec2(i) |
---|
233 | dust_ec(i) = tau_1 * dust_ec1(i) + tau_2 * dust_ec2(i) |
---|
234 | ENDDO |
---|
235 | ! |
---|
236 | !JE IF (test_vent.EQ.(6*2)) THEN |
---|
237 | !JE PRINT *,'6 hrs interval reached' |
---|
238 | !JE print *,'day in read_emission, test_vent = ',julien, test_vent |
---|
239 | !JE ENDIF |
---|
240 | !JE |
---|
241 | !JE test_vent=test_vent+1 |
---|
242 | !JE IF (test_vent.EQ.(6*2)) test_vent=0 !on remet a zero ttes les 6 h |
---|
243 | ! JE |
---|
244 | jH_vent = jH_vent + pdtphys / (24. * 3600.) |
---|
245 | test_vent = test_vent + 1 |
---|
246 | IF (jH_vent>(vent_resol) / 24.) THEN |
---|
247 | test_vent = 0 |
---|
248 | jH_vent = jH_init |
---|
249 | ENDIF |
---|
250 | ! PRINT*,'JE test_vent,test_vent1,jH_vent ', test_vent,test_vent1 |
---|
251 | ! . ,jH_vent |
---|
252 | ! endJEi |
---|
253 | ! |
---|
254 | avgdryrate = 300. / 365. * pdtphys / 86400. |
---|
255 | ! |
---|
256 | DO i = 1, klon |
---|
257 | ! |
---|
258 | IF (cly(i)<9990..AND.wth(i)<9990.) THEN |
---|
259 | zprecipinsoil(i) = zprecipinsoil(i) + & |
---|
260 | (pmflxr(i, 1) + pmflxs(i, 1) + prfl(i, 1) + psfl(i, 1)) * pdtphys |
---|
261 | ! |
---|
262 | clyfac = MIN(16., cly(i) * 0.4 + 8.) ![mm] max amount of water hold in top soil |
---|
263 | drying = avgdryrate * exp(0.03905491 * & |
---|
264 | exp(0.17446 * (t_seri(i, 1) - 273.15))) ! [mm] |
---|
265 | zprecipinsoil(i) = min(max(0., zprecipinsoil(i) - drying), clyfac) ! [mm] |
---|
266 | ENDIF |
---|
267 | ! zprecipinsoil(i)=0.0 ! Temporarely introduced to reproduce obelix result |
---|
268 | ENDDO |
---|
269 | |
---|
270 | ! print *,'cly = ',sum(cly),maxval(cly),minval(cly) |
---|
271 | ! print *,'wth = ',sum(wth),maxval(wth),minval(wth) |
---|
272 | ! print *,'t_seri = ',sum(t_seri),maxval(t_seri),minval(t_seri) |
---|
273 | ! print *,'precipinsoil = ',sum(zprecipinsoil),maxval(zprecipinsoil) |
---|
274 | ! . ,minval(zprecipinsoil) |
---|
275 | icount = 0 |
---|
276 | DO i = 1, klon |
---|
277 | IF (cly(i)>=9990..OR.wth(i)>=9990..OR. & |
---|
278 | t_seri(i, 1)<=273.15.OR.zprecipinsoil(i)>1.e-8) THEN |
---|
279 | dust_ec(i) = 0.0 ! commented out for test dustemtest |
---|
280 | ! print *,'Dust emissions surpressed at grid = ',i |
---|
281 | ! icount=icount+1 |
---|
282 | ENDIF |
---|
283 | ENDDO |
---|
284 | ! |
---|
285 | print *, 'Total N of grids with surpressed emission = ', icount |
---|
286 | print *, 'dust_ec = ', SUM(dust_ec), MINVAL(dust_ec), & |
---|
287 | MAXVAL(dust_ec) |
---|
288 | !nhl Transitory scaling of desert dust emissions |
---|
289 | |
---|
290 | !nhl DO i=1, klon |
---|
291 | !nhl dust_ec(i)=dust_ec(i)/2. |
---|
292 | !nhl ENDDO |
---|
293 | |
---|
294 | !-saving precipitation field to be read in next simulation |
---|
295 | |
---|
296 | IF (lafinphy) THEN |
---|
297 | ! |
---|
298 | CALL gather(zprecipinsoil, zprecipinsoil_glo) |
---|
299 | !$OMP MASTER |
---|
300 | IF (is_mpi_root .AND. is_omp_root) THEN |
---|
301 | |
---|
302 | OPEN(53, file = 'newprecipinsoil.dat', & |
---|
303 | status = 'unknown', form = 'formatted') |
---|
304 | WRITE(53, '(G18.10)') (zprecipinsoil_glo(i), i = 1, klon_glo) |
---|
305 | CLOSE(53) |
---|
306 | ENDIF |
---|
307 | !$OMP END MASTER |
---|
308 | !$OMP BARRIER |
---|
309 | ! |
---|
310 | ENDIF |
---|
311 | ! |
---|
312 | !*********************************************************************** |
---|
313 | ! SEA SALT EMISSIONS |
---|
314 | !*********************************************************************** |
---|
315 | ! |
---|
316 | DO i = 1, klon |
---|
317 | pct_ocean(i) = pctsrf(i, is_oce) |
---|
318 | ENDDO |
---|
319 | |
---|
320 | print *, 'IS_OCE = ', is_oce |
---|
321 | CALL seasalt(v10m_ec, u10m_ec, pct_ocean, lmt_sea_salt) !mgSeaSalt/cm2/s |
---|
322 | ! print *,'SUM, MAX & MIN Sea Salt = ',SUM(lmt_sea_salt), |
---|
323 | ! . MAXVAL(lmt_sea_salt),MINVAL(lmt_sea_salt) |
---|
324 | ! |
---|
325 | !*********************************************************************** |
---|
326 | ! SULFUR & CARBON EMISSIONS |
---|
327 | !*********************************************************************** |
---|
328 | ! |
---|
329 | |
---|
330 | IF (test_day==0) THEN |
---|
331 | print *, 'Computing SULFATE emissions for day : ', iday, julien, & |
---|
332 | step_vent |
---|
333 | CALL condsurfs_new(iday, edgar, flag_dms, & |
---|
334 | lmt_so2ff_l, lmt_so2ff_h, lmt_so2nff, & |
---|
335 | lmt_so2bb_l, lmt_so2bb_h, lmt_so2ba, & |
---|
336 | lmt_so2volc_cont, lmt_altvolc_cont, & |
---|
337 | lmt_so2volc_expl, lmt_altvolc_expl, & |
---|
338 | lmt_dmsbio, lmt_h2sbio, lmt_dms, lmt_dmsconc) |
---|
339 | print *, 'Computing CARBON emissions for day : ', iday, julien, & |
---|
340 | step_vent |
---|
341 | CALL condsurfc_new(iday, & |
---|
342 | lmt_bcff, lmt_bcnff, lmt_bcbb_l, lmt_bcbb_h, & |
---|
343 | lmt_bcba, lmt_omff, lmt_omnff, lmt_ombb_l, & |
---|
344 | lmt_ombb_h, lmt_omnat, lmt_omba) |
---|
345 | print *, 'IDAY = ', iday |
---|
346 | iday = iday + 1 |
---|
347 | print *, 'BCBB_L emissions :', SUM(lmt_bcbb_l), MAXVAL(lmt_bcbb_l) & |
---|
348 | , MINVAL(lmt_bcbb_l) |
---|
349 | print *, 'BCBB_H emissions :', SUM(lmt_bcbb_h), MAXVAL(lmt_bcbb_h) & |
---|
350 | , MINVAL(lmt_bcbb_h) |
---|
351 | ENDIF |
---|
352 | |
---|
353 | !JE test_day=test_day+1 |
---|
354 | !JE IF (test_day.EQ.(24*2.)) THEN |
---|
355 | !JE test_day=0 !on remet a zero ttes les 24 h |
---|
356 | !JE print *,'LAST TIME STEP OF DAY ',julien |
---|
357 | !JE ENDIF |
---|
358 | |
---|
359 | jH_day = jH_day + pdtphys / (24. * 3600.) |
---|
360 | test_day = test_day + 1 |
---|
361 | IF (jH_day>(day_resol) / 24.) THEN |
---|
362 | print *, 'LAST TIME STEP OF DAY ', julien |
---|
363 | test_day = 0 |
---|
364 | jH_day = jH_init |
---|
365 | ENDIF |
---|
366 | ! PRINT*,'test_day,test_day1',test_day,test_day1 |
---|
367 | |
---|
368 | END SUBROUTINE read_newemissions |
---|