Context Navigation

calchim.dbleprec @ 21

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

Creation de repertoires:

chantiers : pour communiquer sur nos projets de modifs
documentation : pour stocker les docs

Ajout de:

libf/phytitan : physique de Titan
libf/chimtitan: chimie de Titan
libf/phyvenus : physique de Venus

File size: 16.9 KB

Line
1	SUBROUTINE calchim(ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim,
2	. ctemp,cplay,cplev,
3	.
4	.
5	. dqyc)
6	implicit none
7
8	c-------------------------------------------------
9	c
10	c Introduction d'une routine chimique
11	c
12	c Auteur: S. Lebonnois, 01/2000 \| 09/2003
13	c adaptation pour Titan 3D: 02/2009
14	c
15	c-------------------------------------------------
16	c
17
18	#include "dimensions.h"
19	#include "dimphy.h"
20	#include "clesphys.h"
21	#include "paramet.h"
22	#include "YOMCST.h"
23
24	#include "titan_for.h"
25	!!! doit etre en accord avec titan.h
26	#include "aerprod.h"
27
28	c Arguments
29	c ---------
30
31	INTEGER ny ! nb de composes (nqmax-nmicro)
32	REAL qy_c(jjm+1,klev,NC) ! Especes chimiques apres adv.+diss.
33	character*10 nomqy_c(NC+1) ! Noms des especes chimiques
34	REAL declin_rad,ls_rad ! declinaison et long solaire en radians
35	REAL dtchim ! pas de temps chimie
36	REAL ctemp(jjm+1,klev) ! Temperature
37	REAL cplay(jjm+1,klev) ! pression (Pa)
38	REAL cplev(jjm+1,klev) ! pression intercouches (Pa)
39
40	REAL dqyc(jjm+1,klev,nqmx) ! Tendances especes chimiques (nqmx, mais en fait NC...)
41
42	c Local variables :
43	c -----------------
44	c variables envoyees dans la chimie: double precision
45
46	integer i,j,l,ic
47	REAL temp_c(klev),press_c(klev) ! T,p(mbar) a 1 lat donnee
48	REAL declin_c ! declinaison en degres
49	REAL cqy(klev,NC) ! frac mol qui seront modifiees
50	REAL surfhaze(klev)
51	REAL cprodaer(klev,4),cmaer(klev,4),ccsn(klev,4),ccsh(klev,4)
52	REAL rinter(klev),nb(klev)
53	REAL a,b
54	character str11,str22
55	integer ntotftop
56	parameter (ntotftop=50)
57	integer nftop(ntotftop),isaisonflux
58	REAL fluxtop(NC),latit,tabletmp(ntotftop),factflux
59	character*10 nomsftop(ntotftop+1)
60	character*20 formt1,formt2,formt3
61
62	c variables locales initialisees au premier appel
63
64	LOGICAL firstcal
65	DATA firstcal/.true./
66	SAVE firstcal
67
68	REAL mass(NC),duree
69	REAL tablefluxtop(NC,jjm+1,5)
70	REAL botCH4
71	DATA botCH4/0.05/
72	REAL krpd(15,ND+1,klev,jjm+1),krate(klev,NR)
73	integer reactif(5,NR),nom_prod(NC),nom_perte(NC)
74	integer prod(200,NC),perte(2,200,NC)
75	SAVE mass,tablefluxtop,botCH4
76	SAVE krpd,krate
77	SAVE reactif,nom_prod,nom_perte,prod,perte
78
79	c-----------------------------------------------------------------------
80	c***********************************************************************
81	c
82	c Initialisations :
83	c ----------------
84
85	duree = dtchim ! passage en real*4 pour appel a routines C
86
87	IF (firstcal) THEN
88
89	print*,'CHIMIE, premier appel'
90
91	c ************************************
92	c Au premier appel, initialisation de toutes les variables
93	c pour les routines de la chimie.
94	c ************************************
95
96	c Verification dimension verticale: coherence titan_for.h et klev
97	c --------------------------------
98
99	if (klev.ne.NLEV) then
100	print*,'PROBLEME de coherence dans la dimension verticale:'
101	. ,klev,NLEV
102	STOP
103	endif
104
105	c Verification du nombre de composes: coherence titan_for.h et nqmax-nmicro
106	c ----------------------------------
107
108	if (ny.ne.NC) then
109	print*,'PROBLEME de coherence dans le nombre de composes:'
110	. ,ny,NC
111	STOP
112	endif
113
114	c calcul de temp_c, densites et press_c a l'equateur:
115	c --------------------------------------------------
116
117	print*,'pression, densites et temp a l equateur (chimie):'
118	print*,'level, press_c, nb, temp_c'
119	DO l=1,klev
120	c temp_c (K):
121	temp_c(l) = ctemp(jjm/2+1,l)
122	c press_c (mbar):
123	press_c(l) = cplay(jjm/2+1,l)/100.
124	c nb (cm-3):
125	nb(l) = 1.e-4press_c(l) / (RKBOLtemp_c(l))
126	print*,l,press_c(l),nb(l),temp_c(l)
127	ENDDO
128
129	c caracteristiques des composes:
130	c -----------------------------
131	mass = 0.0
132	call comp(nomqy_c,mass)
133	print*,' Mass'
134	do ic=1,NC
135	print*,nomqy_c(ic),mass(ic)
136	enddo
137
138
139	c FLUX AU SOMMET: DEPENDANT DE LA LATITUDE ET DE LA SAISON...
140
141	c flux dans la derniere couche:(issu du modele 1D, a 500 km)
142	c -----------------------------
143	c !! en cm-2.s-1 !!
144	c ces valeurs sont converties en cm-3.s-1 (dni/dt) dans a couche 54
145	c
146	c Lecture des tables de flux en fonction lat et saison
147	c ----------------------------------------------------
148
149	WRITE(str2(1:2),'(i2.2)') ntotftop
150	formt1 = '(A10,'//str2//'(3X,A10))'
151	formt2 = '(F12.6,'//str2//'(2X,F13.6))'
152	formt3 = '(F13.6,'//str2//'(2X,F13.6))'
153
154	do j=1,jjp1
155	do ic=1,NC
156	do l=1,5
157	tablefluxtop(ic,j,l) = 0.
158	enddo
159	enddo
160	enddo
161
162	do l=1,4
163	WRITE(str1(1:1),'(i1.1)') l
164	c hx1 -> Ls=RPI/2 / hx4 -> Ls=0
165	c open(10,file="../../INPUT/FLUXTOP/flux500.hs"//str1)
166	open(10,file="FLUXTOP/flux500.hs"//str1)
167	c on remet l=1 et 5 -> Ls=0 / l=2 -> Ls=RPI/2 ...
168
169	read(10,formt1) nomsftop
170	do j=1,ntotftop
171	do i=1,10 !justification a gauche...
172	if (nomsftop(j+1)(i:i).ne.' ') then
173	nomsftop(j) = nomsftop(j+1)(i:)
174	goto 100
175	endif
176	enddo
177	100 continue
178	c print*,nomsftop(j)
179	nftop(j) = 0
180	do i=1,NC
181	if (nomqy_c(i).eq.nomsftop(j)) then
182	nftop(j) = i
183	endif
184	enddo
185	c if (nftop(j).ne.0) print*,nomsftop(j),nomqy_c(nftop(j))
186	enddo
187
188
189	c lecture des flux. Les formats sont un peu alambiques...
190	c a revoir eventuellement
191	do j=1,jjp1/2+1
192	read(10,formt2) latit,tabletmp
193	do i=1,ntotftop
194	if (nftop(i).ne.0) then
195	tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
196	endif
197	enddo
198	enddo
199	do j=jjp1/2+2,jjp1
200	read(10,formt3) latit,tabletmp
201	do i=1,ntotftop
202	if (nftop(i).ne.0) then
203	tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
204	endif
205	enddo
206	enddo
207
208	enddo ! l
209
210	do j=1,jjp1
211	do ic=1,NC
212	tablefluxtop(ic,j,1) = tablefluxtop(ic,j,5)
213	enddo
214	enddo
215
216	c Stockage des composes utilises dans la prod d'aerosols
217	c (aerprod=1) et dans H-> H2 (htoh2=1): utilaer
218	c ! decalage de 1 car utilise dans le c !
219	c ------------------------------------------
220	c + modifs du flux en fonction des obs UVIS et CIRS: C2H2,C2H6,HCN
221	c + modifs des flux qui presentent un probleme evident: C3H8 et C4H10
222	c ------------------------------------------
223
224	do ic=1,NC
225
226	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
227	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
228	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229	c if (nomqy_c(ic).eq."CH4") then
230	c do l=1,llm
231	c do j=1,ip1jmp1
232	c if (qy_c(j,l,ic).le.0.015) qy_c(j,l,ic) = 0.015
233	c enddo
234	c enddo
235	c endif
236	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
237
238	if (nomqy_c(ic).eq."C4H2") then
239	utilaer(10) = ic-1
240	endif
241	if (nomqy_c(ic).eq."HCN") then
242	utilaer(6) = ic-1
243	do j=1,jjp1
244	do i=1,5
245	tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 6.
246	enddo
247	enddo
248	endif
249	if (nomqy_c(ic).eq."HC3N") then
250	utilaer(7) = ic-1
251	endif
252	if (nomqy_c(ic).eq."NCCN") then
253	utilaer(14) = ic-1
254	endif
255	if (nomqy_c(ic).eq."CH3CN") then
256	utilaer(15) = ic-1
257	utilaer(16) = ic-1 ! si pas C2H3CN, CH3CN utilise, mais reac. nulle
258	endif
259	if (nomqy_c(ic).eq."H") then
260	utilaer(1) = ic-1
261	endif
262	if (nomqy_c(ic).eq."H2") then
263	utilaer(2) = ic-1
264	endif
265	if (nomqy_c(ic).eq."C2H2") then
266	utilaer(3) = ic-1
267	do j=1,jjp1
268	do i=1,5
269	c tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 7.3
270	tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 10.
271	enddo
272	enddo
273	endif
274	if (nomqy_c(ic).eq."AC6H6") then
275	utilaer(13) = ic-1
276	endif
277	if (nomqy_c(ic).eq."C2H3CN") then
278	utilaer(16) = ic-1
279	endif
280	if (nomqy_c(ic).eq."C2") then
281	utilaer(4) = ic-1
282	endif
283	if (nomqy_c(ic).eq."C2H") then
284	utilaer(5) = ic-1
285	endif
286	if (nomqy_c(ic).eq."C3N") then
287	utilaer(8) = ic-1
288	endif
289	if (nomqy_c(ic).eq."H2CN") then
290	utilaer(9) = ic-1
291	endif
292	if (nomqy_c(ic).eq."C4H3") then
293	utilaer(11) = ic-1
294	endif
295	if (nomqy_c(ic).eq."AC6H5") then
296	utilaer(12) = ic-1
297	endif
298
299	if (nomqy_c(ic).eq."C2H6") then
300	do j=1,jjp1
301	do i=1,5
302	c tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 640.
303	tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 1200.
304	enddo
305	enddo
306	endif
307	if (nomqy_c(ic).eq."C3H8") then
308	do j=1,jjp1
309	do i=1,5
310	c tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
311	tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-3000.)
312	enddo
313	enddo
314	endif
315	if (nomqy_c(ic).eq."C4H10") then
316	do j=1,jjp1
317	do i=1,5
318	tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
319	enddo
320	enddo
321	endif
322
323	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
324	c if ((nomqy_c(ic).eq."HC3N").or.
325	c $ (nomqy_c(ic).eq."C3N")) then
326	c DO j=1,ip1jmp1
327	c do l=1,34 ! p>~ 1 mbar
328	c qy_c(j,l,ic) = 1.e-30
329	c enddo
330	c ENDDO
331	c endif
332	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
333
334	enddo
335
336	c taux de photodissociations:
337	c --------------------------
338	call disso(krpd,jjp1)
339
340	c reactions chimiques:
341	c -------------------
342	call chimie(nomqy_c,nb,temp_c,krate,reactif,
343	. nom_perte,nom_prod,perte,prod)
344	c print*,'nom_prod, nom_perte:'
345	c do ic=1,NC
346	c print*,nom_prod(ic),nom_perte(ic)
347	c enddo
348	c print*,'premieres prod, perte(1:reaction,2:compagnon):'
349	c do ic=1,NC
350	c print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic)
351	c enddo
352
353	c l = klev-3
354	c print*,'krate a p=',press_c(l),' reactifs et produits:'
355	c do ic=1,ND+1
356	c print,ic,krpd(7,ic,l,4)nb(l)," ",
357	c . nomqy_c(reactif(1,ic)+1),
358	c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
359	c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
360	c enddo
361	c do ic=ND+2,NR
362	c print*,ic,krate(l,ic)," ",
363	c . nomqy_c(reactif(1,ic)+1),
364	c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
365	c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
366	c enddo
367
368	ENDIF ! premier appel
369
370	c***********************************************************************
371	c-----------------------------------------------------------------------
372
373	c calcul declin_c (en degres)
374	c ---------------------------
375
376	declin_c = declin_rad*180./RPI
377	c print*,'declinaison en degre=',declin_c
378
379	c-----------------------------------------------------------------------
380	c
381	c calcul du facteur d'interpolation entre deux saisons pour flux
382	c --------------------------------------------------------------
383
384	isaisonflux = int(ls_rad*2./RPI)+1
385	if (isaisonflux.eq.5) then
386	isaisonflux = 1
387	endif
388	factflux = (sin(ls_rad)-sin((isaisonflux-1)*RPI/2.))/
389	. (sin(isaisonfluxRPI/2.)-sin((isaisonflux-1)RPI/2.))
390
391	c***********************************************************************
392	c***********************************************************************
393	c
394	c BOUCLE SUR LES LATITUDES
395	c
396	DO j=1,jjp1
397
398	c***********************************************************************
399	c***********************************************************************
400
401	c-----------------------------------------------------------------------
402	c
403	c Temperature, pression (mbar), densite (cm-3)
404	c -------------------------------------------
405
406	DO l=1,klev
407	c temp_c (K):
408	temp_c(l) = ctemp(j,l)
409	c press_c (mbar):
410	press_c(l) = cplay(j,l)/100.
411	c nb (cm-3):
412	nb(l) = 1.e-4press_c(l) / (RKBOLtemp_c(l))
413	ENDDO
414
415	c-----------------------------------------------------------------------
416	c
417	c Distances radiales (intercouches, en km)
418	c ----------------------------------------
419
420	rinter(1) = RA/1000.
421	do l=2,klev
422	c A REVOIR: g doit varier avec r !
423	rinter(l) = rinter(l-1) +
424	. (cplev(j,l-1)-cplev(j,l))/cplay(j,l)(RDtemp_c(l)/RG)/1000.
425	c print*,rinter(l)
426	enddo
427
428	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
429	c FLUX AU TOP: interpolation en saison et
430	c conversion en cm-3 s-1 dans la couche klev-1 (NLEV-2 dans le C)
431	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
432	do ic=1,NC
433	fluxtop(ic) = tablefluxtop(ic,j,isaisonflux) *(1-factflux)
434	. + tablefluxtop(ic,j,isaisonflux+1)* factflux
435	fluxtop(ic) = fluxtop(ic)/((rinter(klev)-rinter(klev-1))*1e5)
436	enddo
437
438	c TEST: sortie de l'un des flux avec saveLs et factflux
439	c dans un fichier special pour tracer avec gnuplot
440	c if (j.eq.2) then
441	c open(unit=11,file="flux_80N.txt",status='old',position='append')
442	c write(11,) ls_rad180./RPI, factflux,
443	c . ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
444	c write(11,) ls_rad180./RPI, factflux,
445	c . fluxtop(utilaer(3)+1)((rinter(llm)-rinter(llm-1))1e5), !C2H2
446	c . fluxtop(utilaer(6)+1)((rinter(llm)-rinter(llm-1))1e5) !HCN
447	c close(11)
448	c endif
449	c if (j.eq.17) then
450	c open(unit=11,file="flux_eq.txt",status='old',position='append')
451	c write(11,) ls_rad180./RPI, factflux,
452	c . ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
453	c write(11,) ls_rad180./RPI, factflux,
454	c . fluxtop(utilaer(3)+1)((rinter(llm)-rinter(llm-1))1e5), !C2H2
455	c . fluxtop(utilaer(6)+1)((rinter(llm)-rinter(llm-1))1e5) !HCN
456	c close(11)
457	c endif
458	c FIN TEST: sortie de l'un des flux avec ls et factflux
459
460	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
461
462	if (firstcal.and.(j.eq.1)) then
463	print*,'Alt, densites et temp au pole (chimie):'
464	print*,'level, z_bas, nb, temp_c'
465	do l=1,klev
466	print*,l,rinter(l)-RA/1000.,nb(l),temp_c(l)
467	enddo
468	endif
469
470	if (firstcal.and.(j.eq.jjm/2)) then
471	c print*,'g,mugaz'
472	c print*,g,mugaz
473	print*,'Alt, densites et temp a l equateur (chimie):'
474	print*,'level, z_bas, nb, temp_c'
475	do l=1,klev
476	print*,l,rinter(l)-RA/1000.,nb(l),temp_c(l)
477	enddo
478	endif
479
480	c-----------------------------------------------------------------------
481	c
482	c composition
483	c ------------
484
485	do ic=1,NC
486	do l=1,klev
487	cqy(l,ic) = qy_c(j,l,ic)
488	enddo
489	enddo
490
491	c-----------------------------------------------------------------------
492	c
493	c total haze area (um2/cm3)
494	c -------------------------
495
496	if (htoh2.eq.1) then
497	do l=1,klev
498	! ATTENTION, INVERSION PAR RAPPORT A pg2.F !!!
499	surfhaze(l) = psurfhaze(j,klev+1-l)
500	c if (j.eq.25)
501	c . print*,'psurfhaze en um2/cm3:',surfhaze(l)
502	enddo
503	endif
504
505	c-----------------------------------------------------------------------
506	c
507	c Appel de chimietitan
508	c --------------------
509
510	call gptitan(jjp1,rinter,temp_c,nb,
511	$ nomqy_c,cqy,fluxtop,
512	$ declin_c,duree,(j-1),mass,
513	$ botCH4,krpd,krate,reactif,
514	$ nom_prod,nom_perte,prod,perte,
515	$ aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,
516	$ htoh2,surfhaze)
517
518	c if ( j.eq.jjm/2 )
519	c $ print*,cqy(1,1),cqy(klev,1),cqy(1,2),cqy(klev,2)
520	c if ( j.eq.jjm/2 )
521	c $ print*,qy_c(j,1,1),qy_c(j,klev,1),qy_c(j,1,2),qy_c(j,klev,2)
522
523	c stop
524
525	c Tendances composition
526	c ---------------------
527
528	do ic=1,NC
529	do l=1,klev
530	dqyc(j,l,ic) = (cqy(l,ic) - qy_c(j,l,ic))/dtchim ! en /s
531	enddo
532	enddo
533
534	c-----------------------------------------------------------------------
535	c
536	c production aer
537	c --------------
538
539	if (aerprod.eq.1) then
540
541	do ic=1,4
542	do l=1,klev
543	prodaer(j,l,ic) = cprodaer(l,ic)
544	maer(j,l,ic) = cmaer(l,ic)
545	csn(j,l,ic) = ccsn(l,ic)
546	csh(j,l,ic) = ccsh(l,ic)
547	enddo
548	enddo
549
550	endif
551
552	c***********************************************************************
553	c***********************************************************************
554	c
555	c FIN: BOUCLE SUR LES LATITUDES
556	c
557	ENDDO
558
559	c***********************************************************************
560	c***********************************************************************
561
562
563	firstcal = .false.
564	RETURN
565	END

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

Download in other formats:

Original Format