Context Navigation

phytrac.F @ 398

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

SLebonnois: correction de bugs dans la physique Titan:

effg.F : Z doit etre en km, donc conversion
optc*_1pt_2.F : On utilise cfffv11 et plus optfrac Du coup, les fichiers input testag* ne sont plus necessaires.
phytrac.F : passage de la tendance aerosols en intensif dans tous les cas

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

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

Download in other formats:

Original Format