Context Navigation

photochemistry_asis.F90 @ 2595

Last change on this file since 2595 was 2553, checked in by yjaziri, 3 years ago

Generic GCM:

Chemistry: one input file instead of three for the chemical network
Automatic detection of mono/bi-molecular or quadratic reactions

Property svn:executable set to *

File size: 57.0 KB

Line
1	!*****************************************************************
2	!
3	! Photochemical routine
4	!
5	! Author: Franck Lefevre
6	! Benjamin Charnay
7	! Yassin Jaziri
8	! ------
9	!
10	! Version: 27/05/2016
11	! update 06/03/2021 generic tracer/network + photolysis online (Yassin Jaziri)
12	!
13	!*****************************************************************
14
15	subroutine photochemistry_asis(nlayer, ngrid, &
16	ig, lswitch, zycol, sza, fractcol, ptimestep, press, &
17	alt, temp, dens, zmmean, dist_sol, surfdust1d, &
18	surfice1d, tau, iter,zdtchim)
19
20	use callkeys_mod
21	use comcstfi_mod, only: r,cpp,avocado,pi
22	use tracer_h
23	use types_asis
24	use chimiedata_h
25	use photolysis_mod
26
27	implicit none
28
29	!===================================================================
30	! inputs:
31	!===================================================================
32
33	integer, intent(in) :: nlayer ! number of atmospheric layers
34	integer, intent(in) :: ngrid ! number of atmospheric columns
35
36	integer :: ig ! grid point index
37
38	real :: sza ! solar zenith angle (deg)
39	real :: fractcol ! day fraction
40	real :: ptimestep ! physics timestep (s)
41	real :: press(nlayer) ! pressure (hpa)
42	real :: alt(nlayer) ! altitude (km)
43	real :: temp(nlayer) ! temperature (k)
44	real :: dens(nlayer) ! density (cm-3)
45	real :: zmmean(nlayer) ! mean molar mass (g/mole)
46	real :: dist_sol ! sun distance (au)
47	real :: surfdust1d(nlayer) ! dust surface area (cm2/cm3)
48	real :: surfice1d(nlayer) ! ice surface area (cm2/cm3)
49	real :: tau ! optical depth at 7 hpa
50
51	!===================================================================
52	! input/output:
53	!===================================================================
54
55	real :: zycol(nlayer,nesp) ! chemical species volume mixing ratio
56
57	!===================================================================
58	! output:
59	!===================================================================
60
61	integer :: iter(nlayer) ! iteration counter
62	real :: zdtchim(nlayer) ! temperature modification by ozone
63
64	!===================================================================
65	! local:
66	!===================================================================
67
68	integer :: phychemrat ! (physical timestep)/(nominal chemical timestep)
69	integer :: ilev, iesp, iphot, iw
70	integer :: lswitch
71
72	logical, save :: firstcall = .true.
73
74	real :: stimestep ! standard timestep for the chemistry (s)
75	real :: ctimestep ! real timestep for the chemistry (s)
76	real :: dt_guess ! first-guess timestep (s)
77	real :: dt_corrected ! corrected timestep (s)
78	real :: dt_min = 1. ! minimum allowed timestep (s)
79	real :: dtg ! correction factor for the timestep (s)
80	real :: j(nlayer,nd) ! interpolated photolysis rates (s-1)
81	real :: time ! internal time (between 0 and ptimestep, in s)
82	real :: rho(nlayer) ! mass density (kg/m-3)
83
84	real, dimension(nlayer,nesp) :: rm ! mixing ratios
85	real (kind = 8), dimension(nesp) :: cold ! number densities at previous timestep (molecule.cm-3)
86	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities at current timestep (molecule.cm-3)
87	real (kind = 8), dimension(nesp) :: cnew ! number densities at next timestep (molecule.cm-3)
88
89	! reaction rates
90
91	real (kind = 8), allocatable, save :: v_phot(:,:)
92	real (kind = 8), allocatable, save :: v_3(:,:)
93	real (kind = 8), allocatable, save :: v_4(:,:)
94	real (kind = 8), allocatable, save :: e_phot(:,:) ! photolysis rates by energie (J.mol-1.s-1)
95
96	integer, save :: nreact ! number of reactions in reactions files
97
98	! matrix
99
100	real (kind = 8), dimension(nesp,nesp) :: mat, mat1
101	integer, dimension(nesp) :: indx
102	integer :: code
103
104	! production and loss terms (for first-guess solution only)
105
106	real (kind = 8), dimension(nesp) :: prod, loss
107
108	!===================================================================
109	! initialisation of the reaction indexes
110	!===================================================================
111
112	if (firstcall) then
113	print*,'photochemistry: initialize indexes'
114	call indice(nreact)
115	allocate(v_phot(nlayer,nb_phot_max))
116	allocate(v_3(nlayer,nb_reaction_3_max))
117	allocate(v_4(nlayer,nb_reaction_4_max))
118	allocate(v_phot_3d(ngrid,nlayer,nb_phot_hv_max))
119	allocate(e_phot(nlayer,nb_phot_max))
120	v_phot(:,:) = 0.
121	v_3(:,:) = 0.
122	v_4(:,:) = 0.
123	v_phot_3d(:,:,:) = 0.
124	e_phot(:,:) = 0.
125
126	! Need to be done after subroutine indice
127	if (jonline) then
128	print*,'calchim: Read UV absorption cross-sections'
129	! Jonline cannot run without photolysis
130	if (nb_phot_hv_max == 0) then
131	print*,'Jonline cannot run without photolysis'
132	print*,'set jonline=.false. in callphys.def'
133	print*,'or set photolysis reactions in monoreact file'
134	call abort
135	end if
136	call init_photolysis(nlayer,nreact) ! on-line photolysis
137	allocate(albedo_snow_chim(nw))
138	allocate(albedo_co2_ice_chim(nw))
139	! Step 1 : Initialisation of the Spectral Albedos.
140	DO iw=1,nw
141	albedo_snow_chim(iw)=albedosnow
142	albedo_co2_ice_chim(iw)=albedoco2ice
143
144	! Spectral Snow Albedo Calculation.
145	call albedo_snow_calc(wc(iw)*1.0e-3, &
146	albedo_snow_chim(iw), &
147	albedosnow)
148
149	ENDDO
150	end if
151
152	firstcall = .false.
153	end if
154
155	!===================================================================
156	! initialisation of mixing ratios and densities
157	!===================================================================
158
159	call gcmtochim(nlayer, zycol, lswitch, nesp, dens, rm, c)
160
161	!===================================================================
162	! interpolation of photolysis rates in the lookup table
163	!===================================================================
164
165	if (jonline) then
166	if (sza <= 95.) then
167	call photolysis_online(nlayer, alt, press, temp, zmmean, rm, &
168	tau, sza, dist_sol, v_phot, e_phot, ig, ngrid, nreact)
169	if (ngrid.eq.1) then
170	do iphot = 1,nb_phot_hv_max
171	v_phot(:,iphot) = v_phot(:,iphot)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
172	e_phot(:,iphot) = e_phot(:,iphot)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
173	end do
174	elseif(diurnal .eqv. .false.) then
175	do iphot = 1,nb_phot_hv_max
176	v_phot(:,iphot) = v_phot(:,iphot) * fractcol
177	e_phot(:,iphot) = e_phot(:,iphot) * fractcol
178	end do
179	endif
180	else ! night
181	v_phot(:,:) = 0.
182	e_phot(:,:) = 0.
183	end if
184	! save photolysis for output
185	v_phot_3d(ig,:,:) = v_phot(:,1:nb_phot_hv_max)
186	else if (nb_phot_hv_max /= 0) then
187	call photolysis_asis(nlayer, ngrid, lswitch, press, temp, sza, fractcol,tau, zmmean, dist_sol, &
188	rm(:,indexchim('co2')), rm(:,indexchim('o3')), rm(:,indexchim('ch4')), v_phot)
189	! save photolysis for output
190	v_phot_3d(ig,:,:) = v_phot(:,1:nb_phot_hv_max)
191	end if
192
193	!===================================================================
194	! reaction rates
195	!===================================================================
196
197	call reactionrates(nlayer, nreact, c, lswitch, dens, &
198	press, temp, zmmean, sza, surfdust1d, surfice1d, v_phot, v_3, v_4)
199
200	zdtchim(:) = 0.
201	rho(:) = (press(:)1e2)/(rtemp(:))
202
203	!===================================================================
204	! temperature modification by photolysis reaction
205	!===================================================================
206
207	if (photoheat .and. jonline) then ! for now just with jonline
208
209	do iphot = 1,nb_phot_hv_max
210	zdtchim(:nlayer-1) = zdtchim(:nlayer-1) + e_phot(:nlayer-1,iphot)c(:nlayer-1,indexchim(trim(jlabel(iphot,2))))/(cpp(rho(:nlayer-1)1e-6)avocado)
211	end do
212
213	else
214
215	!! o + o2 -> o3
216	!!dE = 24 ! kcal.mol-1
217	!!zdtchim(:) = zdtchim(:) + 4.18424e3v_4(:,1)c(:,indexchim('o'))c(:,indexchim('o2'))1e6/(cpprho*avocado)
218	!
219	!! o3 -> o2 + o1d
220	!! E(250nm) = 480 kJ.mol-1
221	!j_o3_o1d = 3
222	!zdtchim(:) = zdtchim(:) + 480e3v_phot(:,j_o3_o1d)c(:,indexchim('o3'))/(cpp(rho1e-6)*avocado)
223	!
224	!! o3 -> o2 + o
225	!! E(350nm) = 343 kJ.mol-1
226	!j_o3_o = 4
227	!zdtchim(:) = zdtchim(:) + 343e3v_phot(:,j_o3_o)c(:,indexchim('o3'))/(cpp(rho1e-6)*avocado)
228
229	zdtchim(:) = 0.
230
231	end if
232
233	!===================================================================
234	! stimestep : standard chemical timestep (s)
235	! ctimestep : real chemical timestep (s),
236	! taking into account the physical timestep
237	!===================================================================
238
239	stimestep = 600. ! standard value : 10 mn
240
241	phychemrat = nint(ptimestep/stimestep)
242	phychemrat = 1
243
244	ctimestep = ptimestep/real(phychemrat)
245
246	!===================================================================
247	! loop over levels
248	!===================================================================
249
250	do ilev = 1,lswitch - 1
251
252	! initializations
253
254	time = 0.
255	iter(ilev) = 0
256	dt_guess = ctimestep
257	cold(:) = c(ilev,:)
258
259	! internal loop for the chemistry
260
261	do while (time < ptimestep)
262
263	iter(ilev) = iter(ilev) + 1 ! iteration counter
264
265	! first-guess: fill matrix
266
267	call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
268
269	! adaptative evaluation of the sub time step
270
271	call define_dt(nesp, dt_corrected, dt_guess, ctimestep, cold(:), c(ilev,:), &
272	mat1, prod, loss, dens(ilev))
273
274	if (time + dt_corrected > ptimestep) then
275	dt_corrected = ptimestep - time
276	end if
277
278	! if (dt_corrected /= dt_guess) then ! the timestep has been modified
279
280	! form the matrix identity + mat*dt_corrected
281
282	mat(:,:) = mat1(:,:)*dt_corrected
283	do iesp = 1,nesp
284	mat(iesp,iesp) = 1. + mat(iesp,iesp)
285	end do
286
287	! solve the linear system M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
288
289	cnew(:) = c(ilev,:)
290
291	#ifdef LAPACK
292	call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
293	#else
294	write(,) "photochemistry_asis error, missing LAPACK routine dgesv"
295	call abort
296	#endif
297
298	! end if
299
300	! eliminate small values
301
302	where (cnew(:)/dens(ilev) < 1.e-30)
303	cnew(:) = 0.
304	end where
305
306	! update concentrations
307
308	cold(:) = c(ilev,:)
309	c(ilev,:) = cnew(:)
310	cnew(:) = 0.
311
312	! increment internal time
313
314	time = time + dt_corrected
315	dt_guess = dt_corrected ! first-guess timestep for next iteration
316
317	end do ! while (time < ptimestep)
318
319	end do ! ilev
320
321	!===================================================================
322	! save chemical species for the gcm
323	!===================================================================
324
325	call chimtogcm(nlayer, zycol, lswitch, nesp, dens, c)
326
327	contains
328
329	!================================================================
330
331	subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, &
332	prod, loss, dens)
333
334	!================================================================
335	! iterative evaluation of the appropriate time step dtnew
336	! according to curvature criterion based on
337	! e = 2 Rtol [r Cn+1 -(1-r) Cn + Cn-1 ]/[(1+r) Cn]
338	! with r = (tn - tn-1)/(tn+1 - tn)
339	!================================================================
340
341	implicit none
342
343	! input
344
345	integer :: nesp ! number of species in the chemistry
346
347	real :: dtold, ctimestep
348	real (kind = 8), dimension(nesp) :: cold, ccur
349	real (kind = 8), dimension(nesp,nesp) :: mat1
350	real (kind = 8), dimension(nesp) :: prod, loss
351	real :: dens
352
353	! output
354
355	real :: dtnew
356
357	! local
358
359	real (kind = 8), dimension(nesp) :: cnew
360	real (kind = 8), dimension(nesp,nesp) :: mat
361	real (kind = 8) :: atol, ratio, e, es, coef
362
363	integer :: code, iesp, iter
364	integer, dimension(nesp) :: indx
365
366	real :: dttest
367
368	! parameters
369
370	real (kind = 8), parameter :: dtmin = 10. ! minimum time step (s)
371	real (kind = 8), parameter :: vmrtol = 1.e-11 ! absolute tolerance on vmr
372	real (kind = 8), parameter :: rtol = 1./0.05 ! 1/rtol recommended value : 0.1-0.02
373	integer, parameter :: niter = 3 ! number of iterations
374	real (kind = 8), parameter :: coefmax = 2.
375	real (kind = 8), parameter :: coefmin = 0.1
376	logical :: fast_guess = .true.
377
378	dttest = dtold ! dttest = dtold = dt_guess
379
380	atol = vmrtol*dens ! absolute tolerance in molecule.cm-3
381
382	do iter = 1,niter
383
384	if (fast_guess) then
385
386	! first guess : fast semi-implicit method
387
388	do iesp = 1, nesp
389	cnew(iesp) = (ccur(iesp) + prod(iesp)dttest)/(1. + loss(iesp)dttest)
390	end do
391
392	else
393
394	! first guess : form the matrix identity + mat*dt_guess
395
396	mat(:,:) = mat1(:,:)*dttest
397	do iesp = 1,nesp
398	mat(iesp,iesp) = 1. + mat(iesp,iesp)
399	end do
400
401	! form right-hand side (RHS) of the system
402
403	cnew(:) = ccur(:)
404
405	! solve the linear system M*Cn+1 = Cn (RHS in cnew, then replaced by solution)
406
407	#ifdef LAPACK
408	call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
409	#else
410	write(,) "photochemistry_asis error, missing LAPACK routine dgesv"
411	call abort
412	#endif
413
414	end if
415
416	! ratio old/new subtimestep
417
418	ratio = dtold/dttest
419
420	! e : local error indicatocitr
421
422	e = 0.
423
424	do iesp = 1,nesp
425	es = 2.abs((ratiocnew(iesp) - (1. + ratio)*ccur(iesp) + cold(iesp)) &
426	/(1. + ratio)/max(ccur(iesp),atol))
427
428	if (es > e) then
429	e = es
430	end if
431	end do
432	e = rtol*e
433
434	! timestep correction
435
436	if (e <= 0.) then
437	coef = coefmax
438	else
439	coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
440	end if
441
442	dttest = max(dtmin,dttest*coef)
443	dttest = min(ctimestep,dttest)
444
445	end do ! iter
446
447	! new timestep
448
449	dtnew = dttest
450
451	end subroutine define_dt
452
453
454	!======================================================================
455
456	subroutine reactionrates(nlayer, nreact, c, &
457	lswitch, dens, press, t, zmmean, sza, &
458	surfdust1d, surfice1d, &
459	v_phot, v_3, v_4)
460
461	!================================================================
462	! compute reaction rates !
463	!----------------------------------------------------------------
464	! reaction type array !
465	!----------------------------------------------------------------
466	! A + B --> C + D bimolecular v_4 !
467	! A + A --> B + C quadratic v_3 !
468	! A + C --> B + C quenching v_phot !
469	! A + ice --> B + C heterogeneous v_phot !
470	!================================================================
471
472	use comcstfi_mod
473	use types_asis
474	use pfunc
475	use tracer_h
476	use chimiedata_h
477
478	implicit none
479
480	!----------------------------------------------------------------------
481	! input
482	!----------------------------------------------------------------------
483
484	integer, intent(in) :: nlayer ! number of atmospheric layers
485	integer, intent(in) :: nreact ! number of reactions in reactions files
486	integer :: lswitch ! interface level between lower
487	! atmosphere and thermosphere chemistries
488	real, dimension(nlayer) :: dens ! total number density (molecule.cm-3)
489	real, dimension(nlayer) :: press ! pressure (hPa)
490	real, dimension(nlayer) :: t ! temperature (K)
491	real, dimension(nlayer) :: zmmean ! mean molar mass (g/mole)
492	real :: sza ! solar zenith angle (deg)
493	real, dimension(nlayer) :: surfdust1d ! dust surface area (cm2.cm-3)
494	real, dimension(nlayer) :: surfice1d ! ice surface area (cm2.cm-3)
495	real (kind = 8), dimension(nlayer,nesp) :: c ! species number density (molecule.cm-3)
496
497	!----------------------------------------------------------------------
498	! output
499	!----------------------------------------------------------------------
500
501	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
502	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
503	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
504
505	!----------------------------------------------------------------------
506	! local
507	!----------------------------------------------------------------------
508
509	logical,save :: firstcall = .true.
510	integer :: ilev, ireact
511	integer :: nb_phot, nb_reaction_3, nb_reaction_4
512	integer :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
513	real, dimension(nlayer) :: ratec ! rate coefficient
514
515	!----------------------------------------------------------------------
516	! initialisation
517	!----------------------------------------------------------------------
518
519	nb_phot = nb_phot_hv_max
520	nb_reaction_3 = 0
521	nb_reaction_4 = 0
522
523	nb_hv = 0
524	nb_pfunc1 = 0
525	nb_pfunc2 = 0
526	nb_pfunc3 = 0
527
528	!----------------------------------------------------------------------
529	! reactions
530	!----------------------------------------------------------------------
531
532	do ireact = 1,nreact
533	if (reactions(ireact)%rtype==0) then
534	nb_hv = nb_hv + 1
535	cycle
536	end if
537
538	! get right coefficient rate function
539	if (reactions(ireact)%rtype==1) then
540	nb_pfunc1 = nb_pfunc1 + 1
541	if (reactions(ireact)%third_body==0) then !! third_body check
542	ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
543	else
544	ratec = pfunc1(nlayer,t,c(:,reactions(ireact)%third_body),pfunc1_param(nb_pfunc1))
545	end if
546	else if (reactions(ireact)%rtype==2) then
547	nb_pfunc2 = nb_pfunc2 + 1
548	if (reactions(ireact)%third_body==0) then !! third_body check
549	ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
550	else
551	ratec = pfunc2(nlayer,t,c(:,reactions(ireact)%third_body),pfunc2_param(nb_pfunc2))
552	end if
553	else if (reactions(ireact)%rtype==3) then
554	nb_pfunc3 = nb_pfunc3 + 1
555	if (reactions(ireact)%third_body==0) then !! third_body check
556	ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
557	else
558	ratec = pfunc3(nlayer,t,c(:,reactions(ireact)%third_body),pfunc3_param(nb_pfunc3))
559	end if
560	else
561	print*, 'Error in reactionrates: wrong index coefficient rate vphot'
562	print*, 'rtype(ireact) = ',reactions(ireact)%rtype
563	print*, 'It should be 1 or 2 '
564	print*, 'Please check the reaction ',ireact
565	if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
566	call abort
567	end if
568
569	! fill the right type index
570	if (reactions(ireact)%vtype=="v4") then
571	nb_reaction_4 = nb_reaction_4 + 1
572	v_4(:,nb_reaction_4) = ratec(:)
573	if (reactions(ireact)%three_prod) then ! three_prod check
574	nb_reaction_4 = nb_reaction_4 + 1
575	v_4(:,nb_reaction_4) = ratec(:)
576	end if
577	else if (reactions(ireact)%vtype=="v3") then
578	nb_reaction_3 = nb_reaction_3 + 1
579	v_3(:,nb_reaction_3) = ratec(:)
580	if (reactions(ireact)%three_prod) then ! three_prod check
581	nb_reaction_3 = nb_reaction_3 + 1
582	v_3(:,nb_reaction_3) = ratec(:)
583	end if
584	else if (reactions(ireact)%vtype=="vphot") then
585	nb_phot = nb_phot + 1
586	v_phot(:,nb_phot) = ratec(:)
587	if (reactions(ireact)%three_prod) then ! three_prod check
588	nb_phot = nb_phot + 1
589	v_phot(:,nb_phot) = ratec(:)
590	end if
591	else
592	print*, 'Error in reactionrates: wrong type coefficient rate'
593	print*, 'vtype(ireact) = ',reactions(ireact)%vtype
594	print*, 'It should be v4, v3 or vphot '
595	print*, 'Please check the reaction ',ireact
596	if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
597	call abort
598	end if
599
600	end do
601
602	call reactionrates_HC(nlayer,nesp,dens,t,press,zmmean,sza,c,v_phot,v_4,v_3,nb_phot,nb_reaction_4,nb_reaction_3)
603
604	!===========================================================
605	! check dimensions
606	!===========================================================
607
608	if (firstcall) then
609	if ((nb_phot /= nb_phot_max) .or. &
610	(nb_reaction_3 /= nb_reaction_3_max) .or. &
611	(nb_reaction_4 /= nb_reaction_4_max)) then
612	print*, 'nb_phot = ', nb_phot
613	print*, 'nb_reaction_4 = ', nb_reaction_4
614	print*, 'nb_reaction_3 = ', nb_reaction_3
615	print*, 'wrong dimensions in reactionrates'
616	print*, 'nb_phot_max = ', nb_phot_max
617	print*, 'nb_reaction_4_max = ', nb_reaction_4_max
618	print*, 'nb_reaction_3_max = ', nb_reaction_3_max
619	print*, '------ hard coding reaction ------'
620	print*, 'nphot_hard_coding = ', nphot_hard_coding
621	print*, 'n4_hard_coding = ', n4_hard_coding
622	print*, 'n3_hard_coding = ', n3_hard_coding
623	call abort
624	end if
625	if ((nb_hv /= nb_hv_max) .or. &
626	(nb_pfunc1 /= nb_pfunc1_max) .or. &
627	(nb_pfunc2 /= nb_pfunc2_max) .or. &
628	(nb_pfunc3 /= nb_pfunc3_max)) then
629	print*, 'nb_hv = ', nb_hv
630	print*, 'nb_pfunc1 = ', nb_pfunc1
631	print*, 'nb_pfunc2 = ', nb_pfunc2
632	print*, 'nb_pfunc3 = ', nb_pfunc3
633	print*, 'wrong dimensions in reactionrates'
634	print*, 'nb_hv_max = ', nb_hv_max
635	print*, 'nb_pfunc1_max = ', nb_pfunc1_max
636	print*, 'nb_pfunc2_max = ', nb_pfunc2_max
637	print*, 'nb_pfunc3_max = ', nb_pfunc3_max
638	call abort
639	end if
640	firstcall = .false.
641	end if
642
643	end subroutine reactionrates
644
645	!======================================================================
646
647	subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)
648
649	!======================================================================
650	! filling of the jacobian matrix
651	!======================================================================
652
653	use types_asis
654	use chimiedata_h
655
656	implicit none
657
658	! input
659
660	integer :: ilev ! level index
661	integer :: nesp ! number of species in the chemistry
662	integer, intent(in) :: nlayer ! number of atmospheric layers
663
664	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities
665	real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot
666	real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
667	real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4
668
669	! output
670
671	real (kind = 8), dimension(nesp,nesp), intent(out) :: mat ! matrix
672	real (kind = 8), dimension(nesp), intent(out) :: prod, loss
673
674	! local
675
676	integer :: iesp
677	integer :: ind_phot_2,ind_phot_4,ind_phot_6
678	integer :: ind_3_2,ind_3_4,ind_3_6
679	integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
680	integer :: iphot,i3,i4
681
682	real(kind = jprb) :: eps, eps_4 ! implicit/explicit coefficient
683
684	! initialisations
685
686	mat(:,:) = 0.
687	prod(:) = 0.
688	loss(:) = 0.
689
690	! photodissociations
691	! or reactions a + c -> b + c
692	! or reactions a + ice -> b + c
693
694	do iphot = 1,nb_phot_max
695
696	ind_phot_2 = indice_phot(iphot)%z2
697	ind_phot_4 = indice_phot(iphot)%z4
698	ind_phot_6 = indice_phot(iphot)%z6
699
700	mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
701	mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
702	mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
703
704	loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
705	prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3v_phot(ilev,iphot)c(ilev,ind_phot_2)
706	prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5v_phot(ilev,iphot)c(ilev,ind_phot_2)
707
708	end do
709
710	! reactions a + a -> b + c
711
712	do i3 = 1,nb_reaction_3_max
713
714	ind_3_2 = indice_3(i3)%z2
715	ind_3_4 = indice_3(i3)%z4
716	ind_3_6 = indice_3(i3)%z6
717
718	mat(ind_3_2,ind_3_2) = mat(ind_3_2,ind_3_2) + indice_3(i3)%z1v_3(ilev,i3)c(ilev,ind_3_2)
719	mat(ind_3_4,ind_3_2) = mat(ind_3_4,ind_3_2) - indice_3(i3)%z3v_3(ilev,i3)c(ilev,ind_3_2)
720	mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5v_3(ilev,i3)c(ilev,ind_3_2)
721
722	loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1v_3(ilev,i3)c(ilev,ind_3_2)
723	prod(ind_3_4) = prod(ind_3_4) + indice_3(i3)%z3v_3(ilev,i3)c(ilev,ind_3_2)*c(ilev,ind_3_2)
724	prod(ind_3_6) = prod(ind_3_6) + indice_3(i3)%z5v_3(ilev,i3)c(ilev,ind_3_2)*c(ilev,ind_3_2)
725
726	end do
727
728	! reactions a + b -> c + d
729
730	eps = 1.d-10
731
732	do i4 = 1,nb_reaction_4_max
733
734	ind_4_2 = indice_4(i4)%z2
735	ind_4_4 = indice_4(i4)%z4
736	ind_4_6 = indice_4(i4)%z6
737	ind_4_8 = indice_4(i4)%z8
738
739	eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
740	eps_4 = min(eps_4,1.0_jprb)
741
742	mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1v_4(ilev,i4)(1. - eps_4)*c(ilev,ind_4_4)
743	mat(ind_4_2,ind_4_4) = mat(ind_4_2,ind_4_4) + indice_4(i4)%z1v_4(ilev,i4)eps_4*c(ilev,ind_4_2)
744	mat(ind_4_4,ind_4_2) = mat(ind_4_4,ind_4_2) + indice_4(i4)%z3v_4(ilev,i4)(1. - eps_4)*c(ilev,ind_4_4)
745	mat(ind_4_4,ind_4_4) = mat(ind_4_4,ind_4_4) + indice_4(i4)%z3v_4(ilev,i4)eps_4*c(ilev,ind_4_2)
746	mat(ind_4_6,ind_4_2) = mat(ind_4_6,ind_4_2) - indice_4(i4)%z5v_4(ilev,i4)(1. - eps_4)*c(ilev,ind_4_4)
747	mat(ind_4_6,ind_4_4) = mat(ind_4_6,ind_4_4) - indice_4(i4)%z5v_4(ilev,i4)eps_4*c(ilev,ind_4_2)
748	mat(ind_4_8,ind_4_2) = mat(ind_4_8,ind_4_2) - indice_4(i4)%z7v_4(ilev,i4)(1. - eps_4)*c(ilev,ind_4_4)
749	mat(ind_4_8,ind_4_4) = mat(ind_4_8,ind_4_4) - indice_4(i4)%z7v_4(ilev,i4)eps_4*c(ilev,ind_4_2)
750
751	loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1v_4(ilev,i4)c(ilev,ind_4_4)
752	loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3v_4(ilev,i4)c(ilev,ind_4_2)
753	prod(ind_4_6) = prod(ind_4_6) + indice_4(i4)%z5v_4(ilev,i4)c(ilev,ind_4_2)*c(ilev,ind_4_4)
754	prod(ind_4_8) = prod(ind_4_8) + indice_4(i4)%z7v_4(ilev,i4)c(ilev,ind_4_2)*c(ilev,ind_4_4)
755
756	end do
757
758	end subroutine fill_matrix
759
760	!================================================================
761
762	subroutine indice(nreact)
763
764	!================================================================
765	! set the "indice" arrays used to fill the jacobian matrix !
766	!----------------------------------------------------------------
767	! reaction type !
768	!----------------------------------------------------------------
769	! A + hv --> B + C photolysis indice_phot !
770	! A + B --> C + D bimolecular indice_4 !
771	! A + A --> B + C quadratic indice_3 !
772	! A + C --> B + C quenching indice_phot !
773	! A + ice --> B + C heterogeneous indice_phot !
774	!================================================================
775
776	use types_asis
777	use datafile_mod
778	use ioipsl_getin_p_mod, only: getin_p
779	use tracer_h, only: nesp
780	use chimiedata_h
781	use callkeys_mod, only: jonline
782
783	implicit none
784
785	! output
786
787	integer, intent(out) :: nreact ! number of reactions in reactions files
788
789	! local
790
791	integer :: nb_phot, nb_reaction_3, nb_reaction_4
792	integer :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
793	integer :: iq, ireact
794
795	character(len = 128) :: reactfile ! reactions table file name
796	integer :: nbq ! number of species in reactions file
797	integer :: nlines ! number of lines in reactions file
798	integer :: pnreact ! number of reaction in photochemical reactions file
799	integer :: bnreact ! number of reaction in bimolecular reactions file
800	integer :: qnreact ! number of reaction in quadratic reactions file
801	integer :: specheck(nesp) ! to count the species of traceur.def in reactions file
802	integer :: specheckr(nesp) ! to count the reactant species of traceur.def in reactions file
803	integer :: specheckp(nesp) ! to count the product species of traceur.def in reactions file
804
805	nb_phot = 0
806	nb_reaction_3 = 0
807	nb_reaction_4 = 0
808	nb_phot_hv_max = 0
809
810	nb_hv = 0
811	nb_pfunc1 = 0
812	nb_pfunc2 = 0
813	nb_pfunc3 = 0
814
815	!===========================================================
816	! Read chemical reactions
817	!===========================================================
818
819	! Get number of reactions
820	nlines = 0
821	nreact = 0
822	pnreact = 0
823	bnreact = 0
824	qnreact = 0
825
826	write(,) "Read reaction file"
827	reactfile = "reactfile" ! default
828	call getin_p("reactfile",reactfile) ! default path
829	write(,) " reactfile = ",trim(reactfile)
830
831	call count_react(reactfile,nlines,nreact,nb_phot,nb_reaction_4,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
832
833	!!!!!!!!!!! Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!!!!!
834	nb_phot = nb_phot + nphot_hard_coding
835	nb_reaction_4 = nb_reaction_4 + n4_hard_coding
836	nb_reaction_3 = nb_reaction_3 + n3_hard_coding
837	!!!!!!!!!!! END Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!
838
839	write(,)'number of reaction in reaction files is = ',nreact
840	print*, 'nb_phot = ', nb_phot
841	print*, 'nb_reaction_4 = ', nb_reaction_4
842	print*, 'nb_reaction_3 = ', nb_reaction_3
843	print, '***************'
844	print*, 'nb_hv = ', nb_hv
845	print*, 'nb_pfunc1 = ', nb_pfunc1
846	print*, 'nb_pfunc2 = ', nb_pfunc2
847	print*, 'nb_pfunc3 = ', nb_pfunc3
848	nb_phot_max = nb_phot
849	nb_reaction_4_max = nb_reaction_4
850	nb_reaction_3_max = nb_reaction_3
851	nb_hv_max = nb_hv
852	nb_pfunc1_max = nb_pfunc1
853	nb_pfunc2_max = nb_pfunc2
854	nb_pfunc3_max = nb_pfunc3
855
856	! Allocate
857	allocate(indice_phot(nb_phot_max))
858	allocate(indice_4(nb_reaction_4_max))
859	allocate(indice_3(nb_reaction_3_max))
860	allocate(reactions(nreact))
861	allocate(jlabel(nb_phot_hv_max,2))
862	allocate(jparamline(nb_hv_max))
863	allocate(pfunc1_param(nb_pfunc1_max))
864	allocate(pfunc2_param(nb_pfunc2_max))
865	allocate(pfunc3_param(nb_pfunc3_max))
866
867	! Get reactants, products and number of species in reactfile
868	! inititialize variables
869	nbq = 0
870	specheck(:) = 0
871	specheckr(:) = 0
872	specheckp(:) = 0
873	reactions(:) = reaction('',-1,0,.false.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
874	pfunc1_param(:) = rtype1(0.,0.,0.,0.,0.)
875	pfunc2_param(:) = rtype2(0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
876	pfunc3_param(:) = rtype3(0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
877	nb_phot = 0
878	nb_reaction_4 = 0
879	nb_reaction_3 = 0
880	jlabel(:,:) = ''
881	jparamline(:) = ''
882
883	call get_react(reactfile,nlines,nreact,specheck,specheckr,specheckp,nbq,nb_phot,nb_reaction_4,nb_reaction_3)
884
885	! rewrite jlabel corretly for output
886	do ireact=1,nb_phot_hv_max
887	if (reactions(ireact)%three_prod) then
888	jlabel(ireact+1:nb_phot_hv_max,2) = jlabel(ireact:nb_phot_hv_max-1,2)
889	jlabel(ireact+1:nb_phot_hv_max,1) = jlabel(ireact:nb_phot_hv_max-1,1)
890	jlabel(ireact,1) = trim(jlabel(ireact,1))//'_a'
891	jlabel(ireact+1,1) = trim(jlabel(ireact+1,1))//'_b'
892	end if
893	end do
894
895	!!!!!!!!!!! Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!!!!!
896	call indice_HC(nb_phot,nb_reaction_4,nb_reaction_3)
897	!!!!!!!!!!! END Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!
898
899	write(,)'number of species in reaction files is = ',nbq
900	write(,)'species in reaction files are:'
901
902	do iq=1,nesp
903	if (specheck(iq)==1) print*, chemnoms(iq)
904	end do
905
906	!===========================================================
907	! check species only destroy or product
908	!===========================================================
909
910	do iq=1,nesp
911	if (specheckr(iq)/=specheckp(iq)) then
912	if (specheckr(iq)==0 .and. specheckp(iq)==1) then
913	print*, 'WARNING: ', chemnoms(iq),' is only product'
914	else if (specheckr(iq)==1 .and. specheckp(iq)==0) then
915	print*, 'WARNING: ', chemnoms(iq),' is only destroy'
916	else
917	print*, 'Error in indice: bad value in specheckr or specheckp'
918	call abort
919	end if
920	end if
921	end do
922
923	!===========================================================
924	! check stochiometry
925	!===========================================================
926
927	! If you found a way
928
929	!===========================================================
930	! check dimensions
931	!===========================================================
932
933	if (jonline) then
934	nd = nb_hv_max
935	else if (nb_phot_hv_max /= 0) then
936	print*,'calchim: Read photolysis lookup table'
937	call read_phototable
938	end if
939
940	if ((nb_phot /= nb_phot_max) .or. &
941	(nb_reaction_3 /= nb_reaction_3_max) .or. &
942	(nb_reaction_4 /= nb_reaction_4_max) .or. &
943	(nd /= nb_hv_max)) then
944	print*, 'nb_phot = ', nb_phot
945	print*, 'nb_reaction_4 = ', nb_reaction_4
946	print*, 'nb_reaction_3 = ', nb_reaction_3
947	print*, 'nd = ', nd
948	print*, 'wrong dimensions in indice'
949	print*, 'nb_phot_max = ', nb_phot_max
950	print*, 'nb_reaction_4_max = ', nb_reaction_4_max
951	print*, 'nb_reaction_3_max = ', nb_reaction_3_max
952	print*, 'nb_phot_hv_max = ', nb_phot_hv_max
953	print*, 'nb_hv_max = ', nb_hv_max
954	call abort
955	end if
956
957	end subroutine indice
958
959	!*****************************************************************
960
961	subroutine gcmtochim(nlayer, zycol, lswitch, nesp, dens, rm, c)
962
963	!*****************************************************************
964
965	use callkeys_mod
966
967	implicit none
968
969
970	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
971	! input:
972	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
973
974	integer, intent(in) :: nlayer ! number of atmospheric layers
975	integer, intent(in) :: nesp ! number of species in the chemistry
976	integer :: lswitch ! interface level between chemistries
977
978	real :: zycol(nlayer,nesp) ! volume mixing ratios in the gcm
979	real :: dens(nlayer) ! total number density (molecule.cm-3)
980
981	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
982	! output:
983	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
984
985	real, dimension(nlayer,nesp) :: rm ! volume mixing ratios
986	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities (molecule.cm-3)
987
988	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
989	! local:
990	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
991
992	integer :: l, iesp
993
994	rm(:,:) = 0.
995
996	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
997	! initialise mixing ratios
998	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
999
1000	do l = 1,lswitch-1
1001	rm(l,:) = zycol(l,:)
1002	end do
1003
1004	where (rm(:,:) < 1.e-30)
1005	rm(:,:) = 0.
1006	end where
1007
1008	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1009	! initialise number densities
1010	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1011
1012	do iesp = 1,nesp
1013	do l = 1,lswitch-1
1014	c(l,iesp) = rm(l,iesp)*dens(l)
1015	end do
1016	end do
1017
1018	end subroutine gcmtochim
1019
1020	!*****************************************************************
1021
1022	subroutine chimtogcm(nlayer, zycol, lswitch, nesp, dens, c)
1023
1024
1025	!*****************************************************************
1026
1027	use callkeys_mod
1028
1029	implicit none
1030
1031
1032	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1033	! inputs:
1034	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1035
1036	integer, intent(in) :: nlayer ! number of atmospheric layers
1037	integer, intent(in) :: nesp ! number of species in the chemistry
1038	integer :: lswitch ! interface level between chemistries
1039
1040	real :: dens(nlayer) ! total number density (molecule.cm-3)
1041	real (kind = 8), dimension(nlayer,nesp) :: c ! number densities (molecule.cm-3)
1042
1043	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1044	! output:
1045	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1046
1047	real zycol(nlayer,nesp) ! volume mixing ratios in the gcm
1048
1049	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1050	! local:
1051	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1052
1053	integer l
1054
1055	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1056	! save mixing ratios for the gcm
1057	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1058
1059	do l = 1,lswitch-1
1060	zycol(l,:) = c(l,:)/dens(l)
1061	end do
1062
1063	end subroutine chimtogcm
1064
1065	!*****************************************************************
1066
1067	subroutine split_str(line,words,n,nmax)
1068
1069	!*****************************************************************
1070
1071	implicit none
1072	character(*), intent(in) :: line
1073	integer, intent(in) :: nmax
1074	character(*), intent(out) :: words(nmax)
1075	integer, intent(out) :: n ! number of words at the end
1076	integer :: ios
1077	character(100) :: buf(100) ! large buffer!
1078
1079	n = 0
1080	do
1081	n = n + 1
1082	if (n>nmax) then
1083	print*,'Error in split_str: to much words'
1084	print*,'limit = ',nmax
1085	print*,'change it, if you want, with increasing nmax'
1086	print*,'line is:',line
1087	call abort
1088	end if
1089	read(line,*,iostat=ios) buf(1:n) ! use list-directed input
1090	if (ios==0) then
1091	words(1:n) = buf(1:n) ! if success, copy to the original array
1092	else
1093	n = n-1
1094	exit ! if all the words are obtained, finish
1095	endif
1096	enddo
1097	end subroutine split_str
1098
1099	!*****************************************************************
1100
1101	subroutine find_vtype(reactline,vtype)
1102
1103	!*****************************************************************
1104
1105	use tracer_h, only: nesp
1106	use chimiedata_h, only: indexchim
1107
1108	implicit none
1109
1110	character(len = 50), intent(in) :: reactline ! all reactants of one reaction
1111	character(*), intent(inout) :: vtype ! "v3", "v4" or "vphot"
1112
1113	integer :: nwr ! number of reactant for a reaction
1114	integer,parameter :: nmax = 5 ! number max of reactants or products
1115	character(len = 24) :: words(nmax) ! to get words in reactants and products line
1116	integer :: stochio(nesp) ! stoichiometry of reaction
1117	integer :: iword
1118
1119	! init
1120	stochio(:) = 0
1121
1122	! split reactant
1123	call split_str(reactline,words,nwr,nmax)
1124
1125	! set stochio
1126	do iword = 1,nwr
1127	if (trim(words(iword))=="M") exit
1128	if (trim(words(iword))/="hv" &
1129	.and. trim(words(iword))/="dummy") stochio(indexchim(words(iword))) = stochio(indexchim(words(iword))) + 1
1130	end do
1131
1132	! set vtype
1133	if (sum(stochio)==1) then
1134	vtype = "vphot"
1135	else if (sum(stochio)==2) then
1136	if (any(stochio==2)) then
1137	vtype = "v3"
1138	else
1139	vtype = "v4"
1140	end if
1141	else if (sum(stochio)==3) then
1142	if (any(stochio==2)) then
1143	vtype = "v4"
1144	else
1145	print*,'Error in photochemistry_asis (find_vtype):'
1146	print*,'3 different reactants not implemented'
1147	call abort
1148	end if
1149	else
1150	print*,'Error in photochemistry_asis (find_vtype):'
1151	print*,'more than 2 different reactants not implemented'
1152	call abort
1153	end if
1154
1155	end subroutine find_vtype
1156
1157	!*****************************************************************
1158
1159	subroutine count_react(reactfile,nlines,nreact,nb_phot,nb_reaction_4,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)
1160
1161	!*****************************************************************
1162
1163	use types_asis, only : nb_phot_hv_max
1164
1165	implicit none
1166	character(*), intent(in) :: reactfile ! name of the file to read
1167	integer, intent(inout) :: nlines ! number of lines in the file
1168	integer, intent(out) :: nreact ! real number of reaction
1169	integer, intent(inout) :: nb_phot ! dimension of "vphot" reaction
1170	integer, intent(inout) :: nb_reaction_4 ! dimension of "v4" reaction
1171	integer, intent(inout) :: nb_reaction_3 ! dimension of "v3" reaction
1172	integer, intent(inout) :: nb_hv ! number of typerate 0 reaction
1173	integer, intent(inout) :: nb_pfunc1 ! number of typerate 1 reaction
1174	integer, intent(inout) :: nb_pfunc2 ! number of typerate 2 reaction
1175	integer, intent(inout) :: nb_pfunc3 ! number of typerate 3 reaction
1176
1177	! local
1178	character(len = 50) :: reactline ! all reactants of one reaction
1179	character(len = 50) :: prodline ! all produts of one reaction
1180	integer :: typerate ! get the type of the rate reaction coefficient (pfunc_i)
1181	character(5) :: vtype ! "v3", "v4" or "vphot"
1182	integer :: nwp ! number of products for a reaction
1183	integer,parameter :: nmax = 5 ! number max of reactants or products
1184	character(len = 24) :: words(nmax) ! to get words in reactants and products line
1185	integer :: ierr
1186
1187	nreact = 0
1188	vtype = ''
1189
1190	open(402, form = 'formatted', status = 'old', &
1191	file ='chemnetwork/'//trim(reactfile),iostat=ierr)
1192
1193	if (ierr /= 0) THEN
1194	write(,)'Error : cannot open chemical reaction file : chemnetwork/'//trim(reactfile)
1195	write(,)'It should be in your simulation directory in chemnetwork directory'
1196	write(,)' You can change the input chemical reactions file name in'
1197	write(,)' callphys.def with:'
1198	write(,)' monoreact=filename or bimolreact=filename or quadrareact=filename'
1199	write(,)' depending on what chemical reaction type it is'
1200	call abort
1201	end if
1202
1203	do
1204	read(402,'(A,A,I2)',iostat=ierr) reactline, prodline, typerate
1205	if (ierr<0) exit
1206	! count the number of lines
1207	nlines = nlines + 1
1208	if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
1209	! count the number of reaction
1210	nreact = nreact + 1
1211
1212	call find_vtype(reactline,vtype)
1213	call split_str(prodline,words,nwp,nmax)
1214
1215	! count the dimension of each rate reaction coefficient type
1216	if (trim(vtype)=="vphot") then
1217	nb_phot = nb_phot + 1
1218	! check three product reaction
1219	if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
1220	.and. trim(words(1))/=trim(words(3)) &
1221	.and. trim(words(2))/=trim(words(3))) nb_phot = nb_phot + 1
1222	else if (trim(vtype)=="v4") then
1223	nb_reaction_4 = nb_reaction_4 + 1
1224	! check three product reaction
1225	if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
1226	.and. trim(words(1))/=trim(words(3)) &
1227	.and. trim(words(2))/=trim(words(3))) nb_reaction_4 = nb_reaction_4 + 1
1228	else if (trim(vtype)=="v3") then
1229	nb_reaction_3 = nb_reaction_3 + 1
1230	! check three product reaction
1231	if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
1232	.and. trim(words(1))/=trim(words(3)) &
1233	.and. trim(words(2))/=trim(words(3))) nb_reaction_3 = nb_reaction_3 + 1
1234	else
1235	print*,'Error in photochemistry_asis (count_react):'
1236	print*,'vtype not found'
1237	call abort
1238	end if
1239
1240	! count the number of each rate reaction coefficient type
1241	if (typerate==0) then
1242	nb_hv = nb_hv + 1
1243	nb_phot_hv_max = nb_phot_hv_max + 1
1244	if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
1245	.and. trim(words(1))/=trim(words(3)) &
1246	.and. trim(words(2))/=trim(words(3))) nb_phot_hv_max = nb_phot_hv_max + 1
1247	else if (typerate==1) then
1248	nb_pfunc1 = nb_pfunc1 + 1
1249	else if (typerate==2) then
1250	nb_pfunc2 = nb_pfunc2 + 1
1251	else if (typerate==3) then
1252	nb_pfunc3 = nb_pfunc3 + 1
1253	else
1254	print*, 'Error in indice: wrong index coefficient rate line ',nlines
1255	print*, 'in file : chemnetwork/'//trim(reactfile)
1256	print*, 'It should be 0 for photolysis reations and 1 or 2 for the others'
1257	print*, 'And not : ', typerate
1258	call abort
1259	end if
1260
1261	end if
1262
1263	end do
1264
1265	close(402)
1266
1267	end subroutine count_react
1268
1269	!*****************************************************************
1270
1271	subroutine get_react(reactfile,nlines,nreact,specheck,specheckr,specheckp, &
1272	nbq,nb_phot,nb_reaction_4,nb_reaction_3)
1273
1274	!*****************************************************************
1275
1276	use types_asis
1277	use tracer_h
1278	use chimiedata_h, only: indexchim
1279	use callkeys_mod, only: jonline
1280
1281	implicit none
1282	character(*), intent(in) :: reactfile ! name of the file to read
1283	integer, intent(in) :: nlines ! number of lines in the file
1284	integer, intent(in) :: nreact ! real number of reaction in the file
1285	integer, intent(out) :: nb_phot ! number of "vphot" calculation reactions
1286	integer, intent(out) :: nb_reaction_4 ! number of "vphot" calculation reactions
1287	integer, intent(out) :: nb_reaction_3 ! number of "vphot" calculation reactions
1288	integer, intent(inout) :: specheck(nesp) ! to count the species of traceur.def in reactions file
1289	integer, intent(inout) :: specheckr(nesp) ! to count the reactant species of traceur.def in reactions file
1290	integer, intent(inout) :: specheckp(nesp) ! to count the product species of traceur.def in reactions file
1291	integer, intent(inout) :: nbq ! number of species in reactions file
1292
1293	! local
1294	character(len = 50) :: reactline ! all reactants of one reaction
1295	character(len = 50) :: prodline ! all produts of one reaction
1296	integer :: nwr ! number of reactants for a reaction
1297	integer :: nwp ! number of products for a reaction
1298	integer,parameter :: nmax = 5 ! number max of reactants or products
1299	character(len = 24) :: words(nmax) ! to get words in reactants and products line
1300	real :: nindice(2*nmax) ! stoichiometry of species (for indice variables)
1301	integer :: iindice(2*nmax) ! indice of species (for indice variables)
1302	character(len = 500) :: paramline ! line of reactions parameters
1303	integer :: stochio(nesp) ! stoichiometry of reaction
1304	integer :: ierr, ilines, ireact, i_dummy, iw, ispe, iloc
1305	integer :: nb_phot_hv, nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
1306
1307	i_dummy = 1
1308	ireact = 0
1309	nb_phot = 0
1310	nb_phot_hv = 0
1311	nb_hv = 0
1312	nb_pfunc1 = 0
1313	nb_pfunc2 = 0
1314	nb_pfunc3 = 0
1315
1316	open(402, form = 'formatted', status = 'old', &
1317	file ='chemnetwork/'//trim(reactfile),iostat=ierr)
1318
1319	if (ierr /= 0) THEN
1320	write(,)'Error : cannot open chemical reaction file : chemnetwork/'//trim(reactfile)
1321	write(,)'It should be in your simulation directory in chemnetwork directory'
1322	write(,)' You can change the input chemical reactions file name in'
1323	write(,)' callphys.def with:'
1324	write(,)' monoreact=filename or bimolreact=filename or quadrareact=filename'
1325	write(,)' depending on what chemical reaction type it is'
1326	call abort
1327	end if
1328
1329	do ilines=1,nlines
1330	paramline = ''
1331
1332	read(402,'(A,A,A)') reactline, prodline, paramline
1333
1334	! continue only if it's not a comment line
1335	if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
1336
1337	! increment number of reactions
1338	ireact = ireact + 1
1339
1340	! fill reaction vtype
1341	call find_vtype(reactline,reactions(ireact)%vtype)
1342	if (trim(reactions(ireact)%vtype)=="v4") then
1343	nb_reaction_4 = nb_reaction_4 + 1
1344	else if (trim(reactions(ireact)%vtype)=="v3") then
1345	nb_reaction_3 = nb_reaction_3 + 1
1346	else if (trim(reactions(ireact)%vtype)=="vphot") then
1347	nb_phot = nb_phot + 1
1348	else
1349	print*,'Error in photochemistry_asis (get_react):'
1350	print*,'vtype not found'
1351	call abort
1352	end if
1353
1354	! fill reaction rate type and parameters
1355	if (trim(paramline)=='') then
1356	print*, 'Error in reactfile: where are the parameters - line',ilines
1357	call abort
1358	else
1359	read(paramline,*) reactions(ireact)%rtype
1360	if (reactions(ireact)%rtype==1) then
1361	nb_pfunc1 = nb_pfunc1 + 1
1362	read(paramline,*) reactions(ireact)%rtype, pfunc1_param(nb_pfunc1)
1363	else if (reactions(ireact)%rtype==0) then
1364	nb_hv = nb_hv + 1
1365	nb_phot_hv = nb_phot_hv + 1
1366	if (jonline) then
1367	read(paramline,'(I5,A,A)') reactions(ireact)%rtype, jlabel(nb_hv,1), jparamline(nb_hv)
1368	else
1369	read(paramline,*) reactions(ireact)%rtype, jlabel(nb_hv,1)
1370	end if
1371	else if (reactions(ireact)%rtype==2) then
1372	nb_pfunc2 = nb_pfunc2 + 1
1373	read(paramline,*) reactions(ireact)%rtype, pfunc2_param(nb_pfunc2)
1374	else if (reactions(ireact)%rtype==3) then
1375	nb_pfunc3 = nb_pfunc3 + 1
1376	read(paramline,*) reactions(ireact)%rtype, pfunc3_param(nb_pfunc3)
1377	end if
1378	end if
1379
1380	! WARNING: the implementation is limited to 3 different reactants (for now only 2) and 3 different products
1381	nindice(:) = 0.0 ! 3 first indice for reactants, 3 following for products
1382	iindice(:) = i_dummy ! 3 first indice for reactants, 3 following for products
1383
1384	!-----------------!
1385	!--- reactants ---!
1386	!-----------------!
1387
1388	! init for reactant
1389	stochio(:) = 0
1390	! split reactant
1391	call split_str(reactline,words,nwr,nmax)
1392	! set species that are photolysed
1393	if (reactions(ireact)%rtype==0) jlabel(nb_hv,2) = words(1)
1394	! find reactant stochio
1395	do iw = 1,nwr
1396	! fill third body index
1397	if (trim(words(iw))=="M") then
1398	if (iw==nwr) then
1399	exit
1400	else if (iw==nwr-1) then
1401	reactions(ireact)%third_body = indexchim(words(iw+1))
1402	exit
1403	else
1404	print*, 'Error in reactfile: just only one specie can be after M corresponding to the third body - line',ilines
1405	call abort
1406	end if
1407	end if
1408	! count stochio
1409	if (trim(words(iw))/="hv" &
1410	.and. trim(words(iw))/="dummy") stochio(indexchim(words(iw))) = stochio(indexchim(words(iw))) + 1
1411	end do
1412	! get reaction reactant stochio and indice
1413	iloc = 0
1414	do ispe = 1,nesp
1415	if (stochio(ispe)==0) cycle
1416	iloc = iloc + 1
1417	nindice(iloc) = stochio(ispe)
1418	iindice(iloc) = ispe
1419	! check up: to count the species used in 'reactfile'
1420	if (specheck(ispe)==0) then
1421	specheckr(ispe) = 1
1422	specheck(ispe) = 1
1423	nbq = nbq + 1
1424	else if (specheckr(ispe)==0) then
1425	specheckr(ispe) = 1
1426	end if
1427	end do
1428	! fill reaction reactant stochio and indice
1429	reactions(ireact)%ir1 = nindice(1)
1430	reactions(ireact)%r1 = iindice(1)
1431	reactions(ireact)%ir2 = nindice(2)
1432	reactions(ireact)%r2 = iindice(2)
1433	reactions(ireact)%ir3 = nindice(3)
1434	reactions(ireact)%r3 = iindice(3)
1435
1436	!----------------!
1437	!--- products ---!
1438	!----------------!
1439
1440	! init for products
1441	stochio(:) = 0
1442	! split products
1443	call split_str(prodline,words,nwp,nmax)
1444	! find products stochio
1445	do iw = 1,nwp
1446	stochio(indexchim(words(iw))) = stochio(indexchim(words(iw))) + 1
1447	end do
1448	! get reaction product stochio and indice
1449	iloc = 3
1450	do ispe = 1,nesp
1451	if (stochio(ispe)==0) cycle
1452	iloc = iloc + 1
1453	nindice(iloc) = stochio(ispe)
1454	iindice(iloc) = ispe
1455	! check up: to count the species used in 'reactfile'
1456	if (specheck(ispe)==0) then
1457	specheckp(ispe) = 1
1458	specheck(ispe) = 1
1459	nbq = nbq + 1
1460	else if (specheckp(ispe)==0) then
1461	specheckp(ispe) = 1
1462	end if
1463	end do
1464	! fill reaction product stochio and indice
1465	reactions(ireact)%ip1 = nindice(4)
1466	reactions(ireact)%p1 = iindice(4)
1467	reactions(ireact)%ip2 = nindice(5)
1468	reactions(ireact)%p2 = iindice(5)
1469	reactions(ireact)%ip3 = nindice(6)
1470	reactions(ireact)%p3 = iindice(6)
1471	! set reaction three prod to true with checking in position 6 if there is three different products
1472	if (nindice(6)/=0) reactions(ireact)%three_prod = .true.
1473
1474
1475	!-------------------------!
1476	!--- fill indice array ---!
1477	!-------------------------!
1478
1479	if (trim(reactions(ireact)%vtype)=="v4") then
1480	if (nindice(6)==0) then ! reaction with 1 or 2 products
1481	! z4spec(ir1,r1,ir2,r2,ip1,p1,ip2,p2)
1482	indice_4(nb_reaction_4) = z4spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(4), iindice(4), nindice(5), iindice(5))
1483	else ! reaction with 3 different products
1484	indice_4(nb_reaction_4) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(4), iindice(4), nindice(5)/2., iindice(5))
1485	nb_reaction_4 = nb_reaction_4 + 1
1486	indice_4(nb_reaction_4) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(6), iindice(6), nindice(5)/2., iindice(5))
1487	end if
1488	else if (trim(reactions(ireact)%vtype)=="v3") then
1489	if (nindice(6)==0) then ! reaction with 1 or 2 products
1490	! z3spec(ir1,r1,ip1,p1,ip2,p2)
1491	indice_3(nb_reaction_3) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
1492	else ! reaction with 3 different products
1493	indice_3(nb_reaction_3) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
1494	nb_reaction_3 = nb_reaction_3 + 1
1495	indice_3(nb_reaction_3) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
1496	end if
1497	else if (trim(reactions(ireact)%vtype)=="vphot") then
1498	if (reactions(ireact)%rtype==0) then
1499	if (nindice(6)==0) then ! reaction with 1 or 2 products
1500	! z3spec(ir1,r1,ip1,p1,ip2,p2)
1501	indice_phot(nb_phot_hv) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
1502	else ! reaction with 3 different products
1503	indice_phot(nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
1504	nb_phot = nb_phot + 1
1505	nb_phot_hv = nb_phot_hv + 1
1506	indice_phot(nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
1507	end if
1508	else
1509	if (nindice(6)==0) then ! reaction with 1 or 2 products
1510	! z3spec(ir1,r1,ip1,p1,ip2,p2)
1511	indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
1512	else ! reaction with 3 different products
1513	indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
1514	nb_phot = nb_phot + 1
1515	indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
1516	end if
1517	end if
1518	else
1519	print*,'Error in photochemistry_asis (get_react):'
1520	print*,'vtype',reactions(ireact)%vtype,' not implemented'
1521	call abort
1522	end if
1523
1524	end if ! end if comment line
1525
1526	end do ! end loop on lines
1527
1528	close(402)
1529
1530	end subroutine get_react
1531
1532	end subroutine photochemistry_asis

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

Download in other formats:

Original Format