Context Navigation

calchim.F90 @ 1941

Last change on this file since 1941 was 1928, checked in by jvatant, 8 years ago
Forget prints in r1927 ! --JVO
File size: 22.5 KB

Line
1	SUBROUTINE calchim(ngrid,qy_c,declin,dtchim, &
2	ctemp,cpphi,cplay,cplev,czlay,czlev,dqyc)
3
4	!---------------------------------------------------------------------------------------------------------
5	!
6	! Purpose : Interface subroutine to photochemical C model for Titan GCM.
7	! -------
8	! The subroutine computes the chemical processes for a single vertical column.
9	!
10	! - Only tendencies are returned.
11	! - With moyzon_ch=.true. and input vectors zonally averaged
12	! the calculation is done only once per lat. band
13	!
14	! Authors: + S. Lebonnois : 01/2000 \| 09/2003
15	! ------- adaptation for Titan 3D : 02/2009
16	! adaptation for // : 04/2013
17	! extension chemistry up to 1300km : 10/2013
18	!
19	! + J. Vatant d'Ollone
20	! + 02/17 - adaptation for the issu du generic
21	! + 01/18 - 03/18 - Major transformations :
22	! - Upper chemistry fields are now stored in startfi
23	! and defined on a pressure grid from Vervack profile
24	! - Altitudes above GCM top are computed using correct g(z)
25	! to be coherent with obs and with dz=10km for UV processes
26	! but this should be even better if all physiq "under" was doing this !
27	! - These modifs enables to run chemistry with others resolution than 32x48x55 !
28	! - Only the actinic fluxes are still read in a 49-lat input but interp. on lat grid
29	! - Chemistry can still be done in 2D
30	! -> Calcul. once per band lat and put same tendency in all longi.
31	! Check for negs in physiq_mod.
32	! -> If procs sharing a lat band, no problem, the calcul will just be done twice.
33	! -> Will not work with Dynamico, where the chemistry will have to be done in 3D.
34	! ( and there'll be work to do to get rid of averaged fields )
35	!
36	! + STILL TO DO : Replug the interaction with haze (cf titan.old) -> to see with JB.
37	!---------------------------------------------------------------------------------------------------------
38
39	! --------------------------------------------------------------------
40	! Structure :
41	! -----------
42	! 0. Declarations
43	! I. Init and firstcall
44	! 1. Read and store Vervack profile
45	! 2. Compute planetar averaged atm. properties
46	! 3. Init compound caracteristics
47	! 4. Init photodissociations rates from actinic fluxes
48	! 5. Init chemical reactions
49	! 6. Init eddy diffusion coeff
50	! II. Loop on latitudes/grid-points
51	! 0. Check on 2D chemistry
52	! 1. Compute atm. properties at grid points
53	! 2. Interpolate photodissociation rates at lat,alt,dec
54	! 3. Read composition
55	! 4. Call main solver gptitan C routine
56	! 5. Calculate output tendencies on advected tracers
57	! 6. Update upper chemistry fields
58	! 0bis. If 2D chemsitry, don't recalculate if needed
59	! -----------------------------------------------------------------
60
61	USE comchem_h
62	USE dimphy
63	USE datafile_mod, ONLY: datadir
64	USE comcstfi_mod, ONLY: g, rad, pi, r, kbol
65	USE geometry_mod, ONLY: latitude
66	USE logic_mod, ONLY: moyzon_ch
67	USE moyzon_mod, ONLY: tmoy, playmoy, phimoy, zlaymoy, zlevmoy
68
69	IMPLICIT NONE
70
71	! ------------------------------------------
72	! ********* 0. Declarations ***********
73	! ------------------------------------------
74
75	! Arguments
76	! ---------
77
78	INTEGER, INTENT(IN) :: ngrid ! Number of atmospheric columns.
79	REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(IN) :: qy_c ! Chemical species on GCM layers after adv.+diss. (mol/mol).
80	REAL*8, INTENT(IN) :: declin ! Solar declination (rad).
81	REAL*8, INTENT(IN) :: dtchim ! Chemistry timsetep (s).
82	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: ctemp ! Mid-layer temperature (K).
83	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: cpphi ! Mid-layer geopotential (m2.s-2).
84	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: cplay ! Mid-layer pressure (Pa).
85	REAL*8, DIMENSION(ngrid,klev+1), INTENT(IN) :: cplev ! Inter-layer pressure (Pa).
86	REAL*8, DIMENSION(ngrid,klev), INTENT(IN) :: czlay ! Mid-layer altitude (m).
87	REAL*8, DIMENSION(ngrid,klev+1), INTENT(IN) :: czlev ! Inter-layer altitude (m).
88
89	REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(OUT) :: dqyc ! Chemical species tendencies on GCM layers (mol/mol/s).
90
91	! Local variables :
92	! -----------------
93
94	INTEGER :: i , l, ic, ig, igm1
95
96	INTEGER :: dec, idec, ipres, ialt, klat
97
98	REAL*8 :: declin_c ! Declination (deg).
99	REAL*8 :: factp, factalt, factdec, factlat, krpddec, krpddecp1, krpddecm1
100	REAL*8 :: temp1, logp
101
102	! Variables sent into chemistry module (must be in double precision)
103	! ------------------------------------------------------------------
104
105	REAL*8, DIMENSION(nlaykim_tot) :: temp_c ! Temperature (K).
106	REAL*8, DIMENSION(nlaykim_tot) :: press_c ! Pressure (Pa).
107	REAL*8, DIMENSION(nlaykim_tot) :: phi_c ! Geopotential (m2.s-2) - actually not sent in chem. module but used to compute alts.
108	REAL*8, DIMENSION(nlaykim_tot) :: nb ! Density (cm-3).
109
110	REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy ! Chemical species in whole column (mol/mol) sent to chem. module.
111	REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy0 ! " " " " " " before modifs.
112
113	REAL*8 :: surfhaze(nlaykim_tot)
114	REAL*8 :: cprodaer(nlaykim_tot,4), cmaer(nlaykim_tot,4)
115	REAL*8 :: ccsn(nlaykim_tot,4), ccsh(nlaykim_tot,4)
116
117	REAL*8, DIMENSION(nlaykim_tot) :: rmil ! Mid-layer distance (km) to planetographic center.
118	REAL*8, DIMENSION(nlaykim_tot) :: rinter ! Inter-layer distance (km) to planetographic center (RA grid in chem. module).
119	! NB : rinter is on nlaykim_tot too, we don't care of the uppermost layer upper boundary altitude.
120
121	! Saved variables initialized at firstcall
122	! ----------------------------------------
123
124	LOGICAL, SAVE :: firstcall = .TRUE.
125	!$OMP THREADPRIVATE(firstcall)
126
127	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: kedd ! Eddy mixing coefficient for upper chemistry (cm^2.s-1)
128	!$OMP THREADPRIVATE(kedd)
129
130	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: mass ! Molar mass of the compunds (g.mol-1)
131	REAL*8, DIMENSION(:,:), ALLOCATABLE, SAVE :: md
132	!$OMP THREADPRIVATE(mass,md)
133
134	REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: r1d, ct1d, p1d, t1d ! Vervack profile
135	! JVO 18 : No threadprivate for those as they'll be read in tcp.ver by master
136
137	REAL*8, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: krpd ! Photodissociations rate table
138	REAL*8, DIMENSION(:,:) , ALLOCATABLE, SAVE :: krate ! Reactions rate ( photo + chem )
139	!$OMP THREADPRIVATE(krpd,krate)
140
141	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: nom_prod, nom_perte
142	INTEGER, DIMENSION(:,:), ALLOCATABLE, SAVE :: reactif
143	INTEGER, DIMENSION(:,:), ALLOCATABLE, SAVE :: prod
144	INTEGER, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: perte
145	!$OMP THREADPRIVATE(nom_prod,nom_perte,reactif,prod,perte)
146
147	! TEMPORARY : Dummy parameters without microphysics
148	! Here to keep the whole stuff running without modif chem. module
149	! ---------------------------------------------------------------
150
151	INTEGER :: utilaer(16)
152	INTEGER :: aerprod = 0
153	INTEGER :: htoh2 = 0
154
155	! -----------------------------------------------------------------------
156	! *************** I. Initialisations and Firstcall ******************
157	! -----------------------------------------------------------------------
158
159	IF (firstcall) THEN
160
161	PRINT*, 'CHIMIE, premier appel'
162
163	call check(nlaykim_tot,klev-15,nlrt_kim,nkim)
164
165	ALLOCATE(r1d(131))
166	ALLOCATE(ct1d(131))
167	ALLOCATE(p1d(131))
168	ALLOCATE(t1d(131))
169
170	ALLOCATE(md(nlaykim_tot,nkim))
171	ALLOCATE(mass(nkim))
172
173	ALLOCATE(kedd(nlaykim_tot))
174
175	ALLOCATE(krate(nlaykim_tot,nr_kim))
176	ALLOCATE(krpd(15,nd_kim+1,nlrt_kim,nlat_actfluxes))
177
178	ALLOCATE(nom_prod(nkim))
179	ALLOCATE(nom_perte(nkim))
180	ALLOCATE(reactif(5,nr_kim))
181	ALLOCATE(prod(200,nkim))
182	ALLOCATE(perte(2,200,nkim))
183
184	! 1. Read Vervack profile "tcp.ver", once for all
185	! -----------------------------------------------
186
187	!$OMP MASTER
188	OPEN(11,FILE=TRIM(datadir)//'/tcp.ver',STATUS='old')
189	READ(11,*)
190	DO i=1,131
191	READ(11,*) r1d(i), t1d(i), ct1d(i), p1d(i)
192	! For debug :
193	! -----------
194	! PRINT*, "tcp.ver", r1d(i), t1d(i), ct1d(i), p1d(i)
195	ENDDO
196	CLOSE(11)
197	!$OMP END MASTER
198	!$OMP BARRIER
199
200	! 2. Calculation of temp_c, densities and altitudes in planetary average
201	! ----------------------------------------------------------------------
202
203	! JVO18 : altitudes calculated in firstcall are just for diag, as I set kedd in pressure grid
204
205	! a. For GCM layers we just copy-paste
206
207	!$OMP MASTER
208	PRINT*,'Init chemistry : pressure, density, temperature ... :'
209	PRINT*,'level, rmil (km), press_c (mbar), nb (cm-3), temp_c(K)'
210	!$OMP END MASTER
211	!$OMP BARRIER
212
213	IF (ngrid.NE.1) THEN
214	DO l=1,klev
215	rinter(l) = (zlevmoy(l)+rad)/1000.0 ! km
216	rmil(l) = (zlaymoy(l)+rad)/1000.0 ! km
217	temp_c(l) = tmoy(l) ! K
218	phi_c(l) = phimoy(l) ! m2.s-2
219	press_c(l) = playmoy(l)/100. ! mbar
220	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
221	!$OMP MASTER
222	PRINT*, l, rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
223	!$OMP END MASTER
224	!$OMP BARRIER
225	ENDDO
226	rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000.
227	ELSE
228	DO l=1,klev
229	rinter(l) = (czlev(1,l)+rad)/1000.0 ! km
230	rmil(l) = (czlay(1,l)+rad)/1000.0 ! km
231	temp_c(l) = ctemp(1,l) ! K
232	phi_c(l) = cpphi(1,l) ! m2.s-2
233	press_c(l) = cplay(1,l)/100. ! mbar
234	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
235	PRINT*, l, rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
236	ENDDO
237	rinter(klev+1)=(czlev(1,klev+1)+rad)/1000.
238	ENDIF ! if ngrid.ne.1
239
240	! b. Extension in upper atmosphere with Vervack profile
241	! NB : Maybe the transition klev/klev+1 is harsh ...
242
243	ipres=1
244	DO l=klev+1,nlaykim_tot
245	press_c(l) = preskim(l-klev) / 100.0
246	DO i=ipres,130
247	IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
248	ipres=i
249	ENDIF
250	ENDDO
251	factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
252
253	nb(l) = exp( log(ct1d(ipres))(1-factp) + log(ct1d(ipres+1)) factp )
254	temp_c(l) = t1d(ipres)(1-factp) + t1d(ipres+1)factp
255	ENDDO
256
257	! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
258	! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
259
260	DO l=klev+1,nlaykim_tot
261
262	! Compute geopotential on the upper grid, to have correct altitude
263	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
264	! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
265	! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
266
267	temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
268	phi_c(l) = phi_c(l-1) + rtemp1log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
269
270	rmil(l) = ( gradrad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
271
272	!$OMP MASTER
273	PRINT*, l , rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
274	!$OMP END MASTER
275	!$OMP BARRIER
276	ENDDO
277
278	DO l=klev+2,nlaykim_tot
279	rinter(l) = 0.5*(rmil(l-1) + rmil(l))
280	ENDDO
281
282	! 3. Compounds caracteristics
283	! ---------------------------
284	mass(:) = 0.0
285	call comp(nomqy_c,nb,temp_c,mass,md)
286	PRINT*,' Mass'
287	DO ic=1,nkim
288	PRINT*, nomqy_c(ic), mass(ic)
289	ENDDO
290
291	! 4. Photodissociation rates
292	! --------------------------
293	call disso(krpd,nlat_actfluxes)
294
295	! 5. Init. chemical reactions with planetary average T prof.
296	! ----------------------------------------------------------
297
298	! NB : Chemical reactions rate are assumed to be constant within the T range of Titan's atm
299	! so we fill their krate once for all but krate for photodiss will be filled at each timestep
300
301	call chimie(nomqy_c,nb,temp_c,krate,reactif, &
302	nom_perte,nom_prod,perte,prod)
303
304	! 6. Eddy mixing coefficients (constant with time and space)
305	! ----------------------------------------------------------
306
307	kedd(:) = 1.e3 ! Default value =/= zero
308
309	! NB : Eddy coeffs from 1D models (e.g. Lavvas et al 08) in altitude but they're rather linked to pressure
310	! Below GCM top we have dynamic mixing and for levs < nld=klev-15 the chem. solver ignores diffusion
311
312	! First calculate kedd for upper chemistry layers
313	DO l=klev+1,nlaykim_tot
314	logp=-log10(press_c(l))
315	! 1E6 at 300 km ~ 10-1 mbar
316	IF ( logp.ge.1.0 .and. logp.le.4.0 ) THEN
317	kedd(l) = 10.*(6.0+1.3(logp-1.0)/3.0)
318	! 2E7 at 600 km ~ 10-4 mbar
319	ELSEIF ( logp.gt.4.0 .and. logp.le.6.0 ) THEN
320	kedd(l) = 10.*(7.3+0.35(logp-4.0))
321	! 1E8 above 900 km ~ 10-6 mbar
322	ELSEIF ( logp.gt.6.0 ) THEN
323	kedd(l) = 1.e8
324	ENDIF
325	ENDDO
326
327	! Then adjust 15 last layers profile fading to default value depending on kedd(ptop)
328	DO l=klev-15,klev
329	temp1 = ( log10(press_c(l)/press_c(klev-15)) ) / ( log10(press_c(klev+1)/press_c(klev-15)) )
330	kedd(l) = 10.*( 3.0 + log10(kedd(klev+1)/1.e3) temp1 )
331	ENDDO
332
333	firstcall = .FALSE.
334	ENDIF ! firstcall
335
336	declin_c = declin*180./pi
337
338	! -----------------------------------------------------------------------
339	! ********************* II. Loop on latitudes ***********************
340	! -----------------------------------------------------------------------
341
342	DO ig=1,ngrid
343
344	IF (ig.eq.1) THEN
345	igm1=1
346	ELSE
347	igm1=ig-1
348	ENDIF
349
350	! If 2D chemistry, trick to do the calculation only once per latitude band within the chunk
351	! NB1 : Will be obsolete with DYNAMICO, the chemistry will necessarly be 3D
352	! NB2 : Test of same latitude with dlat=0.5 : I think that if you run sims better than 1/10th degree then
353	! either it's with Dynamico and doesn't apply OR it is more than enough in terms of "preco / calc time" !
354	! -------------------------------------------------------------------------------------------------------
355
356	IF ( ( moyzon_ch .AND. ( ig.EQ.1 .OR. (ABS(latitude(ig)-latitude(igm1)).GT.0.1*pi/180.)) ) .OR. (.NOT. moyzon_ch) ) THEN
357
358	! 1. Compute altitude for the grid point with hydrostat. equilib.
359	! ---------------------------------------------------------------
360
361	! a. For GCM layers we just copy-paste
362
363	DO l=1,klev
364	rinter(l) = (czlev(ig,l)+rad)/1000.0 ! km
365	rmil(l) = (czlay(ig,l)+rad)/1000.0 ! km
366	temp_c(l) = ctemp(ig,l) ! K
367	phi_c(l) = cpphi(ig,l) ! m2.s-2
368	press_c(l) = cplay(ig,l)/100. ! mbar
369	nb(l) = 1.e-4press_c(l) / (kboltemp_c(l)) ! cm-3
370	ENDDO
371	rinter(klev+1)=(czlev(ig,klev+1)+rad)/1000.
372
373	! b. Extension in upper atmosphere with Vervack profile
374	! NB : The transition klev/klev+1 is maybe harsh if we have correct g above but not under
375
376	ipres=1
377	DO l=klev+1,nlaykim_tot
378	press_c(l) = preskim(l-klev) / 100.0
379	DO i=ipres,130
380	IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
381	ipres=i
382	ENDIF
383	ENDDO
384	factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
385
386	nb(l) = exp( log(ct1d(ipres))(1-factp) + log(ct1d(ipres+1)) factp )
387	temp_c(l) = t1d(ipres)(1-factp) + t1d(ipres+1)factp
388	ENDDO
389
390	! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
391	! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
392
393	DO l=klev+1,nlaykim_tot
394
395	! Compute geopotential on the upper grid, to have correct altitude
396	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
397	! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
398	! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
399
400	temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
401	phi_c(l) = phi_c(l-1) + rtemp1log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
402
403	rmil(l) = ( gradrad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
404	ENDDO
405
406	DO l=klev+2,nlaykim_tot
407	rinter(l) = 0.5*(rmil(l-1) + rmil(l))
408	ENDDO
409
410	! 2. From krpd, compute krate for dissociations (declination-latitude-altitude interpolation)
411	! -------------------------------------------------------------------------------------------
412
413	! a. Calculate declination dependence
414
415	if ((declin_c*10+267).lt.14.) then
416	idec = 0
417	dec = 0
418	else
419	if ((declin_c*10+267).gt.520.) then
420	idec = 14
421	dec = 534
422	else
423	idec = 1
424	dec = 27
425	do while( (declin_c*10+267).ge.real(dec+20) )
426	dec = dec+40
427	idec = idec+1
428	enddo
429	endif
430	endif
431	if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
432	factdec = ( declin_c - (dec-267)/10. ) / 4.
433	else
434	factdec = ( declin_c - (dec-267)/10. ) / 2.7
435	endif
436
437	! b. Calculate klat for interpolation on fixed latitudes of actinic fluxes input
438
439	klat=1
440	DO i=1,nlat_actfluxes
441	IF (latitude(ig).LT.lat_actfluxes(i)) klat=i
442	ENDDO
443	IF (klat==nlat_actfluxes) THEN ! avoid rounding problems
444	klat = nlat_actfluxes-1
445	factlat = 1.0
446	ELSE
447	factlat = (latitude(ig)-lat_actfluxes(klat))/(lat_actfluxes(klat+1)-lat_actfluxes(klat))
448	ENDIF
449
450	! c. Altitude loop
451
452	DO l=1,nlaykim_tot
453
454	! Calculate ialt for interpolation in altitude (krpd every 2 km)
455	ialt = int((rmil(l)-rad/1000.)/2.)+1
456	factalt = (rmil(l)-rad/1000.)/2.-(ialt-1)
457
458	! Altitude can go above top limit of UV levels - in this case we keep the 1310km top fluxes
459	IF (ialt.GT.nlrt_kim-1) THEN
460	ialt = nlrt_kim-1 ! avoid out-of-bound array
461	factalt = 1.0
462	ENDIF
463
464	DO i=1,nd_kim+1 ! nd_kim+1 is dissociation of N2 by GCR
465
466	krpddec = ( krpd(idec+1,i,ialt ,klat) * (1.0-factalt) &
467	+ krpd(idec+1,i,ialt+1,klat) * factalt ) * (1.0-factlat) &
468	* ( krpd(idec+1,i,ialt ,klat+1) * (1.0-factalt) &
469	+ krpd(idec+1,i,ialt+1,klat+1) * factalt ) * factlat
470
471	if ( factdec.lt.0. ) then
472	krpddecm1 = ( krpd(idec ,i,ialt ,klat) * (1.0-factalt) &
473	+ krpd(idec ,i,ialt+1,klat) * factalt ) * (1.0-factlat) &
474	* ( krpd(idec ,i,ialt ,klat+1) * (1.0-factalt) &
475	+ krpd(idec ,i,ialt+1,klat+1) * factalt ) * factlat
476	krate(l,i) = krpddecm1 * abs(factdec) + krpddec * ( 1.0 + factdec)
477	else if ( factdec.gt.0. ) then
478	krpddecp1 = ( krpd(idec+2,i,ialt ,klat) * (1.0-factalt) &
479	+ krpd(idec+2,i,ialt+1,klat) * factalt ) * (1.0-factlat) &
480	* ( krpd(idec+2,i,ialt ,klat+1) * (1.0-factalt) &
481	+ krpd(idec+2,i,ialt+1,klat+1) * factalt ) * factlat
482	krate(l,i) = krpddecp1 * factdec + krpddec * ( 1.0 - factdec)
483	else if ( factdec.eq.0. ) then
484	krate(l,i) = krpddec
485	endif
486
487	ENDDO ! i=1,nd_kim+1
488	ENDDO ! l=1,nlaykim_tot
489
490	! 3. Read composition
491	! -------------------
492
493	DO ic=1,nkim
494	DO l=1,klev
495	cqy(l,ic) = qy_c(ig,l,ic) ! advected tracers for the GCM part converted to molar frac.
496	ENDDO
497
498	DO l=1,nlaykim_up
499	cqy(klev+l,ic) = ykim_up(ic,ig,l) ! ykim_up for the upper atm.
500	ENDDO
501	ENDDO
502
503	cqy0(:,:) = cqy(:,:) ! Stores compo. before modifs
504
505	! 4. Call main Titan chemistry C routine
506	! --------------------------------------
507
508	call gptitan(rinter,temp_c,nb, &
509	nomqy_c,cqy, &
510	dtchim,latitude(ig)*180./pi,mass,md, &
511	kedd,botCH4,krate,reactif, &
512	nom_prod,nom_perte,prod,perte, &
513	aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh, &
514	htoh2,surfhaze)
515
516	! 5. Calculates tendencies on composition for advected tracers
517	! ------------------------------------------------------------
518	DO ic=1,nkim
519	DO l=1,klev
520	dqyc(ig,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! (mol/mol/s)
521	ENDDO
522	ENDDO
523
524	! 6. Update ykim_up
525	! -----------------
526	DO ic=1,nkim
527	DO l=1,nlaykim_up
528	ykim_up(ic,ig,l) = cqy(klev+l,ic)
529	ENDDO
530	ENDDO
531	! NB: The full vertical composition grid will be created only for the outputs
532
533
534	ELSE ! In 2D chemistry, if following grid point at same latitude, same zonal mean so don't do calculations again !
535	dqyc(ig,:,:) = dqyc(igm1,:,:) ! will be put back in 3D with longitudinal variations assuming same relative tendencies within a lat band
536	ykim_up(:,ig,:) = ykim_up(:,igm1,:) ! no horizontal mixing in upper layers -> no longitudinal variations
537	ENDIF
538
539	ENDDO
540
541	END SUBROUTINE calchim

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

Download in other formats:

Original Format