Context Navigation

muphys3D.F @ 201

Last change on this file since 201 was 176, checked in by slebonnois, 14 years ago
SLebonnois: petite correction pour pouvoir tourner avec clouds=0 sans etre embetes par ntype
File size: 15.9 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+1,3) ! tx prod glace (puit/source) ! ind klev+1 = surface ?
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	real ttq(ngrid,klev,nmicro,2)
137	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(ihor,klev+1-j,1)=solesp(ihor,klev+1-j,1) +
354	& (c1i(j,i)-q(IHOR,j,i+2*nrad)/dzb(j))
355	solesp(ihor,klev+1-j,2)=solesp(ihor,klev+1-j,2) +
356	& (c2i(j,i)-q(IHOR,j,i+3*nrad)/dzb(j))
357	solesp(ihor,klev+1-j,3)=solesp(ihor,klev+1-j,3) +
358	& (c3i(j,i)-q(IHOR,j,i+4*nrad)/dzb(j))
359	enddo
360	enddo
361
362	* concerne les gouttes (tableaux cn,c1,c2):
363
364	c do j=1,klev
365	c do i=1,nrad
366	c tttq(ihor,klev-j+1,i,1) = c(j,i)*dzb(j)
367	c tttq(ihor,klev-j+1,i+nrad,1) = cni(j,i)*dzb(j)
368	c tttq(ihor,klev-j+1,i+2nrad,1) = c1i(j,i)dzb(j)
369	c tttq(ihor,klev-j+1,i+3nrad,1) = c2i(j,i)dzb(j)
370	c tttq(ihor,klev-j+1,i+4nrad,1) = c3i(j,i)dzb(j)
371	c enddo
372	c enddo
373
374	c call begintime(tt0)
375	call snuages(ngrid,cni,c1i,c2i,c3i,c,ddt,
376	& klev,nrad,ihor,
377	& flxesp_i,tau_drop,prec)
378	c call endtime(tt0,tt1)
379	c ttsclds=ttsclds+tt1
380
381	c do j=1,klev
382	c do i=1,nrad
383	c tttq(ihor,klev-j+1,i,2) = c(j,i)*dzb(j)
384	c tttq(ihor,klev-j+1,i+nrad,2) = cni(j,i)*dzb(j)
385	c tttq(ihor,klev-j+1,i+2nrad,2) = c1i(j,i)dzb(j)
386	c tttq(ihor,klev-j+1,i+3nrad,2) = c2i(j,i)dzb(j)
387	c tttq(ihor,klev-j+1,i+4nrad,2) = c3i(j,i)dzb(j)
388	c enddo
389	c enddo
390	ENDIF ! flag nuages :)
391
392
393	* on recompose le tableau de traceurs ici.
394	*--------------------------------------------------
395
396	do i=1,nrad
397	do j=1,klev
398	q(IHOR,j,i) = c(j,i)*dzb(j) ! nombre aerosols /m^2
399	if (clouds.eq.1) then
400	q(IHOR,j,i+ nrad) = cni(j,i)*dzb(j) ! nombre noyaux /m^2
401	q(IHOR,j,i+2nrad) = c1i(j,i)dzb(j) ! concentration volume glace/m^2
402	q(IHOR,j,i+3nrad) = c2i(j,i)dzb(j) ! concentration volume glace/m^2
403	q(IHOR,j,i+4nrad) = c3i(j,i)dzb(j) ! concentration volume glace/m^2
404	endif
405	enddo
406	enddo
407	if (clouds.eq.1) then
408	do j=1,klev
409	gaz1(IHOR,klev-j+1) = gazc1(j) ! fraction molaire
410	gaz2(IHOR,klev-j+1) = gazc2(j) ! fraction molaire
411	gaz3(IHOR,klev-j+1) = gazc3(j) ! fraction molaire
412	enddo
413	endif
414
415	ENDDO ! Fin de la boucle IHOR
416
417	102 CONTINUE ! la premiere fois, c'est une boucle vide!
418
419
420	c***************************************************************
421	c FIN: on renvoie les nouvelles valeurs q(t+dt)=q(t) + dq(t)
422	c
423	c Pour les aerosols, les noyaux, les glaces et les vapeurs modifiees...
424	c
425	c***************************************************************
426
427	do n=1,ngrid
428	do i=1,nmicro
429	do j=1,klev ! j de 1 a 54
430	tq(n,j,i) = q(n,klev+1-j,i)
431	enddo
432	enddo
433	enddo
434
435
436	c do j=1,15
437	cc CH4 -- cnuages
438	c write(210,'(I4,3(ES24.16,1X))') j,
439	c & sum(ttq(40,j,2nrad+1:3nrad,1)),
440	c & sum(ttq(40,j,2nrad+1:3nrad,2)),
441	c & sum(ttq(40,j,2nrad+1:3nrad,2))-
442	c & sum(ttq(40,j,2nrad+1:3nrad,1))
443	cc C2H6 -- cnuages
444	c write(211,'(I4,3(ES24.16,1X))') j,
445	c & sum(ttq(40,j,3nrad+1:4nrad,1)),
446	c & sum(ttq(40,j,3nrad+1:4nrad,2)),
447	c & sum(ttq(40,j,3nrad+1:4nrad,2))-
448	c & sum(ttq(40,j,3nrad+1:4nrad,1))
449	cc C2H2 -- cnuages
450	c write(212,'(I4,3(ES24.16,1X))') j,
451	c & sum(ttq(40,j,4nrad+1:5nrad,1)),
452	c & sum(ttq(40,j,4nrad+1:5nrad,2)),
453	c & sum(ttq(40,j,4nrad+1:5nrad,2))-
454	c & sum(ttq(40,j,4nrad+1:5nrad,1))
455	c
456	cc CH4 -- snuages
457	c write(310,'(I4,3(ES24.16,1X))') j,
458	c & sum(tttq(40,j,2nrad+1:3nrad,1)),
459	c & sum(tttq(40,j,2nrad+1:3nrad,2)),
460	c & sum(tttq(40,j,2nrad+1:3nrad,2))-
461	c & sum(tttq(40,j,2nrad+1:3nrad,1))
462	cc C2H6 -- snuages
463	c write(311,'(I4,3(ES24.16,1X))') j,
464	c & sum(tttq(40,j,3nrad+1:4nrad,1)),
465	c & sum(tttq(40,j,3nrad+1:4nrad,2)),
466	c & sum(tttq(40,j,3nrad+1:4nrad,2))-
467	c & sum(tttq(40,j,3nrad+1:4nrad,1))
468	cc C2H2 -- snuages
469	c write(312,'(I4,3(ES24.16,1X))') j,
470	c & sum(tttq(40,j,4nrad+1:5nrad,1)),
471	c & sum(tttq(40,j,4nrad+1:5nrad,2)),
472	c & sum(tttq(40,j,4nrad+1:5nrad,2))-
473	c & sum(tttq(40,j,4nrad+1:5nrad,1))
474	c enddo
475	c write(210,*) "NEWLINE"
476	c write(211,*) "NEWLINE"
477	c write(212,*) "NEWLINE"
478	c write(310,*) "NEWLINE"
479	c write(311,*) "NEWLINE"
480	c write(312,*) "NEWLINE"
481
482
483	c do j=1,20
484	c write(410,'(I4,3(ES24.16,1X))') j,
485	c & flxesp_i(40,j,1),flxesp_i(40,j,2),flxesp_i(40,j,3)
486	c write(510,'(I4,3(ES24.16,1X))') j,
487	c & solesp(40,j,1),solesp(40,j,2),solesp(40,j,3)
488	c enddo
489	c write(410,*) "NEWLINE"
490	c write(510,*) "NEWLINE"
491
492	IPREM=1 ! LA PROCHAINE FOIS NE SERA PLUS LA 1ERE
493
494
495	16 return
496
497	end
498
499
500	c---------------------------------------------------------------------

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

Download in other formats:

Original Format