Context Navigation

calchim.F @ 1243

Last change on this file since 1243 was 1126, checked in by slebonnois, 11 years ago
SL: update of Titan photochemical module to include computation of chemistry up to 1300 km
File size: 18.3 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
135	c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
136	do l=klev+1,NLEV
137	rinter(l) = rinter(klev)
138	& + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
139	rmil(l) = rmil(klev)
140	& + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
141	enddo
142
143	c lecture de tcp.ver, une seule fois
144	c remplissage r1d,t1d,ct1d,p1d
145	open(11,file='../../INPUT/tcp.ver',status='old')
146	read(11,*) dummy
147	do i=1,131
148	read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
149	c print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
150	enddo
151	close(11)
152
153	c extension pour klev+1 a NLEV avec tcp.ver
154	c la jonction klev/klev+1 est brutale... Tant pis ?
155	ialt=1
156	do l=klev+1,NLEV
157	zalt=rmil(l)-RA/1000.
158	do i=ialt,130
159	if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
160	ialt=i
161	endif
162	enddo
163	factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
164	press_c(l) = exp( log(p1d(ialt)) * (1-factalt)
165	& + log(p1d(ialt+1)) * factalt )
166	nb(l) = exp( log(ct1d(ialt)) * (1-factalt)
167	& + log(ct1d(ialt+1)) * factalt )
168	temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
169	print*,l,zalt,press_c(l),nb(l),temp_c(l)
170	enddo
171
172	c caracteristiques des composes:
173	c -----------------------------
174	mass(:) = 0.0
175	call comp(nomqy_c,nb,temp_c,mass,md)
176	print*,' Mass'
177	do ic=1,NC
178	print*,nomqy_c(ic),mass(ic)
179	c if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
180	enddo
181
182
183	c Stockage des composes utilises dans la prod d aerosols
184	c (aerprod=1) et dans H-> H2 (htoh2=1): utilaer
185	c ! decalage de 1 car utilise dans le c !
186
187	do ic=1,NC
188
189	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190	c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
191	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
192	c if (nomqy_c(ic).eq."CH4") then
193	c do l=1,llm
194	c do j=1,ip1jmp1
195	c if (qy_c(j,l,ic).le.0.015) qy_c(j,l,ic) = 0.015
196	c enddo
197	c enddo
198	c endif
199	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
200
201	if (nomqy_c(ic).eq."C4H2") then
202	utilaer(10) = ic-1
203	endif
204	if (nomqy_c(ic).eq."HCN") then
205	utilaer(6) = ic-1
206	endif
207	if (nomqy_c(ic).eq."HC3N") then
208	utilaer(7) = ic-1
209	endif
210	if (nomqy_c(ic).eq."NCCN") then
211	utilaer(14) = ic-1
212	endif
213	if (nomqy_c(ic).eq."CH3CN") then
214	utilaer(15) = ic-1
215	utilaer(16) = ic-1 ! si pas C2H3CN, CH3CN utilise, mais reac. nulle
216	endif
217	if (nomqy_c(ic).eq."H") then
218	utilaer(1) = ic-1
219	endif
220	if (nomqy_c(ic).eq."H2") then
221	utilaer(2) = ic-1
222	endif
223	if (nomqy_c(ic).eq."C2H2") then
224	utilaer(3) = ic-1
225	endif
226	if (nomqy_c(ic).eq."AC6H6") then
227	utilaer(13) = ic-1
228	endif
229	if (nomqy_c(ic).eq."C2H3CN") then
230	utilaer(16) = ic-1
231	endif
232	if (nomqy_c(ic).eq."C2") then
233	utilaer(4) = ic-1
234	endif
235	if (nomqy_c(ic).eq."C2H") then
236	utilaer(5) = ic-1
237	endif
238	if (nomqy_c(ic).eq."C3N") then
239	utilaer(8) = ic-1
240	endif
241	if (nomqy_c(ic).eq."H2CN") then
242	utilaer(9) = ic-1
243	endif
244	if (nomqy_c(ic).eq."C4H3") then
245	utilaer(11) = ic-1
246	endif
247	if (nomqy_c(ic).eq."AC6H5") then
248	utilaer(12) = ic-1
249	endif
250
251	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
252	c if ((nomqy_c(ic).eq."HC3N").or.
253	c $ (nomqy_c(ic).eq."C3N")) then
254	c DO j=1,ip1jmp1
255	c do l=1,34 ! p>~ 1 mbar
256	c qy_c(j,l,ic) = 1.e-30
257	c enddo
258	c ENDDO
259	c endif
260	c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
261
262	enddo
263
264	c taux de photodissociations:
265	c --------------------------
266	call disso(krpd,jjp1)
267
268	c reactions chimiques:
269	c -------------------
270	call chimie(nomqy_c,nb,temp_c,krate,reactif,
271	. nom_perte,nom_prod,perte,prod)
272	c print*,'nom_prod, nom_perte:'
273	c do ic=1,NC
274	c print*,nom_prod(ic),nom_perte(ic)
275	c enddo
276	c print*,'premieres prod, perte(1:reaction,2:compagnon):'
277	c do ic=1,NC
278	c print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic)
279	c enddo
280
281	c l = klev-3
282	c print*,'krate a p=',press_c(l),' reactifs et produits:'
283	c do ic=1,ND+1
284	c print,ic,krpd(7,ic,l,4)nb(l)," ",
285	c . nomqy_c(reactif(1,ic)+1),
286	c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
287	c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
288	c enddo
289	c do ic=ND+2,NR
290	c print*,ic,krate(l,ic)," ",
291	c . nomqy_c(reactif(1,ic)+1),
292	c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
293	c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
294	c enddo
295
296
297	c init kedd
298	c ---------
299	c kedd stays constant with time and space
300	c => linked to pressure rather than altitude...
301
302	kedd(:) = 5e2 ! valeur mise par defaut
303	! UNITE ? doit etre ok pour gptitan
304	do l=1,NLEV
305	zalt=rmil(l)-RA/1000. ! z en km
306	if ((zalt.ge.250.).and.(zalt.le.400.)) then
307	kedd(l) = 10.**(3.+(zalt-250.)/50.)
308	! 1E3 at 250 km
309	! 1E6 at 400 km
310	elseif ((zalt.gt.400.).and.(zalt.le.600.)) then
311	kedd(l) = 10.*(6.+1.3(zalt-400.)/200.)
312	! 2E7 at 600 km
313	elseif ((zalt.gt.600.).and.(zalt.le.900.)) then
314	kedd(l) = 10.*(7.3+0.7(zalt-600.)/300.)
315	! 1E8 above 900 km
316	elseif ( zalt.gt.900. ) then
317	kedd(l) = 1.e8
318	endif
319	enddo
320
321	ENDIF ! premier appel
322
323	c***********************************************************************
324	c-----------------------------------------------------------------------
325
326	c calcul declin_c (en degres)
327	c ---------------------------
328
329	declin_c = declin_rad*180./RPI
330	c print*,'declinaison en degre=',declin_c
331
332	c***********************************************************************
333	c***********************************************************************
334	c
335	c BOUCLE SUR LES LATITUDES
336	c
337	DO j=1,nlon
338
339	if (j.eq.1) then
340	jm1=1
341	else
342	jm1=j-1
343	endif
344
345	if((j.eq.1).or.(klat(j).ne.klat(jm1))) then
346
347	c***********************************************************************
348	c***********************************************************************
349
350	c-----------------------------------------------------------------------
351	c
352	c Distance radiale (intercouches, en km)
353	c ----------------------------------------
354
355	do l=1,klev
356	rinter(l) = (RA+czlev(j,l))/1000.
357	rmil(l) = (RA+czlay(j,l))/1000.
358	c print*,rinter(l)
359	enddo
360
361	c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
362	do l=klev+1,NLEV
363	rinter(l) = rinter(klev)
364	& + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
365	rmil(l) = rmil(klev)
366	& + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
367	enddo
368
369	c-----------------------------------------------------------------------
370	c
371	c Temperature, pression (mbar), densite (cm-3)
372	c -------------------------------------------
373
374	DO l=1,klev
375	c temp_c (K):
376	temp_c(l) = ctemp(j,l)
377	c press_c (mbar):
378	press_c(l) = cplay(j,l)/100.
379	c nb (cm-3):
380	nb(l) = 1.e-4press_c(l) / (RKBOLtemp_c(l))
381	ENDDO
382	c extension pour klev+1 a NLEV avec tcp.ver
383	ialt=1
384	do l=klev+1,NLEV
385	zalt=rmil(l)-RA/1000.
386	do i=ialt,130
387	if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
388	ialt=i
389	endif
390	enddo
391	factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
392	press_c(l) = exp( log(p1d(ialt)) * (1-factalt)
393	& + log(p1d(ialt+1)) * factalt )
394	nb(l) = exp( log(ct1d(ialt)) * (1-factalt)
395	& + log(ct1d(ialt+1)) * factalt )
396	temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
397	enddo
398
399	c-----------------------------------------------------------------------
400	c
401	c passage krpd => krate
402	c ---------------------
403	c initialisation krate pour dissociations
404
405	if ((declin_c*10+267).lt.14.) then
406	declinint = 0
407	dec = 0
408	else
409	if ((declin_c*10+267).gt.520.) then
410	declinint = 14
411	dec = 534
412	else
413	declinint = 1
414	dec = 27
415	do while( (declin_c*10+267).ge.real(dec+20) )
416	dec = dec+40
417	declinint = declinint+1
418	enddo
419	endif
420	endif
421	if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
422	factdec = ( declin_c - (dec-267)/10. ) / 4.
423	else
424	factdec = ( declin_c - (dec-267)/10. ) / 2.7
425	endif
426
427	do l=1,NLEV
428
429	c INTERPOL EN ALT POUR k (krpd tous les deux km)
430	ialt = int((rmil(l)-RA/1000.)/2.)+1
431	factalt = (rmil(l)-RA/1000.)/2.-(ialt-1)
432
433	do i=1,ND+1
434	krpddecm1 = krpd(declinint ,i,ialt ,klat(j)) * (1-factalt)
435	& + krpd(declinint ,i,ialt+1,klat(j)) * factalt
436	krpddec = krpd(declinint+1,i,ialt ,klat(j)) * (1-factalt)
437	& + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
438	krpddecp1 = krpd(declinint+2,i,ialt ,klat(j)) * (1-factalt)
439	& + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
440
441	! ND+1 correspond a la dissociation de N2 par les GCR
442	if ( factdec.lt.0. ) then
443	krate(l,i) = krpddecm1 * abs(factdec)
444	& + krpddec * ( 1 + factdec)
445	endif
446	if ( factdec.gt.0. ) then
447	krate(l,i) = krpddecp1 * factdec
448	& + krpddec * ( 1 - factdec)
449	endif
450	if ( factdec.eq.0. ) krate(l,i) = krpddec
451	enddo
452	enddo
453
454	c-----------------------------------------------------------------------
455	c
456	c composition
457	c ------------
458
459	do ic=1,NC
460	do l=1,klev
461	cqy(l,ic) = qy_c(j,l,ic)
462	cqy0(l,ic) = cqy(l,ic)
463	enddo
464	enddo
465
466	c lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
467
468	WRITE(str2,'(i2.2)') klat(j)
469	fich = "comp_"//str2//".dat"
470	inquire (file=fich,exist=okfic)
471	if (okfic) then
472	open(11,file=fich,status='old')
473	c premiere ligne=declin
474	read(11,'(A15)') ficdec
475	write(curdec,'(E15.5)') declin_c
476	c si la declin est la meme: ce fichier a deja ete reecrit
477	c on lit la colonne t-1 au lieu de la colonne t
478	c (cas d une bande de latitude partagee par 2 procs)
479	do ic=1,NC
480	read(11,'(A10)') name
481	if (name.ne.nomqy_c(ic)) then
482	print*,"probleme lecture ",fich
483	print*,name,nomqy_c(ic)
484	stop
485	endif
486	if (ficdec.eq.curdec) then
487	do l=klev+1,NLEV
488	read(11,*) ficalt,cqy(l,ic),newq
489	enddo
490	else
491	do l=klev+1,NLEV
492	read(11,*) ficalt,oldq,cqy(l,ic)
493	enddo
494	endif
495	enddo
496	close(11)
497	else ! le fichier n'est pas la
498	do ic=1,NC
499	do l=klev+1,NLEV
500	cqy(l,ic)=cqy(klev,ic)
501	enddo
502	enddo
503	endif
504	cqy0 = cqy
505
506	c-----------------------------------------------------------------------
507	c
508	c total haze area (um2/cm3)
509	c -------------------------
510
511	if (htoh2.eq.1) then
512	do l=1,klev
513	! ATTENTION, INVERSION PAR RAPPORT A pg2.F !!!
514	surfhaze(l) = psurfhaze(j,klev+1-l)
515	c if (j.eq.25)
516	c . print*,'psurfhaze en um2/cm3:',surfhaze(l)
517	enddo
518	endif
519
520	c-----------------------------------------------------------------------
521	c
522	c Appel de chimietitan
523	c --------------------
524
525	call gptitan(rinter,temp_c,nb,
526	$ nomqy_c,cqy,
527	$ duree,(klat(j)-1),mass,md,
528	$ kedd,botCH4,krate,reactif,
529	$ nom_prod,nom_perte,prod,perte,
530	$ aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,
531	$ htoh2,surfhaze)
532
533	c Tendances composition
534	c ---------------------
535
536	do ic=1,NC
537	do l=1,klev
538	dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! en /s
539	enddo
540	enddo
541
542	c-----------------------------------------------------------------------
543	c
544	c production aer
545	c --------------
546
547	if (aerprod.eq.1) then
548
549	do ic=1,4
550	do l=1,klev
551	prodaer(j,l,ic) = cprodaer(l,ic)
552	maer(j,l,ic) = cmaer(l,ic)
553	csn(j,l,ic) = ccsn(l,ic)
554	csh(j,l,ic) = ccsh(l,ic)
555	enddo
556	enddo
557
558	endif
559
560	c-----------------------------------------------------------------------
561	c
562	c sauvegarde compo sur NLEV
563	c -------------------------
564
565	c dans fichier compo_klat(j) (01 à 49)
566
567	open(11,file=fich,status='replace')
568	c premiere ligne=declin
569	write(11,'(E15.5)') declin_c
570	do ic=1,NC
571	write(11,'(A10)') nomqy_c(ic)
572	do l=klev+1,NLEV
573	write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-RA/1000.,
574	. cqy0(l,ic),cqy(l,ic)
575	enddo
576	enddo
577	close(11)
578
579	c***********************************************************************
580	c***********************************************************************
581
582	c FIN: BOUCLE SUR LES LATITUDES
583
584	else ! same latitude, we don't do calculations again
585	dqyc(j,:,:) = dqyc(jm1,:,:)
586	if (aerprod.eq.1) then
587	prodaer(j,:,:) = prodaer(jm1,:,:)
588	maer(j,:,:) = maer(jm1,:,:)
589	csn(j,:,:) = csn(jm1,:,:)
590	csh(j,:,:) = csh(jm1,:,:)
591	endif
592	endif
593
594	ENDDO
595
596	c***********************************************************************
597	c***********************************************************************
598
599
600	firstcal = .false.
601	RETURN
602	END

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

Download in other formats:

Original Format