Context Navigation

phytrac.F @ 1243

Last change on this file since 1243 was 1126, checked in by slebonnois, 11 years ago
SL: update of Titan photochemical module to include computation of chemistry up to 1300 km
File size: 25.6 KB

Line
1	SUBROUTINE phytrac (firstcall,lastcall,
2	. nqmax,nmicro,ptimestep,appkim,dtkim,
3	. pplev,pplay,delp,ptemp,pmu0,pfract,pdecli,
4	. lonsol,
5	. pu,pv,pzlev,pzlay,ftsol,
6	. tr_seri,qaer,d_tr_mph,d_tr_kim,
7	. fclat,reservoir)
8
9	c======================================================================
10	c S. Lebonnois, mai 2008
11	c
12	c Arguments:
13	c
14	c firstcall----input-L-variable logique indiquant le premier passage
15	c lastcall-----input-L-variable logique indiquant le dernier passage
16	c nqmax--------input-I-nombre de traceurs (total)
17	c nmicro-------input-I-nombre de traceurs microphysiques !! doivent etre toujours en premiers!!
18	c ptimestep----input-R-pas d integration pour la physique (seconde)
19	c appkim-------input-I-appel a la chimie
20	c dtkim--------input-R-pas de temps chimique (seconde)
21	c pplev--------input-R-pression pour chaque inter-couche (en Pa)
22	c pplay--------input-R-pression pour chaque couche (en Pa)
23	c delp---------input-R-epaisseur d'une couche (en Pa)
24	c ptemp--------input-R-temperature (K)
25	c pmu0---------input-R-cos angle zenithal
26	c pfract-------input-R-fractional day
27	c pdecli-------input-R-declinaison en radian
28	c lonsol-------input-R-longitude solaire en radian
29	c pu-----------input-R-vitesse dans la direction X (de O a E) en m/s (1ere couche)
30	c pv-----------input-R-vitesse Y (de S a N) en m/s (1ere couche)
31	c pzlev--------input-R-altitude pour chaque inter-couche (en m)
32	c pzlay--------input-R-altitude pour chaque couche (en m)
33	c ftsol--------input-R-temperature au sol (en K)
34	c tr_seri------input-R-mass mixing ratio traceurs (kg/kg)
35	c d_tr_mph----output-R-tendance microphysique de "qx" (kg/kg/s)
36	c d_tr_kim----output-R-tendance chimique de "qx" (kg/kg/s)
37	c fclat--------output-R-flux de chaleur latente d'evaporation du reservoir CH4 (J/m2/s)
38	c reservoir----outpur-R-un reservoir de surface !!! (m)
39	c======================================================================
40	USE infotrac
41	use dimphy
42	USE common_mod, only: rmcbar,xfbar,ncount,
43	& flxesp_i,tau_drop,tau_aer,solesp,precip,
44	& evapch4,occcld_m,occcld,satch4,satc2h6,satc2h2,rmcloud
45	USE moyzon_mod
46	USE write_field_phy
47	IMPLICIT none
48	#include "dimensions.h"
49	#include "clesphys.h"
50	#include "YOMCST.h"
51	#include "microtab.h"
52	#include "varmuphy.h"
53	#include "itemps.h"
54	#include "logic.h"
55
56	c======================================================================
57	c Variables argument:
58	c
59	LOGICAL firstcall,lastcall
60	INTEGER nqmax,nmicro,nlat,appkim
61	REAL ptimestep,dtkim
62	REAL pplev(klon,klev+1),pplay(klon,klev+1),delp(klon,klev)
63	REAL ptemp(klon,klev)
64	REAL pmu0(klon), pfract(klon), pdecli, lonsol
65	REAL pu(klon),pv(klon)
66	REAL pzlev(klon,klev+1),pzlay(klon,klev)
67	REAL ftsol(klon)
68	REAL tr_seri(klon,klev,nqmax)
69	REAL qaer(klon,klev,nqmax)
70	REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
71	REAL fclat(klon)
72	REAL reservoir(klon)
73
74	c======================================================================
75	c Local variables
76	REAL qaer0(klon,klev,nqmax)
77	REAL prec(klon,5)
78
79	REAL rcloud(klon,klev,nrad),xfrac(klon,klev,4)
80
81	REAL vcl,nuc,xgsn,xmsn,xesn,xasn
82
83
84	ReAL gaz1(klon,klev),gaz2(klon,klev),gaz3(klon,klev)
85
86	REAL socccld
87
88	c grandeurs en moyennes zonales
89	REAL zplev(klon,klev+1),zplay(klon,klev)
90	REAL zzlev(klon,klev+1),zzlay(klon,klev)
91	REAL ztemp(klon,klev), delpbar(klon,klev)
92	real temp_eq(klev),press_eq(klev)
93	REAL qaer0bar(klon,klev,nqmax) ! et non nmicro... Permet nmicro=0.
94	REAL zdqmufi(klon,klev,nqmax)
95	REAL ychim(klon,klev,nqmax-nmicro)
96	c La saturation n est calculee qu une seule fois: sauvegarde qysat
97	c La chimie n est pas calculee tous les pas, il faut donc
98	c sauvegarder les sorties de la chimie
99	REAL,save,allocatable :: qysat(:,:),pdyfi(:,:,:)
100
101	character*10 nomqy(nqmax-nmicro+1)
102	integer i,j,k,l,iq,ig0
103
104	c indice des esp chimiques utilisees dans la microfi
105	integer icldch4,icldc2h6,icldc2h2
106	save icldch4,icldc2h6,icldc2h2
107
108	real fte,ftm,Lvch4
109
110	REAL tmp,ex,kmin,kmax,dqsq
111	REAL dqch4
112
113	c======================================================================
114	c======================================================================
115
116	if (firstcall) then
117	allocate(qysat(klev,nqmax-nmicro),pdyfi(klon,klev,nqmax-nmicro))
118
119	c -------- Quelques verifications au demarrage sur les tailles des tableaux.
120	IF (microfi.ge.1) then
121	c Faire de la microphysique sans traceurs... bon courage !
122	if (nmicro.le.0) then
123	print*,"aLeRtE cRiTiQuE !!!"
124	print*,"Vous faites de la microphysique sans traceurs"
125	print*,"microphysique..."
126	print*,"Je m'arrete et vous laisse reflechir !"
127	stop
128	endif
129	c Nombre de traceurs incompatibles avec la microphysique.
130	if ((nmicro.ne.ntype*nrad).and.(clouds.eq.1)) then
131	print*,"aLeRtE cRiTiQuE !!!"
132	print*,"Nb trac imcompatible avec la microphysique."
133	print,nmicro,ntypenrad
134	stop
135	endif
136	if ((nmicro.ne.nrad).and.(clouds.eq.0)) then
137	print*,"aLeRtE cRiTiQuE !!!"
138	print*,"Nb trac imcompatible avec la microphysique."
139	print*,nmicro,nrad
140	stop
141	endif
142	ENDIF
143
144	endif ! firstcall
145
146	c RAZ des sorties : les moyennes se font directement dans IOIPSL :
147	c
148	flxesp_i(:,:,:) = 0.
149	tau_drop(:,:) = 0.
150	tau_aer(:,:,:) = 0.
151	solesp(:,:,:) = 0.
152	precip(:,:) = 0. ! c'est uniquement une sortie en um/s
153	c
154	prec(:,:) = 0. ! c'est la variable temporaire des precipitions de la microfi
155	! prec est en m (metre precipitable)
156
157	c-----------------------------------------------------------------------
158	c convertion moyennes zonales et changement d unites pour microphy
159	c ---------------------------------
160
161	c print*,'CONVERSION 2D ET CHANGEMENT UNITES (PHYTRAC)'
162
163	c -------------------
164	c Gestion de la temperature et de la pression :
165	c Utilisation des moyennes zonales:
166
167	c soit la chimie est active, soit la microphysique se fait en 2D.
168	IF (chimi.or.microfi.eq.1) THEN
169	zplev(:,:) = zplevbar(:,:)
170	zplay(:,:) = zplaybar(:,:)
171	zzlev(:,:) = zzlevbar(:,:)
172	zzlay(:,:) = zzlaybar(:,:)
173	ztemp(:,:) = ztfibar(:,:)
174	ychim = 0.0
175	ENDIF
176
177	c Si la microphysique est faite en 2D:
178	IF (microfi.eq.1) THEN
179	DO l=1,llm
180	DO i = 1, klon
181	delpbar(i,l) = zplevbar(i,l) - zplevbar(i,l+1)
182	ENDDO
183	ENDDO
184	c Traceurs microphysiques: passage en extensif: n/kg --> n/m^2
185	DO iq=1,nmicro
186	qaer(:,:,iq) = zqfibar(:,:,iq)*delpbar(:,:)/RG
187	DO l=1,klev
188	DO i = 1, klon
189	if (qaer(i,l,iq).lt.0.) then
190	print*,"NEGS ICI ICI !!!!",qaer(i,l,iq),i,l,iq
191	qaer(i,l,iq)=0.
192	c stop
193	endif
194	if (delpbar(i,l).lt.0.) then
195	print*,"NEGS DELP ICI !!!!",i,l,iq,delpbar(i,l)
196	stop
197	endif
198	ENDDO
199	ENDDO
200	qaer0(:,:,iq)= tr_seri(:,:,iq)*delp(:,:)/RG
201	qaer0bar(:,:,iq) = qaer(:,:,iq)
202	ENDDO
203	ENDIF
204
205	c Si la microphysique est faite en 3D:
206	IF (microfi.eq.2) THEN
207	zplev(:,:) = pplev(:,:)
208	zplay(:,:) = pplay(:,:)
209	zzlev(:,:) = pzlev(:,:)
210	zzlay(:,:) = pzlay(:,:)
211	ztemp(:,:) = ptemp(:,:)
212	c Traceurs microphysiques: passage en extensif: n/kg --> n/m^2
213	DO iq=1,nmicro
214	qaer(:,:,iq) = tr_seri(:,:,iq)*delp(:,:)/RG
215	qaer0(:,:,iq)= tr_seri(:,:,iq)*delp(:,:)/RG
216	ENDDO
217	ENDIF
218
219	do l=1,llm
220	temp_eq = tmoy
221	press_eq = playmoy/100. ! en mbar
222	enddo
223
224	c -------------------
225	c Extraction des gaz pour les nuages
226	IF ((microfi.ge.1).and.(clouds.eq.1)) THEN
227
228	c recuperation des indices des gaz qui nous interesse
229	if (firstcall) then
230	icldch4=-1
231	icldc2h6=-1
232	icldc2h2=-1
233	do i=1,nqmax
234	if (tname(i).eq."CH4") then
235	icldch4=i
236	c ich4=i
237	elseif (tname(i).eq."C2H6") then
238	icldc2h6=i
239	elseif (tname(i).eq."C2H2") then
240	icldc2h2=i
241	endif
242	enddo
243	if (icldch4 .eq.-1 .or.
244	& icldc2h6.eq.-1 .or.
245	& icldc2h2.eq.-1 ) then
246	print*, "Sacrebleu !!!"
247	print*, "Vous voulez faire des nuages sans gaz."
248	print*, "Mais vous etes inconscient. Je vais m'arreter la"
249	print*, "pour vous laisser reflechir au probleme"
250	STOP
251	endif
252	endif ! firstcall
253
254	c Saturation et fraction molaire CLOUD
255	c Calcul des saturations pour les esp chimique de la muphy des nuages.
256	c On le fait ici pour les sortir dans physiq.F sans avoir a surcharger la routine.
257	c Elles passent ensuite dans un common pour passer dans les I/O.
258
259	DO l=1,llm
260	DO i = 1, klon
261	call ch4sat(ptemp(i,l),pplay(i,l),tmp) !tmp en kg/kg !
262	satch4(i,l) = tr_seri(i,l,icldch4)/(tmp*28./16.)
263
264	call c2h6sat(ptemp(i,l),pplay(i,l),tmp)
265	satc2h6(i,l) =tr_seri(i,l,icldc2h6)/(tmp*28./30.)
266
267	call c2h2sat(ptemp(i,l),pplay(i,l),tmp)
268	satc2h2(i,l) =tr_seri(i,l,icldc2h2)/(tmp*28./26.)
269	ENDDO
270	ENDDO
271
272	c Copie des gaz (en 3D) <== UNIQUEMENT SI ON FAIT DES NUAGES
273	if (moyzon_mu) then
274	gaz1(:,:) = zqfibar(:,:,icldch4)
275	gaz2(:,:) = zqfibar(:,:,icldc2h6)
276	gaz3(:,:) = zqfibar(:,:,icldc2h2)
277	else
278	gaz1(:,:) = tr_seri(:,:,icldch4)
279	gaz2(:,:) = tr_seri(:,:,icldc2h6)
280	gaz3(:,:) = tr_seri(:,:,icldc2h2)
281	endif
282
283	endif ! microfi.ge.1 + clouds.eq.1
284	c -------------------
285
286	c AUTRES TRACEURS
287
288	if (nqmax.gt.nmicro) then
289	do iq=nmicro+1,nqmax
290	if (moyzon_ch) then
291	ychim(:,:,iq-nmicro) = zqfibar(:,:,iq)
292	else
293	ychim(:,:,iq-nmicro) = tr_seri(:,:,iq)
294	endif
295	nomqy(iq-nmicro) = tname(iq)
296	c print*,iq-nmicro,nomqy(iq-nmicro)
297	enddo
298	nomqy(nqmax-nmicro+1) = "HV"
299	endif
300
301	c-----------------------------------------------------------------------
302	c initialisation des qysat au premier appel:
303	c ---------------------------------
304
305	c!! ATTENTION, qysat pris uniquement a l'equateur
306	c!! justifie puisque dans cette region, les var de t et p sont faibles...
307
308	if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
309	call inicondens(nqmax-nmicro,press_eq,temp_eq,nomqy,qysat)
310	endif
311
312	c-----------------------------------------------------------------------
313	c Appel de la microphysique en 2D/3D !!!!!!
314	c --------------------------
315
316	IF(firstcall) THEN
317	print*,'MICROPHYSIQUE ',MICROFI
318	ENDIF
319
320	c call begintime(tt0)
321	IF (MICROFI.eq.0) THEN
322	c PAS DE MICROPHYSIQUE :
323	IF (firstcall) THEN
324	print*,'MICROPHYSIQUE OFF-LINE',MICROFI
325	ENDIF
326	ELSE
327	zdqmufi = 0. ! ne sert que pour chimi pour condensation
328	call muphys(klon,
329	& zplev,zplay,zzlev,zzlay,
330	& ztemp,qaer,gaz1,gaz2,gaz3,
331	& nmicro,ptimestep,
332	& pmu0,pfract,
333	c -------- sorties diagnostiques
334	& flxesp_i,
335	& tau_drop,tau_aer,
336	& solesp,prec)
337
338	c NOTES :
339	c Ici toutes nos sorties sont des champs 3D...(meme les diagnostiques)
340	c On a rien a faire mis a part copier les dq dans les d_tr
341
342	ENDIF
343	c call endtime(tt0,tt1)
344	c ttmuphys=ttmuphys+tt1
345
346	c-----------------------------------------------------------------------
347	c Gestion des sources
348	c -------------
349	c
350	IF (clouds.eq.1) THEN
351	IF (microfi.eq.1) THEN
352	c On repasse les gaz en 3D si on a fait de la microphysique en 2D
353	gaz1(:,:)=gaz1(:,:)*tr_seri(:,:,icldch4)/zqfibar(:,:,icldch4)
354	gaz2(:,:)=gaz2(:,:)*tr_seri(:,:,icldc2h6)/zqfibar(:,:,icldc2h6)
355	gaz3(:,:)=gaz3(:,:)*tr_seri(:,:,icldc2h2)/zqfibar(:,:,icldc2h2)
356	ENDIF
357	c Mise a jour du reservoir de CH4 (ie : seul le CH4 remplit le reservoir)
358	DO i=1,klon
359	reservoir(i) = reservoir(i)+prec(i,1)
360	ENDDO
361
362	CALL sources(klon,klev,ptimestep,z0,
363	& pu,pv,pplev,pzlay,pzlev,
364	& gaz1,gaz2,gaz3,
365	& ftsol,evapch4,reservoir)
366
367	ENDIF
368	c-----------------------------------------------------------------------
369	c Condensation
370	c -------------
371
372	IF ((chimi).and.(nqmax.gt.nmicro)) then
373
374	c tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
375	c print*,'Condensation'
376
377	do iq=1,nqmax-nmicro
378	do l=1,llm
379	do j=1,klon
380	if (ychim(j,l,iq).gt.qysat(l,iq)) then
381	zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
382	. / ptimestep ! / dt
383	endif
384	enddo
385	enddo
386	enddo
387
388	ENDIF
389
390	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
391	c eventuellement, modif initiale de la compo
392	c
393	c tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
394	c
395	c if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
396	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
397	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
398	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
399	c do iq=1,nqmax-nmicro
400	c if (nomqy(iq).eq."CH4") then
401	c do l=1,llm
402	c do j=1,klon
403	c if (ychim(j,l,iq).le.0.015) then
404	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
405	c . / ptimestep ! / dt
406	c endif
407	c enddo
408	c enddo
409	c endif
410	c enddo
411	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
412	c
413	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
414	c!!!remise de C2H2 a 1.e-5 max !!!!!!!!!!!!!
415	c!!!remise de C2H6 a 3.e-5 max !!!!!!!!!!!!!
416	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
417	c do iq=1,nqmax-nmicro
418	c if (nomqy(iq).eq."C2H2") then
419	c do l=1,llm
420	c do j=1,klon
421	c if (ychim(j,l,iq).gt.1.e-5) then
422	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
423	c . / ptimestep ! / dt
424	c endif
425	c enddo
426	c enddo
427	c endif
428	c if (nomqy(iq).eq."C2H6") then
429	c do l=1,llm
430	c do j=1,klon
431	c if (ychim(j,l,iq).gt.3.e-5) then
432	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
433	c . / ptimestep ! / dt
434	c endif
435	c enddo
436	c enddo
437	c endif
438	c enddo
439	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
440	c endif
441	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
442
443	c ----- commentaire de fin (mise a jour des profil de fraction molaire)
444
445	c-----------------------------------------------------------------------
446	c Appel de la chimie
447	c --------------------------
448
449	if((appkim.eq.1).and.(chimi)) then
450	print*,'On passe dans la CHIMIE'
451
452	c do iq=1,nqmax-nmicro
453	c if (nomqy(iq).eq."C2H2") then
454	c print*,"C2H2top=",ychim(:,klev,iq)
455	c endif
456	c enddo
457
458	c Appel Chimie
459	c ------------
460	CALL calchim(klon,nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
461	. ztemp,zplay,zplev,zzlay,zzlev,
462	. pdyfi)
463	c ychim ne doit pas etre modifie, pdyfi en /s
464
465	endif
466
467	c-----------------------------------------------------------------------
468	c retour des tendances vers 3D
469	c ---------------------------------
470
471	c call WriteField_phy('phytrac_qaer01',qaer(:,:,1),klev)
472
473	c===============================
474	c TRACEURS MICROPHYSIQUES
475	c
476	c===============================
477	c ---> pas de microphysique
478	IF (microfi.eq.0) THEN
479	DO iq=1,nmicro
480	d_tr_mph(:,:,iq)=0.
481	ENDDO
482	ENDIF
483	c===============================
484	c ---> microphysique 2D
485	c IF (microfi.eq.1) THEN
486	c DO iq=1,nmicro
487	c DO l=1,llm
488	c DO i=1,klon
489	c on repasse le champ de traceurs en 3D (pas les tendances)
490	c qaer est ce qui entre dans muphy, donc la moyenne zonale
491	c qaer0 est la valeur initiale du champ
492	c qaer0bar est la moyenne zonale initiale
493	c la variation relative pour une bande de latitude est donc (qaer/qaer0bar)
494	c la nouvelle valeur en un point (3D) est donc qaer0*(qaer/qaer0bar)
495	c et la tendance: qaer0*(qaer/qaer0bar)-qaer0
496	c un petit patch :
497	c Si la moyenne zonale au depart est "nulle" :
498	c On a quand meme le droit de produire des traceurs dans la cellule.
499	c On considere donc que la valeur de sortie 3D correspond a la valeur de sortie 2D.
500	c Cela permet aussi entre autre d eviter les NaN pour les traceurs des nuages !
501	c (au dessus de la tropo pas de nuages donc qaer(nrad+1:ntype*nrad) = 0 !!!)
502	c IF (qaer0bar(i,l,iq).gt.1e-100) THEN
503	c qaer(i,l,iq) = qaer0(i,l,iq) *
504	c & qaer(i,l,iq)/qaer0bar(i,l,iq)
505	c ENDIF
506	c La tendance correspond a (qaer-qaer0)/ptimestep
507	c d_tr_mph(i,l,iq) = (qaer(i,l,iq)-qaer0(i,l,iq))/
508	c & ptimestep
509	c ENDDO
510	c ENDDO
511	c ENDDO
512	c ---> microphysique 3D
513	c ELSEIF(microfi.gt.1) THEN
514	c DO iq=1,nmicro
515	c d_tr_mph(:,:,iq)=(qaer(:,:,iq)-qaer0(:,:,iq))/ptimestep
516	c ENDDO
517	c ENDIF ! microfi
518
519	c DO iq=1,nmicro
520	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
521	c d_tr_mph(:,:,iq) = d_tr_mph(:,:,iq)*RG/delp(:,:)
522	c ENDDO
523
524	c===============================
525	c TOUT CE QUI EST AU-DESSUS NE MARCHE PAS: PLEIN DE NEGS...
526	c CA MARCHE EN 3D, MAIS PAS EN MOY ZON...
527	c===============================
528
529	c ---> microphysique 2D
530	IF (microfi.eq.1) THEN
531	DO iq=1,nmicro
532	DO l=1,llm
533	DO i=1,klon
534	c ici, qaer correspond a la moy zonale modifiee par la microphys.
535	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
536	c en mettant ceci:
537	d_tr_mph(i,l,iq) = (qaer(i,l,iq)*RG/delpbar(i,l)
538	& -qaer0(i,l,iq)*RG/delp(i,l))/
539	& ptimestep
540	c on remplace le champ 3D initial par la valeur modifiee de sa moyenne zonale
541	c => on remet un champ uniforme en zonal...
542	ENDDO
543	ENDDO
544	ENDDO
545	c ---> microphysique 3D
546	ELSEIF(microfi.gt.1) THEN
547	DO iq=1,nmicro
548	d_tr_mph(:,:,iq)=(qaer(:,:,iq)-qaer0(:,:,iq))/ptimestep
549	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
550	d_tr_mph(:,:,iq) = d_tr_mph(:,:,iq)*RG/delp(:,:)
551	ENDDO
552	ENDIF ! microfi
553
554	c===============================
555
556	c AUTRES TRACEURS
557
558	if ((chimi).and.(nqmax.gt.nmicro)) then
559	c on passe de pdyfi (tendance chimique en /s calculee quand chimie appelee)
560	c a d_tr_kim (tendance chimique 3D en /s, passee a physiq)
561	c et de zdqmufi a d_tr_mph (tendance condensation 3D en /s passee a physiq)
562
563	DO iq=nmicro+1,nqmax
564	d_tr_kim(:,:,iq) = pdyfi(:,:,iq-nmicro)
565	& *tr_seri(:,:,iq)/ychim(:,:,iq-nmicro)
566	d_tr_mph(:,:,iq) = zdqmufi(:,:,iq)
567	& *tr_seri(:,:,iq)/ychim(:,:,iq-nmicro)
568	ENDDO
569
570	endif ! chimi
571
572	c--------------------------------------------------
573	c CONDENSATION VIA MICROFI
574	c----------------------
575	c La microphysique avec nuages doit se faire obligatoirement en 3D. (FAUX ACTUELLEMENT)
576	c Rien n empeche de faire la chimie en 2D. Cependant pour prendre en compte la
577	c condensation due a la microfi (en 3D) on recalcule la tendance finale pour
578	c les especes concernees (CH4, C2H6 pour le moment).
579	IF (microfi.ge.1.and.clouds.eq.1) THEN
580	c condensation CH4
581	d_tr_mph(:,:,icldch4) =(gaz1(:,:)-tr_seri(:,:,icldch4))
582	& /ptimestep
583	c condensation C2H6
584	d_tr_mph(:,:,icldc2h6)=(gaz2(:,:)-tr_seri(:,:,icldc2h6))
585	& /ptimestep
586	c condensation C2H2
587	d_tr_mph(:,:,icldc2h2)=(gaz3(:,:)-tr_seri(:,:,icldc2h2))
588	& /ptimestep
589	ENDIF
590
591	c--------------------------------------------------
592	c MISE A JOUR CH4 : (pour refixer la fraction
593	c molaire)
594	c--------------------------------------------------
595	c IF (firstcall) THEN
596	c do i=1,klon
597	c do j=1,llm
598	c call ch4sat(ptemp(i,j),pplay(i,j),tmp) !tmp en kg/kg !
599	c tmp=0.950.85tmp*28./16.
600	c if (pplay(i,j).lt.20000.) then
601	c dqch4 = 1.4e-2
602	c else
603	c dqch4 = tmp
604	c endif
605	c d_tr_mph(i,j,icldch4)=(-tr_seri(i,j,icldch4)+dqch4)/
606	c & ptimestep
607	c enddo
608	c enddo
609	c
610	c ENDIF
611
612	c--------------------------------------------------
613	c CONVERSION PRECIPITATION :
614	c en microns/secondes
615	c--------------------------------------------------
616	precip = prec * 1.e6 / ptimestep
617
618
619	c--------------------------------------------------
620	c CALCUL DU FLUX DE CHALEUR LATENTE D EVAPORATION
621	c DU METHANE
622	c--------------------------------------------------
623	IF (clouds.eq.1) THEN
624	DO i=1,klon
625	fte= (1.-ftsol(i)/305.5)
626	ftm= (1.-ftsol(i)/190.5)
627	if(ftm.le.1.e-3) ftm=1.e-3
628	if(fte.le.1.e-3) fte=1.e-3
629	Lvch4 =8.314190.4
630	& (7.08ftm0.354+10.951.1e-2ftm*0.456)
631	& /mch4
632	! evapch4 en m3/m2 {ok}
633	! 425 en kg/m3
634	! Lv en J/kg {ok}
635	! ptimestep en s {ok}
636	fclat(i)=(evapch4(i)Lvch4rhoi_ch4) ! en J/m2/s
637	ENDDO
638	ENDIF
639
640	c--------------------------------------------------
641	c GESTION DES RAYONS DE GOUTTES POUR TR
642	c--------------------------------------------------
643	IF (clouds.eq.1) THEN
644
645	c Calcul du rayon des gouttes par bin ...
646	c----------------------------------------
647	DO i=1,klon
648	DO j=1,klev
649	DO iq=1,nrad
650	* Rayon minimum selon la quantite de noyaux
651	IF (qaer(i,j,iq+nrad) .le. 1.e-5) THEN
652	rcloud(i,j,iq) = 1.e-10
653	ELSE
654	rcloud(i,j,iq)=
655	& ((qaer(i,j,iq+2*nrad)/qaer(i,j,iq+nrad)+
656	& qaer(i,j,iq+3*nrad)/qaer(i,j,iq+nrad) +
657	& v_e(iq))0.75/RPI)*(1./3.)
658	ENDIF
659	ENDDO
660	ENDDO
661	ENDDO
662
663	c .... et de leur rayon moyen total (tt bins confondu)
664	c------------------------------------------------------
665	DO i=1,klon
666	socccld=0.
667	DO j=klev,1,-1 !de haut en bas pour le calcul des opacites
668	vcl=0.
669	nuc=0.
670	xgsn=0.
671	xmsn=0.
672	xesn=0.
673	xasn=0.
674	DO iq=1,nrad
675	vcl=vcl+qaer(i,j,iq+2*nrad)+
676	& qaer(i,j,iq+3*nrad)+
677	& qaer(i,j,iq+4*nrad)+
678	& v_e(iq)*qaer(i,j,iq+nrad) ! volume des gouttes
679	nuc=nuc+qaer(i,j,iq+nrad) ! nombre de noyaux
680	xgsn=xgsn+qaer(i,j,iq+nrad)*v_e(iq) ! volume de noyaux
681	xmsn=xmsn+qaer(i,j,iq+2*nrad) ! volume de methane
682	xesn=xesn+qaer(i,j,iq+3*nrad) ! volume d' ethane
683	xasn=xasn+qaer(i,j,iq+4*nrad) ! volume d' acethylene
684	ENDDO
685	IF (nuc .le. 1.e-5) THEN
686	rmcloud(i,j)=1.0e-10
687	xfrac(i,j,:)=0.
688	ELSE
689	IF(xgsn/vcl.lt.0. .or. xgsn/vcl.gt.1.001)
690	& print*, 'PB AVEC XFRAC:', i,j,xgsn,vcl
691	rmcloud(i,j)= ! rayon moyen des gouttes
692	& (vcl/nuc0.75/RPI)*(1./3.)
693	xfrac(i,j,1)=xgsn/vcl ! fraction volumique noyau/goutte
694	xfrac(i,j,2)=xmsn/vcl ! fraction volumique CH4/goutte
695	xfrac(i,j,3)=xesn/vcl ! fraction volumique C2H6/goutte
696	xfrac(i,j,4)=xasn/vcl ! fraction volumique C2H2/goutte
697	c calcul du rayon moyen (moyenne temporelle)
698	rmcbar(i,j)=rmcbar(i,j)+rmcloud(i,j)
699	xfbar(i,j,:)=xfbar(i,j,:)+xfrac(i,j,:)
700	ncount(i,j) = ncount(i,j)+1
701	ENDIF
702	socccld=socccld+RPI(rmcloud(i,j)2.)nuc
703	occcld(i,j)=socccld
704	ENDDO
705	ENDDO
706	c
707	c OCCCLD
708	c Calcul le nombre d occurence d un nuage
709	c d opacite comprise en kmin et kmax
710	c k kmin kmax
711	c 1 0.0000000 0.10000000
712	c 2 0.10000000 0.17782794
713	c 3 0.17782794 0.31622776
714	c 4 0.31622776 0.56234139
715	c 5 0.56234139 1.0000000
716	c 6 1.0000000 1.7782795
717	c 7 1.7782795 3.1622777
718	c 8 3.1622777 5.6234136
719	c 9 5.6234136 10.000000
720	c 10 10.000000 17.782795
721	c 11 17.782795 31.622778
722	c 12 31.622778 100000.00
723	c
724	c mise a zero de occld_m
725	occcld_m=0.
726	DO i=1,klon
727	DO j=1,klev
728	DO k=1,12
729	ex=10.**(0.25)
730	kmin=0.
731	kmax=1.e5
732	if(k.ne.1) kmin=0.1ex*(k-2)
733	if(k.ne.12) kmax=0.1ex*(k-1)
734	if(occcld(i,j).ge.kmin .and. occcld(i,j).lt.kmax)
735	& occcld_m(i,j,k)=1.
736	ENDDO
737	ENDDO
738	ENDDO
739	ENDIF ! fin condition clouds => pas besoin de calculer des rayons
740
741	RETURN
742	END
743

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

Download in other formats:

Original Format