Context Navigation

new_cloud_sedim.F @ 3556

Last change on this file since 3556 was 3323, checked in by flefevre, 8 months ago
1) revised cloud microphysical parameters (this changes the results) 2) the number of modes is passed as an argument to cloud routines (this does not change the results) 3) some cosmetics to chemparam_mod.F90 (this does not change the results)
Property svn:executable set to ``*
File size: 14.7 KB

Line
1	subroutine new_cloud_sedim(nbr_mode, n_lon, n_lev, ptimestep,
2	$ pmidlay, pbndlay, pt, pq,
3	$ d_tr_chem, pdqsed,
4	$ nq, F_sed)
5
6	USE ioipsl
7	USE dimphy
8	USE chemparam_mod
9	IMPLICIT NONE
10
11	c-----------------------------------------------------------------------
12	c declarations:
13	c -------------
14	#include "YOMCST.h"
15	c#include "dimphys.h"
16	c#include "comcstfi.h"
17	c#include "tracer.h"
18	c#include "callkeys.h"
19	c
20	c arguments:
21	c ----------
22
23	INTEGER n_lon ! number of horizontal grid points
24	INTEGER n_lev ! number of atmospheric layers
25	integer nbr_mode
26	REAL ptimestep ! physics time step (s)
27	REAL pmidlay(n_lon,n_lev) ! pressure at middle layers (Pa)
28	REAL pt(n_lon,n_lev) ! temperature at mid-layer (l)
29	REAL pbndlay(n_lon,n_lev+1) ! pressure at layer boundaries
30
31	c Traceurs :
32	integer nq ! number of tracers
33	real pq(n_lon,n_lev,nq) ! tracers (kg/kg)
34	real pdqsed(n_lon,n_lev,2) ! tendency due to sedimentation (kg/kg)
35	real d_tr_chem(n_lon,n_lev,nq)! tendency due to chemistry and clouds (kg/kg)
36
37	c local:
38	c ------
39	integer imode
40	integer ig
41	integer iq
42	integer l
43
44	real zlev(n_lon,n_lev+1) ! altitude at layer boundaries
45	real zlay(n_lon,n_lev) ! altitude at the midlle layer
46	real zqi_wv(n_lon,n_lev) ! to locally store H2O tracer
47	real zqi_sa(n_lon,n_lev) ! to locally store H2SO4 tracer
48	real m_lay (n_lon,n_lev) ! Layer Pressure over gravity (Dp/g == kg.m-2)
49	real wq(n_lon,n_lev+1) ! displaced tracer mass (kg.m-2)
50
51	c Physical constant
52	c ~~~~~~~~~~~~~~~~~
53	c Gas molecular viscosity (N.s.m-2)
54	c real,parameter :: visc=1.e-5 ! CO2
55	REAL :: VISCOSITY_CO2
56	c Effective gas molecular radius (m)
57	real,parameter :: molrad=2.2e-10 ! CO2
58
59	REAL :: qmass
60	c masse volumique du coeur (kg.m-3)
61	c ATTENTION ! DOIT ETRE COHERENT AVEC new_cloud_venus !
62	REAL, PARAMETER :: rho_core = 2500.0
63
64	REAL, DIMENSION(n_lon,n_lev+1) :: wgt_SA ! Fraction of H2SO4 in droplet local
65
66	c Stokes speed and sedimentation flux variable
67	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
68
69	REAL :: A1,A2,A3,A4, ! coeff du DL du Flux de sedimentation
70	+ D_stokes, ! coeff de la vitesse de Stokes
71	+ Rp_DL, ! "Point" du DL
72	+ l_mean, ! libre parcours moyen (m)
73	+ a,b_exp,c ! coeff du calcul du Flux de sedimentation
74	REAL, DIMENSION(n_lon,n_lev+1) ::
75	+ F_sed ! Flux de sedimentation (kg.m-2.s-1 puis en output kg.m-2)
76
77	REAL :: R_mode0 ! Rayon mode 0 (m), rayon le plus frequent
78
79	! PRINT*,'RHO_DROPLET new_cloud_sedim.F'
80	! PRINT*,'rho_droplet',rho_droplet(16,21)
81	! PRINT*,'T',pt(16,21),'WSA',WH2SO4(16,21)
82
83	c-----------------------------------------------------------------------
84	c 1. Initialization
85	c -----------------
86
87	! update water vapour and sulfuric acid mixing ratios
88
89	zqi_wv(:,:) = pq(:,:,i_h2oliq) + d_tr_chem(:,:,i_h2oliq)*ptimestep
90	zqi_sa(:,:) = pq(:,:,i_h2so4liq)
91	$ + d_tr_chem(:,:,i_h2so4liq)*ptimestep
92
93	wgt_SA(:,1:n_lev) = wh2so4(:,:)
94
95	c Init F_sed
96	F_sed(:,:) = 0.0E+0
97
98	c Au niveau top+1 , tout égal a 0
99	wgt_SA(:,n_lev+1) = 0.0E+0
100
101	c Computing the different layer properties
102	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
103	c m_lay (kg.m-2)
104	c Ici g=8.87, conflit pour g entre #include "YOMCST.h"
105	c et #include "comcstfi.h"
106
107	do l=1,n_lev
108	do ig=1, n_lon
109	m_lay(ig,l)=(pbndlay(ig,l) - pbndlay(ig,l+1)) /8.87E+0
110	IF (m_lay(ig,l).LE.0.0) THEN
111	PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
112	PRINT*,'!!!! m_lay <= 0 !!!!'
113	PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
114	ENDIF
115	end do
116	end do
117
118	c Computing sedimentation for droplet "tracer"
119	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
120	c pbndlay(:,51)=0 (en parallèle c'est sûr), ne pas l'utiliser pour Fse
121
122	c Sedimentation pour une gouttelette mode 3 de type J. Cimino, 1982, Icarus
123	c c.a.d 97% radius due a solide 3% radius acide sulfurique
124	DO imode=1, nbr_mode - 1
125	DO l = cloudmin, cloudmax
126	DO ig=1,n_lon
127
128	c RD=1000.RNAVORKBOL/RMD avec RMD=43.44 Masse molaire atm venus en g.mol-1
129	D_stokes=((rho_droplet(ig,l)-pmidlay(ig,l)/(RD*pt(ig,l))))
130	& (2./9.)(RG/VISCOSITY_CO2(pt(ig,l)))
131
132	l_mean=(pt(ig,l)/pmidlay(ig,l))*
133	& (0.707R/(4.RPI* molradmolrad RNAVO))
134
135	R_mode0=R_MEDIAN(ig,l,imode)*
136	& EXP(-LOG(STDDEV(ig,l,imode))**2.)
137	IF ((l_mean/(R_mode0)).GT.10.) THEN
138	Rp_DL=R_MEDIAN(ig,l,imode)*
139	& EXP(3.LOG(STDDEV(ig,l,imode))*2.)
140	ELSE
141	Rp_DL=R_MEDIAN(ig,l,imode)*
142	& EXP(4.LOG(STDDEV(ig,l,imode))*2.)
143	ENDIF
144
145	a=1.246*l_mean
146
147	c=0.87/l_mean
148
149	b_exp=0.42l_meanEXP(-c*Rp_DL)
150
151	A1=a+b_exp(1.+cRp_DL
152	& +0.5(Rp_DLc)**2
153	& +1./6.(Rp_DLc)**3)
154
155	A2=1.-b_exp*(c
156	& +Rp_DLc*2
157	& +0.5(Rp_DL2)(c**3))
158
159	A3=0.5b_exp(c*2+Rp_DLc**3)
160
161	A4=-b_exp1./6.c**3
162
163	c Addition des Flux de tous les modes presents
164	F_sed(ig,l)=F_sed(ig,l)+(rho_droplet(ig,l)4./3.RPI*
165	& NBRTOT(ig,l,imode)1.0E6D_stokes*(
166	& A1R_MEDIAN(ig,l,imode)*4
167	& EXP(8.0LOG(STDDEV(ig,l,imode))**2.)
168	& +A2R_MEDIAN(ig,l,imode)*5
169	& EXP(12.5LOG(STDDEV(ig,l,imode))**2.)
170	& +A3R_MEDIAN(ig,l,imode)*6
171	& EXP(18.0LOG(STDDEV(ig,l,imode))**2.)
172	& +A4R_MEDIAN(ig,l,imode)*7
173	& EXP(24.5LOG(STDDEV(ig,l,imode))**2.)))
174
175	c PRINT*,' APRES dTime: F_sed=',F_sed(ig,l), ig, l
176
177	IF (F_sed(ig,l).GT.m_lay(ig,l)) THEN
178	PRINT*,'==============================================='
179	PRINT*,'WARNING On a epuise la couche', ig, l
180	PRINT*,'On epuise pas une couche avec une espèce
181	& minoritaire, c est pas bien maaaaaal'
182	PRINT*,'Water',zqi_wv(ig,l),'Sulfuric Acid',zqi_sa(ig,l)
183	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
184	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
185	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
186	& rho_droplet(ig,l)
187	PRINT*,'Ntot',NBRTOT(ig,l,:)
188	PRINT*,'StdDev',STDDEV(ig,l,:),'Rmed',R_MEDIAN(ig,l,:)
189	PRINT*,'K_MASS',K_MASS(ig,l,:)
190	PRINT*,'WSA',WH2SO4(ig,l),'RHO',rho_droplet(ig,l)
191
192	c ELSE
193	c
194	c PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
195	c PRINT*,'WARNING On a PAS epuise la couche', ig, l
196	c PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
197	c PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
198	c PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
199	c & rho_droplet(ig,l)(ig,l)
200	c PRINT,'Ntot',NBRTOT(ig,l),'Ntot m3',NBRTOT(ig,l)1.0e6
201	c PRINT*,'StdDev',STDDEV(ig,l),'Rmed',R_MEDIAN(ig,l)
202	STOP
203	ENDIF
204
205	IF (F_sed(ig,l).LT.0.0d0) THEN
206	PRINT*,"F_sed est négatif !!!"
207	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
208	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
209	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Pmid',pmidlay(ig,l)
210	PRINT*,'Temp',pt(ig,l),'Rho',
211	& rho_droplet(ig,l)
212	PRINT*,'Ntot',NBRTOT(ig,l,imode),'Ntot m3',
213	& NBRTOT(ig,l,imode)*1.0e6
214	PRINT*,'StdDev',STDDEV(ig,l,imode),'Rmed',
215	& R_MEDIAN(ig,l,imode)
216	PRINT*,'A1',A1,'A2',A2
217	PRINT*,'A3',A1,'A4',A2
218	PRINT*,'D_stokes',D_stokes
219	STOP
220	ENDIF
221
222	ENDDO
223
224	c ELSE
225	c F_sed(:,l)=0.0d0
226	c ENDIF
227
228	ENDDO
229	ENDDO
230
231	c****************************************************************
232	c On calcule le F_sed du mode 3 + coeff*(Fsed1 + Fsed2)
233	c****************************************************************
234	DO l = cloudmin, cloudmax
235	DO ig=1,n_lon
236
237	c calcul de qmass
238	qmass=(rho_coreqrad*3)/
239	& (rho_coreqrad3+rho_droplet(ig,l)(1.-qrad**3))
240
241	c RD=1000.RNAVORKBOL/RMD avec RMD=43.44 Masse molaire atm venus en g.mol-1
242	D_stokes=(((qmassrho_core+(1.-qmass)rho_droplet(ig,l))
243	& -pmidlay(ig,l)/(RD*pt(ig,l))))
244	& (2./9.)(RG/VISCOSITY_CO2(pt(ig,l)))
245
246	l_mean=(pt(ig,l)/pmidlay(ig,l))*
247	& (0.707R/(4.RPI* molradmolrad RNAVO))
248
249	R_mode0=R_MEDIAN(ig,l,3)*
250	& EXP(-LOG(STDDEV(ig,l,3))**2.)
251	IF ((l_mean/(R_mode0)).GT.10.) THEN
252	Rp_DL=R_MEDIAN(ig,l,3)*
253	& EXP(3.LOG(STDDEV(ig,l,3))*2.)
254	ELSE
255	Rp_DL=R_MEDIAN(ig,l,3)*
256	& EXP(4.LOG(STDDEV(ig,l,3))*2.)
257	ENDIF
258
259	a=1.246*l_mean
260
261	c=0.87/l_mean
262
263	b_exp=0.42l_meanEXP(-c*Rp_DL)
264
265	A1=a+b_exp(1.+cRp_DL
266	& +0.5(Rp_DLc)**2
267	& +1./6.(Rp_DLc)**3)
268
269	A2=1.-b_exp*(c
270	& +Rp_DLc*2
271	& +0.5(Rp_DL2)(c**3))
272
273	A3=0.5b_exp(c*2+Rp_DLc**3)
274
275	A4=-b_exp1./6.c**3
276
277	c Addition des Flux de tous les modes presents
278	F_sed(ig,l)=F_sed(ig,l)
279	& +((1.-qmass)/(1.-qmassK_MASS(ig,l,3)))(
280	& (qmassrho_core+(1.-qmass)rho_droplet(ig,l))4./3.RPI*
281	& NBRTOT(ig,l,3)1.0E6D_stokes*(
282	& A1R_MEDIAN(ig,l,3)*4
283	& EXP(8.0LOG(STDDEV(ig,l,3))**2.)
284	& +A2R_MEDIAN(ig,l,3)*5
285	& EXP(12.5LOG(STDDEV(ig,l,3))**2.)
286	& +A3R_MEDIAN(ig,l,3)*6
287	& EXP(18.0LOG(STDDEV(ig,l,3))**2.)
288	& +A4R_MEDIAN(ig,l,3)*7
289	& EXP(24.5LOG(STDDEV(ig,l,3))**2.)))
290
291	c PRINT*,' APRES dTime: F_sed=',F_sed(ig,l), ig, l
292
293	IF (F_sed(ig,l).GT.m_lay(ig,l)) THEN
294	PRINT*,'==============================================='
295	PRINT*,'WARNING On a epuise la couche', ig, l
296	PRINT*,'On epuise pas une couche avec une espèce
297	& minoritaire, c est pas bien maaaaaal'
298	PRINT*,'Water',zqi_wv(ig,l),'Sulfuric Acid',zqi_sa(ig,l)
299	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
300	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
301	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
302	& rho_droplet(ig,l)
303	PRINT*,'Ntot',NBRTOT(ig,l,:)
304	PRINT*,'StdDev',STDDEV(ig,l,:),'Rmed',R_MEDIAN(ig,l,:)
305	PRINT*,'K_MASS',K_MASS(ig,l,:)
306	PRINT*,'WSA',WH2SO4(ig,l),'RHO',rho_droplet(ig,l)
307
308	c ELSE
309	c
310	c PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
311	c PRINT*,'WARNING On a PAS epuise la couche', ig, l
312	c PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
313	c PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
314	c PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
315	c & rho_droplet(ig,l)(ig,l)
316	c PRINT,'Ntot',NBRTOT(ig,l),'Ntot m3',NBRTOT(ig,l)1.0e6
317	c PRINT*,'StdDev',STDDEV(ig,l),'Rmed',R_MEDIAN(ig,l)
318	STOP
319	ENDIF
320
321	IF (F_sed(ig,l).LT.0.0d0) THEN
322	PRINT*,"F_sed est négatif !!!"
323	PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
324	PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
325	PRINT*,'Pbnd top',pbndlay(ig,l+1),'Pmid',pmidlay(ig,l)
326	PRINT*,'Temp',pt(ig,l),'Rho',
327	& rho_droplet(ig,l)
328	PRINT*,'Ntot',NBRTOT(ig,l,imode),'Ntot m3',
329	& NBRTOT(ig,l,imode)*1.0e6
330	PRINT*,'StdDev',STDDEV(ig,l,imode),'Rmed',
331	& R_MEDIAN(ig,l,imode)
332	PRINT*,'A1',A1,'A2',A2
333	PRINT*,'A3',A1,'A4',A2
334	PRINT*,'D_stokes',D_stokes
335	STOP
336	ENDIF
337
338	ENDDO
339
340	c ELSE
341	c F_sed(:,l)=0.0d0
342	c ENDIF
343
344	ENDDO
345
346	c Passage du Flux au Flux pour un pas de temps (== kg.m-2)
347
348	F_sed(:,:) = F_sed(:,:)*ptimestep
349
350	!=========================================================
351	! compute tendency due to sedimentation
352	!=========================================================
353
354	! h2so4
355
356	do l = 1,n_lev
357	do ig = 1,n_lon
358	zqi_sa(ig,l) = zqi_sa(ig,l)
359	$ + (F_sed(ig,l+1)*wgt_SA(ig,l+1)
360	$ - F_sed(ig,l)*wgt_SA(ig,l))/m_lay(ig,l)
361	! if (zqi_sa(ig,l) < 0.) THEN
362	! print*,'STOP sedim on epuise tout le H2SO4l present'
363	! print*,'point ',ig,'level ',l
364	! print*,'zqi_sa = ', zqi_sa(ig,l)
365	! STOP
366	! zqi_sa(ig,l) = 0.
367	! end if
368	zqi_sa(ig,l) = max(zqi_sa(ig,l), 0.)
369	pdqsed(ig,l,1) = zqi_sa(ig,l) - pq(ig,l,i_h2so4liq)
370	end do
371	end do
372
373	! h2o
374
375	do l = 1, n_lev
376	do ig=1,n_lon
377	zqi_wv(ig,l) = zqi_wv(ig,l)
378	$ + (F_sed(ig,l+1)*(1. - wgt_SA(ig,l+1))
379	& - F_sed(ig,l)*(1. - wgt_SA(ig,l)))
380	& /m_lay(ig,l)
381	! if (zqi_wv(ig,l) < 0.) THEN
382	! print*,'STOP sedim on epuise tout le H2Ol present'
383	! print*,'point ',ig,'level ',l
384	! print*,'zqi_wv = ', zqi_wv(ig,l)
385	! STOP
386	! zqi_wv(ig,l) = 0.
387	! end if
388	zqi_wv(ig,l) = max(zqi_wv(ig,l), 0.)
389	pdqsed(ig,l,2) = zqi_wv(ig,l) - pq(ig,l,i_h2oliq)
390	end do
391	end do
392
393	c Save output file in 1D model
394	c ============================
395	c IF (n_lon .EQ. 1) THEN
396	c PRINT*,'Save output sedim'
397	c DO l = 1, n_lev
398	c DO ig=1,n_lon
399	c WRITE(77,"(i4,','11(e15.8,','))") l,pdqsed(ig,l),zqi(ig,l),
400	c & (WH2SO4(ig,l)*pq(ig,l,i_h2so4liq)+
401	c & (1.-WH2SO4(ig,l))*pq(ig,l,i_h2oliq)),
402	c & pq(ig,l,i_h2so4liq),pq(ig,l,i_h2oliq)
403	c ENDDO
404	c ENDDO
405	c ENDIF
406
407	RETURN
408	END
409
410	*******************************************************************************
411	REAL FUNCTION VISCOSITY_CO2(temp)
412	c Aurélien Stolzenbach 2015
413	c Calcul de la viscosité dynamique du CO2 80°K -> 300°K
414	c Viscosité dynamique en Pa.s
415	c Source: Johnston & Grilly (1942)
416
417	c température en °K
418	REAL, INTENT(IN) :: temp
419
420	REAL :: denom, numer
421
422	c Calcul de la viscosité dynamique grâce à la formule de Jones (Lennard-Jones (1924))
423
424	numer = 200.**(2.27/4.27)-0.435
425	denom = temp**(2.27/4.27)-0.435
426
427	VISCOSITY_CO2 = (numer/denom)1015.(temp/200.)**(3./2.)
428
429	c convertion de Poises*1e7 -> Pa.s
430	VISCOSITY_CO2 = VISCOSITY_CO2*1.e-8
431
432	END FUNCTION VISCOSITY_CO2
433	*******************************************************************************
434
435

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

Download in other formats:

Original Format