Context Navigation

new_cloud_sedim.F @ 3094

Last change on this file since 3094 was 2464, checked in by slebonnois, 4 years ago
SL: major update related to extension to 250 km. New EUV heating as used in Martian GCM, debug of a few routines, adding He, clean up, updating mmean regularly, modification of the order of processes, add the possibility to use Hedin composition in rad transfer
Property svn:executable set to ``*
File size: 15.0 KB

Line
1	subroutine new_cloud_sedim(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	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
36	c local:
37	c ------
38	integer imode
39	integer ig
40	integer iq
41	integer l
42
43	real zlev(n_lon,n_lev+1) ! altitude at layer boundaries
44	real zlay(n_lon,n_lev) ! altitude at the midlle layer
45	real zqi_wv(n_lon,n_lev) ! to locally store H2O tracer
46	real zqi_sa(n_lon,n_lev) ! to locally store H2SO4 tracer
47	real m_lay (n_lon,n_lev) ! Layer Pressure over gravity (Dp/g == kg.m-2)
48	real wq(n_lon,n_lev+1) ! displaced tracer mass (kg.m-2)
49
50	c Physical constant
51	c ~~~~~~~~~~~~~~~~~
52	c Gas molecular viscosity (N.s.m-2)
53	c real,parameter :: visc=1.e-5 ! CO2
54	REAL :: VISCOSITY_CO2
55	c Effective gas molecular radius (m)
56	real,parameter :: molrad=2.2e-10 ! CO2
57
58	c Ratio radius shell model du mode 3
59	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
60	c Ce ratio correspond aux mesures effectuées par J. Cimino (1982), Icarus
61	c Fixer ce parametre a 0 revient a une gouttelette pure en liquide acide sulfurique
62	c ATTENTION ! DOIT ETRE COHERENT AVEC new_cloud_venus !
63	REAL, PARAMETER :: qrad = 0.97
64	REAL :: qmass
65	c masse volumique du coeur (kg.m-3)
66	c ATTENTION ! DOIT ETRE COHERENT AVEC new_cloud_venus !
67	REAL, PARAMETER :: rho_core = 2500.0
68
69	REAL, DIMENSION(n_lon,n_lev+1) :: wgt_SA ! Fraction of H2SO4 in droplet local
70
71	c Stokes speed and sedimentation flux variable
72	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
73
74	REAL :: A1,A2,A3,A4, ! coeff du DL du Flux de sedimentation
75	+ D_stokes, ! coeff de la vitesse de Stokes
76	+ Rp_DL, ! "Point" du DL
77	+ l_mean, ! libre parcours moyen (m)
78	+ a,b_exp,c ! coeff du calcul du Flux de sedimentation
79	REAL, DIMENSION(n_lon,n_lev+1) ::
80	+ F_sed ! Flux de sedimentation (kg.m-2.s-1 puis en output kg.m-2)
81
82	REAL :: R_mode0 ! Rayon mode 0 (m), rayon le plus frequent
83
84	! PRINT*,'RHO_DROPLET new_cloud_sedim.F'
85	! PRINT*,'rho_droplet',rho_droplet(16,21)
86	! PRINT*,'T',pt(16,21),'WSA',WH2SO4(16,21)
87
88	c-----------------------------------------------------------------------
89	c 1. Initialization
90	c -----------------
91
92	! update water vapour and sulfuric acid mixing ratios
93
94	zqi_wv(:,:) = pq(:,:,i_h2oliq) + d_tr_chem(:,:,i_h2oliq)*ptimestep
95	zqi_sa(:,:) = pq(:,:,i_h2so4liq)
96	$ + d_tr_chem(:,:,i_h2so4liq)*ptimestep
97
98	wgt_SA(:,1:n_lev) = 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
351	c Passage du Flux au Flux pour un pas de temps (== kg.m-2)
352
353	F_sed(:,:) = F_sed(:,:)*ptimestep
354
355	!=========================================================
356	! compute tendency due to sedimentation
357	!=========================================================
358
359	! h2so4
360
361	do l = 1,n_lev
362	do ig = 1,n_lon
363	zqi_sa(ig,l) = zqi_sa(ig,l)
364	$ + (F_sed(ig,l+1)*wgt_SA(ig,l+1)
365	$ - F_sed(ig,l)*wgt_SA(ig,l))/m_lay(ig,l)
366	! if (zqi_sa(ig,l) < 0.) THEN
367	! print*,'STOP sedim on epuise tout le H2SO4l present'
368	! print*,'point ',ig,'level ',l
369	! print*,'zqi_sa = ', zqi_sa(ig,l)
370	! STOP
371	! zqi_sa(ig,l) = 0.
372	! end if
373	zqi_sa(ig,l) = max(zqi_sa(ig,l), 0.)
374	pdqsed(ig,l,1) = zqi_sa(ig,l) - pq(ig,l,i_h2so4liq)
375	end do
376	end do
377
378	! h2o
379
380	do l = 1, n_lev
381	do ig=1,n_lon
382	zqi_wv(ig,l) = zqi_wv(ig,l)
383	$ + (F_sed(ig,l+1)*(1. - wgt_SA(ig,l+1))
384	& - F_sed(ig,l)*(1. - wgt_SA(ig,l)))
385	& /m_lay(ig,l)
386	! if (zqi_wv(ig,l) < 0.) THEN
387	! print*,'STOP sedim on epuise tout le H2Ol present'
388	! print*,'point ',ig,'level ',l
389	! print*,'zqi_wv = ', zqi_wv(ig,l)
390	! STOP
391	! zqi_wv(ig,l) = 0.
392	! end if
393	zqi_wv(ig,l) = max(zqi_wv(ig,l), 0.)
394	pdqsed(ig,l,2) = zqi_wv(ig,l) - pq(ig,l,i_h2oliq)
395	end do
396	end do
397
398	c Save output file in 1D model
399	c ============================
400	c IF (n_lon .EQ. 1) THEN
401	c PRINT*,'Save output sedim'
402	c DO l = 1, n_lev
403	c DO ig=1,n_lon
404	c WRITE(77,"(i4,','11(e15.8,','))") l,pdqsed(ig,l),zqi(ig,l),
405	c & (WH2SO4(ig,l)*pq(ig,l,i_h2so4liq)+
406	c & (1.-WH2SO4(ig,l))*pq(ig,l,i_h2oliq)),
407	c & pq(ig,l,i_h2so4liq),pq(ig,l,i_h2oliq)
408	c ENDDO
409	c ENDDO
410	c ENDIF
411
412	RETURN
413	END
414
415	*******************************************************************************
416	REAL FUNCTION VISCOSITY_CO2(temp)
417	c Aurélien Stolzenbach 2015
418	c Calcul de la viscosité dynamique du CO2 80°K -> 300°K
419	c Viscosité dynamique en Pa.s
420	c Source: Johnston & Grilly (1942)
421
422	c température en °K
423	REAL, INTENT(IN) :: temp
424
425	REAL :: denom, numer
426
427	c Calcul de la viscosité dynamique grâce à la formule de Jones (Lennard-Jones (1924))
428
429	numer = 200.**(2.27/4.27)-0.435
430	denom = temp**(2.27/4.27)-0.435
431
432	VISCOSITY_CO2 = (numer/denom)1015.(temp/200.)**(3./2.)
433
434	c convertion de Poises*1e7 -> Pa.s
435	VISCOSITY_CO2 = VISCOSITY_CO2*1.e-8
436
437	END FUNCTION VISCOSITY_CO2
438	*******************************************************************************
439
440

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

Download in other formats:

Original Format