Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lsc_scav_orig.F90 @ 5434

Last change on this file since 5434 was 2175, checked in by jescribano, 10 years ago
SPLA code included for first time
File size: 9.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format