source: LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/traccoag_mod.F90 @ 5411

Last change on this file since 5411 was 5220, checked in by abarral, 3 months ago

Correct errors induced in r5218 r5219
  • Property svn:keywords set to Id
File size: 19.3 KB
Line 
1MODULE traccoag_mod
2
3  ! This module calculates the concentration of aerosol particles in certain size bins
4  ! considering coagulation and sedimentation.
5  USE lmdz_abort_physic, ONLY: abort_physic
6CONTAINS
7
8  SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
9          presnivs, xlat, xlon, pphis, pphi, &
10          t_seri, pplay, paprs, sh, rh, tr_seri)
11
12    USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
13            budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, &
14            R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode, reff_sulfate
15    USE dimphy
16    USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
17    USE aerophys
18    USE lmdz_geometry, ONLY: cell_area, boundslat
19    USE lmdz_grid_phy
20    USE lmdz_phys_mpi_data, ONLY: is_mpi_root
21    USE lmdz_phys_para, ONLY: gather, scatter
22    USE phys_cal_mod, ONLY: year_len, year_cur, mth_cur, day_cur, hour
23    USE sulfate_aer_mod
24    USE phys_local_var_mod, ONLY: stratomask
25    USE lmdz_yomcst
26    USE lmdz_print_control, ONLY: lunout
27    USE strataer_local_var_mod
28
29    IMPLICIT NONE
30
31    ! Input argument
32    !---------------
33    REAL, INTENT(IN) :: pdtphys    ! Pas d'integration pour la physique (seconde)
34    REAL, INTENT(IN) :: gmtime     ! Heure courante
35    LOGICAL, INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
36    INTEGER, INTENT(IN) :: julien     ! Jour julien
37
38    REAL, DIMENSION(klev), INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
39    REAL, DIMENSION(klon), INTENT(IN) :: xlat    ! latitudes pour chaque point
40    REAL, DIMENSION(klon), INTENT(IN) :: xlon    ! longitudes pour chaque point
41    REAL, DIMENSION(klon), INTENT(IN) :: pphis   ! geopotentiel du sol
42    REAL, DIMENSION(klon, klev), INTENT(IN) :: pphi    ! geopotentiel de chaque couche
43
44    REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri  ! Temperature
45    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay   ! pression pour le mileu de chaque couche (en Pa)
46    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
47    REAL, DIMENSION(klon, klev), INTENT(IN) :: sh      ! humidite specifique
48    REAL, DIMENSION(klon, klev), INTENT(IN) :: rh      ! humidite relative
49
50    ! Output argument
51    !----------------
52    REAL, DIMENSION(klon, klev, nbtr), INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
53
54    ! Local variables
55    !----------------
56    REAL :: m_aer_emiss_vol_daily ! daily injection mass emission
57    REAL :: m_aer               ! aerosol mass
58    INTEGER :: it, k, i, j, ilon, ilev, itime, i_int, ieru
59    LOGICAL, DIMENSION(klon, klev) :: is_strato           ! true = above tropopause, false = below
60    REAL, DIMENSION(klon, klev) :: m_air_gridbox       ! mass of air in every grid box [kg]
61    REAL, DIMENSION(klon_glo, klev, nbtr) :: tr_seri_glo         ! Concentration Traceur [U/KgA]
62    REAL, DIMENSION(klev + 1) :: altLMDz             ! altitude of layer interfaces in m
63    REAL, DIMENSION(klev) :: f_lay_emiss         ! fraction of emission for every vertical layer
64    REAL :: f_lay_sum           ! sum of layer emission fractions
65    REAL :: alt                 ! altitude for integral calculation
66    INTEGER, PARAMETER :: n_int_alt = 10        ! number of subintervals for integration over Gaussian emission profile
67    REAL, DIMENSION(nbtr_bin) :: r_bin               ! particle radius in size bin [m]
68    REAL, DIMENSION(nbtr_bin) :: r_lower             ! particle radius at lower bin boundary [m]
69    REAL, DIMENSION(nbtr_bin) :: r_upper             ! particle radius at upper bin boundary [m]
70    REAL, DIMENSION(nbtr_bin) :: m_part_dry          ! mass of one dry particle in size bin [kg]
71    REAL :: zrho                ! Density of air [kg/m3]
72    REAL :: zdz                 ! thickness of atm. model layer in m
73    REAL, DIMENSION(klev) :: zdm                 ! mass of atm. model layer in kg
74    REAL, DIMENSION(klon, klev) :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
75    REAL :: emission            ! emission
76    REAL :: theta_min, theta_max ! for SAI computation between two latitudes
77    REAL :: dlat_loc
78    REAL :: latmin, latmax, lonmin, lonmax ! lat/lon min/max for injection
79    REAL :: sigma_alt, altemiss ! injection altitude + sigma for distrib
80    REAL :: pdt, stretchlong     ! physic timestep, stretch emission over one day
81
82    INTEGER :: injdur_sai          ! injection duration for SAI case [days]
83    INTEGER :: yr, is_bissext
84    REAL :: samoment2, samoment3! 2nd and 3rd order moments of size distribution
85
86    IF (is_mpi_root .AND. flag_verbose_strataer) THEN
87      WRITE(lunout, *) 'in traccoag: date from phys_cal_mod =', year_cur, '-', mth_cur, '-', day_cur, '-', hour
88      WRITE(lunout, *) 'IN traccoag flag_emit: ', flag_emit
89    ENDIF
90
91    !   radius [m]
92    DO it = 1, nbtr_bin
93      r_bin(it) = mdw(it) / 2.
94    ENDDO
95
96    !--set boundaries of size bins
97    DO it = 1, nbtr_bin
98      IF (it==1) THEN
99        r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it))
100        r_lower(it) = r_bin(it)**2. / r_upper(it)
101      ELSEIF (it==nbtr_bin) THEN
102        r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1))
103        r_upper(it) = r_bin(it)**2. / r_lower(it)
104      ELSE
105        r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1))
106        r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it))
107      ENDIF
108    ENDDO
109
110    IF (debutphy .AND. is_mpi_root) THEN
111      DO it = 1, nbtr_bin
112        WRITE(lunout, *) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
113      ENDDO
114    ENDIF
115
116    !--initialising LOGICAL is_strato from stratomask
117    is_strato(:, :) = .FALSE.
118    WHERE (stratomask>0.5) is_strato = .TRUE.
119
120    IF(flag_new_strat_compo) THEN
121      IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo
122      ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
123      CALL stracomp_kelvin(sh, t_seri, pplay)
124    ELSE
125      IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo
126      ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
127      ! H2SO4 mass fraction in aerosol (%)
128      CALL stracomp(sh, t_seri, pplay)
129
130      ! aerosol density (gr/cm3)
131      CALL denh2sa(t_seri)
132
133      ! compute factor for converting dry to wet radius (for every grid box)
134      f_r_wet(:, :) = (dens_aer_dry / (DENSO4(:, :) * 1000.) / (R2SO4(:, :) / 100.))**(1. / 3.)
135    ENDIF
136
137    !--calculate mass of air in every grid box
138    DO ilon = 1, klon
139      DO ilev = 1, klev
140        m_air_gridbox(ilon, ilev) = (paprs(ilon, ilev) - paprs(ilon, ilev + 1)) / RG * cell_area(ilon)
141      ENDDO
142    ENDDO
143
144    !--initialise emission diagnostics
145    IF (nErupt > 0 .AND. (flag_emit == 1 .OR. flag_emit == 4)) budg_emi(:, 1) = 0.0
146    budg_emi_ocs(:) = 0.0
147    budg_emi_so2(:) = 0.0
148    budg_emi_h2so4(:) = 0.0
149    budg_emi_part(:) = 0.0
150
151    !--sulfur emission, depending on chosen scenario (flag_emit)
152    SELECT CASE(flag_emit)
153
154    CASE(0) ! background aerosol
155      ! do nothing (no emission)
156
157    CASE(1) ! volcanic eruption
158      !--only emit on day of eruption
159      ! stretch emission over one day of Pinatubo eruption
160      DO ieru = 1, nErupt
161        IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
162                day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru) + injdur)) THEN
163
164          ! daily injection mass emission
165          m_aer = m_aer_emiss_vol(ieru, 1) / (REAL(injdur) * REAL(ponde_lonlat_vol(ieru)))
166          !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
167          m_aer = m_aer * (mSO2mol / mSatom)
168
169          WRITE(lunout, *) 'IN traccoag m_aer_emiss_vol(ieru)=', m_aer_emiss_vol(ieru, 1), &
170                  'ponde_lonlat_vol(ieru)=', ponde_lonlat_vol(ieru), '(injdur*ponde_lonlat_vol(ieru))', &
171                  (injdur * ponde_lonlat_vol(ieru)), 'm_aer_emiss_vol_daily=', m_aer, 'ieru=', ieru
172          WRITE(lunout, *) 'IN traccoag, dlon=', dlon
173
174          latmin = xlat_min_vol(ieru)
175          latmax = xlat_max_vol(ieru)
176          lonmin = xlon_min_vol(ieru)
177          lonmax = xlon_max_vol(ieru)
178          altemiss = altemiss_vol(ieru)
179          sigma_alt = sigma_alt_vol(ieru)
180          pdt = pdtphys
181          ! stretch emission over one day of eruption
182          stretchlong = 1.
183
184          CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, tr_seri, &
185                  m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, &
186                  stretchlong, 1, 0)
187
188        ENDIF ! emission period
189      ENDDO ! eruption number
190
191    CASE(2) ! stratospheric aerosol injections (SAI)
192
193      ! Computing duration of SAI in days...
194      ! ... starting from 0...
195      injdur_sai = 0
196      ! ... then adding whole years from first to (n-1)th...
197      DO yr = year_emit_sai_start, year_emit_sai_end - 1
198        ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
199        is_bissext = (MOD(yr, 4)==0) .AND. (MOD(yr, 100) /= 0 .OR. MOD(yr, 400) == 0)
200        injdur_sai = injdur_sai + 365 + is_bissext
201      ENDDO
202      ! ... then subtracting part of the first year where no injection yet...
203      is_bissext = (MOD(year_emit_sai_start, 4)==0) .AND. (MOD(year_emit_sai_start, 100) /= 0 .OR. MOD(year_emit_sai_start, 400) == 0)
204      SELECT CASE(mth_emit_sai_start)
205      CASE(2)
206        injdur_sai = injdur_sai - 31
207      CASE(3)
208        injdur_sai = injdur_sai - 31 - 28 - is_bissext
209      CASE(4)
210        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31
211      CASE(5)
212        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30
213      CASE(6)
214        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31
215      CASE(7)
216        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30
217      CASE(8)
218        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31
219      CASE(9)
220        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31
221      CASE(10)
222        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30
223      CASE(11)
224        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31
225      CASE(12)
226        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31 - 30
227      END SELECT
228      injdur_sai = injdur_sai - day_emit_sai_start + 1
229      ! ... then adding part of the n-th year
230      is_bissext = (MOD(year_emit_sai_end, 4)==0) .AND. (MOD(year_emit_sai_end, 100) /= 0 .OR. MOD(year_emit_sai_end, 400) == 0)
231      SELECT CASE(mth_emit_sai_end)
232      CASE(2)
233        injdur_sai = injdur_sai + 31
234      CASE(3)
235        injdur_sai = injdur_sai + 31 + 28 + is_bissext
236      CASE(4)
237        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31
238      CASE(5)
239        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30
240      CASE(6)
241        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31
242      CASE(7)
243        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30
244      CASE(8)
245        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31
246      CASE(9)
247        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31
248      CASE(10)
249        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30
250      CASE(11)
251        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31
252      CASE(12)
253        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30
254      END SELECT
255      injdur_sai = injdur_sai + day_emit_sai_end
256      ! A security: are SAI dates of injection consistent?
257      IF (injdur_sai <= 0) THEN
258        CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.', 1)
259      ENDIF
260      ! Injection in itself
261      IF ((year_emit_sai_start <= year_cur) &
262              .AND. (year_cur <= year_emit_sai_end) &
263              .AND. (mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur) &
264              .AND. (mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end) &
265              .AND. (day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur) &
266              .AND. (day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end)) THEN
267
268        m_aer = m_aer_emiss_sai
269        !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
270        m_aer = m_aer * (mSO2mol / mSatom)
271
272        latmin = xlat_sai
273        latmax = xlat_sai
274        lonmin = xlon_sai
275        lonmax = xlon_sai
276        altemiss = altemiss_sai
277        sigma_alt = sigma_alt_sai
278        pdt = 0.
279        ! stretch emission over whole year (360d)
280        stretchlong = FLOAT(year_len)
281
282        CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, m_air_gridbox, tr_seri, &
283                m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, &
284                stretchlong, 1, 0)
285
286        budg_emi_so2(:) = budg_emi(:, 1) * mSatom / mSO2mol
287      ENDIF ! Condition over injection dates
288
289    CASE(3) ! --- SAI injection over a single band of longitude and between
290      !     lat_min and lat_max
291
292      m_aer = m_aer_emiss_sai
293      !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
294      m_aer = m_aer * (mSO2mol / mSatom)
295
296      latmin = xlat_min_sai
297      latmax = xlat_max_sai
298      lonmin = xlon_sai
299      lonmax = xlon_sai
300      altemiss = altemiss_sai
301      sigma_alt = sigma_alt_sai
302      pdt = 0.
303      ! stretch emission over whole year (360d)
304      stretchlong = FLOAT(year_len)
305
306      CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, m_air_gridbox, tr_seri, &
307              m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, &
308              stretchlong, 1, 0)
309
310      budg_emi_so2(:) = budg_emi(:, 1) * mSatom / mSO2mol
311
312    END SELECT ! emission scenario (flag_emit)
313
314    !--read background concentrations of OCS and SO2 and lifetimes from input file
315    !--update the variables defined in phys_local_var_mod
316    CALL interp_sulf_input(debutphy, pdtphys, paprs, tr_seri)
317
318    !--convert OCS to SO2 in the stratosphere
319    CALL ocs_to_so2(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato)
320
321    !--convert SO2 to H2SO4
322    CALL so2_to_h2so4(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato)
323
324    !--common routine for nucleation and condensation/evaporation with adaptive timestep
325    CALL micphy_tstep(pdtphys, tr_seri, t_seri, pplay, paprs, rh, is_strato)
326
327    !--CALL coagulation routine
328    CALL coagulate(pdtphys, mdw, tr_seri, t_seri, pplay, dens_aer, is_strato)
329
330    !--CALL sedimentation routine
331    CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
332
333    !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
334    surf_PM25_sulf(:) = 0.0
335    DO i = 1, klon
336      DO it = 1, nbtr_bin
337        IF (mdw(it) < 2.5e-6) THEN
338          !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
339          !assume that particles consist of ammonium sulfate at the surface (132g/mol)
340          !and are dry at T = 20 deg. C and 50 perc. humidity
341          surf_PM25_sulf(i) = surf_PM25_sulf(i) + tr_seri(i, 1, it + nbtr_sulgas) &
342                  * 132. / 98. * dens_aer_dry * 4. / 3. * RPI * (mdw(it) / 2.)**3 &
343                  * pplay(i, 1) / t_seri(i, 1) / RD * 1.e9
344        ENDIF
345      ENDDO
346    ENDDO
347
348    !--compute
349    !     sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr)
350    !     SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET)
351    !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr)
352    !     nd_mode: DRY(?) particle concentration in different modes (part/m3)
353    sulfmmr(:, :) = 0.0
354    SAD_sulfate(:, :) = 0.0
355    sulfmmr_mode(:, :, :) = 0.0
356    nd_mode(:, :, :) = 0.0
357    reff_sulfate(:, :) = 0.0
358
359    DO i = 1, klon
360      DO j = 1, klev
361        samoment2 = 0.0
362        samoment3 = 0.0
363        DO it = 1, nbtr_bin
364          !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
365          !assume that particles consist of ammonium sulfate at the surface (132g/mol)
366          !and are dry at T = 20 deg. C and 50 perc. humidity
367
368          !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr)
369          !     equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it
370          sulfmmr_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) &        ! [DRY part/kgA in bin it]
371                  * (4. / 3.) * RPI * (mdw(it) / 2.)**3.   &                   ! [mdw: dry diameter in m]
372                  * dens_aer_dry                                       ! [dry aerosol mass density in kg/m3]
373
374          !     sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio)
375          !     equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA
376          sulfmmr(i, j) = sulfmmr(i, j) + sulfmmr_mode(i, j, it)
377
378          !     nd_mode: particle concentration in different modes (DRY part/m3)
379          nd_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) &             ! [DRY part/kgA in bin it]
380                  * pplay(i, j) / t_seri(i, j) / RD                           ! [air mass concentration in kg air /m3A]
381
382          IF(flag_new_strat_compo) THEN
383            !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
384            SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
385                    * 4. * RPI * (mdw(it) * f_r_wetB(i, j, it) / 2.)**2. &       ! [WET SA of part it in m2]
386                    * 1.e-2                                              ! conversion from m2/m3 to cm2/cm3A!    samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3)
387            samoment2 = samoment2 + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
388                    & * (mdw(it) * f_r_wetB(i, j, it) / 2.)**2.                     ! [WET SA of part it in m2]
389            !    samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3)
390            samoment3 = samoment3 + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
391                    & * (mdw(it) * f_r_wetB(i, j, it) / 2.)**3.                     ! [WET SA of part it in m2]
392          ELSE
393            !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
394            SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
395                    * 4. * RPI * (mdw(it) * f_r_wet(i, j) / 2.)**2. &           ! [WET SA of part it in m2]
396                    * 1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
397            !    samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3)
398            samoment2 = samoment2 + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
399                    & * (mdw(it) * f_r_wet(i, j) / 2.)**2.                          ! [WET SA of part it in m2]
400            !    samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3)
401            samoment3 = samoment3 + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
402                    & * (mdw(it) * f_r_wet(i, j) / 2.)**3.                          ! [WET SA of part it in m2]
403
404          END IF
405        END DO!     reff_sulfate: effective radius of WET sulfate aerosols (cm)
406        reff_sulfate(i, j) = (samoment3 / samoment2) &
407                & * 1.e2                                              ! conversion from m to cm
408      END DO
409    END DO
410  END SUBROUTINE traccoag
411
412END MODULE traccoag_mod
Note: See TracBrowser for help on using the repository browser.