Context Navigation

new_cloud_sedim.F @ 2200

Last change on this file since 2200 was 2200, checked in by flefevre, 5 years ago
correction de bugs dans le supercycling de la chimie
Property svn:executable set to ``*
File size: 15.9 KB

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

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

Download in other formats:

Original Format