source: LMDZ6/trunk/libf/phylmd/Dust/read_newemissions.F @ 3786

Last change on this file since 3786 was 3786, checked in by asima, 4 years ago

Makes LMDZ-SPLA work again.
Minimal - but not minor - modifications required by going from LMDZ5 to LMDZ6.
For the time being, SPLA input files (including correction coefficient files for aerosols emissions) are only available for the 128x88 grid zoomed on N Africa, used by Jeronimo Escribano and Binta Diallo.
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)
61c as      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.