Context Navigation

phytrac.F @ 175

Last change on this file since 175 was 175, checked in by slebonnois, 14 years ago

S.LEBONNOIS:

Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
Quelques ajouts dans la doc.

File size: 31.2 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	IMPLICIT none
43	#include "dimensions.h"
44	#include "clesphys.h"
45	#include "YOMCST.h"
46	#include "microtab.h"
47	#include "varmuphy.h"
48	#include "diagmuphy.h"
49	#include "itemps.h"
50
51	c======================================================================
52	c Variables argument:
53	c
54	LOGICAL firstcall,lastcall
55	INTEGER nqmax,nmicro,nlat,appkim
56	REAL ptimestep,dtkim
57	REAL pplev(klon,klev+1),pplay(klon,klev+1),delp(klon,klev)
58	REAL ptemp(klon,klev)
59	REAL pmu0(klon), pfract(klon), pdecli, lonsol
60	REAL pu(klon),pv(klon)
61	REAL pzlev(klon,klev+1),pzlay(klon,klev)
62	REAL ftsol(klon)
63	REAL tr_seri(klon,klev,nqmax)
64	REAL qaer(klon,klev,nqmax)
65	REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
66	REAL fclat(klon)
67	REAL reservoir(klon)
68
69	c======================================================================
70	c Local variables
71	REAL qaer0(klon,klev,nqmax)
72
73	c ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
74	INTEGER ngrid,NLAYER
75	PARAMETER (ngrid=(jjm-1)*iim+2) ! = klon
76	PARAMETER (NLAYER=llm) ! = klev
77	* common relatifs au nuages
78	real rmcbar(ngrid,NLAYER),xfbar(ngrid,NLAYER,4)
79	integer ncount(ngrid,NLAYER)
80	COMMON/rnuabar/ncount,rmcbar,xfbar
81
82	REAL rcloud(klon,klev,nrad),xfrac(klon,klev,4)
83
84	REAL vcl,nuc,xgsn,xmsn,xesn,xasn
85
86
87	ReAL gaz1(klon,klev),gaz2(klon,klev),gaz3(klon,klev)
88
89	REAL socccld
90
91	c grandeurs en moyennes zonales
92	REAL zplev(jjm+1,klev+1),zplay(jjm+1,klev),ztsol(jjm+1)
93	REAL zzlev(jjm+1,klev+1),zzlay(jjm+1,klev)
94	REAL ztemp(jjm+1,klev),zmu0(jjm+1),zfract(jjm+1)
95	real temp_eq(klev),press_eq(klev)
96	REAL zqaer(jjm+1,klev,nqmax) ! et non nmicro... Permet nmicro=0.
97	REAL zqaer0(jjm+1,klev,nqmax)
98	REAL zdqmufi(jjm+1,klev,nqmax)
99	REAL ychim(jjm+1,klev,nqmax-nmicro)
100	REAL zgaz1(jjm+1,klev),zgaz2(jjm+1,klev),zgaz3(jjm+1,klev)
101	REAL zgaz10(jjm+1,klev),zgaz20(jjm+1,klev),zgaz30(jjm+1,klev)
102	c La saturation n est calculee qu une seule fois: sauvegarde qysat
103	c La chimie n est pas calculee tous les pas, il faut donc
104	c sauvegarder les sorties de la chimie
105	REAL,save,allocatable :: qysat(:,:),pdyfi(:,:,:)
106
107	character*10 nomqy(nqmax-nmicro+1)
108	integer i,j,k,l,iq,ig0
109
110	REAL zprec(jjm+1,5),zsolesp(jjm+1,klev+1,3),
111	& zflxesp_i(jjm+1,klev,3)
112	REAL ztau_drop(jjm+1,klev),ztau_aer(jjm+1,klev,nrad)
113	c
114	c indice des esp chimiques utilisees dans la microfi
115	integer icldch4,icldc2h6,icldc2h2
116	save icldch4,icldc2h6,icldc2h2
117
118	real fte,ftm,Lvch4
119
120	REAL tmp,ex,kmin,kmax,dqsq
121	REAL dqch4
122	c REAL ch4(jjm+1),ch4b(jjm+1),dch4(jjm+1),ch4c(jjm+1,llm)
123	c integer ich4
124	c common/ch4ind/ich4
125
126	c======================================================================
127	c======================================================================
128
129	if (firstcall) then
130	allocate(qysat(klev,nqmax-nmicro),pdyfi(jjm+1,klev,nqmax-nmicro))
131
132	c -------- Quelques verifications au demarrage sur les tailles des tableaux.
133	IF (microfi.ge.1) then
134	c Faire de la microphysique sans traceurs... bon courage !
135	if (nmicro.le.0) then
136	print*,"aLeRtE cRiTiQuE !!!"
137	print*,"Vous faites de la microphysique sans traceurs"
138	print*,"microphysique..."
139	print*,"Je m'arrete et vous laisse reflechir !"
140	stop
141	endif
142	c Nombre de traceurs incompatibles avec la microphysique.
143	if (nmicro.ne.ntype*nrad) then
144	print*,"aLeRtE cRiTiQuE !!!"
145	print*,"Nb trac imcompatible avec la microphysique."
146	print,nmicro,ntypenrad
147	stop
148	endif
149	ENDIF
150
151	endif
152
153	c RAZ des sorties : les moyennes se font directement dans IOIPSL :
154	c
155	flxesp_i(:,:,:) = 0.
156	tau_drop(:,:) = 0.
157	tau_aer(:,:,:) = 0.
158	solesp(:,:,:) = 0.
159	prec(:,:) = 0.
160
161	c-----------------------------------------------------------------------
162	c convertion moyennes zonales et changement d unites pour microphy
163	c ---------------------------------
164
165	c print*,'CONVERSION 2D ET CHANGEMENT UNITES (PHYTRAC)'
166
167	c -------------------
168	c Gestion de la temperature et de la pression :
169	c soit la chimie est active, soit la microphysique se fait en 2D.
170	IF (chimi.or.microfi.eq.1) THEN
171
172	zplev = 0.0
173	zplay = 0.0
174	zzlev = 0.0
175	zzlay = 0.0
176	ztemp = 0.0
177	zqaer = 0.0
178	ychim = 0.0
179	zmu0 = 0.0
180	zfract= 0.0
181	zgaz1 = 0.0
182	zgaz2 = 0.0
183	zgaz3 = 0.0
184	zprec = 0.0
185	zflxesp_i = 0.0
186	ztau_drop = 0.0
187	ztau_aer = 0.0
188	zsolesp = 0.0
189
190	do l=1,llm+1
191	zplev(1,l) = pplev(1,l)
192	zzlev(1,l) = pzlev(1,l)
193	do j=2,jjm
194	ig0=1+(j-2)*iim
195	do i=1,iim
196	zplev(j,l) = zplev(j,l) + pplev(ig0+i,l)/iim
197	zzlev(j,l) = zzlev(j,l) + pzlev(ig0+i,l)/iim
198	enddo
199	enddo
200	zplev(jjm+1,l) = pplev(klon,l)
201	zzlev(jjm+1,l) = pzlev(klon,l)
202	enddo
203
204	do l=1,llm
205	ztemp(1,l) = ptemp(1,l)
206	zplay(1,l) = pplay(1,l)
207	zzlay(1,l) = pzlay(1,l)
208	do j=2,jjm
209	ig0=1+(j-2)*iim
210	do i=1,iim
211	ztemp(j,l) = ztemp(j,l) + ptemp(ig0+i,l)/iim
212	zplay(j,l) = zplay(j,l) + pplay(ig0+i,l)/iim
213	zzlay(j,l) = zzlay(j,l) + pzlay(ig0+i,l)/iim
214	enddo
215	enddo
216	ztemp(jjm+1,l) = ptemp(klon,l)
217	zplay(jjm+1,l) = pplay(klon,l)
218	zzlay(jjm+1,l) = pzlay(klon,l)
219	temp_eq = ztemp((jjm+1)/2,:)
220	press_eq = zplay((jjm+1)/2,:)/100. ! en mbar
221	enddo
222
223	ENDIF ! chimi or microfi=1
224
225	c -----------------------------
226	c Gestion des variables de la microphysique :
227	c
228	c -------------------
229	if (microfi.ge.1) then
230
231	c Traceurs microphysiques: passage en extensif: n/kg --> n/m^2 (2D ou 3D passage obligatoire)
232	DO iq=1,nmicro
233	c print*,tname(iq)
234	DO l=1,llm
235	DO i = 1, klon
236	qaer(i,l,iq) = tr_seri(i,l,iq)*delp(i,l)/RG
237	ENDDO
238	ENDDO
239	ENDDO
240	c copie du tableau de traceur :
241	qaer0(:,:,:)=qaer(:,:,:)
242	c
243	c -------------------
244	c Extraction des gaz pour les nuages
245	c
246	c recuperation des indices des gaz qui nous interesse
247	if (firstcall) then
248	if (clouds.eq.1) then
249	icldch4=-1
250	icldc2h6=-1
251	icldc2h2=-1
252	do i=1,nqmax
253	if (tname(i).eq."CH4") then
254	icldch4=i
255	c ich4=i
256	elseif (tname(i).eq."C2H6") then
257	icldc2h6=i
258	elseif (tname(i).eq."C2H2") then
259	icldc2h2=i
260	endif
261	enddo
262	if (icldch4 .eq.-1 .or.
263	& icldc2h6.eq.-1 .or.
264	& icldc2h2.eq.-1 ) then
265	print*, "Sacrebleu !!!"
266	print*, "Vous voulez faire des nuages sans gaz."
267	print*, "Mais vous etes inconscient. Je vais m'arreter la"
268	print*, "pour vous laisser reflechir au probleme"
269	STOP
270	endif
271	endif ! clouds=1
272	endif ! firstcall
273
274	c Saturation et fraction molaire CLOUD
275	c Calcul des saturations pour les esp chimique de la muphy des nuages.
276	c On le fait ici pour les sortir dans physiq.F sans avoir a surcharger la routine.
277	c Elles passent ensuite dans un common pour passer dans les I/O.
278	c
279	c-------------------------------------------
280	IF (clouds.eq.1) THEN
281	DO l=1,llm
282	DO i = 1, klon
283	call ch4sat(ptemp(i,l),pplay(i,l),tmp) !tmp en kg/kg !
284	satch4(i,l) = tr_seri(i,l,icldch4)/(tmp*28./16.)
285
286	call c2h6sat(ptemp(i,l),pplay(i,l),tmp)
287	satc2h6(i,l) =tr_seri(i,l,icldc2h6)/(tmp*28./30.)
288
289	call c2h2sat(ptemp(i,l),pplay(i,l),tmp)
290	satc2h2(i,l) =tr_seri(i,l,icldc2h2)/(tmp*28./26.)
291
292	ENDDO
293	ENDDO
294
295	c Copie des gaz (en 3D) <== UNIQUEMENT SI ON FAIT DES NUAGES
296	gaz1(:,:) = tr_seri(:,:,icldch4)
297	gaz2(:,:) = tr_seri(:,:,icldc2h6)
298	gaz3(:,:) = tr_seri(:,:,icldc2h2)
299
300	ENDIF ! clouds=1
301
302	endif ! microfi.ge.1
303
304	c -------------------
305	c Si microfi = 1 on est en 2D :
306	c conversion des inputs de muphys
307	IF (microfi.eq.1) THEN
308
309	zmu0(1) = pmu0(1)
310	zfract(1) = pfract(1)
311	do j=2,jjm
312	ig0=1+(j-2)*iim
313	do i=1,iim
314	zmu0(j) = zmu0(j) + pmu0(ig0+i)/iim
315	zfract(j) = zfract(j) + pfract(ig0+i)/iim
316	enddo
317	enddo
318	zmu0(jjm+1) = pmu0(klon)
319	zfract(jjm+1) = pfract(klon)
320	c
321	c traceurs 3D --> 2D
322	c
323	do iq=1,nqmax
324	do l=1,llm
325	zqaer(1,l,iq) = qaer(1,l,iq)
326	do j=2,jjm
327	ig0=1+(j-2)*iim
328	do i=1,iim
329	zqaer(j,l,iq) = zqaer(j,l,iq) + qaer(ig0+i,l,iq)/iim
330	enddo
331	enddo
332	zqaer(jjm+1,l,iq) = qaer(klon,l,iq)
333	enddo
334	enddo
335	c copie du tableau de traceur
336	zqaer0(:,:,:) = zqaer(:,:,:)
337	c
338	c gaz 3D --> 2D <=== UNIQUEMENT SI ON FAIT DES NUAGES.
339	c
340	if (clouds.eq.1) then
341	do l=1,llm
342	zgaz1(1,l) = gaz1(1,l)
343	zgaz2(1,l) = gaz2(1,l)
344	zgaz3(1,l) = gaz3(1,l)
345	do j=2,jjm
346	ig0=1+(j-2)*iim
347	do i=1,iim
348	zgaz1(j,l) = zgaz1(j,l) + gaz1(ig0+i,l)/iim
349	zgaz2(j,l) = zgaz2(j,l) + gaz2(ig0+i,l)/iim
350	zgaz3(j,l) = zgaz3(j,l) + gaz3(ig0+i,l)/iim
351	enddo
352	enddo
353	zgaz1(jjm+1,l) = gaz1(klon,l)
354	zgaz2(jjm+1,l) = gaz2(klon,l)
355	zgaz3(jjm+1,l) = gaz3(klon,l)
356	enddo
357
358	zgaz10=zgaz1
359	zgaz20=zgaz2
360	zgaz30=zgaz3
361	endif ! clouds=1
362
363	endif ! microfi=1
364
365	c AUTRES TRACEURS
366
367	if (nqmax.gt.nmicro) then
368	do iq=nmicro+1,nqmax
369	do l=1,llm
370	ychim(1,l,iq-nmicro) = tr_seri(1,l,iq)
371	do j=2,jjm
372	ig0=1+(j-2)*iim
373	do i=1,iim
374	ychim(j,l,iq-nmicro) = ychim(j,l,iq-nmicro)
375	. + tr_seri(ig0+i,l,iq)/iim
376	enddo
377	enddo
378	ychim(jjm+1,l,iq-nmicro) = tr_seri(klon,l,iq)
379	enddo
380	nomqy(iq-nmicro) = tname(iq)
381	c print*,iq-nmicro,nomqy(iq-nmicro)
382	enddo
383	nomqy(nqmax-nmicro+1) = "HV"
384	endif
385
386	c-----------------------------------------------------------------------
387	c initialisation des qysat au premier appel:
388	c ---------------------------------
389
390	c!! ATTENTION, qysat pris uniquement a l'equateur
391	c!! justifie puisque dans cette region, les var de t et p sont faibles...
392
393	if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
394	call inicondens(nqmax-nmicro,press_eq,temp_eq,nomqy,qysat)
395	endif
396
397	c-----------------------------------------------------------------------
398	c Appel de la microphysique en 2D/3D !!!!!!
399	c --------------------------
400
401	IF(firstcall) THEN
402	print*,'MICROPHYSIQUE ',MICROFI
403	ENDIF
404
405	c call begintime(tt0)
406	IF (MICROFI.eq.0) THEN
407	c PAS DE MICROPHYSIQUE :
408	c On appelle juste rdf pour creer la grille de rayons.
409	IF (firstcall) THEN
410	print*,'MICROPHYSIQUE OFF-LINE',MICROFI
411	call rdf()
412	ENDIF
413	c NOTES :
414	c L'appel de rdf ne sert a rien ici mis a part pour le TR. Si cet
415	c appel a deja lieu dans le TR inutile de le refaire ici.
416	c Je ne sais pas exactement comment marche les modules en F90
417	c Mais je recopie les valeurs du common/part/ de rdf pour
418	c les mettre dans un common interne a la microphysique (voir varmuphy.h)
419	c DONC J'AI BESOIN D'AVOIR ACCES A L'ANCIEN COMMON !!!
420	c
421	ELSEIF (MICROFI.eq.1) THEN
422	c MICROPHYSIQUE 2D :
423	c Les input/output comportent le prefixe z pour 2D :)
424	zdqmufi = 0. ! ne sert que pour chimi pour condensation
425	call muphys(jjm+1,
426	& zplev,zplay,zzlev,zzlay,
427	& ztemp,zqaer,zgaz1,zgaz2,zgaz3,
428	& nmicro,ptimestep,
429	& zmu0,zfract,
430	c -------- sorties diagnostiques
431	& zflxesp_i,
432	& ztau_drop,ztau_aer,
433	& zsolesp,zprec)
434	ELSE
435	c MICROPHYSIQUE 3D :
436	c Les input sont des champs 3D directement !
437	call muphys(klon,
438	& pplev,pplay,pzlev,pzlay,
439	& ptemp,qaer,gaz1,gaz2,gaz3,
440	& nmicro,ptimestep,
441	& pmu0,pfract,
442	c ------ sorties diagnostiques
443	& flxesp_i,
444	& tau_drop,tau_aer,
445	& solesp,prec)
446	c
447	c NOTES :
448	c Ici toutes nos sorties sont des champs 3D...(meme les diagnostiques)
449	c On a rien a faire mis a part copier les dq dans les d_tr
450	c
451	ENDIF
452	c call endtime(tt0,tt1)
453	c ttmuphys=ttmuphys+tt1
454
455	c-----------------------------------------------------------------------
456	c Mise a jour des sorties de muphys
457	c -------------
458	c En 2D on copie les sorties de muphys de la grille LATxALT
459	c sur la grille complete.
460	IF (microfi.eq.1) THEN
461	c precipitations
462	DO l=1,5
463	prec(1,l) = zprec(1,l)
464	ig0 = 2
465	DO j=2,jjm
466	DO i = 1, iim
467	prec(ig0,l) = zprec(j,l)
468	ig0 = ig0 + 1
469	ENDDO
470	ENDDO
471	prec(ig0,l) = zprec(jjm+1,l)
472	ENDDO
473	c taux sedimentation
474	DO l=1,llm
475	c taux sed goutte
476	IF (clouds.eq.1) THEN
477	tau_drop(1,l) = ztau_drop(1,l)
478	ig0 = 2
479	DO j=2,jjm
480	DO i = 1, iim
481	tau_drop(ig0,l) = ztau_drop(j,l)
482	ig0 = ig0 + 1
483	ENDDO
484	ENDDO
485	tau_drop(ig0,l) = ztau_drop(jjm+1,l)
486	ENDIF
487	c taux sed aer
488	DO iq=1,nrad
489	tau_aer(1,l,iq) = ztau_aer(1,l,iq)
490	ig0 = 2
491	DO j=2,jjm
492	DO i = 1, iim
493	tau_aer(ig0,l,iq) = ztau_aer(j,l,iq)
494	ig0 = ig0 + 1
495	ENDDO
496	ENDDO
497	tau_aer(ig0,l,iq) = ztau_aer(jjm+1,l,iq)
498	ENDDO
499	ENDDO
500	c flux glace / production glace
501	IF (clouds.eq.1) THEN
502	DO iq=1,3
503	DO l=1,llm+1
504	if (l.le.llm) flxesp_i(1,l,iq) = zflxesp_i(1,l,iq)
505	solesp(1,l,iq) = zsolesp(1,l,iq)
506	ig0 = 2
507	DO j=2,jjm
508	DO i = 1, iim
509	if (l.le.llm) flxesp_i(ig0,l,iq)=zflxesp_i(1,l,iq)
510	solesp(ig0,l,iq) = zsolesp(1,l,iq)
511	ig0 = ig0 + 1
512	ENDDO
513	ENDDO
514	if (l.le.llm) flxesp_i(ig0,l,iq)=zflxesp_i(jjm+1,l,iq)
515	solesp(ig0,l,iq) = zsolesp(jjm+1,l,iq)
516	ENDDO
517	ENDDO
518	ENDIF
519	ENDIF
520
521	c-----------------------------------------------------------------------
522	c Gestion des sources
523	c -------------
524	c
525	IF (clouds.eq.1) THEN
526	IF (microfi.eq.1) THEN
527	c On repasse les gaz en 3D si on a fait de la microphysique en 2D
528	DO l=1,llm
529	gaz1(1,l) = zgaz1(1,l)
530	gaz2(1,l) = zgaz2(1,l)
531	gaz3(1,l) = zgaz3(1,l)
532	ig0 = 2
533	DO j=2,jjm
534	DO i = 1, iim
535	gaz1(ig0,l) = zgaz1(j,l)* gaz1(ig0,l) /zgaz10(j,l)
536	gaz2(ig0,l) = zgaz2(j,l)* gaz2(ig0,l) /zgaz20(j,l)
537	gaz3(ig0,l) = zgaz3(j,l)* gaz3(ig0,l) /zgaz30(j,l)
538	ig0 = ig0 + 1
539	ENDDO
540	ENDDO
541	gaz1(ig0,l) = zgaz1(jjm+1,l)
542	gaz2(ig0,l) = zgaz2(jjm+1,l)
543	gaz3(ig0,l) = zgaz3(jjm+1,l)
544	ENDDO
545	ENDIF
546	c Mise a jour du reservoir de CH4 (ie : seul le CH4 remplit le reservoir)
547	DO i=1,klon
548	reservoir(i) = reservoir(i)+prec(i,1)
549	ENDDO
550	c Calcul des sources :
551	c ch4=0.
552	c ch4(1) = gaz1(1,1)
553	c do j=2,jjm
554	c ig0=1+(j-2)*iim
555	c do i=1,iim
556	c ch4(j)= ch4(j) + gaz1(ig0+i,1)/iim
557	c enddo
558	c enddo
559	c ch4(jjm+1) = gaz1(ig0,1)
560
561	CALL sources(klon,klev,ptimestep,z0,
562	& pu,pv,pplev,pzlay,pzlev,
563	& gaz1,gaz2,gaz3,
564	& ftsol,solesp,reservoir)
565
566	c ch4b=0.
567	c ch4b(1) = gaz1(1,1)
568	c do j=2,jjm
569	c ig0=1+(j-2)*iim
570	c do i=1,iim
571	c ch4b(j)= ch4b(j) + gaz1(ig0+i,1)/iim
572	c enddo
573	c enddo
574	c ch4b(jjm+1) = gaz1(ig0,1)
575	c do j=1,jjm+1
576	c write(499,*) j,ch4(j),ch4b(j)
577	c enddo
578	c write(499,*) ""
579	ENDIF
580	c-----------------------------------------------------------------------
581	c Condensation
582	c -------------
583
584	IF ((chimi).and.(nqmax.gt.nmicro)) then
585
586	c tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
587	c print*,'Condensation'
588
589	do iq=1,nqmax-nmicro
590	do l=1,llm
591	do j=1,jjm+1
592	if (ychim(j,l,iq).gt.qysat(l,iq)) then
593	zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
594	. / ptimestep ! / dt
595	endif
596	enddo
597	enddo
598	enddo
599
600	ENDIF
601
602	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
603	c eventuellement, modif initiale de la compo
604	c
605	c tendance (en /s) passee sur zdqmufi(nmicro+1 a nqmax)
606	c
607	c if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
608	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
609	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
610	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
611	c do iq=1,nqmax-nmicro
612	c if (nomqy(iq).eq."CH4") then
613	c do l=1,llm
614	c do j=1,jjm+1
615	c if (ychim(j,l,iq).le.0.015) then
616	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
617	c . / ptimestep ! / dt
618	c endif
619	c enddo
620	c enddo
621	c endif
622	c enddo
623	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
624	c
625	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
626	c!!!remise de C2H2 a 1.e-5 max !!!!!!!!!!!!!
627	c!!!remise de C2H6 a 3.e-5 max !!!!!!!!!!!!!
628	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
629	c do iq=1,nqmax-nmicro
630	c if (nomqy(iq).eq."C2H2") then
631	c do l=1,llm
632	c do j=1,jjm+1
633	c if (ychim(j,l,iq).gt.1.e-5) then
634	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
635	c . / ptimestep ! / dt
636	c endif
637	c enddo
638	c enddo
639	c endif
640	c if (nomqy(iq).eq."C2H6") then
641	c do l=1,llm
642	c do j=1,jjm+1
643	c if (ychim(j,l,iq).gt.3.e-5) then
644	c zdqmufi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
645	c . / ptimestep ! / dt
646	c endif
647	c enddo
648	c enddo
649	c endif
650	c enddo
651	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
652	c endif
653	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
654
655	c ----- commentaire de fin (mise a jour des profil de fraction molaire)
656
657	c-----------------------------------------------------------------------
658	c Appel de la chimie
659	c --------------------------
660
661	if((appkim.eq.1).and.(chimi)) then
662	print*,'On passe dans la CHIMIE'
663
664	c do iq=1,nqmax-nmicro
665	c if (nomqy(iq).eq."C2H2") then
666	c print*,"C2H2top=",ychim(:,klev,iq)
667	c endif
668	c enddo
669
670	c Appel Chimie
671	c ------------
672	CALL calchim(nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
673	. ztemp,zplay,zplev,
674	. pdyfi)
675	c ychim ne doit pas etre modifie, pdyfi en /s
676
677	endif
678
679	c-----------------------------------------------------------------------
680	c retour des tendances vers 3D
681	c ---------------------------------
682
683	c TRACEURS MICROPHYSIQUES
684	c
685	c ---> pas de microphysique
686	IF (microfi.eq.0) THEN
687	DO iq=1,nmicro
688	d_tr_mph(:,:,iq)=0.
689	ENDDO
690	ENDIF
691	c ---> microphysique 2D
692	IF (microfi.eq.1) THEN
693	c on repasse le champ de traceurs en 3D (pas les tendances)
694	DO iq=1,nmicro
695	DO l=1,llm
696	qaer(1,l,iq) = zqaer(1,l,iq)
697	ig0 = 2
698	DO j=2,jjm
699	DO i = 1, iim
700	c un petit patch :
701	c Si la moyenne zonale au depart est "nulle" :
702	c On a quand meme le droit de produire des traceurs dans la cellule.
703	c On considere donc que la valeur de sortie 3D correspond a la valeur de sortie 2D.
704	c Cela permet aussi entre autre d'eviter les NaN pour les traceurs des nuages !
705	c (au dessus de la tropo pas de nuages donc qaer(nrad+1:ntype*nrad) = 0 !!!)
706	IF (zqaer0(j,l,iq).lt.1e-100) THEN
707	qaer(ig0,l,iq) = zqaer(j,l,iq)
708	ELSE
709	qaer(ig0,l,iq) = zqaer(j,l,iq) *
710	& qaer0(ig0,l,iq)/zqaer0(j,l,iq)
711	ENDIF
712	ig0 = ig0 + 1
713	ENDDO
714	ENDDO
715	qaer(ig0,l,iq) = zqaer(jjm+1,l,iq)
716	ENDDO
717	ENDDO
718	c La tendances correspond a (qaer-qaer0)/ptimestep
719	DO iq=1,nmicro
720	DO i=1,klon
721	DO l=1,llm
722	d_tr_mph(i,l,iq) = (qaer(i,l,iq)-qaer0(i,l,iq))/
723	& ptimestep
724	ENDDO
725	ENDDO
726	ENDDO
727	c ---> microphysique 3D
728	ELSEIF(microfi.gt.1) THEN
729	DO iq=1,nmicro
730	DO l=1,llm
731	DO i = 1, klon
732	d_tr_mph(i,l,iq)=(qaer(i,l,iq)-qaer0(i,l,iq))/ptimestep
733	ENDDO
734	ENDDO
735	ENDDO
736
737	do iq=1,nmicro
738	DO l=1,llm
739	DO i = 1, klon
740	c Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
741	d_tr_mph(i,l,iq) = d_tr_mph(i,l,iq)*RG/delp(i,l)
742	ENDDO
743	ENDDO
744	enddo
745
746	ENDIF ! microfi
747
748	c AUTRES TRACEURS
749
750	if ((chimi).and.(nqmax.gt.nmicro)) then
751	c on passe de pdyfi (tendance chimique en /s calculee quand chimie appelee)
752	c a d_tr_kim (tendance chimique 3D en /s, passee a physiq)
753	c et de zdqmufi a d_tr_mph (tendance condensation 3D en /s passee a physiq)
754
755	DO iq=nmicro+1,nqmax
756	DO l=1,llm
757	d_tr_kim(1,l,iq) = pdyfi(1,l,iq-nmicro)
758	d_tr_mph(1,l,iq) = zdqmufi(1,l,iq)
759	ig0 = 2
760	DO j=2,jjm
761	DO i = 1, iim
762	d_tr_kim(ig0,l,iq) = pdyfi(j,l,iq-nmicro)
763	& *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
764	d_tr_mph(ig0,l,iq) = zdqmufi(j,l,iq)
765	& *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
766	ig0 = ig0 + 1
767	ENDDO
768	ENDDO
769	d_tr_kim(ig0,l,iq) = pdyfi(jjm+1,l,iq-nmicro)
770	d_tr_mph(ig0,l,iq) = zdqmufi(jjm+1,l,iq)
771	ENDDO
772	ENDDO
773
774	endif ! chimi
775
776	c--------------------------------------------------
777	c CONDENSATION VIA MICROFI
778	c----------------------
779	c La microphysique avec nuages doit se faire obligatoirement en 3D.
780	c Rien n'empeche de faire la chimie en 2D. Cependant pour prendre en compte la
781	c condensation due a la microfi (en 3D) on recalcule la tendance finale pour
782	c les especes concernees (CH4, C2H6 pour le moment).
783	IF (microfi.ge.1.and.clouds.eq.1) THEN
784	DO i=1,klon
785	DO l=1,klev
786	c condensation CH4
787	d_tr_mph(i,l,icldch4)=(gaz1(i,l)-tr_seri(i,l,icldch4))
788	& /ptimestep
789	c condensation C2H6
790	d_tr_mph(i,l,icldc2h6)=(gaz2(i,l)-tr_seri(i,l,icldc2h6))
791	& /ptimestep
792	c condensation C2H2
793	d_tr_mph(i,l,icldc2h2)=(gaz3(i,l)-tr_seri(i,l,icldc2h2))
794	& /ptimestep
795	ENDDO
796	ENDDO
797	ENDIF
798	c ch4c=0.
799	c do l=1,llm
800	c ch4c(1,l) = tr_seri(1,l,icldch4)
801	c do j=2,jjm
802	c ig0=1+(j-2)*iim
803	c do i=1,iim
804	c ch4c(j,l)= ch4c(j,l)+tr_seri(ig0+i,l,icldch4)/iim
805	c enddo
806	c enddo
807	c ch4c(jjm+1,l) = tr_seri(klon,l,icldch4)
808	c enddo
809	c do l=1,llm
810	c write(500,*) pplay(25,l),ch4c(25,l)
811	c enddo
812	c write(500,*) ""
813
814
815	c--------------------------------------------------
816	c MISE A JOUR CH4 : (pour refixer la fraction
817	c molaire)
818	c--------------------------------------------------
819	c IF (firstcall) THEN
820	c do i=1,klon
821	c do j=1,llm
822	c call ch4sat(ptemp(i,j),pplay(i,j),tmp) !tmp en kg/kg !
823	c tmp=0.950.85tmp*28./16.
824	c if (pplay(i,j).lt.20000.) then
825	c dqch4 = 1.4e-2
826	c else
827	c dqch4 = tmp
828	c endif
829	c d_tr_mph(i,j,icldch4)=(-tr_seri(i,j,icldch4)+dqch4)/
830	c & ptimestep
831	c enddo
832	c enddo
833	c
834	c ENDIF
835
836	c--------------------------------------------------
837	c CALCUL DU FLUX DE CHALEUR LATENTE D'EVAPORATION
838	c DU METHANE
839	c--------------------------------------------------
840	IF (clouds.eq.1) THEN
841	DO i=1,klon
842	fte= (1.-ftsol(i)/305.5)
843	ftm= (1.-ftsol(i)/190.5)
844	if(ftm.le.1.e-3) ftm=1.e-3
845	if(fte.le.1.e-3) fte=1.e-3
846	Lvch4 =8.314190.4
847	& (7.08ftm0.354+10.951.1e-2ftm*0.456)
848	& /mch4
849	! solcxhy en m3/m2 {ok}
850	! 425 en kg/m3
851	! Lv en J/kg {ok}
852	! ptimestep en s {ok}
853	fclat(i)=(solesp(i,klev+1,1)Lvch4rhoi_ch4) ! en J/m2/s
854	ENDDO
855	ENDIF
856
857	c--------------------------------------------------
858	c GESTION DES RAYONS DE GOUTTES POUR TR
859	c--------------------------------------------------
860	IF (clouds.eq.1) THEN
861
862	c Calcul du rayon des gouttes par bin ...
863	c----------------------------------------
864	DO i=1,klon
865	DO j=1,klev
866	DO iq=1,nrad
867	* Rayon minimum selon la quantité de noyaux
868	IF (qaer(i,j,iq+nrad) .le. 1.e-5) THEN
869	rcloud(i,j,iq) = 1.e-10
870	ELSE
871	rcloud(i,j,iq)=
872	& ((qaer(i,j,iq+2*nrad)/qaer(i,j,iq+nrad)+
873	& qaer(i,j,iq+3*nrad)/qaer(i,j,iq+nrad) +
874	& v_e(iq))0.75/3.1415926353)*(1./3.)
875	ENDIF
876	ENDDO
877	ENDDO
878	ENDDO
879
880	c .... et de leur rayon moyen total (tt bins confondu)
881	c------------------------------------------------------
882	DO i=1,klon
883	socccld=0.
884	DO j=klev,1,-1 !de haut en bas pour le calcul des opacites
885	vcl=0.
886	nuc=0.
887	xgsn=0.
888	xmsn=0.
889	xesn=0.
890	xasn=0.
891	DO iq=1,nrad
892	vcl=vcl+qaer(i,j,iq+2*nrad)+
893	& qaer(i,j,iq+3*nrad)+
894	& qaer(i,j,iq+4*nrad)+
895	& v_e(iq)*qaer(i,j,iq+nrad) ! volume des gouttes
896	nuc=nuc+qaer(i,j,iq+nrad) ! nombre de noyaux
897	xgsn=xgsn+qaer(i,j,iq+nrad)*v_e(iq) ! volume de noyaux
898	xmsn=xmsn+qaer(i,j,iq+2*nrad) ! volume de methane
899	xesn=xesn+qaer(i,j,iq+3*nrad) ! volume d' ethane
900	xasn=xasn+qaer(i,j,iq+4*nrad) ! volume d' acethylene
901	ENDDO
902	IF (nuc .le. 1.e-5) THEN
903	rmcloud(i,j)=1.0e-10
904	xfrac(i,j,:)=0.
905	ELSE
906	IF(xgsn/vcl.lt.0. .or. xgsn/vcl.gt.1.001)
907	& print*, 'PB AVEC XFRAC:', i,j,xgsn,vcl
908	rmcloud(i,j)= ! rayon moyen des gouttes
909	& (vcl/nuc0.75/3.1415926353)*(1./3.)
910	xfrac(i,j,1)=xgsn/vcl ! fraction volumique noyau/goutte
911	xfrac(i,j,2)=xmsn/vcl ! fraction volumique CH4/goutte
912	xfrac(i,j,3)=xesn/vcl ! fraction volumique C2H6/goutte
913	xfrac(i,j,4)=xasn/vcl ! fraction volumique C2H2/goutte
914	c calcul du rayon moyen (moyenne temporelle)
915	rmcbar(i,j)=rmcbar(i,j)+rmcloud(i,j)
916	xfbar(i,j,:)=xfbar(i,j,:)+xfrac(i,j,:)
917	ncount(i,j) = ncount(i,j)+1
918	ENDIF
919	socccld=socccld+3.1415926(rmcloud(i,j)2.)nuc
920	occcld(i,j)=socccld
921	ENDDO
922	ENDDO
923	c
924	c OCCCLD
925	c Calcul le nombre d'occurence d'un nuage
926	c d'opacité comprise en kmin et kmax
927	c k kmin kmax
928	c 1 0.0000000 0.10000000
929	c 2 0.10000000 0.17782794
930	c 3 0.17782794 0.31622776
931	c 4 0.31622776 0.56234139
932	c 5 0.56234139 1.0000000
933	c 6 1.0000000 1.7782795
934	c 7 1.7782795 3.1622777
935	c 8 3.1622777 5.6234136
936	c 9 5.6234136 10.000000
937	c 10 10.000000 17.782795
938	c 11 17.782795 31.622778
939	c 12 31.622778 100000.00
940	c
941	DO i=1,klon
942	DO j=1,klev
943	DO k=1,12
944	ex=10.**(0.25)
945	kmin=0.
946	kmax=1.e5
947	if(k.ne.1) kmin=0.1ex*(k-2)
948	if(k.ne.12) kmax=0.1ex*(k-1)
949	if(occcld(i,j).ge.kmin .and. occcld(i,j).lt.kmax)
950	& occcld_m(i,j,k)=1.
951	ENDDO
952	ENDDO
953	ENDDO
954	ENDIF ! fin condition clouds => pas besoin de calculer des rayons
955
956	RETURN
957	END

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

Download in other formats:

Original Format