Context Navigation

calchim.F @ 537

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

SLebonnois: modification de makelmdz et create_make_gcm pour pouvoir
compiler la chimie titan. Pas de raison que ca gene les autres.
Dans cette version, les compilations de Venus et Titan fonctionnent.

Phytitan: modifications pour pouvoir compiler correctement.
Il ne manque plus que physiq.F a faire.

File size: 16.9 KB

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

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

Download in other formats:

Original Format