Context Navigation

photochemistry_asis.F90 @ 3567

Last change on this file since 3567 was 3312, checked in by yjaziri, 8 months ago

Generic PCM:

Photochemistry: correction for using diurnal equal true in 1D
No correction factor with diurnal equal true for photolysis rate
+ writediagspecUV change name to flux_surf for clarity

Property svn:executable set to *

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

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

Download in other formats:

Original Format