Context Navigation

new_cloud_sedim.F @ 1543

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

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

Property svn:executable set to *

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

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

Download in other formats:

Original Format