Context Navigation

calchim.F @ 1379

Last change on this file since 1379 was 1379, checked in by slebonnois, 10 years ago
SL: correction of a bug in Titan's chemistry
File size: 18.5 KB

Line
1	SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim,
2	. ctemp,cplay,cplev,czlay,czlev,
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 adaptation pour // : 04/2013
12	c extension chimie jusqu a 1300km : 10/2013
13	c
14	c-------------------------------------------------
15	c
16	use dimphy
17	use common_mod, only:utilaer,maer,prodaer,csn,csh,psurfhaze,
18	. NLEV,NLRT,NC,ND,NR
19	USE comgeomphy, only: rlatd
20	USE moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
21	implicit none
22	#include "dimensions.h"
23	#include "clesphys.h"
24	#include "paramet.h"
25	#include "YOMCST.h"
26
27	c Arguments
28	c ---------
29
30	INTEGER nlon ! nb of horiz points
31	INTEGER ny ! nb de composes (nqmax-nmicro)
32	REAL qy_c(nlon,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(nlon,klev) ! Temperature
37	REAL cplay(nlon,klev) ! pression (Pa)
38	REAL cplev(nlon,klev+1) ! pression intercouches (Pa)
39	REAL czlay(nlon,klev) ! altitude (m)
40	REAL czlev(nlon,klev+1) ! altitude intercouches (m)
41
42	REAL dqyc(nlon,klev,NC) ! Tendances especes chimiques
43
44	c Local variables :
45	c -----------------
46
47	integer i,j,l,ic,jm1
48
49	c variables envoyees dans la chimie: double precision
50
51	REAL temp_c(NLEV)
52	REAL press_c(NLEV) ! T,p(mbar) a 1 lat donnee
53	REAL cqy(NLEV,NC) ! frac mol qui seront modifiees
54	REAL cqy0(NLEV,NC) ! frac mol avant modif
55	REAL surfhaze(NLEV)
56	REAL cprodaer(NLEV,4),cmaer(NLEV,4)
57	REAL ccsn(NLEV,4),ccsh(NLEV,4)
58	! rmil: milieu de couche, grille pour K,D,p,ct,T,y
59	! rinter: intercouche (grille RA dans la chimie)
60	REAL rmil(NLEV),rinter(NLEV),nb(NLEV)
61	REAL,save :: kedd(NLEV)
62
63	REAL a,b
64	character str11,str22
65	REAL latit
66	character*20 formt1,formt2,formt3
67
68	c variables locales initialisees au premier appel
69
70	LOGICAL firstcal
71	DATA firstcal/.true./
72	SAVE firstcal
73
74	integer dec,declinint,ialt
75	REAL declin_c ! declinaison en degres
76	real factalt,factdec,krpddec,krpddecp1,krpddecm1
77	real duree
78	REAL,save :: mass(NC)
79	REAL,save,allocatable :: md(:,:)
80	REAL,save :: botCH4
81	DATA botCH4/0.05/
82	real,save :: r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
83	REAL,save,allocatable :: krpd(:,:,:,:),krate(:,:)
84	integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC)
85	integer,save :: prod(200,NC),perte(2,200,NC)
86	character dummy30,fich7,ficdec15,curdec15,name*10
87	real ficalt,oldq,newq,zalt
88	logical okfic
89
90	c-----------------------------------------------------------------------
91	c***********************************************************************
92	c
93	c Initialisations :
94	c ----------------
95
96	duree = dtchim ! passage en real*4 pour appel a routines C
97
98	IF (firstcal) THEN
99
100	print*,'CHIMIE, premier appel'
101
102	c ************************************
103	c Au premier appel, initialisation de toutes les variables
104	c pour les routines de la chimie.
105	c ************************************
106
107	allocate(krpd(15,ND+1,NLRT,jjp1),krate(NLEV,NR),md(NLEV,NC))
108
109	c Verification du nombre de composes: coherence common_mod et nqmax-nmicro
110	c ----------------------------------
111
112	if (ny.ne.NC) then
113	print*,'PROBLEME de coherence dans le nombre de composes:'
114	. ,ny,NC
115	STOP
116	endif
117
118	c calcul de temp_c, densites et press_c en moyenne planetaire :
119	c --------------------------------------------------------------
120
121	print*,'pression, densites et temp (init chimie):'
122	print*,'level, press_c, nb, temp_c'
123	DO l=1,klev
124	rinter(l) = (zlevmoy(l)+RA)/1000.
125	rmil(l) = (zlaymoy(l)+RA)/1000.
126	c temp_c (K):
127	temp_c(l) = tmoy(l)
128	c press_c (mbar):
129	press_c(l) = playmoy(l)/100.
130	c nb (cm-3):
131	nb(l) = 1.e-4press_c(l) / (RKBOLtemp_c(l))
132	print*,l,rmil(l)-RA/1000.,press_c(l),nb(l),temp_c(l)
133	ENDDO
134	rinter(klev+1)=(zlevmoy(klev+1)+RA)/1000.
135
136	c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
137	do l=klev+2,NLEV
138	rinter(l) = rinter(klev+1)
139	& + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1)
140	rmil(l-1) = (rinter(l-1)+rinter(l))/2.
141	enddo
142	rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2.
143
144	c lecture de tcp.ver, une seule fois
145	c remplissage r1d,t1d,ct1d,p1d
146	open(11,file='../../INPUT/tcp.ver',status='old')
147	read(11,*) dummy
148	do i=1,131
149	read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
150	c print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
151	enddo
152	close(11)
153
154	c extension pour klev+1 a NLEV avec tcp.ver
155	c la jonction klev/klev+1 est brutale... Tant pis ?
156	ialt=1
157	do l=klev+1,NLEV
158	zalt=rmil(l)-RA/1000.
159	do i=ialt,130
160	if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
161	ialt=i
162	endif
163	enddo
164	factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
165	press_c(l) = exp( log(p1d(ialt)) * (1-factalt)
166	& + log(p1d(ialt+1)) * factalt )
167	nb(l) = exp( log(ct1d(ialt)) * (1-factalt)
168	& + log(ct1d(ialt+1)) * factalt )
169	temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
170	print*,l,zalt,press_c(l),nb(l),temp_c(l)
171	enddo
172
173	c caracteristiques des composes:
174	c -----------------------------
175	mass(:) = 0.0
176	call comp(nomqy_c,nb,temp_c,mass,md)
177	print*,' Mass'
178	do ic=1,NC
179	print*,nomqy_c(ic),mass(ic)
180	c if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
181	enddo
182
183
184	c Stockage des composes utilises dans la prod d aerosols
185	c (aerprod=1) et dans H-> H2 (htoh2=1): utilaer
186	c ! decalage de 1 car utilise dans le c !
187
188	do ic=1,NC
189
190	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
192	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
193	c if (nomqy_c(ic).eq."CH4") then
194	c do l=1,llm
195	c do j=1,ip1jmp1
196	c if (qy_c(j,l,ic).le.0.015) qy_c(j,l,ic) = 0.015
197	c enddo
198	c enddo
199	c endif
200	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
201
202	if (nomqy_c(ic).eq."C4H2") then
203	utilaer(10) = ic-1
204	endif
205	if (nomqy_c(ic).eq."HCN") then
206	utilaer(6) = ic-1
207	endif
208	if (nomqy_c(ic).eq."HC3N") then
209	utilaer(7) = ic-1
210	endif
211	if (nomqy_c(ic).eq."NCCN") then
212	utilaer(14) = ic-1
213	endif
214	if (nomqy_c(ic).eq."CH3CN") then
215	utilaer(15) = ic-1
216	utilaer(16) = ic-1 ! si pas C2H3CN, CH3CN utilise, mais reac. nulle
217	endif
218	if (nomqy_c(ic).eq."H") then
219	utilaer(1) = ic-1
220	endif
221	if (nomqy_c(ic).eq."H2") then
222	utilaer(2) = ic-1
223	endif
224	if (nomqy_c(ic).eq."C2H2") then
225	utilaer(3) = ic-1
226	endif
227	if (nomqy_c(ic).eq."AC6H6") then
228	utilaer(13) = ic-1
229	endif
230	if (nomqy_c(ic).eq."C2H3CN") then
231	utilaer(16) = ic-1
232	endif
233	if (nomqy_c(ic).eq."C2") then
234	utilaer(4) = ic-1
235	endif
236	if (nomqy_c(ic).eq."C2H") then
237	utilaer(5) = ic-1
238	endif
239	if (nomqy_c(ic).eq."C3N") then
240	utilaer(8) = ic-1
241	endif
242	if (nomqy_c(ic).eq."H2CN") then
243	utilaer(9) = ic-1
244	endif
245	if (nomqy_c(ic).eq."C4H3") then
246	utilaer(11) = ic-1
247	endif
248	if (nomqy_c(ic).eq."AC6H5") then
249	utilaer(12) = ic-1
250	endif
251
252	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
253	c if ((nomqy_c(ic).eq."HC3N").or.
254	c $ (nomqy_c(ic).eq."C3N")) then
255	c DO j=1,ip1jmp1
256	c do l=1,34 ! p>~ 1 mbar
257	c qy_c(j,l,ic) = 1.e-30
258	c enddo
259	c ENDDO
260	c endif
261	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
262
263	enddo
264
265	c taux de photodissociations:
266	c --------------------------
267	call disso(krpd,jjp1)
268
269	c reactions chimiques:
270	c -------------------
271	call chimie(nomqy_c,nb,temp_c,krate,reactif,
272	. nom_perte,nom_prod,perte,prod)
273	c print*,'nom_prod, nom_perte:'
274	c do ic=1,NC
275	c print*,nom_prod(ic),nom_perte(ic)
276	c enddo
277	c print*,'premieres prod, perte(1:reaction,2:compagnon):'
278	c do ic=1,NC
279	c print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic)
280	c enddo
281
282	c l = klev-3
283	c print*,'krate a p=',press_c(l),' reactifs et produits:'
284	c do ic=1,ND+1
285	c print,ic,krpd(7,ic,l,4)nb(l)," ",
286	c . nomqy_c(reactif(1,ic)+1),
287	c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
288	c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
289	c enddo
290	c do ic=ND+2,NR
291	c print*,ic,krate(l,ic)," ",
292	c . nomqy_c(reactif(1,ic)+1),
293	c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
294	c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
295	c enddo
296
297
298	c init kedd
299	c ---------
300	c kedd stays constant with time and space
301	c => linked to pressure rather than altitude...
302
303	kedd(:) = 5e2 ! valeur mise par defaut
304	! UNITE ? doit etre ok pour gptitan
305	do l=1,NLEV
306	zalt=rmil(l)-RA/1000. ! z en km
307	if ((zalt.ge.250.).and.(zalt.le.400.)) then
308	kedd(l) = 10.**(3.+(zalt-250.)/50.)
309	! 1E3 at 250 km
310	! 1E6 at 400 km
311	elseif ((zalt.gt.400.).and.(zalt.le.600.)) then
312	kedd(l) = 10.*(6.+1.3(zalt-400.)/200.)
313	! 2E7 at 600 km
314	elseif ((zalt.gt.600.).and.(zalt.le.900.)) then
315	kedd(l) = 10.*(7.3+0.7(zalt-600.)/300.)
316	! 1E8 above 900 km
317	elseif ( zalt.gt.900. ) then
318	kedd(l) = 1.e8
319	endif
320	enddo
321
322	ENDIF ! premier appel
323
324	c***********************************************************************
325	c-----------------------------------------------------------------------
326
327	c calcul declin_c (en degres)
328	c ---------------------------
329
330	declin_c = declin_rad*180./RPI
331	c print*,'declinaison en degre=',declin_c
332
333	c***********************************************************************
334	c***********************************************************************
335	c
336	c BOUCLE SUR LES LATITUDES
337	c
338	DO j=1,nlon
339
340	if (j.eq.1) then
341	jm1=1
342	else
343	jm1=j-1
344	endif
345
346	if((j.eq.1).or.(klat(j).ne.klat(jm1))) then
347
348	c***********************************************************************
349	c***********************************************************************
350
351	c-----------------------------------------------------------------------
352	c
353	c Distance radiale (intercouches, en km)
354	c ----------------------------------------
355
356	do l=1,klev
357	rinter(l) = (RA+czlev(j,l))/1000.
358	rmil(l) = (RA+czlay(j,l))/1000.
359	c print*,rinter(l)
360	enddo
361	rinter(klev+1)=(RA+czlev(j,klev+1))/1000.
362
363	c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
364	do l=klev+2,NLEV
365	rinter(l) = rinter(klev+1)
366	& + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1)
367	rmil(l-1) = (rinter(l-1)+rinter(l))/2.
368	enddo
369	rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2.
370
371	c-----------------------------------------------------------------------
372	c
373	c Temperature, pression (mbar), densite (cm-3)
374	c -------------------------------------------
375
376	DO l=1,klev
377	c temp_c (K):
378	temp_c(l) = ctemp(j,l)
379	c press_c (mbar):
380	press_c(l) = cplay(j,l)/100.
381	c nb (cm-3):
382	nb(l) = 1.e-4press_c(l) / (RKBOLtemp_c(l))
383	ENDDO
384	c extension pour klev+1 a NLEV avec tcp.ver
385	ialt=1
386	do l=klev+1,NLEV
387	zalt=rmil(l)-RA/1000.
388	do i=ialt,130
389	if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
390	ialt=i
391	endif
392	enddo
393	factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
394	press_c(l) = exp( log(p1d(ialt)) * (1-factalt)
395	& + log(p1d(ialt+1)) * factalt )
396	nb(l) = exp( log(ct1d(ialt)) * (1-factalt)
397	& + log(ct1d(ialt+1)) * factalt )
398	temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
399	enddo
400
401	c-----------------------------------------------------------------------
402	c
403	c passage krpd => krate
404	c ---------------------
405	c initialisation krate pour dissociations
406
407	if ((declin_c*10+267).lt.14.) then
408	declinint = 0
409	dec = 0
410	else
411	if ((declin_c*10+267).gt.520.) then
412	declinint = 14
413	dec = 534
414	else
415	declinint = 1
416	dec = 27
417	do while( (declin_c*10+267).ge.real(dec+20) )
418	dec = dec+40
419	declinint = declinint+1
420	enddo
421	endif
422	endif
423	if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
424	factdec = ( declin_c - (dec-267)/10. ) / 4.
425	else
426	factdec = ( declin_c - (dec-267)/10. ) / 2.7
427	endif
428
429	do l=1,NLEV
430
431	c INTERPOL EN ALT POUR k (krpd tous les deux km)
432	ialt = int((rmil(l)-RA/1000.)/2.)+1
433	factalt = (rmil(l)-RA/1000.)/2.-(ialt-1)
434
435	do i=1,ND+1
436	krpddecm1 = krpd(declinint ,i,ialt ,klat(j)) * (1-factalt)
437	& + krpd(declinint ,i,ialt+1,klat(j)) * factalt
438	krpddec = krpd(declinint+1,i,ialt ,klat(j)) * (1-factalt)
439	& + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
440	krpddecp1 = krpd(declinint+2,i,ialt ,klat(j)) * (1-factalt)
441	& + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
442
443	! ND+1 correspond a la dissociation de N2 par les GCR
444	if ( factdec.lt.0. ) then
445	krate(l,i) = krpddecm1 * abs(factdec)
446	& + krpddec * ( 1 + factdec)
447	endif
448	if ( factdec.gt.0. ) then
449	krate(l,i) = krpddecp1 * factdec
450	& + krpddec * ( 1 - factdec)
451	endif
452	if ( factdec.eq.0. ) krate(l,i) = krpddec
453	enddo
454	enddo
455
456	c-----------------------------------------------------------------------
457	c
458	c composition
459	c ------------
460
461	do ic=1,NC
462	do l=1,klev
463	cqy(l,ic) = qy_c(j,l,ic)
464	cqy0(l,ic) = cqy(l,ic)
465	enddo
466	enddo
467
468	c lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
469
470	WRITE(str2,'(i2.2)') klat(j)
471	fich = "comp_"//str2//".dat"
472	inquire (file=fich,exist=okfic)
473	if (okfic) then
474	open(11,file=fich,status='old')
475	c premiere ligne=declin
476	read(11,'(A15)') ficdec
477	write(curdec,'(E15.5)') declin_c
478	c si la declin est la meme: ce fichier a deja ete reecrit
479	c on lit la colonne t-1 au lieu de la colonne t
480	c (cas d une bande de latitude partagee par 2 procs)
481	do ic=1,NC
482	read(11,'(A10)') name
483	if (name.ne.nomqy_c(ic)) then
484	print*,"probleme lecture ",fich
485	print*,name,nomqy_c(ic)
486	stop
487	endif
488	if (ficdec.eq.curdec) then
489	do l=klev+1,NLEV
490	read(11,*) ficalt,cqy(l,ic),newq
491	enddo
492	else
493	do l=klev+1,NLEV
494	read(11,*) ficalt,oldq,cqy(l,ic)
495	enddo
496	endif
497	enddo
498	close(11)
499	else ! le fichier n'est pas la
500	do ic=1,NC
501	do l=klev+1,NLEV
502	cqy(l,ic)=cqy(klev,ic)
503	enddo
504	enddo
505	endif
506	cqy0 = cqy
507
508	c-----------------------------------------------------------------------
509	c
510	c total haze area (um2/cm3)
511	c -------------------------
512
513	if (htoh2.eq.1) then
514	do l=1,klev
515	! ATTENTION, INVERSION PAR RAPPORT A pg2.F !!!
516	surfhaze(l) = psurfhaze(j,klev+1-l)
517	c if (j.eq.25)
518	c . print*,'psurfhaze en um2/cm3:',surfhaze(l)
519	enddo
520	endif
521
522	c-----------------------------------------------------------------------
523	c
524	c Appel de chimietitan
525	c --------------------
526
527	call gptitan(rinter,temp_c,nb,
528	$ nomqy_c,cqy,
529	$ duree,(klat(j)-1),mass,md,
530	$ kedd,botCH4,krate,reactif,
531	$ nom_prod,nom_perte,prod,perte,
532	$ aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,
533	$ htoh2,surfhaze)
534
535	c Tendances composition
536	c ---------------------
537
538	do ic=1,NC
539	do l=1,klev
540	dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! en /s
541	enddo
542	enddo
543
544	c-----------------------------------------------------------------------
545	c
546	c production aer
547	c --------------
548
549	if (aerprod.eq.1) then
550
551	do ic=1,4
552	do l=1,klev
553	prodaer(j,l,ic) = cprodaer(l,ic)
554	maer(j,l,ic) = cmaer(l,ic)
555	csn(j,l,ic) = ccsn(l,ic)
556	csh(j,l,ic) = ccsh(l,ic)
557	enddo
558	enddo
559
560	endif
561
562	c-----------------------------------------------------------------------
563	c
564	c sauvegarde compo sur NLEV
565	c -------------------------
566
567	c dans fichier compo_klat(j) (01 à 49)
568
569	open(11,file=fich,status='replace')
570	c premiere ligne=declin
571	write(11,'(E15.5)') declin_c
572	do ic=1,NC
573	write(11,'(A10)') nomqy_c(ic)
574	do l=klev+1,NLEV
575	write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-RA/1000.,
576	. cqy0(l,ic),cqy(l,ic)
577	enddo
578	enddo
579	close(11)
580
581	c***********************************************************************
582	c***********************************************************************
583
584	c FIN: BOUCLE SUR LES LATITUDES
585
586	else ! same latitude, we don't do calculations again
587	dqyc(j,:,:) = dqyc(jm1,:,:)
588	if (aerprod.eq.1) then
589	prodaer(j,:,:) = prodaer(jm1,:,:)
590	maer(j,:,:) = maer(jm1,:,:)
591	csn(j,:,:) = csn(jm1,:,:)
592	csh(j,:,:) = csh(jm1,:,:)
593	endif
594	endif
595
596	ENDDO
597
598	c***********************************************************************
599	c***********************************************************************
600
601
602	firstcal = .false.
603	RETURN
604	END

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

Download in other formats:

Original Format