source: LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.f90 @ 5435

Last change on this file since 5435 was 5367, checked in by abarral, 3 weeks ago

(WIP) Turn implicit into explicit declarations

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