Context Navigation

new_cloud_sedim.F @ 1687

Last change on this file since 1687 was 1687, checked in by slebonnois, 8 years ago
SL: bugs corrections outputs/clouds_AStol/upper_atmosphere/start2archive
Property svn:executable set to ``*
File size: 16.1 KB

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

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

Download in other formats:

Original Format