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

Last change on this file since 5267 was 5246, checked in by abarral, 4 days ago

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

File size: 13.2 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 mod_grid_phy_lmdz
23  USE mod_phys_lmdz_para
24
25  IMPLICIT NONE
26
27
28  INCLUDE "dimensions.h"
29   ! INCLUDE 'dimphy.h'
30  INCLUDE 'paramet.h'
31  INCLUDE 'chem.h'
32  INCLUDE 'chem_spla.h'
33   ! 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
47  !
48  !   Emissions:
49  !   ---------
50  !
51  !---------------------------- 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
54  ! je      REAL u10m_ec1(klon), v10m_ec1(klon), dust_ec1(klon)
55  ! 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  ! as      REAL u10m_nc(iip1,jjp1), v10m_nc(iip1,jjp1)
62  REAL :: u10m_ec(klon), v10m_ec(klon), dust_ec(klon)
63  !      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
70  ! je     SAVE u10m_ec2, v10m_ec2, dust_ec2
71  ! je      SAVE u10m_ec1, v10m_ec1, dust_ec1   ! Added on titane
72  ! je      SAVE cly, wth, zprecipinsoil        ! Added on titane
73  ! SAVE cly, wth, zprecipinsoil, u10m_ec2, v10m_ec2, dust_ec2
74  !------------------------- 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
80  !------------------------ 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
87  !------------------------- 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)
104  !
105  !  Lessivage
106  !  ---------
107  !
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
112  !
113  !  Variable interne
114  !  ----------------
115  !
116  INTEGER :: icount
117  REAL :: tau_1, tau_2
118  REAL :: max_flux, min_flux
119  INTRINSIC MIN, MAX
120  !
121  ! JE: Changes due to new pdtphys in new physics.
122  !  REAL windintime ! time in hours of the wind input files resolution
123  !  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
134  ! je
135  ! allocate if necessary
136  !
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))
148  ! end je nov2013
149  !
150  !***********************************************************************
151  ! DUST EMISSIONS
152  !***********************************************************************
153  !
154  IF (debutphy) THEN
155  !---Fields are read only at the beginning of the period
156  !--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)
165  ! 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)
172  ! 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
200  !
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
220  ! ubicacion original
221  !  test_vent=test_vent+1
222  !  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
237  !
238  !JE      IF (test_vent.EQ.(6*2)) THEN
239  !JE        PRINT *,'6 hrs interval reached'
240  !JE        print *,'day in read_emission, test_vent = ',julien, test_vent
241  !JE      ENDIF
242  !JE
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
245  ! 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
254  ! endJEi
255  !
256  avgdryrate=300./365.*pdtphys/86400.
257  !
258  DO i=1, klon
259  !
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
263  !
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
286  !
287  print *,'Total N of grids with surpressed emission = ',icount
288  print *,'dust_ec = ',SUM(dust_ec),MINVAL(dust_ec), &
289        MAXVAL(dust_ec)
290  !nhl Transitory scaling of desert dust emissions
291
292  !nhl      DO i=1, klon
293  !nhl         dust_ec(i)=dust_ec(i)/2.
294  !nhl      ENDDO
295
296  !-saving precipitation field to be read in next simulation
297
298  IF (lafinphy) THEN
299  !
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
311  !
312  ENDIF
313  !
314  !***********************************************************************
315  ! SEA SALT EMISSIONS
316  !***********************************************************************
317  !
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)
326  !
327  !***********************************************************************
328  ! SULFUR & CARBON EMISSIONS
329  !***********************************************************************
330  !
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          lmt_so2ff_l, lmt_so2ff_h, lmt_so2nff, &
337          lmt_so2bb_l, lmt_so2bb_h, lmt_so2ba, &
338          lmt_so2volc_cont, lmt_altvolc_cont, &
339          lmt_so2volc_expl, lmt_altvolc_expl, &
340          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          lmt_bcff,lmt_bcnff,lmt_bcbb_l,lmt_bcbb_h, &
345          lmt_bcba,lmt_omff,lmt_omnff,lmt_ombb_l, &
346          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
371END SUBROUTINE read_newemissions
Note: See TracBrowser for help on using the repository browser.