source: LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90 @ 4744

Last change on this file since 4744 was 4601, checked in by dcugnet, 16 months ago

StratAer? commit, N. Lebas
  • Property svn:keywords set to Id
File size: 14.6 KB
RevLine 
[3526]1!
2! $Id: traccoag_mod.F90 4601 2023-06-30 22:07:30Z jyg $
3!
[2690]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, &
[2752]13       t_seri, pplay, paprs, sh, rh, tr_seri)
[2690]14
[2752]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
[2690]17
18    USE dimphy
[4293]19    USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
[2690]20    USE aerophys
[3526]21    USE geometry_mod, ONLY : cell_area, boundslat
[2690]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
[4601]25    USE phys_cal_mod, ONLY : year_len, year_cur, mth_cur, day_cur, hour
[2690]26    USE sulfate_aer_mod
27    USE phys_local_var_mod, ONLY: stratomask
28    USE YOMCST
[3526]29    USE print_control_mod, ONLY: lunout
[4601]30    USE strataer_local_var_mod
31   
[2690]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!----------------
[3526]59    REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
[4601]60    REAL                                   :: m_aer               ! aerosol mass
[3526]61    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
[2690]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
[2699]68    REAL                                   :: alt                 ! altitude for integral calculation
[2690]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
[2752]76    REAL,DIMENSION(klev)                   :: zdm                 ! mass of atm. model layer in kg
[2690]77    REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
[2752]78    REAL                                   :: emission            ! emission
[3526]79    REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
80    REAL                                   :: dlat_loc
[4601]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   
[4513]85    INTEGER                                :: injdur_sai          ! injection duration for SAI case [days]
86    INTEGER                                :: yr, is_bissext
[2690]87
[4601]88    IF (is_mpi_root .AND. flag_verbose_strataer) THEN
[3526]89       WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
[4601]90       WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit
[2690]91    ENDIF
[3526]92   
[2690]93    DO it=1, nbtr_bin
94      r_bin(it)=mdw(it)/2.
95    ENDDO
96
97!--set boundaries of size bins
98    DO it=1, nbtr_bin
99    IF (it.EQ.1) THEN
100      r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
101      r_lower(it)=r_bin(it)**2./r_upper(it)
102    ELSEIF (it.EQ.nbtr_bin) THEN
103      r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
104      r_upper(it)=r_bin(it)**2./r_lower(it)
105    ELSE
106      r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
107      r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
108    ENDIF
109    ENDDO
110
111    IF (debutphy .and. is_mpi_root) THEN
112      DO it=1, nbtr_bin
[3526]113        WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
[2690]114      ENDDO
115    ENDIF
116
117!--initialising logical is_strato from stratomask
118    is_strato(:,:)=.FALSE.
[2695]119    WHERE (stratomask.GT.0.5) is_strato=.TRUE.
[2690]120
121! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
122! H2SO4 mass fraction in aerosol (%)
123    CALL stracomp(sh,t_seri,pplay)
124
125! aerosol density (gr/cm3)
126    CALL denh2sa(t_seri)
127
128! compute factor for converting dry to wet radius (for every grid box)
129    f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
130
131!--calculate mass of air in every grid box
132    DO ilon=1, klon
[4601]133       DO ilev=1, klev
134          m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
135       ENDDO
[2690]136    ENDDO
[4601]137   
[2752]138!--initialise emission diagnostics
[4601]139    budg_emi(:,1)=0.0
[2752]140    budg_emi_ocs(:)=0.0
141    budg_emi_so2(:)=0.0
142    budg_emi_h2so4(:)=0.0
143    budg_emi_part(:)=0.0
144
[4601]145!--sulfur emission, depending on chosen scenario (flag_emit)
146    SELECT CASE(flag_emit)
[2690]147
148    CASE(0) ! background aerosol
149      ! do nothing (no emission)
150
151    CASE(1) ! volcanic eruption
152      !--only emit on day of eruption
153      ! stretch emission over one day of Pinatubo eruption
[3526]154       DO ieru=1, nErupt
155          IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
156               day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
[4601]157
158             ! daily injection mass emission
159             m_aer=m_aer_emiss_vol(ieru,1)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
160             !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
161             m_aer=m_aer*(mSO2mol/mSatom)
162             
163             WRITE(lunout,*) 'IN traccoag m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru,1), &
[3526]164                  'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
[4601]165                  (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer,'ieru=',ieru
[3526]166             WRITE(lunout,*) 'IN traccoag, dlon=',dlon
[4601]167             
168             latmin=xlat_min_vol(ieru)
169             latmax=xlat_max_vol(ieru)
170             lonmin=xlon_min_vol(ieru)
171             lonmax=xlon_max_vol(ieru)
172             altemiss = altemiss_vol(ieru)
173             sigma_alt = sigma_alt_vol(ieru)
174             pdt=pdtphys
175             ! stretch emission over one day of eruption
176             stretchlong = 1.
177             
178             CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
179                  m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
180                  stretchlong,1,0)
181             
[3526]182          ENDIF ! emission period
183       ENDDO ! eruption number
184       
[2690]185    CASE(2) ! stratospheric aerosol injections (SAI)
186!
[4513]187     ! Computing duration of SAI in days...
188     ! ... starting from 0...
189     injdur_sai = 0
190     ! ... then adding whole years from first to (n-1)th...
191     DO yr = year_emit_sai_start, year_emit_sai_end-1
192       ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
193       is_bissext = (MOD(yr,4)==0) .AND. (MOD(yr,100) /= 0 .OR. MOD(yr,400) == 0)
194       injdur_sai = injdur_sai+365+is_bissext
195     ENDDO
196     ! ... then subtracting part of the first year where no injection yet...
197     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)
198     SELECT CASE(mth_emit_sai_start)
199     CASE(2)
200        injdur_sai = injdur_sai-31
201     CASE(3)
202        injdur_sai = injdur_sai-31-28-is_bissext
203     CASE(4)
204        injdur_sai = injdur_sai-31-28-is_bissext-31
205     CASE(5)
206        injdur_sai = injdur_sai-31-28-is_bissext-31-30
207     CASE(6)
208        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31
209     CASE(7)
210        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30
211     CASE(8)
212        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31
213     CASE(9)
214        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31
215     CASE(10)
216        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30
217     CASE(11)
218        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31
219     CASE(12)
220        injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31-30
221     END SELECT
222     injdur_sai = injdur_sai-day_emit_sai_start+1
223     ! ... then adding part of the n-th year
224     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)
225     SELECT CASE(mth_emit_sai_end)
226     CASE(2)
227        injdur_sai = injdur_sai+31
228     CASE(3)
229        injdur_sai = injdur_sai+31+28+is_bissext
230     CASE(4)
231        injdur_sai = injdur_sai+31+28+is_bissext+31
232     CASE(5)
233        injdur_sai = injdur_sai+31+28+is_bissext+31+30
234     CASE(6)
235        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31
236     CASE(7)
237        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30
238     CASE(8)
239        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31
240     CASE(9)
241        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31
242     CASE(10)
243        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30
244     CASE(11)
245        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31
246     CASE(12)
247        injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31+30
248     END SELECT
249     injdur_sai = injdur_sai+day_emit_sai_end
250     ! A security: are SAI dates of injection consistent?
251     IF (injdur_sai <= 0) THEN
252        CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.',1)
253     ENDIF
254     ! Injection in itself
255     IF (( year_emit_sai_start <= year_cur ) &
256        .AND. ( year_cur <= year_emit_sai_end ) &
257        .AND. ( mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur ) &
258        .AND. ( mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end ) &
259        .AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) &
260        .AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN
[4601]261       
262       m_aer=m_aer_emiss_sai
263       !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
264       m_aer=m_aer*(mSO2mol/mSatom)
265       
266       latmin=xlat_sai
267       latmax=xlat_sai
268       lonmin=xlon_sai
269       lonmax=xlon_sai
270       altemiss = altemiss_sai
271       sigma_alt = sigma_alt_sai
272       pdt=0.
273       ! stretch emission over whole year (360d)
274       stretchlong=FLOAT(year_len)
275       
276       CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
277            m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
278            stretchlong,1,0)
279       
280       budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
[4513]281     ENDIF ! Condition over injection dates
282
[3526]283    CASE(3) ! --- SAI injection over a single band of longitude and between
284            !     lat_min and lat_max
285
[4601]286       m_aer=m_aer_emiss_sai
287       !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
288       m_aer=m_aer*(mSO2mol/mSatom)
[3526]289
[4601]290       latmin=xlat_min_sai
291       latmax=xlat_max_sai
292       lonmin=xlon_sai
293       lonmax=xlon_sai
294       altemiss = altemiss_sai
295       sigma_alt = sigma_alt_sai
296       pdt=0.
297       ! stretch emission over whole year (360d)
298       stretchlong=FLOAT(year_len)
[3526]299
[4601]300       CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
301            m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
302            stretchlong,1,0)
[3526]303
[4601]304       budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
305       
306    END SELECT ! emission scenario (flag_emit)
[3526]307
[2690]308!--read background concentrations of OCS and SO2 and lifetimes from input file
[2695]309!--update the variables defined in phys_local_var_mod
310    CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
[2690]311
312!--convert OCS to SO2 in the stratosphere
[2752]313    CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]314
315!--convert SO2 to H2SO4
[2752]316    CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]317
318!--common routine for nucleation and condensation/evaporation with adaptive timestep
319    CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
320
321!--call coagulation routine
322    CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
323
[4601]324!--call sedimentation routine
[2690]325    CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
326
327!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
328    surf_PM25_sulf(:)=0.0
329    DO i=1,klon
330      DO it=1, nbtr_bin
331        IF (mdw(it) .LT. 2.5e-6) THEN
332          !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
[3526]333          !assume that particles consist of ammonium sulfate at the surface (132g/mol)
334          !and are dry at T = 20 deg. C and 50 perc. humidity
[2690]335          surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
336                           & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &
[2752]337                           & *pplay(i,1)/t_seri(i,1)/RD*1.e9
[2690]338        ENDIF
339      ENDDO
340    ENDDO
[4601]341   
[2690]342  END SUBROUTINE traccoag
343
344END MODULE traccoag_mod
Note: See TracBrowser for help on using the repository browser.