Context Navigation

phytrac.F @ 1000

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

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

Download in other formats:

Original Format