source: LMDZ6/trunk/libf/phylmd/Dust/read_newemissions.f90 @ 5326

Last change on this file since 5326 was 5326, checked in by abarral, 24 hours ago

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