Context Navigation

photochemistry_asis.F90 @ 3893

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

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

Download in other formats:

Original Format