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, 11 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.