Context Navigation

new_cloud_sedim.F @ 1530

Last change on this file since 1530 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

Property svn:executable set to *

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

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

Download in other formats:

Original Format