source: LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/read_newemissions.f90 @ 5119

Last change on this file since 5119 was 5119, checked in by abarral, 4 months ago

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way

File size: 13.5 KB
Line 
1! Routine to read the emissions of the different species
2!
3SUBROUTINE 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
368END SUBROUTINE read_newemissions
Note: See TracBrowser for help on using the repository browser.