My Project
 All Classes Files Functions Variables Macros
lsc_scav.F90
Go to the documentation of this file.
1 !$Id $
2 
3 SUBROUTINE lsc_scav(pdtime,it,iflag_lscav,oliq,flxr,flxs,rneb,beta_fisrt, &
4  beta_v1,pplay,paprs,t,tr_seri,d_tr_insc, &
5  d_tr_bcscav,d_tr_evap,qprls)
6  USE ioipsl
7  USE dimphy
8  USE mod_grid_phy_lmdz
9  USE mod_phys_lmdz_para
10  USE traclmdz_mod
11  USE infotrac,ONLY : nbtr
12  USE comgeomphy
13  USE iophy
14  IMPLICIT NONE
15 !=====================================================================
16 ! Objet : depot humide (lessivage et evaporation) de traceurs
17 ! Inspired by routines of Olivier Boucher (mars 1998)
18 ! author R. Pilon 10 octobre 2012
19 ! last modification 16/01/2013 (reformulation partie evaporation)
20 !=====================================================================
21 
22  include "dimensions.h"
23  include "chem.h"
24  include "YOMCST.h"
25  include "YOECUMF.h"
26 
27  REAL,INTENT(IN) :: pdtime ! time step (s)
28  INTEGER,INTENT(IN) :: it ! tracer number
29  INTEGER,INTENT(IN) :: iflag_lscav ! LS scavenging param:
30 ! 3=Reddy_Boucher2004, 4=3+RPilon.
31  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxr ! flux precipitant de pluie
32  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxs ! flux precipitant de neige
33  REAL,INTENT(IN) :: oliq ! contenu en eau liquide dans le nuage (kg/kg)
34  REAL,DIMENSION(klon,klev),INTENT(IN) :: rneb
35  REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression
36  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression
37  REAL,DIMENSION(klon,klev),INTENT(IN) :: t ! temperature
38 ! tracers
39  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! q de traceur
40  REAL,DIMENSION(klon,klev),INTENT(IN) :: beta_fisrt ! taux de conversion de l'eau cond
41  REAL,DIMENSION(klon,klev),INTENT(OUT) :: beta_v1 ! -- (originale version)
42  REAL,DIMENSION(klon) :: his_dh ! tendance de traceur integre verticalement
43  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: d_tr_insc ! tendance du traceur
44  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: d_tr_bcscav ! tendance de traceur
45  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: d_tr_evap
46  REAL,DIMENSION(klon,nbtr),INTENT(OUT) :: qprls !jyg: concentration tra dans pluie LS a la surf.
47  REAL :: dxin,dxev ! tendance temporaire de traceur incloud
48  REAL,DIMENSION(klon,klev) :: dxbc ! tendance temporaire de traceur bc
49 
50 
51 ! variables locales
52  LOGICAL,SAVE :: debut=.true.
53 !$OMP THREADPRIVATE(debut)
54 !
55  REAL,PARAMETER :: henry=1.4 ! constante de Henry en mol/l/atm ~1.4 for gases
56  REAL :: henry_t ! constante de Henry a T t (mol/l/atm)
57  REAL,PARAMETER :: kk=2900. ! coefficient de dependence en T (K)
58  REAL :: f_a ! rapport de la phase aqueuse a la phase gazeuse
59  REAL :: beta ! taux de conversion de l'eau en pluie
60 
61  INTEGER :: i, k
62  REAL,DIMENSION(klon,klev) :: scav ! water liquid content / fraction aqueuse du constituant
63  REAL,DIMENSION(klon,klev) :: zrho
64  REAL,DIMENSION(klon,klev) :: zdz
65  REAL,DIMENSION(klon,klev) :: zmass ! layer mass
66 
67  REAL :: frac_ev ! cste pour la reevaporation : dropplet shrinking
68 ! frac_ev = frac_gas ou frac_aer
69  REAL,PARAMETER :: frac_gas=1.0 ! cste pour la reevaporation pour les gaz
70  REAL :: frac_aer ! cste pour la reevaporation pour les particules
71  REAL,DIMENSION(klon,klev) :: deltap ! P(i+1)-P(i)
72  REAL,DIMENSION(klon,klev) :: beta_ev ! dP/P(i+1)
73 
74 ! 101.325 m3/l x Pa/atm
75 ! R Pa.m3/mol/K
76 ! cste de dissolution pour le depot humide
77  REAL,SAVE :: frac_fine_scav
78  REAL,SAVE :: frac_coar_scav
79 !$OMP THREADPRIVATE(frac_fine_scav, frac_coar_scav)
80 
81 ! below-cloud scav variables
82 ! aerosol : alpha_r=0.001, gas 0.001 (Pruppacher & Klett 1967)
83  REAL,SAVE :: alpha_r ! coefficient d'impaction pour la pluie
84  REAL,SAVE :: alpha_s ! coefficient d'impaction pour la neige
85  REAL,SAVE :: r_r ! mean raindrop radius (m)
86  REAL,SAVE :: r_s ! mean snow crystal radius (m)
87 !$OMP THREADPRIVATE(alpha_r, alpha_s, R_r, R_s)
88  REAL :: pr, ps, ice, water
89  real :: conserv
90 !
91 !!!!!!!!!!!!!!!!!!!! choix lessivage !!!!!!!!!!!!!!!!!!!!!!!!
92 !! logical,save :: inscav_fisrt
93 !!! $OMP THREADPRIVATE(inscav_first)
94 !
95 !!!!!!!!!!!!!!!!!!!!!!!!!!!
96  IF (debut) THEN
97 !
98 ! inscav_fisrt=.true.
99 ! call getin('inscav_fisrt',inscav_fisrt)
100 ! if(inscav_fisrt) then
101 ! print*,'beta from fisrtilp.F90, beta = (z_cond - z_oliq)/z_cond, inscav_fisrt=',inscav_fisrt
102 ! else
103 ! print*,'beta from Reddy and Bocuher 2004 (original version), inscav_fisrt=',inscav_fisrt
104 ! endif
105 !
106  alpha_r=0.001 ! coefficient d'impaction pour la pluie
107  alpha_s=0.01 ! coefficient d'impaction pour la neige
108  r_r=0.001 ! mean raindrop radius (m)
109  r_s=0.001 ! mean snow crystal radius (m)
110  frac_fine_scav=0.7
111  frac_coar_scav=0.7
112 ! frac_aer=0.5 ~ droplet size shrinks by evap
113  frac_aer=0.5
114 !
115  OPEN(99,file='lsc_scav_param.data',status='old', &
116  form='formatted',err=9999)
117  READ(99,*,end=9998) alpha_r
118  READ(99,*,end=9998) alpha_s
119  READ(99,*,end=9998) r_r
120  READ(99,*,end=9998) r_s
121  READ(99,*,end=9998) frac_fine_scav
122  READ(99,*,end=9998) frac_coar_scav
123  READ(99,*,end=9998) frac_aer
124 9998 Continue
125  CLOSE(99)
126 9999 Continue
127 
128  print*,'alpha_r',alpha_r
129  print*,'alpha_s',alpha_s
130  print*,'R_r',r_r
131  print*,'R_s',r_s
132  print*,'frac_fine_scav',frac_fine_scav
133  print*,'frac_coar_scav',frac_coar_scav
134  print*,'frac_aer ev',frac_aer
135 
136 !
137  ENDIF !(debut)
138 !!!!!!!!!!!!!!!!!!!!!!!!!!!
139 !
140 ! initialization
141  dxin=0.
142  dxev=0.
143  beta_ev=0.
144 
145  DO i=1,klon
146  his_dh(i)=0.
147  ENDDO
148 
149  DO k=1,klev
150  DO i=1, klon
151  dxbc(i,k)=0.
152  beta_v1(i,k)=0.
153  deltap(i,k)=0.
154  ENDDO
155  ENDDO
156 
157  DO k=1,klev
158  DO i=1, klon
159  d_tr_insc(i,k,it)=0.
160  d_tr_bcscav(i,k,it)=0.
161  d_tr_evap(i,k,it)=0.
162  ENDDO
163  ENDDO
164 
165 ! pressure and size of the layer
166  DO k=klev-1, 1, -1
167  DO i=1, klon
168  zrho(i,k)=pplay(i,k)/t(i,k)/rd
169  zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/rg
170  zmass(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
171  ENDDO
172  ENDDO
173 
174  IF (it.gt.1) THEN ! aerosol
175  frac_ev=frac_aer
176  ELSE ! gas
177  frac_ev=frac_gas
178  ENDIF
179 
180  IF(it.gt.1) then ! aerosol
181  DO k=1, klev
182  DO i=1, klon
183  scav(i,k)=frac_fine_scav
184  ENDDO
185  ENDDO
186  ELSE ! gas
187  DO k=1, klev
188  DO i=1, klon
189  henry_t=henry*exp(-kk*(1./298.-1./t(i,k))) ! mol/l/atm
190  f_a=henry_t/101.325*r*t(i,k)*oliq*zrho(i,k)/rho_water
191  scav(i,k)=f_a/(1.+f_a)
192  ENDDO
193  ENDDO
194  ENDIF
195 
196  DO k=klev-1, 1, -1
197  DO i=1, klon
198 ! incloud scavenging
199 ! if(inscav_fisrt) then
200  if (iflag_lscav .eq. 4) then
201  beta=beta_fisrt(i,k)*rneb(i,k)
202  else
203  beta=flxr(i,k)-flxr(i,k+1)+flxs(i,k)-flxs(i,k+1)
204 ! beta=beta/zdz(i,k)/oliq/zrho(i,k)
205  beta=beta/zmass(i,k)/oliq
206  beta=max(0.,beta)
207  endif ! (iflag_lscav .eq. 4)
208  beta_v1(i,k)=beta !! for output
209 !
210  dxin=tr_seri(i,k,it)*(exp(-scav(i,k)*beta*pdtime)-1.)
211 ! his_dh(i)=his_dh(i)-dxin*zrho(i,k)*zdz(i,k)/pdtime ! kg/m2/s
212  his_dh(i)=his_dh(i)-dxin*zmass(i,k)/pdtime ! kg/m2/s
213  d_tr_insc(i,k,it)=dxin
214 
215 ! below-cloud impaction
216  IF(it.eq.1) then
217  d_tr_bcscav(i,k,it)=0.
218  ELSE
219  pr=0.5*(flxr(i,k)+flxr(i,k+1))
220  ps=0.5*(flxs(i,k)+flxs(i,k+1))
221  water=pr*alpha_r/r_r/rho_water
222  ice=ps*alpha_s/r_s/rho_ice
223  dxbc(i,k)=-3./4.*tr_seri(i,k,it)*pdtime*(water+ice)
224 ! add tracers from below cloud scav in his_dh
225  his_dh(i)=his_dh(i)-dxbc(i,k)*zmass(i,k)/pdtime ! kg/m2/s
226  d_tr_bcscav(i,k,it)=dxbc(i,k)
227  ENDIF
228 
229 ! reevaporation
230  deltap(i,k)=flxr(i,k+1)+flxs(i,k+1)-flxr(i,k)-flxs(i,k)
231  deltap(i,k)=max(deltap(i,k),0.)
232 
233  if(flxr(i,k+1)+flxs(i,k+1).gt.1.e-16) then
234  beta_ev(i,k)=deltap(i,k)/(flxr(i,k+1)+flxs(i,k+1))
235  else
236  beta_ev(i,k)=0.
237  endif
238 
239  beta_ev(i,k)=max(min(1.,beta_ev(i,k)),0.)
240 
241 !jyg
242 
243  if(abs(1-(1-frac_ev)*beta_ev(i,k)).gt.1.e-16) then
244 ! remove tracers from precipitation owing to release by evaporation in his_dh
245 ! dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
246  dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zmass(i,k)) &
247  /(1 -(1-frac_ev)*beta_ev(i,k))
248  his_dh(i)=his_dh(i)*(1 - frac_ev*beta_ev(i,k) / (1 -(1-frac_ev)*beta_ev(i,k)))
249  else
250 ! dxev=his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k))
251  dxev=his_dh(i) *pdtime/(zmass(i,k))
252  his_dh(i)=0.
253  endif
254 ! print*, k, 'beta_ev',beta_ev
255 ! remove tracers from precipitation owing to release by evaporation in his_dh
256 !! dxev=frac_ev*deltaP(i,k)*pdtime * his_dh(i) /(zrho(i,k)*zdz(i,k))
257 !rplmd
258 !! dxev=frac_ev*deltaP(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
259 !! /max(flxr(i,k)+flxs(i,k),1.e-16)
260 
261 
262  d_tr_evap(i,k,it)=dxev
263 !! tendency is further added in phytrac x = x + dx
264  ENDDO !! do i
265  ENDDO !! do k
266 
267 !jyg (20130114)
268  DO i = 1,klon
269  qprls(i,it) = his_dh(i)/max(flxr(i,1)+flxs(i,1),1.e-16)
270  ENDDO
271 !jyg end
272 
273 
274 ! test de conservation
275  conserv=0.
276 ! DO k= klev,1,-1
277 ! DO i=1, klon
278 ! conserv=conserv+d_tr_insc(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG &
279 ! +d_tr_bcscav(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG &
280 ! +d_tr_evap(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG
281 ! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)'),&
282 ! k,'lsc conserv ',conserv,'insc',d_tr_insc(i,k,it),'bc',d_tr_bcscav(i,k,it),'ev',d_tr_evap(i,k,it)
283 ! ENDDO
284 ! ENDDO
285 
286 END SUBROUTINE lsc_scav