source: LMDZ5/trunk/libf/phylmd/Dust/read_newemissions.F @ 4947

Last change on this file since 4947 was 2630, checked in by fhourdin, 8 years ago

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