source: LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/lmdz_read_newemissions.f90 @ 5182

Last change on this file since 5182 was 5160, checked in by abarral, 3 months ago

Put .h into modules

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