Context Navigation

muphys3D.F @ 3303

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

Venus and Titan GCMs:

Adaptation wrt previous changes for Titan and Venus where
longitude and latitude arrays (in phycommon/geometry_mod) were overwritten
with values from startphy.nc files, where values are given in degrees.
For the sake of homegeneity with other physics package, revert to "default"
behaviour: longitude/latitude are in radians and longitude_deg/latitude_deg
are in degrees.
Also added checking in phyetat0 that the longitude/latitude read in the
restartphy.nc files match the ones provided by the dynamics.

File size: 17.3 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	c use radcommon_h, only : volume,rayon,vrat,drayon,dvolume
67	USE geometry_mod, only: latitude ! in radians
68
69	IMPLICIT NONE
70	#include "dimensions.h"
71	#include "microtab.h"
72	#include "varmuphy.h"
73	#include "clesphys.h"
74	#include "itemps.h"
75
76	integer ngrid
77
78	integer iq,nmicro
79	real ptimestep
80	real pdpsrf(ngrid)
81
82	c a la place de radcommon_h:
83	common/part/vaer,raer,vrat,draer,dvaer
84	real vaer(nrad),raer(nrad),vrat,
85	& draer(nrad),dvaer(nrad)
86
87	c*************************************
88	c declaration des variables internes *
89	c*************************************
90
91	c sources *
92	c------------------*
93
94	REAL plev(ngrid,klev+1)
95	REAL play(ngrid,klev)
96	REAL zlev(ngrid,klev+1)
97	REAL zlay(ngrid,klev)
98	REAL pu(ngrid),pv(ngrid)
99	REAL pmu0(ngrid),pfract(ngrid)
100	REAL tpt(ngrid,klev)
101	REAL tq(ngrid,klev,nmicro),
102	& gaz1(ngrid,klev),
103	& gaz2(ngrid,klev),
104	& gaz3(ngrid,klev)
105
106	c OUTPUT !!!!! c
107	c note : gaz1,...,gazN sont aussi des outputs, ils ont modifié tout au long de muphys.
108	c--------------------c
109	REAL pdq(ngrid,klev,nmicro)
110	REAL flxesp_i(ngrid,klev,3) ! flx esp GLACE
111	REAL solesp(ngrid,klev,3) ! tx prod glace (puit/source)
112	REAL tau_drop(ngrid,klev)
113	REAL tau_aer(ngrid,klev,nrad)
114	REAL prec(ngrid,5)
115
116	c LOCAL
117	c--------------------c
118	real q(ngrid,klev,nmicro)
119	REAL taused(klev,nrad)
120	integer jsup,jinf,h,jalt,ihor,k,im1
121
122	c microphysique *
123	c---------------------*
124	real c(klev,nrad), cni(klev,nrad)
125	real c1i(klev,nrad),c2i(klev,nrad),c3i(klev,nrad)
126	real gazc1(klev),gazc2(klev),gazc3(klev)
127	real ddt
128
129	real vcl,nuc,r,xgsn,xmsn
130	real zz,effg,xlog,rapport
131
132	integer IPREM,i,n,j,l
133	integer ibid
134	save IPREM
135	data IPREM/0/
136
137
138	c real ttq(ngrid,klev,nmicro,2)
139	c real tttq(ngrid,klev,nmicro,2)
140
141
142	c **************************
143	c INITIALISATION DE TABLEAUX
144	c **************************
145	c A NE FAIRE QU'UNE FOIS
146	c **************************
147	c
148	IF (IPREM.eq.0) THEN
149
150	IF (ngrid.ne.klon) THEN
151	print*,"aLeRte :"
152	print*,"microfi, mais ngrid.ne.klon"
153	print*,ngrid,klon
154	stop "je m'arrete... (muphys3D)"
155	ENDIF
156
157	c initialisation des constantes de la microphysique :
158	c ----------------------------------------------
159	call inimphycst()
160
161
162	c initialisation des c(z,r), c1i(z,r), c2i(z,r)
163	c ----------------------------------------------
164
165	do i=1,nmicro
166	do n=1,ngrid
167	do j=klev,1,-1
168	q(n,klev+1-j,i)=tq(n,j,i) ! glaces...
169	if(i.le.2nrad) q(n,klev+1-j,i)=tq(n,j,i)tcorrect ! noyaux+aerosol
170	pdq(n,j,i)=0.
171	enddo
172	enddo
173	enddo
174	c ici, les tableaux definissant la structure des aerosols sont
175	c remplis: rf,df(nmicro),r(nmicro),v(nmicro)......
176	call rdf()
177	c ici on recopie la grille dans un common specifique a la microfi...
178	c v_e = volume
179	c r_e = rayon
180	c vrat_e = vrat
181	c dr_e = drayon
182	c dv_e = dvolume
183	v_e = vaer
184	r_e = raer
185	vrat_e = vrat
186	dr_e = draer
187	dv_e = dvaer
188	c
189	ELSE
190
191	c les tq() doivent etre en nombre d'aerosols / cases
192
193	do j=1,klev ! j de 1 a 119
194	do n=1,ngrid
195	do i=1,nmicro
196	q(n,j,i)=tq(n,klev-j+1,i)
197	pdq(n,j,i)=0.
198	enddo
199	enddo
200	enddo
201
202	ENDIF ! FIN IPREM
203
204
205	c-----------------------------------------------------
206	c !! La premiere fois, on ne passe pas par
207	c !! q--->c et par pg3.F
208	c !! on passe directement au remplissage c-->q
209	IF (IPREM.eq.0) goto 102
210	c-----------------------------------------------------
211
212	c****************************************
213	c *
214	c ADAPTATION GCM > micro *
215	c *
216	c****************************************
217
218
219	c correpondance des couches / sens GCM > microphysique
220	c-----------------------------------------------------
221
222	c***************************************************************
223	do IHOR=1,NGRID ! GRANDE BOUCLE HORIZONTALE / SEPARATION DES COLONNES
224
225	if (IHOR.eq.1) then
226	im1=1
227	else
228	im1=IHOR-1
229	endif
230
231	c***************************************************************
232	c On refait les calculs si on est au premier point
233	c OU si on change de latitude
234	c OU si on calcule la microfi en 3D
235	c***************************************************************
236	if((IHOR.eq.1)
237	& .or.(latitude(IHOR).ne.latitude(im1))
238	& .or.(microfi.eq.2)) then
239	c***************************************************************
240
241	c Ici, on initialise la grille verticale et les
242	c variables communes aux routines de microphysique.
243	c*******************************************************
244	call inimuphy(ihor,plev(ihor,:),play(ihor,:),
245	& zlev(ihor,:),zlay(ihor,:),
246	& tpt(ihor,:))
247
248	c Ici, on scinde les tableaux aerosols et glaces
249	c*******************************************************
250
251	if (clouds.eq.0) then
252	if(nrad .ne. nmicro) then
253	print*,"aLeRte : nrad != nmicro"
254	print*,'nmicro= ',nmicro
255	print*,'nrad= ',nrad
256	stop "je m'arrete..."
257	endif
258	else
259	if(nrad .ne. nmicro/ntype) then
260	print*,"aLeRte : nrad != nmicro/ntype"
261	print*,'nmicro= ',nmicro
262	print*,'ntype=',ntype
263	print*,'nmicro/ntype= ',nmicro/ntype
264	print*,'nrad= ',nrad
265	stop "je m'arrete..."
266	endif
267	endif
268
269	do i=1,nrad
270	do j=1,klev
271	c(j,i) =q(IHOR,j,i )/dzb(j) ! concentration aerosols/m^3
272	if (clouds.eq.1) then
273	cni(j,i)=q(IHOR,j,i+ nrad)/dzb(j) ! concentration noyaux /m^3
274	c1i(j,i)=q(IHOR,j,i+2*nrad)/dzb(j) ! concentration volume glace/m^3
275	c2i(j,i)=q(IHOR,j,i+3*nrad)/dzb(j) ! concentration volume glace /m^3
276	c3i(j,i)=q(IHOR,j,i+4*nrad)/dzb(j) ! concentration volume glace /m^3
277	endif
278	enddo
279	enddo
280	if (clouds.eq.1) then
281	do j=1,klev
282	gazc1(j) =gaz1(IHOR,klev-j+1) ! fraction molaire CH4
283	gazc2(j) =gaz2(IHOR,klev-j+1) ! fraction molaire C2H6
284	gazc3(j) =gaz3(IHOR,klev-j+1) ! fraction molaire C2H2
285	enddo
286	endif
287
288	c****************************************
289	c
290	c FIN DE LA PREPARATION:
291	c
292	c
293	c ON APPELLE LES MODELES MICROPHYSIQUES
294	c
295	c - brume (coagulation + sedimentation)
296	c - nuages (nucleation + condensation)
297	c - sedimentation nuages
298	c
299	c
300	c****************************************
301
302
303	do j=1,klev
304	solesp(ihor,klev+1-j,:) = 0.
305	do i=1,nrad
306	tau_aer(ihor,klev+1-j,i)=0.
307	enddo
308	enddo
309
310	ddt=ptimestep
311
312	* concerne les aerosols (tableau c):
313	c call begintime(tt0)
314	call brume(ngrid,c,ddt,klev,nrad,taused,ihor,
315	& pmu0(ihor),pfract(ihor),
316	& prec)
317	c call endtime(tt0,tt1)
318	c tthaze=tthaze+tt1
319
320	* concerne aerosols + gouttes (tableaux c,cn,c1,c2):
321
322	do j=1,klev
323	do i=1,nrad
324	tau_aer(ihor,klev+1-j,i)=tau_aer(ihor,klev+1-j,i)
325	& +taused(j,i)
326	enddo
327	enddo
328
329	IF (clouds.eq.1) THEN
330
331	c do j=1,klev
332	c do i=1,nrad
333	c ttq(ihor,klev-j+1,i,1) = c(j,i)*dzb(j)
334	c ttq(ihor,klev-j+1,i+nrad,1) = cni(j,i)*dzb(j)
335	c ttq(ihor,klev-j+1,i+2nrad,1) = c1i(j,i)dzb(j)
336	c ttq(ihor,klev-j+1,i+3nrad,1) = c2i(j,i)dzb(j)
337	c ttq(ihor,klev-j+1,i+4nrad,1) = c3i(j,i)dzb(j)
338	c enddo
339	c enddo
340
341	c call begintime(tt0)
342	call cnuages(c,c1i,c2i,c3i,cni,
343	& gazc1,gazc2,gazc3,ddt)
344	c call endtime(tt0,tt1)
345	c ttcclds=ttcclds+tt1
346
347	c verification des valeurs de c,cni,c1i,c2i et c3i.
348	c Lorsque l'on vide completement une case, on peut avoir des chiffres negatifs :s
349	do j=1,klev
350	do i=1,nrad
351	c(j,i)= MAX(c(j,i),0.)
352	cni(j,i)= MAX(cni(j,i),0.)
353	c1i(j,i)= MAX(c1i(j,i),0.)
354	c2i(j,i)= MAX(c2i(j,i),0.)
355	c3i(j,i)= MAX(c3i(j,i),0.)
356	c ttq(ihor,klev-j+1,i,2) = c(j,i)*dzb(j)
357	c ttq(ihor,klev-j+1,i+nrad,2) = cni(j,i)*dzb(j)
358	c ttq(ihor,klev-j+1,i+2nrad,2) = c1i(j,i)dzb(j)
359	c ttq(ihor,klev-j+1,i+3nrad,2) = c2i(j,i)dzb(j)
360	c ttq(ihor,klev-j+1,i+4nrad,2) = c3i(j,i)dzb(j)
361	enddo
362	enddo
363
364
365	do j=1,klev
366	do i=1,nrad
367	! solesp en m3/m3 pour passer en m3/m2 il faut faire :
368	! (c1i(j,i)dzb(j) -q(IHOR,j,i+2nrad))
369	solesp(ihor,klev+1-j,1)=solesp(ihor,klev+1-j,1) +
370	& (c1i(j,i)-q(IHOR,j,i+2*nrad)/dzb(j))
371	solesp(ihor,klev+1-j,2)=solesp(ihor,klev+1-j,2) +
372	& (c2i(j,i)-q(IHOR,j,i+3*nrad)/dzb(j))
373	solesp(ihor,klev+1-j,3)=solesp(ihor,klev+1-j,3) +
374	& (c3i(j,i)-q(IHOR,j,i+4*nrad)/dzb(j))
375	enddo
376	enddo
377
378	* concerne les gouttes (tableaux cn,c1,c2):
379
380	c do j=1,klev
381	c do i=1,nrad
382	c tttq(ihor,klev-j+1,i,1) = c(j,i)*dzb(j)
383	c tttq(ihor,klev-j+1,i+nrad,1) = cni(j,i)*dzb(j)
384	c tttq(ihor,klev-j+1,i+2nrad,1) = c1i(j,i)dzb(j)
385	c tttq(ihor,klev-j+1,i+3nrad,1) = c2i(j,i)dzb(j)
386	c tttq(ihor,klev-j+1,i+4nrad,1) = c3i(j,i)dzb(j)
387	c enddo
388	c enddo
389
390	c call begintime(tt0)
391	call snuages(ngrid,cni,c1i,c2i,c3i,c,ddt,
392	& klev,nrad,ihor,
393	& flxesp_i,tau_drop,prec)
394	c call endtime(tt0,tt1)
395	c ttsclds=ttsclds+tt1
396
397	c do j=1,klev
398	c do i=1,nrad
399	c tttq(ihor,klev-j+1,i,2) = c(j,i)*dzb(j)
400	c tttq(ihor,klev-j+1,i+nrad,2) = cni(j,i)*dzb(j)
401	c tttq(ihor,klev-j+1,i+2nrad,2) = c1i(j,i)dzb(j)
402	c tttq(ihor,klev-j+1,i+3nrad,2) = c2i(j,i)dzb(j)
403	c tttq(ihor,klev-j+1,i+4nrad,2) = c3i(j,i)dzb(j)
404	c enddo
405	c enddo
406	ENDIF ! flag nuages :)
407
408
409	* on recompose le tableau de traceurs ici.
410	*--------------------------------------------------
411
412	do i=1,nrad
413	do j=1,klev
414	q(IHOR,j,i) = c(j,i)*dzb(j) ! nombre aerosols /m^2
415	if (clouds.eq.1) then
416	q(IHOR,j,i+ nrad) = cni(j,i)*dzb(j) ! nombre noyaux /m^2
417	q(IHOR,j,i+2nrad) = c1i(j,i)dzb(j) ! concentration volume glace/m^2
418	q(IHOR,j,i+3nrad) = c2i(j,i)dzb(j) ! concentration volume glace/m^2
419	q(IHOR,j,i+4nrad) = c3i(j,i)dzb(j) ! concentration volume glace/m^2
420	endif
421	enddo
422	enddo
423	if (clouds.eq.1) then
424	do j=1,klev
425	gaz1(IHOR,klev-j+1) = gazc1(j) ! fraction molaire
426	gaz2(IHOR,klev-j+1) = gazc2(j) ! fraction molaire
427	gaz3(IHOR,klev-j+1) = gazc3(j) ! fraction molaire
428	enddo
429	endif
430
431	c***************************************************************
432	else ! same latitude, we don't do calculations again
433	q(ihor,:,:) = q(im1,:,:)
434	tau_aer(ihor,:,:) = tau_aer(im1,:,:)
435	prec(ihor,:) = prec(im1,:)
436	if (clouds.eq.1) then
437	solesp(ihor,:,:) = solesp(im1,:,:)
438	flxesp_i(ihor,:,:) = flxesp_i(im1,:,:)
439	tau_drop(ihor,:) = tau_drop(im1,:)
440	gaz1(ihor,:) = gaz1(im1,:)
441	gaz2(ihor,:) = gaz2(im1,:)
442	gaz3(ihor,:) = gaz3(im1,:)
443	endif
444	endif
445
446	ENDDO ! Fin de la boucle IHOR
447	c***************************************************************
448
449	102 CONTINUE ! la premiere fois, c'est une boucle vide!
450
451
452	c***************************************************************
453	c FIN: on renvoie les nouvelles valeurs q(t+dt)=q(t) + dq(t)
454	c
455	c Pour les aerosols, les noyaux, les glaces et les vapeurs modifiees...
456	c
457	c***************************************************************
458
459	do n=1,ngrid
460	do i=1,nmicro
461	do j=1,klev ! j de 1 a 54
462	tq(n,j,i) = q(n,klev+1-j,i)
463	enddo
464	enddo
465	enddo
466
467
468	c do j=1,15
469	cc CH4 -- cnuages
470	c write(210,'(I4,3(ES24.16,1X))') j,
471	c & sum(ttq(40,j,2nrad+1:3nrad,1)),
472	c & sum(ttq(40,j,2nrad+1:3nrad,2)),
473	c & sum(ttq(40,j,2nrad+1:3nrad,2))-
474	c & sum(ttq(40,j,2nrad+1:3nrad,1))
475	cc C2H6 -- cnuages
476	c write(211,'(I4,3(ES24.16,1X))') j,
477	c & sum(ttq(40,j,3nrad+1:4nrad,1)),
478	c & sum(ttq(40,j,3nrad+1:4nrad,2)),
479	c & sum(ttq(40,j,3nrad+1:4nrad,2))-
480	c & sum(ttq(40,j,3nrad+1:4nrad,1))
481	cc C2H2 -- cnuages
482	c write(212,'(I4,3(ES24.16,1X))') j,
483	c & sum(ttq(40,j,4nrad+1:5nrad,1)),
484	c & sum(ttq(40,j,4nrad+1:5nrad,2)),
485	c & sum(ttq(40,j,4nrad+1:5nrad,2))-
486	c & sum(ttq(40,j,4nrad+1:5nrad,1))
487	c
488	cc CH4 -- snuages
489	c write(310,'(I4,3(ES24.16,1X))') j,
490	c & sum(tttq(40,j,2nrad+1:3nrad,1)),
491	c & sum(tttq(40,j,2nrad+1:3nrad,2)),
492	c & sum(tttq(40,j,2nrad+1:3nrad,2))-
493	c & sum(tttq(40,j,2nrad+1:3nrad,1))
494	cc C2H6 -- snuages
495	c write(311,'(I4,3(ES24.16,1X))') j,
496	c & sum(tttq(40,j,3nrad+1:4nrad,1)),
497	c & sum(tttq(40,j,3nrad+1:4nrad,2)),
498	c & sum(tttq(40,j,3nrad+1:4nrad,2))-
499	c & sum(tttq(40,j,3nrad+1:4nrad,1))
500	cc C2H2 -- snuages
501	c write(312,'(I4,3(ES24.16,1X))') j,
502	c & sum(tttq(40,j,4nrad+1:5nrad,1)),
503	c & sum(tttq(40,j,4nrad+1:5nrad,2)),
504	c & sum(tttq(40,j,4nrad+1:5nrad,2))-
505	c & sum(tttq(40,j,4nrad+1:5nrad,1))
506	c enddo
507	c write(210,*) "NEWLINE"
508	c write(211,*) "NEWLINE"
509	c write(212,*) "NEWLINE"
510	c write(310,*) "NEWLINE"
511	c write(311,*) "NEWLINE"
512	c write(312,*) "NEWLINE"
513
514
515	c do j=1,20
516	c write(410,'(I4,3(ES24.16,1X))') j,
517	c & flxesp_i(40,j,1),flxesp_i(40,j,2),flxesp_i(40,j,3)
518	c write(510,'(I4,3(ES24.16,1X))') j,
519	c & solesp(40,j,1),solesp(40,j,2),solesp(40,j,3)
520	c enddo
521	c write(410,*) "NEWLINE"
522	c write(510,*) "NEWLINE"
523
524	IPREM=1 ! LA PROCHAINE FOIS NE SERA PLUS LA 1ERE
525
526
527	16 return
528
529	end
530
531
532	c---------------------------------------------------------------------

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

Download in other formats:

Original Format