Context Navigation

calchim.F90 @ 2365

Last change on this file since 2365 was 2326, checked in by jvatant, 5 years ago
Update Titan reference photochemistry (reaction constants,branching ratios, condensation rates) according to Vuitton et al 2019. --JVO
File size: 22.5 KB

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

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

Download in other formats:

Original Format