source: LMDZ6/branches/LMDZ-COSP/libf/phylmd/StratAer/traccoag_mod.F90 @ 5305

Last change on this file since 5305 was 3677, checked in by oboucher, 4 years ago

Changed the way to initialise nbtr_bin and other dimensions and indices
in the StratAer? module based on infotrac_phy rather than infotrac.

Also added a missing $OMP THREADPRIVATE(nqperes)

  • Property svn:keywords set to Id
File size: 19.4 KB
Line 
1!
2! $Id: traccoag_mod.F90 3677 2020-05-06 15:18:32Z abarral $
3!
4MODULE traccoag_mod
5!
6! This module calculates the concentration of aerosol particles in certain size bins
7! considering coagulation and sedimentation.
8!
9CONTAINS
10
11  SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
12       presnivs, xlat, xlon, pphis, pphi, &
13       t_seri, pplay, paprs, sh, rh, tr_seri)
14
15    USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
16        & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part
17
18    USE dimphy
19    USE infotrac_phy
20    USE aerophys
21    USE geometry_mod, ONLY : cell_area, boundslat
22    USE mod_grid_phy_lmdz
23    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
24    USE mod_phys_lmdz_para, only: gather, scatter
25    USE phys_cal_mod
26    USE sulfate_aer_mod
27    USE phys_local_var_mod, ONLY: stratomask
28    USE YOMCST
29    USE print_control_mod, ONLY: lunout
30    USE strataer_mod
31    USE phys_cal_mod, ONLY : year_len
32
33    IMPLICIT NONE
34
35! Input argument
36!---------------
37    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
38    REAL,INTENT(IN)    :: gmtime     ! Heure courante
39    LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
40    INTEGER,INTENT(IN) :: julien     ! Jour julien
41
42    REAL,DIMENSION(klev),INTENT(IN)        :: presnivs! pressions approximat. des milieux couches (en PA)
43    REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point
44    REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point
45    REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! geopotentiel du sol
46    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel de chaque couche
47
48    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
49    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
50    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
51    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
52    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative   
53
54! Output argument
55!----------------
56    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
57
58! Local variables
59!----------------
60    REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
61    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
62    LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
63    REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
64    REAL,DIMENSION(klon_glo,klev,nbtr)     :: tr_seri_glo         ! Concentration Traceur [U/KgA] 
65    REAL,DIMENSION(klev+1)                 :: altLMDz             ! altitude of layer interfaces in m
66    REAL,DIMENSION(klev)                   :: f_lay_emiss         ! fraction of emission for every vertical layer
67    REAL                                   :: f_lay_sum           ! sum of layer emission fractions
68    REAL                                   :: alt                 ! altitude for integral calculation
69    INTEGER,PARAMETER                      :: n_int_alt=10        ! number of subintervals for integration over Gaussian emission profile
70    REAL,DIMENSION(nbtr_bin)               :: r_bin               ! particle radius in size bin [m]
71    REAL,DIMENSION(nbtr_bin)               :: r_lower             ! particle radius at lower bin boundary [m]
72    REAL,DIMENSION(nbtr_bin)               :: r_upper             ! particle radius at upper bin boundary [m]
73    REAL,DIMENSION(nbtr_bin)               :: m_part_dry          ! mass of one dry particle in size bin [kg]
74    REAL                                   :: zrho                ! Density of air [kg/m3]
75    REAL                                   :: zdz                 ! thickness of atm. model layer in m
76    REAL,DIMENSION(klev)                   :: zdm                 ! mass of atm. model layer in kg
77    REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
78    REAL                                   :: emission            ! emission
79    REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
80    REAL                                   :: dlat_loc
81
82    IF (is_mpi_root) THEN
83       WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
84       WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit
85    ENDIF
86   
87    DO it=1, nbtr_bin
88      r_bin(it)=mdw(it)/2.
89    ENDDO
90
91!--set boundaries of size bins
92    DO it=1, nbtr_bin
93    IF (it.EQ.1) THEN
94      r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
95      r_lower(it)=r_bin(it)**2./r_upper(it)
96    ELSEIF (it.EQ.nbtr_bin) THEN
97      r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
98      r_upper(it)=r_bin(it)**2./r_lower(it)
99    ELSE
100      r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
101      r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
102    ENDIF
103    ENDDO
104
105    IF (debutphy .and. is_mpi_root) THEN
106      DO it=1, nbtr_bin
107        WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
108      ENDDO
109    ENDIF
110
111!--initialising logical is_strato from stratomask
112    is_strato(:,:)=.FALSE.
113    WHERE (stratomask.GT.0.5) is_strato=.TRUE.
114
115! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
116! H2SO4 mass fraction in aerosol (%)
117    CALL stracomp(sh,t_seri,pplay)
118
119! aerosol density (gr/cm3)
120    CALL denh2sa(t_seri)
121
122! compute factor for converting dry to wet radius (for every grid box)
123    f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
124
125!--calculate mass of air in every grid box
126    DO ilon=1, klon
127    DO ilev=1, klev
128      m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
129    ENDDO
130    ENDDO
131
132!    IF (debutphy) THEN
133!      CALL gather(tr_seri, tr_seri_glo)
134!      IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
135!--initialising tracer concentrations to zero
136!        DO it=1, nbtr
137!        tr_seri(:,:,it)=0.0
138!        ENDDO
139!      ENDIF
140!    ENDIF
141
142!--initialise emission diagnostics
143    budg_emi_ocs(:)=0.0
144    budg_emi_so2(:)=0.0
145    budg_emi_h2so4(:)=0.0
146    budg_emi_part(:)=0.0
147
148!--sulfur emission, depending on chosen scenario (flag_sulf_emit)
149    SELECT CASE(flag_sulf_emit)
150
151    CASE(0) ! background aerosol
152      ! do nothing (no emission)
153
154    CASE(1) ! volcanic eruption
155      !--only emit on day of eruption
156      ! stretch emission over one day of Pinatubo eruption
157       DO ieru=1, nErupt
158          IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
159               day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
160             !
161             ! daily injection mass emission - NL
162             m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
163             WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), &
164                  'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
165                  (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru
166             WRITE(lunout,*) 'IN traccoag, dlon=',dlon
167             DO i=1,klon
168                !Pinatubo eruption at 15.14N, 120.35E
169                dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
170                WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc
171                IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. &
172                     xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN
173                   !
174                   WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur
175                   WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily
176                   !         compute altLMDz
177                   altLMDz(:)=0.0
178                   DO k=1, klev
179                      zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
180                      zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
181                      zdz=zdm(k)/zrho                      !thickness of layer in m
182                      altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
183                   ENDDO
184
185                   SELECT CASE(flag_sulf_emit_distrib)
186                   
187                   CASE(0) ! Gaussian distribution
188                   !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
189                   f_lay_sum=0.0
190                   DO k=1, klev
191                      f_lay_emiss(k)=0.0
192                      DO i_int=1, n_int_alt
193                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
194                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* &
195                              &              exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)*   &
196                              &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
197                      ENDDO
198                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
199                   ENDDO
200                   
201                   CASE(1) ! Uniform distribution
202                   ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the
203                   ! height of the injection, centered around altemiss_vol(ieru)
204                   DO k=1, klev
205                      f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- &
206                      & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru))
207                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
208                   ENDDO
209
210                   END SELECT        ! End CASE over flag_sulf_emit_distrib)
211
212                   WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru)
213                   WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss
214                   !correct for step integration error
215                   f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
216                   !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
217                   !vertically distributed emission
218                   DO k=1, klev
219                      ! stretch emission over one day of Pinatubo eruption
220                      emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
221                      tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
222                      budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
223                   ENDDO
224                ENDIF ! emission grid cell
225             ENDDO ! klon loop
226             WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily
227          ENDIF ! emission period
228       ENDDO ! eruption number
229       
230    CASE(2) ! stratospheric aerosol injections (SAI)
231!
232      DO i=1,klon
233!       SAI standard scenario with continuous emission from 1 grid point at the equator
234!       SAI emission on single month
235!       SAI continuous emission o
236        dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
237        IF  ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. &
238          &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
239!
240!         compute altLMDz
241          altLMDz(:)=0.0
242          DO k=1, klev
243            zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
244            zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
245            zdz=zdm(k)/zrho                      !thickness of layer in m
246            altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
247          ENDDO
248
249          SELECT CASE(flag_sulf_emit_distrib)
250
251          CASE(0) ! Gaussian distribution
252          !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
253          f_lay_sum=0.0
254               DO k=1, klev
255                     f_lay_emiss(k)=0.0
256                     DO i_int=1, n_int_alt
257                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
258                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
259                         &              exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   &
260                         &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
261                     ENDDO
262                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
263               ENDDO
264
265          CASE(1) ! Uniform distribution
266          f_lay_sum=0.0
267          ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
268          ! the height of the injection, centered around altemiss_sai
269               DO k=1, klev
270                    f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
271                    & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
272                    f_lay_sum=f_lay_sum+f_lay_emiss(k)
273               ENDDO
274
275          END SELECT ! Gaussian or uniform distribution
276
277          !correct for step integration error
278          f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
279          !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
280          !vertically distributed emission
281          DO k=1, klev
282            ! stretch emission over whole year (360d)
283            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400. 
284            tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
285            budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
286          ENDDO
287
288!          !emission as monodisperse particles with 0.1um dry radius (BIN21)
289!          !vertically distributed emission
290!          DO k=1, klev
291!            ! stretch emission over whole year (360d)
292!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400.
293!            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
294!            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
295!          ENDDO
296        ENDIF ! emission grid cell
297      ENDDO ! klon loop
298
299    CASE(3) ! --- SAI injection over a single band of longitude and between
300            !     lat_min and lat_max
301
302    DO i=1,klon
303!       SAI scenario with continuous emission
304        dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
305        theta_min = max(xlat(i)-dlat_loc,xlat_min_sai)
306        theta_max = min(xlat(i)+dlat_loc,xlat_max_sai)
307        IF  ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. &
308          &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
309!
310!         compute altLMDz
311          altLMDz(:)=0.0
312          DO k=1, klev
313            zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
314            zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
315            zdz=zdm(k)/zrho                      !thickness of layer in m
316            altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
317          ENDDO
318
319          SELECT CASE(flag_sulf_emit_distrib)
320
321          CASE(0) ! Gaussian distribution
322          !compute distribution of emission to vertical model layers (based on
323          !Gaussian peak in altitude)
324          f_lay_sum=0.0
325               DO k=1, klev
326                     f_lay_emiss(k)=0.0
327                     DO i_int=1, n_int_alt
328                         alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
329                         f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
330                         & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   &
331                         & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
332                     ENDDO
333                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
334               ENDDO
335
336          CASE(1) ! Uniform distribution
337          f_lay_sum=0.0
338          ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
339          ! the height of the injection, centered around altemiss_sai
340               DO k=1, klev
341                    f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
342                    & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
343                    f_lay_sum=f_lay_sum+f_lay_emiss(k)
344               ENDDO
345
346          END SELECT ! Gaussian or uniform distribution
347
348          !correct for step integration error
349          f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
350          !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
351          !vertically distributed emission
352          DO k=1, klev
353            ! stretch emission over whole year (360d)
354            emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ &
355                      & FLOAT(year_len)/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ &
356                      & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI))
357            tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
358            budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
359          ENDDO
360
361!          !emission as monodisperse particles with 0.1um dry radius (BIN21)
362!          !vertically distributed emission
363!          DO k=1, klev
364!            ! stretch emission over whole year (360d)
365!            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400
366!            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
367!            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
368!          ENDDO
369        ENDIF ! emission grid cell
370      ENDDO ! klon loop
371
372    END SELECT ! emission scenario (flag_sulf_emit)
373
374!--read background concentrations of OCS and SO2 and lifetimes from input file
375!--update the variables defined in phys_local_var_mod
376    CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
377
378!--convert OCS to SO2 in the stratosphere
379    CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
380
381!--convert SO2 to H2SO4
382    CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
383
384!--common routine for nucleation and condensation/evaporation with adaptive timestep
385    CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
386
387!--call coagulation routine
388    CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
389
390!--call sedimentation routine
391    CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
392
393!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
394    surf_PM25_sulf(:)=0.0
395    DO i=1,klon
396      DO it=1, nbtr_bin
397        IF (mdw(it) .LT. 2.5e-6) THEN
398          !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
399          !assume that particles consist of ammonium sulfate at the surface (132g/mol)
400          !and are dry at T = 20 deg. C and 50 perc. humidity
401          surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
402                           & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &
403                           & *pplay(i,1)/t_seri(i,1)/RD*1.e9
404        ENDIF
405      ENDDO
406    ENDDO
407
408!    CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2')
409!    DO it=1, nbtr_bin
410!      CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ')
411!    ENDDO
412
413  END SUBROUTINE traccoag
414
415END MODULE traccoag_mod
Note: See TracBrowser for help on using the repository browser.