Context Navigation

muphys3D.F @ 855

Last change on this file since 855 was 474, checked in by slebonnois, 13 years ago
Update of Titan physics for clouds.
File size: 16.0 KB

Line
1	subroutine muphys(ngrid,
2	& plev,play,zlev,zlay,
3	& tpt,tq,gaz1,gaz2,gaz3,
4	& nmicro,ptimestep,
5	& pmu0,pfract,
6	* Sorties diagnostiques
7	& flxesp_i,
8	& tau_drop,tau_aer,
9	& solesp,prec)
10
11	c
12	c
13	c CETTE NOUVELLE ROUTINE DE MICROPHYSIQUE GERE
14	c LA MICROPHYSIQUE DES AEROSOLS ET CELLE DES NUAGES
15	c EN UN SEUL APPEL...
16	c
17	c TOUS LES TRACEURS AEROSOL+NOYAUX+2*GLACES SONT CONTENUS DANS
18	c LE TABLEAU TQ ET LEUR TENDANCES DANS TDQ. CES TABLEAUX
19	c SONT COMPOSES DE TROIS PARTIES:
20	c
21	c TQ( 1, nmicro/4 pour les aerosols
22	c + nmicro/4+1, 2*nmicro/4 pour les noyaux
23	c + 2nmicro/4+1, 3nmicro/4 pour la glace 1
24	c + 3*nmicro/4+1, nmicro ) pour la glace 2
25	c
26	c pour les aerosols, les noyaux, la glace 1 et la glace 2. la separation
27	c puis la concatenation de fait juste avant l'appel aux routines
28	c dans lesquels la separation est necessaire.
29	c
30	c _______
31	c \| \| ____
32	c EX: --------> QAER() ----> \| M \| \
33	c / \| I \| \
34	c / \| C \| \
35	c TQ() ----------> QGLACE1() ---> \| R \| ------ TDQ()
36	c \ \| O \| /
37	c \ \| P \| /
38	c --------> QGLACE2() ----> \| H \| ----/
39	c \ \| Y \| /
40	c \ \| \| /
41	c -----> QNOYAUX() \| \| _/
42	c \| \|
43	c -------
44	c
45	c
46	c
47	c DANS LA MICROPHYSIQUE, SE TROUVENT IMBRIQUES LES PROCESSUS SUIVANTS
48	c
49	c
50	c
51	c
52	c
53	c NUCLEATION/SEDIMENTATION ---> n_ethane.F / n_methane.F
54	c
55	c MICROPHYSIQUE AEROSOLS ---> brume.F
56	c
57	c SEDIMENTATION DES GOUTTES ---> snuages.F
58	c
59	c
60	c
61	c
62	c
63	c
64	c------------------------------------------------------
65	use dimphy
66	IMPLICIT NONE
67	#include "dimensions.h"
68	#include "microtab.h"
69	#include "varmuphy.h"
70	#include "clesphys.h"
71	#include "itemps.h"
72
73	integer ngrid
74
75	integer iq,nmicro
76
77	common/part/vaer,raer,vrat,draer,dvaer
78	real vaer(nrad),raer(nrad),vrat,
79	& draer(nrad),dvaer(nrad)
80
81	real ptimestep
82
83	real pdpsrf(ngrid)
84
85	c*************************************
86	c declaration des variables internes *
87	c*************************************
88
89	c sources *
90	c------------------*
91
92	REAL plev(ngrid,klev+1)
93	REAL play(ngrid,klev)
94	REAL zlev(ngrid,klev+1)
95	REAL zlay(ngrid,klev)
96	REAL pu(ngrid),pv(ngrid)
97	REAL pmu0(ngrid),pfract(ngrid)
98	REAL tpt(ngrid,klev)
99	REAL tq(ngrid,klev,nmicro),
100	& gaz1(ngrid,klev),
101	& gaz2(ngrid,klev),
102	& gaz3(ngrid,klev)
103
104	c OUTPUT !!!!! c
105	c note : gaz1,...,gazN sont aussi des outputs, ils ont modifié tout au long de muphys.
106	c--------------------c
107	REAL pdq(ngrid,klev,nmicro)
108	REAL flxesp_i(ngrid,klev,3) ! flx esp GLACE
109	REAL solesp(ngrid,klev,3) ! tx prod glace (puit/source)
110	REAL tau_drop(ngrid,klev)
111	REAL tau_aer(ngrid,klev,nrad)
112	REAL prec(ngrid,5)
113
114	c LOCAL
115	c--------------------c
116	real q(ngrid,klev,nmicro)
117	REAL taused(klev,nrad)
118	integer jsup,jinf,h,jalt,ihor,k
119
120	c microphysique *
121	c---------------------*
122	real c(klev,nrad), cni(klev,nrad)
123	real c1i(klev,nrad),c2i(klev,nrad),c3i(klev,nrad)
124	real gazc1(klev),gazc2(klev),gazc3(klev)
125	real ddt
126
127	real vcl,nuc,r,xgsn,xmsn
128	real zz,effg,xlog,rapport
129
130	integer IPREM,i,n,j,l
131	integer ibid
132	save IPREM
133	data IPREM/0/
134
135
136	c real ttq(ngrid,klev,nmicro,2)
137	c real tttq(ngrid,klev,nmicro,2)
138
139
140	c **************************
141	c INITIALISATION DE TABLEAUX
142	c **************************
143	c A NE FAIRE QU'UNE FOIS
144	c **************************
145	c
146	IF (IPREM.eq.0) THEN
147
148	IF (microfi.eq.1) THEN
149	IF (ngrid.ne.jjm+1) THEN
150	print*,"aLeRte :"
151	print*,"microfi en 2D mais ngrid.ne.jjm+1"
152	print*,ngrid,jjm+1
153	stop "je m'arrete..."
154	ENDIF
155	ELSEIF (microfi.eq.2) THEN
156	IF (ngrid.ne.klon) THEN
157	print*,"aLeRte :"
158	print*,"microfi en 3D mais ngrid.ne.klon"
159	print*,ngrid,klon
160	stop "je m'arrete..."
161	ENDIF
162	ENDIF
163
164	c initialisation des constantes de la microphysique :
165	c ----------------------------------------------
166	call inimphycst()
167
168
169	c initialisation des c(z,r), c1i(z,r), c2i(z,r)
170	c ----------------------------------------------
171
172	do i=1,nmicro
173	do n=1,ngrid
174	do j=klev,1,-1
175	q(n,klev+1-j,i)=tq(n,j,i) ! glaces...
176	if(i.le.2nrad) q(n,klev+1-j,i)=tq(n,j,i)tcorrect ! noyaux+aerosol
177	pdq(n,j,i)=0.
178	enddo
179	enddo
180	enddo
181	c ici, les tableaux definissant la structure des aerosols sont
182	c remplis: rf,df(nmicro),r(nmicro),v(nmicro)......
183	call rdf()
184	c ici on recopie la grille dans un common specifique a la microfi...
185	v_e = vaer
186	r_e = raer
187	vrat_e = vrat
188	dr_e = draer
189	dv_e = dvaer
190	c
191	ELSE
192
193	c les tq() doivent etre en nombre d'aerosols / cases
194
195	do j=1,klev ! j de 1 a 119
196	do n=1,ngrid
197	do i=1,nmicro
198	q(n,j,i)=tq(n,klev-j+1,i)
199	pdq(n,j,i)=0.
200	enddo
201	enddo
202	enddo
203
204	ENDIF ! FIN IPREM
205
206
207	c-----------------------------------------------------
208	c !! La premiere fois, on ne passe pas par
209	c !! q--->c et par pg3.F
210	c !! on passe directement au remplissage c-->q
211	IF (IPREM.eq.0) goto 102
212	c-----------------------------------------------------
213
214	c****************************************
215	c *
216	c ADAPTATION GCM > micro *
217	c *
218	c****************************************
219
220
221	c correpondance des couches / sens GCM > microphysique
222	c-----------------------------------------------------
223	c
224	do IHOR=1,NGRID ! GRANDE BOUCLE HORIZONTALE / SEPARATION DES COLONNES
225
226
227	c Ici, on initialise la grille verticale et les
228	c variables communes aux routines de microphysique.
229	c*******************************************************
230	call inimuphy(ihor,plev(ihor,:),play(ihor,:),
231	& zlev(ihor,:),zlay(ihor,:),
232	& tpt(ihor,:))
233
234	c Ici, on scinde les tableaux aerosols et glaces
235	c*******************************************************
236
237	if (clouds.eq.0) then
238	if(nrad .ne. nmicro) then
239	print*,"aLeRte : nrad != nmicro"
240	print*,'nmicro= ',nmicro
241	print*,'nrad= ',nrad
242	stop "je m'arrete..."
243	endif
244	else
245	if(nrad .ne. nmicro/ntype) then
246	print*,"aLeRte : nrad != nmicro/ntype"
247	print*,'nmicro= ',nmicro
248	print*,'ntype=',ntype
249	print*,'nmicro/ntype= ',nmicro/ntype
250	print*,'nrad= ',nrad
251	stop "je m'arrete..."
252	endif
253	endif
254
255	do i=1,nrad
256	do j=1,klev
257	c(j,i) =q(IHOR,j,i )/dzb(j) ! concentration aerosols/m^3
258	if (clouds.eq.1) then
259	cni(j,i)=q(IHOR,j,i+ nrad)/dzb(j) ! concentration noyaux /m^3
260	c1i(j,i)=q(IHOR,j,i+2*nrad)/dzb(j) ! concentration volume glace/m^3
261	c2i(j,i)=q(IHOR,j,i+3*nrad)/dzb(j) ! concentration volume glace /m^3
262	c3i(j,i)=q(IHOR,j,i+4*nrad)/dzb(j) ! concentration volume glace /m^3
263	endif
264	enddo
265	enddo
266	if (clouds.eq.1) then
267	do j=1,klev
268	gazc1(j) =gaz1(IHOR,klev-j+1) ! fraction molaire CH4
269	gazc2(j) =gaz2(IHOR,klev-j+1) ! fraction molaire C2H6
270	gazc3(j) =gaz3(IHOR,klev-j+1) ! fraction molaire C2H2
271	enddo
272	endif
273
274	c****************************************
275	c
276	c FIN DE LA PREPARATION:
277	c
278	c
279	c ON APPELLE LES MODELES MICROPHYSIQUES
280	c
281	c - brume (coagulation + sedimentation)
282	c - nuages (nucleation + condensation)
283	c - sedimentation nuages
284	c
285	c
286	c****************************************
287
288
289	do j=1,klev
290	solesp(ihor,klev+1-j,:) = 0.
291	do i=1,nrad
292	tau_aer(ihor,klev+1-j,i)=0.
293	enddo
294	enddo
295
296	ddt=ptimestep
297
298	* concerne les aerosols (tableau c):
299	c call begintime(tt0)
300	call brume(ngrid,c,ddt,klev,nrad,taused,ihor,
301	& pmu0(ihor),pfract(ihor),
302	& prec)
303	c call endtime(tt0,tt1)
304	c tthaze=tthaze+tt1
305
306	* concerne aerosols + gouttes (tableaux c,cn,c1,c2):
307
308	do j=1,klev
309	do i=1,nrad
310	tau_aer(ihor,klev+1-j,i)=tau_aer(ihor,klev+1-j,i)
311	& +taused(j,i)
312	enddo
313	enddo
314
315	IF (clouds.eq.1) THEN
316
317	c do j=1,klev
318	c do i=1,nrad
319	c ttq(ihor,klev-j+1,i,1) = c(j,i)*dzb(j)
320	c ttq(ihor,klev-j+1,i+nrad,1) = cni(j,i)*dzb(j)
321	c ttq(ihor,klev-j+1,i+2nrad,1) = c1i(j,i)dzb(j)
322	c ttq(ihor,klev-j+1,i+3nrad,1) = c2i(j,i)dzb(j)
323	c ttq(ihor,klev-j+1,i+4nrad,1) = c3i(j,i)dzb(j)
324	c enddo
325	c enddo
326
327	c call begintime(tt0)
328	call cnuages(c,c1i,c2i,c3i,cni,
329	& gazc1,gazc2,gazc3,ddt)
330	c call endtime(tt0,tt1)
331	c ttcclds=ttcclds+tt1
332
333	c verification des valeurs de c,cni,c1i,c2i et c3i.
334	c Lorsque l'on vide completement une case, on peut avoir des chiffres negatifs :s
335	do j=1,klev
336	do i=1,nrad
337	c(j,i)= MAX(c(j,i),0.)
338	cni(j,i)= MAX(cni(j,i),0.)
339	c1i(j,i)= MAX(c1i(j,i),0.)
340	c2i(j,i)= MAX(c2i(j,i),0.)
341	c3i(j,i)= MAX(c3i(j,i),0.)
342	c ttq(ihor,klev-j+1,i,2) = c(j,i)*dzb(j)
343	c ttq(ihor,klev-j+1,i+nrad,2) = cni(j,i)*dzb(j)
344	c ttq(ihor,klev-j+1,i+2nrad,2) = c1i(j,i)dzb(j)
345	c ttq(ihor,klev-j+1,i+3nrad,2) = c2i(j,i)dzb(j)
346	c ttq(ihor,klev-j+1,i+4nrad,2) = c3i(j,i)dzb(j)
347	enddo
348	enddo
349
350
351	do j=1,klev
352	do i=1,nrad
353	! solesp en m3/m3 pour passer en m3/m2 il faut faire :
354	! (c1i(j,i)dzb(j) -q(IHOR,j,i+2nrad))
355	solesp(ihor,klev+1-j,1)=solesp(ihor,klev+1-j,1) +
356	& (c1i(j,i)-q(IHOR,j,i+2*nrad)/dzb(j))
357	solesp(ihor,klev+1-j,2)=solesp(ihor,klev+1-j,2) +
358	& (c2i(j,i)-q(IHOR,j,i+3*nrad)/dzb(j))
359	solesp(ihor,klev+1-j,3)=solesp(ihor,klev+1-j,3) +
360	& (c3i(j,i)-q(IHOR,j,i+4*nrad)/dzb(j))
361	enddo
362	enddo
363
364	* concerne les gouttes (tableaux cn,c1,c2):
365
366	c do j=1,klev
367	c do i=1,nrad
368	c tttq(ihor,klev-j+1,i,1) = c(j,i)*dzb(j)
369	c tttq(ihor,klev-j+1,i+nrad,1) = cni(j,i)*dzb(j)
370	c tttq(ihor,klev-j+1,i+2nrad,1) = c1i(j,i)dzb(j)
371	c tttq(ihor,klev-j+1,i+3nrad,1) = c2i(j,i)dzb(j)
372	c tttq(ihor,klev-j+1,i+4nrad,1) = c3i(j,i)dzb(j)
373	c enddo
374	c enddo
375
376	c call begintime(tt0)
377	call snuages(ngrid,cni,c1i,c2i,c3i,c,ddt,
378	& klev,nrad,ihor,
379	& flxesp_i,tau_drop,prec)
380	c call endtime(tt0,tt1)
381	c ttsclds=ttsclds+tt1
382
383	c do j=1,klev
384	c do i=1,nrad
385	c tttq(ihor,klev-j+1,i,2) = c(j,i)*dzb(j)
386	c tttq(ihor,klev-j+1,i+nrad,2) = cni(j,i)*dzb(j)
387	c tttq(ihor,klev-j+1,i+2nrad,2) = c1i(j,i)dzb(j)
388	c tttq(ihor,klev-j+1,i+3nrad,2) = c2i(j,i)dzb(j)
389	c tttq(ihor,klev-j+1,i+4nrad,2) = c3i(j,i)dzb(j)
390	c enddo
391	c enddo
392	ENDIF ! flag nuages :)
393
394
395	* on recompose le tableau de traceurs ici.
396	*--------------------------------------------------
397
398	do i=1,nrad
399	do j=1,klev
400	q(IHOR,j,i) = c(j,i)*dzb(j) ! nombre aerosols /m^2
401	if (clouds.eq.1) then
402	q(IHOR,j,i+ nrad) = cni(j,i)*dzb(j) ! nombre noyaux /m^2
403	q(IHOR,j,i+2nrad) = c1i(j,i)dzb(j) ! concentration volume glace/m^2
404	q(IHOR,j,i+3nrad) = c2i(j,i)dzb(j) ! concentration volume glace/m^2
405	q(IHOR,j,i+4nrad) = c3i(j,i)dzb(j) ! concentration volume glace/m^2
406	endif
407	enddo
408	enddo
409	if (clouds.eq.1) then
410	do j=1,klev
411	gaz1(IHOR,klev-j+1) = gazc1(j) ! fraction molaire
412	gaz2(IHOR,klev-j+1) = gazc2(j) ! fraction molaire
413	gaz3(IHOR,klev-j+1) = gazc3(j) ! fraction molaire
414	enddo
415	endif
416
417	ENDDO ! Fin de la boucle IHOR
418
419	102 CONTINUE ! la premiere fois, c'est une boucle vide!
420
421
422	c***************************************************************
423	c FIN: on renvoie les nouvelles valeurs q(t+dt)=q(t) + dq(t)
424	c
425	c Pour les aerosols, les noyaux, les glaces et les vapeurs modifiees...
426	c
427	c***************************************************************
428
429	do n=1,ngrid
430	do i=1,nmicro
431	do j=1,klev ! j de 1 a 54
432	tq(n,j,i) = q(n,klev+1-j,i)
433	enddo
434	enddo
435	enddo
436
437
438	c do j=1,15
439	cc CH4 -- cnuages
440	c write(210,'(I4,3(ES24.16,1X))') j,
441	c & sum(ttq(40,j,2nrad+1:3nrad,1)),
442	c & sum(ttq(40,j,2nrad+1:3nrad,2)),
443	c & sum(ttq(40,j,2nrad+1:3nrad,2))-
444	c & sum(ttq(40,j,2nrad+1:3nrad,1))
445	cc C2H6 -- cnuages
446	c write(211,'(I4,3(ES24.16,1X))') j,
447	c & sum(ttq(40,j,3nrad+1:4nrad,1)),
448	c & sum(ttq(40,j,3nrad+1:4nrad,2)),
449	c & sum(ttq(40,j,3nrad+1:4nrad,2))-
450	c & sum(ttq(40,j,3nrad+1:4nrad,1))
451	cc C2H2 -- cnuages
452	c write(212,'(I4,3(ES24.16,1X))') j,
453	c & sum(ttq(40,j,4nrad+1:5nrad,1)),
454	c & sum(ttq(40,j,4nrad+1:5nrad,2)),
455	c & sum(ttq(40,j,4nrad+1:5nrad,2))-
456	c & sum(ttq(40,j,4nrad+1:5nrad,1))
457	c
458	cc CH4 -- snuages
459	c write(310,'(I4,3(ES24.16,1X))') j,
460	c & sum(tttq(40,j,2nrad+1:3nrad,1)),
461	c & sum(tttq(40,j,2nrad+1:3nrad,2)),
462	c & sum(tttq(40,j,2nrad+1:3nrad,2))-
463	c & sum(tttq(40,j,2nrad+1:3nrad,1))
464	cc C2H6 -- snuages
465	c write(311,'(I4,3(ES24.16,1X))') j,
466	c & sum(tttq(40,j,3nrad+1:4nrad,1)),
467	c & sum(tttq(40,j,3nrad+1:4nrad,2)),
468	c & sum(tttq(40,j,3nrad+1:4nrad,2))-
469	c & sum(tttq(40,j,3nrad+1:4nrad,1))
470	cc C2H2 -- snuages
471	c write(312,'(I4,3(ES24.16,1X))') j,
472	c & sum(tttq(40,j,4nrad+1:5nrad,1)),
473	c & sum(tttq(40,j,4nrad+1:5nrad,2)),
474	c & sum(tttq(40,j,4nrad+1:5nrad,2))-
475	c & sum(tttq(40,j,4nrad+1:5nrad,1))
476	c enddo
477	c write(210,*) "NEWLINE"
478	c write(211,*) "NEWLINE"
479	c write(212,*) "NEWLINE"
480	c write(310,*) "NEWLINE"
481	c write(311,*) "NEWLINE"
482	c write(312,*) "NEWLINE"
483
484
485	c do j=1,20
486	c write(410,'(I4,3(ES24.16,1X))') j,
487	c & flxesp_i(40,j,1),flxesp_i(40,j,2),flxesp_i(40,j,3)
488	c write(510,'(I4,3(ES24.16,1X))') j,
489	c & solesp(40,j,1),solesp(40,j,2),solesp(40,j,3)
490	c enddo
491	c write(410,*) "NEWLINE"
492	c write(510,*) "NEWLINE"
493
494	IPREM=1 ! LA PROCHAINE FOIS NE SERA PLUS LA 1ERE
495
496
497	16 return
498
499	end
500
501
502	c---------------------------------------------------------------------

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

Download in other formats:

Original Format