Context Navigation

photolysis_mod.F90 @ 3567

Last change on this file since 3567 was 3370, checked in by yjaziri, 7 months ago

Generic PCM:

Change directory path for photochem stellar spectra
Just to organize datagcm
Now in datagcm/stellar_spectra/photochem_stellar_spectra/

File size: 32.9 KB

Line
1	module photolysis_mod
2
3	!***********************************************************************
4	!
5	! subject:
6	! --------
7	!
8	! module for photolysis online
9	!
10	! VERSION: Extracted from LMDZ.MARS work of Franck Lefevre (Yassin Jaziri)
11	! April 2019 - Yassin Jaziri add updates generic input (Yassin Jaziri)
12	!
13	!***********************************************************************
14
15	use types_asis
16	implicit none
17
18	! photolysis
19
20	integer, save :: nabs ! number of absorbing gases
21
22	! spectral grid
23
24	integer, save :: nw ! number of spectral intervals (low-res)
25	integer, save :: mopt ! high-res/low-res switch
26
27	real, allocatable, save :: wl(:) ! lower wavelength for each interval (nm)
28	real, allocatable, save :: wc(:) ! center wavelength for each interval (nm)
29	real, allocatable, save :: wu(:) ! upper wavelength for each interval (nm)
30	real, save :: photoheat_lmax ! maximum wavelength until photochemical heating is calculated
31
32	! solar flux
33
34	real, allocatable, save :: fstar1AU(:) ! stellar flux (photon.s-1.nm-1.cm-2) at 1 au
35
36	! cross-sections and quantum yields
37
38	real, allocatable, save :: qy(:,:) ! photodissociation yield
39	real, allocatable, save :: xs(:,:,:) ! absorption cross-section (cm2)
40	real, allocatable, save :: sj(:,:,:) ! general cross-section array by photodissociation yield (cm2)
41	real, allocatable, save :: xs_temp(:,:) ! absorption cross-section temperature (K)
42	integer, allocatable, save :: tdim(:) ! absorption cross-section temperature dimension
43	logical, allocatable, save :: jlabelbis(:) ! check in jlabel if the species is included in more than 1 photolysis
44
45	contains
46
47	subroutine init_photolysis(nlayer,nreact)
48
49	use types_asis, only: reactions
50	use tracer_h
51	use chimiedata_h, only: indexchim
52	use ioipsl_getin_p_mod, only: getin_p
53
54	integer, intent(in) :: nlayer ! number of atmospheric layers
55	integer, intent(in) :: nreact ! number of reactions in reactions files
56	integer :: iphot,tdim_max,i3prod,ij,iw,ilay,ireact
57	integer :: specheck(nesp)
58
59	! initialise on-line photolysis
60
61	! mopt = 1 high-resolution
62	! mopt = 2 low-resolution (recommended for gcm use)
63	! mopt = 3 input file reading grid
64
65	mopt = 3
66
67	! set wavelength grid
68
69	call gridw(nw,wl,wc,wu,mopt)
70
71	allocate(fstar1AU(nw))
72
73	! read and grid solar flux data
74
75	call rdsolarflux(nw,wl,wc,fstar1AU)
76
77	! read maximum wavelength until photochemical heating is calculated
78	write(,) "Maximum wavelength until photochemical heating is calculated"
79	photoheat_lmax=wc(nw-1) ! default value last point of wavelength grid
80	call getin_p("photoheat_lmax",photoheat_lmax)
81	write(,) "photoheat_lmax = ",photoheat_lmax
82
83	! calculate nabs number of absorbing gases
84	allocate(jlabelbis(nb_phot_hv_max))
85	nabs = 0
86	specheck(:) = 0
87	jlabelbis(:) = .true.
88	do iphot=1,nb_phot_hv_max
89	if (specheck(indexchim(trim(jlabel(iphot,2)))) .eq. 0) then
90	specheck(indexchim(trim(jlabel(iphot,2)))) = 1
91	nabs = nabs + 1
92	else
93	jlabelbis(iphot) = .false.
94	endif
95	end do
96
97	! Get temperature dimension and allocate
98	allocate(tdim(nb_hv_max))
99	tdim(:) = 1
100	tdim_max = 1
101	do iphot=1,nb_hv_max
102	read(jparamline(iphot),*) tdim(iphot)
103	if (tdim(iphot) > tdim_max) tdim_max = tdim(iphot)
104	end do
105
106	! allocation
107	allocate(qy(nw,nb_hv_max))
108	allocate(xs(tdim_max,nw,nb_hv_max))
109	allocate(sj(nlayer,nw,nb_phot_hv_max))
110	allocate(xs_temp(tdim_max,nb_hv_max))
111	xs(:,:,:) = 0.
112	xs_temp(:,:) = 0.
113
114	print*, 'WARNING: for all branching photolysis of one specie'
115	print*, ' the cross section files should be the same'
116	print*, ' they are used for optical depths calculation'
117	print*, ' only the quantum yield should be different'
118
119	i3prod = 0
120	ireact = 0
121
122	! read and grid cross-sections
123	do iphot=1,nb_hv_max
124	ireact = ireact + 1
125	call rdxsi(nw,wl,xs(:tdim(iphot),:,iphot),jparamline(iphot),xs_temp(:tdim(iphot),iphot),qy(:,iphot),tdim(iphot),jlabel(iphot+i3prod,2))
126	do while(reactions(ireact)%rtype/=0)
127	ireact = ireact + 1
128	end do
129	if (reactions(ireact)%three_prod) i3prod = i3prod + 1
130	end do
131
132	! init sj for no temperature dependent cross sections (tdim.eq.1)
133	iphot = 0
134	ij = 0
135	ireact = 0
136	do while(iphot<nb_phot_hv_max)
137	ij = ij + 1
138	iphot = iphot + 1
139	ireact = ireact + 1
140	if (tdim(ij).eq.1) then
141	do iw = 1,nw-1
142	do ilay = 1,nlayer
143	sj(ilay,iw,iphot) = xs(1,iw,ij)*qy(iw,ij)
144	end do
145	end do
146	endif
147	do while(reactions(ireact)%rtype/=0)
148	ireact = ireact + 1
149	end do
150	if (reactions(ireact)%three_prod) then
151	iphot = iphot + 1
152	if (tdim(ij).eq.1) then
153	do iw = 1,nw-1
154	do ilay = 1,nlayer
155	sj(ilay,iw,iphot) = sj(ilay,iw,iphot-1)
156	end do
157	end do
158	endif
159	end if
160	enddo
161
162	end subroutine init_photolysis
163
164	!==============================================================================
165
166	subroutine gridw(nw,wl,wc,wu,mopt)
167
168	use datafile_mod, only: datadir
169	use ioipsl_getin_p_mod, only: getin_p
170
171	!==========================================================================
172	! Create the wavelength grid for all interpolations and radiative transfer
173	! calculations. Grid may be irregularly spaced. Wavelengths are in nm.
174	! No gaps are allowed within the wavelength grid.
175	!==========================================================================
176
177	implicit none
178
179	! input
180
181	integer :: mopt ! high-res/low-res switch
182
183	! output
184
185	integer :: nw ! number of wavelength grid points
186	real, allocatable :: wl(:), wc(:), wu(:) ! lower, center, upper wavelength for each interval
187
188	! local
189
190	real :: wincr ! wavelength increment
191	integer :: iw, kw, ierr
192	character(len=200) :: fil
193	character(len=150) :: gridwfile
194
195	! mopt = 1 high-resolution mode (3789 intervals)
196	!
197	! 0-108 nm : 1.0 nm
198	! 108-124 nm : 0.1 nm
199	! 124-175 nm : 0.5 nm
200	! 175-205 nm : 0.01 nm
201	! 205-365 nm : 0.5 nm
202	! 365-850 nm : 5.0 nm
203	!
204	! mopt = 2 low-resolution mode
205	!
206	! 0-60 nm : 6.0 nm
207	! 60-80 nm : 2.0 nm
208	! 80-85 nm : 5.0 nm
209	! 85-117 nm : 2.0 nm
210	! 117-120 nm : 5.0 nm
211	! 120-123 nm : 0.2 nm
212	! 123-163 nm : 5.0 nm
213	! 163-175 nm : 2.0 nm
214	! 175-205 nm : 0.5 nm
215	! 205-245 nm : 5.0 nm
216	! 245-415 nm : 10.0 nm
217	! 415-815 nm : 50.0 nm
218	!
219	! mopt = 3 input file reading grid
220	!
221	! read "gridw.dat" in datadir/cross_sections/
222
223	if (mopt == 1) then ! high-res
224
225	nw = 3789
226	allocate(wl(nw))
227	allocate(wu(nw))
228	allocate(wc(nw))
229	! define wavelength intervals of width 1.0 nm from 0 to 108 nm:
230
231	kw = 0
232	wincr = 1.0
233	do iw = 0, 107
234	kw = kw + 1
235	wl(kw) = real(iw)
236	wu(kw) = wl(kw) + wincr
237	wc(kw) = (wl(kw) + wu(kw))/2.
238	end do
239
240	! define wavelength intervals of width 0.1 nm from 108 to 124 nm:
241
242	wincr = 0.1
243	do iw = 1080, 1239, 1
244	kw = kw + 1
245	wl(kw) = real(iw)/10.
246	wu(kw) = wl(kw) + wincr
247	wc(kw) = (wl(kw) + wu(kw))/2.
248	end do
249
250	! define wavelength intervals of width 0.5 nm from 124 to 175 nm:
251
252	wincr = 0.5
253	do iw = 1240, 1745, 5
254	kw = kw + 1
255	wl(kw) = real(iw)/10.
256	wu(kw) = wl(kw) + wincr
257	wc(kw) = (wl(kw) + wu(kw))/2.
258	end do
259
260	! define wavelength intervals of width 0.01 nm from 175 to 205 nm:
261
262	wincr = 0.01
263	do iw = 17500, 20499, 1
264	kw = kw + 1
265	wl(kw) = real(iw)/100.
266	wu(kw) = wl(kw) + wincr
267	wc(kw) = (wl(kw) + wu(kw))/2.
268	end do
269
270	! define wavelength intervals of width 0.5 nm from 205 to 365 nm:
271
272	wincr = 0.5
273	do iw = 2050, 3645, 5
274	kw = kw + 1
275	wl(kw) = real(iw)/10.
276	wu(kw) = wl(kw) + wincr
277	wc(kw) = (wl(kw) + wu(kw))/2.
278	end do
279
280	! define wavelength intervals of width 5.0 nm from 365 to 855 nm:
281
282	wincr = 5.0
283	do iw = 365, 850, 5
284	kw = kw + 1
285	wl(kw) = real(iw)
286	wu(kw) = wl(kw) + wincr
287	wc(kw) = (wl(kw) + wu(kw))/2.
288	end do
289	wl(kw+1) = wu(kw)
290
291	!============================================================
292
293	else if (mopt == 2) then ! low-res
294
295	nw = 162
296	allocate(wl(nw))
297	allocate(wu(nw))
298	allocate(wc(nw))
299
300	! define wavelength intervals of width 6.0 nm from 0 to 60 nm:
301
302	kw = 0
303	wincr = 6.0
304	DO iw = 0, 54, 6
305	kw = kw + 1
306	wl(kw) = real(iw)
307	wu(kw) = wl(kw) + wincr
308	wc(kw) = (wl(kw) + wu(kw))/2.
309	END DO
310
311	! define wavelength intervals of width 2.0 nm from 60 to 80 nm:
312
313	wincr = 2.0
314	DO iw = 60, 78, 2
315	kw = kw + 1
316	wl(kw) = real(iw)
317	wu(kw) = wl(kw) + wincr
318	wc(kw) = (wl(kw) + wu(kw))/2.
319	END DO
320
321	! define wavelength intervals of width 5.0 nm from 80 to 85 nm:
322
323	wincr = 5.0
324	DO iw = 80, 80
325	kw = kw + 1
326	wl(kw) = real(iw)
327	wu(kw) = wl(kw) + wincr
328	wc(kw) = (wl(kw) + wu(kw))/2.
329	END DO
330
331	! define wavelength intervals of width 2.0 nm from 85 to 117 nm:
332
333	wincr = 2.0
334	DO iw = 85, 115, 2
335	kw = kw + 1
336	wl(kw) = real(iw)
337	wu(kw) = wl(kw) + wincr
338	wc(kw) = (wl(kw) + wu(kw))/2.
339	END DO
340
341	! define wavelength intervals of width 3.0 nm from 117 to 120 nm:
342
343	wincr = 3.0
344	DO iw = 117, 117
345	kw = kw + 1
346	wl(kw) = real(iw)
347	wu(kw) = wl(kw) + wincr
348	wc(kw) = (wl(kw) + wu(kw))/2.
349	END DO
350
351	! define wavelength intervals of width 0.2 nm from 120 to 123 nm:
352
353	wincr = 0.2
354	DO iw = 1200, 1228, 2
355	kw = kw + 1
356	wl(kw) = real(iw)/10.
357	wu(kw) = wl(kw) + wincr
358	wc(kw) = (wl(kw) + wu(kw))/2.
359	ENDDO
360
361	! define wavelength intervals of width 5.0 nm from 123 to 163 nm:
362
363	wincr = 5.0
364	DO iw = 123, 158, 5
365	kw = kw + 1
366	wl(kw) = real(iw)
367	wu(kw) = wl(kw) + wincr
368	wc(kw) = (wl(kw) + wu(kw))/2.
369	ENDDO
370
371	! define wavelength intervals of width 2.0 nm from 163 to 175 nm:
372
373	wincr = 2.0
374	DO iw = 163, 173, 2
375	kw = kw + 1
376	wl(kw) = real(iw)
377	wu(kw) = wl(kw) + wincr
378	wc(kw) = (wl(kw) + wu(kw))/2.
379	ENDDO
380
381	! define wavelength intervals of width 0.5 nm from 175 to 205 nm:
382
383	wincr = 0.5
384	DO iw = 1750, 2045, 5
385	kw = kw + 1
386	wl(kw) = real(iw)/10.
387	wu(kw) = wl(kw) + wincr
388	wc(kw) = (wl(kw) + wu(kw))/2.
389	ENDDO
390
391	! define wavelength intervals of width 5.0 nm from 205 to 245 nm:
392
393	wincr = 5.
394	DO iw = 205, 240, 5
395	kw = kw + 1
396	wl(kw) = real(iw)
397	wu(kw) = wl(kw) + wincr
398	wc(kw) = (wl(kw) + wu(kw))/2.
399	ENDDO
400
401	! define wavelength intervals of width 10.0 nm from 245 to 415 nm:
402
403	wincr = 10.0
404	DO iw = 245, 405, 10
405	kw = kw + 1
406	wl(kw) = real(iw)
407	wu(kw) = wl(kw) + wincr
408	wc(kw) = (wl(kw) + wu(kw))/2.
409	ENDDO
410
411	! define wavelength intervals of width 50.0 nm from 415 to 865 nm:
412
413	wincr = 50.0
414	DO iw = 415, 815, 50
415	kw = kw + 1
416	wl(kw) = real(iw)
417	wu(kw) = wl(kw) + wincr
418	wc(kw) = (wl(kw) + wu(kw))/2.
419	ENDDO
420
421	wl(kw+1) = wu(kw)
422
423	!============================================================
424
425	else if (mopt == 3) then ! input file
426
427	! look for a " gridwfile= ..." option in def files
428	write(,) "Input wavelenght grid files for photolysis online is:"
429	gridwfile = "gridw.dat" ! default
430	call getin_p("gridwfile",gridwfile) ! default path
431	write(,) " gridwfile = ",trim(gridwfile)
432	write(,) 'Please use 1 and only 1 skipping line in ',trim(gridwfile)
433
434	! Opening file
435	fil = trim(datadir)//'/cross_sections/'//gridwfile
436	print*, 'wavelenght grid : ', fil
437	OPEN(UNIT=10,FILE=fil,STATUS='old',iostat=ierr)
438
439	if (ierr /= 0) THEN
440	write(,)'Error : cannot open wavelenght grid file ', trim(gridwfile)
441	write(,)'It should be in :',trim(datadir),'/cross_sections/'
442	write(,)'1) You can change the directory in callphys.def'
443	write(,)' with:'
444	write(,)' datadir=/path/to/the/directory'
445	write(,)'2) You can change the input wavelenght grid file name in'
446	write(,)' callphys.def with:'
447	write(,)' gridwfile=filename'
448	stop
449	end if
450
451	READ(10,*)
452
453	READ(10,*) nw
454	allocate(wl(nw))
455	allocate(wu(nw))
456	allocate(wc(nw))
457
458	kw = 0
459	READ(10,*) wl(1)
460	DO iw = 2, nw
461	READ(10,*) wl(iw)
462	wu(iw-1) = wl(iw)
463	wc(iw-1) = (wl(iw-1) + wu(iw-1))/2.
464	kw = kw + 1
465	ENDDO
466
467	close(10)
468	end if ! mopt
469
470	print*, 'number of spectral intervals : ', kw+1
471
472	end subroutine gridw
473
474	!==============================================================================
475
476	subroutine rdsolarflux(nw,wl,wc,fstar1AU)
477
478	! Read and re-grid solar flux data.
479
480	use datafile_mod, only: datadir
481	use ioipsl_getin_p_mod, only: getin_p
482
483	implicit none
484
485	! input
486
487	integer :: nw ! number of wavelength grid points
488	real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval
489
490	! output
491
492	real, dimension(nw) :: fstar1AU ! stellar flux (photon.s-1.nm-1.cm-2)
493
494	! local
495
496	integer :: iw, nhead, n, i, ierr, kin, naddpnt
497
498	real, parameter :: deltax = 1.e-4
499	real, allocatable :: x1(:), y1(:) ! input solar flux
500	real, dimension(nw) :: yg1 ! gridded solar flux
501
502	character(len=200) :: fil
503	character(len=150) :: stellarflux
504
505	kin = 10 ! input logical unit
506	naddpnt = 4 ! number of adding point
507	nhead = 1 ! number of skipping line at the beggining of the file
508
509	! look for a " stellarflux= ..." option in def files
510	write(,) "Input stellar spectra files for photolysis online is:"
511	stellarflux = "Sun_Claire2012_0.0Ga.txt" ! default
512	call getin_p("stellarflux",stellarflux) ! default path
513	write(,) " stellarflux = ",trim(stellarflux)
514	write(,) 'Please use ',nhead,' and only ',nhead,' skipping line in ',trim(stellarflux)
515
516	! Opening file
517	fil = trim(datadir)//'/stellar_spectra/photochem_stellar_spectra/'//stellarflux
518	print*, 'solar flux : ', fil
519	OPEN(UNIT=kin,FILE=fil,STATUS='old',iostat=ierr)
520
521	if (ierr /= 0) THEN
522	write(,)'Error : cannot open stellarflux file ', trim(stellarflux)
523	write(,)'It should be in :',trim(datadir),'/stellar_spectra/photochem_stellar_spectra/'
524	write(,)'1) You can change the data directory in callphys.def'
525	write(,)' with:'
526	write(,)' datadir=/path/to/the/directory'
527	write(,)'2) You can change the input stellarflux file name in'
528	write(,)' callphys.def with:'
529	write(,)' stellarflux=filename'
530	stop
531	end if
532
533	! Get number of line in the file
534	DO i = 1, nhead
535	READ(kin,*)
536	ENDDO
537	n = 0
538	do
539	read(kin,*,iostat=ierr)
540	if (ierr<0) exit
541	n = n + 1
542	end do
543	rewind(kin)
544	DO i = 1, nhead
545	READ(kin,*)
546	ENDDO
547	allocate(x1(n+naddpnt))
548	allocate(y1(n+naddpnt))
549
550	! Read stellar flux
551	DO i = 1, n
552	READ(kin,*) x1(i), y1(i) !-> x1 nm and y1 in photon.s-1.nm-1.cm-2
553	ENDDO
554	CLOSE (kin)
555
556	! Interpolation on the grid
557	CALL addpnt(x1,y1,n+naddpnt,n,x1(1)*(1.-deltax),0.)
558	CALL addpnt(x1,y1,n+naddpnt-1,n, 0.,0.)
559	CALL addpnt(x1,y1,n+naddpnt-2,n,x1(n)*(1.+deltax),0.)
560	CALL addpnt(x1,y1,n+naddpnt-3,n, 1.e+38,0.)
561	CALL inter2(nw,wl,yg1,n,x1,y1,ierr)
562	IF (ierr .NE. 0) THEN
563	WRITE(,) ierr, fil
564	STOP
565	ENDIF
566
567	! factor to convert to photon.s-1.nm-1.cm-2 :
568	! 5.039e11 = 1.e-41e-9/(hc = 6.62e-342.998e8)
569
570	! fstar1AU need to be in photon.s-1.nm-1.cm-2
571	DO iw = 1, nw-1
572	!f(iw) = yg1(iw)wc(iw)5.039e11 ! If yg1 in w.m-2.nm-1
573	fstar1AU(iw) = yg1(iw)
574	ENDDO
575
576	end subroutine rdsolarflux
577
578	!==============================================================================
579
580	subroutine addpnt ( x, y, ld, n, xnew, ynew )
581
582	!-----------------------------------------------------------------------------*
583	!= PURPOSE: =*
584	!= Add a point <xnew,ynew> to a set of data pairs <x,y>. x must be in =*
585	!= ascending order =*
586	!-----------------------------------------------------------------------------*
587	!= PARAMETERS: =*
588	!= X - REAL vector of length LD, x-coordinates (IO)=*
589	!= Y - REAL vector of length LD, y-values (IO)=*
590	!= LD - INTEGER, dimension of X, Y exactly as declared in the calling (I)=*
591	!= program =*
592	!= N - INTEGER, number of elements in X, Y. On entry, it must be: (IO)=*
593	!= N < LD. On exit, N is incremented by 1. =*
594	!= XNEW - REAL, x-coordinate at which point is to be added (I)=*
595	!= YNEW - REAL, y-value of point to be added (I)=*
596	!-----------------------------------------------------------------------------*
597
598	IMPLICIT NONE
599
600	! calling parameters
601
602	INTEGER ld, n
603	REAL x(ld), y(ld)
604	REAL xnew, ynew
605	INTEGER ierr
606
607	! local variables
608
609	INTEGER insert
610	INTEGER i
611
612	!-----------------------------------------------------------------------
613
614	! initialize error flag
615
616	ierr = 0
617
618	! check n<ld to make sure x will hold another point
619
620	IF (n .GE. ld) THEN
621	WRITE(0,*) '>>> ERROR (ADDPNT) <<< Cannot expand array '
622	WRITE(0,*) ' All elements used.'
623	STOP
624	ENDIF
625
626	insert = 1
627	i = 2
628
629	! check, whether x is already sorted.
630	! also, use this loop to find the point at which xnew needs to be inserted
631	! into vector x, if x is sorted.
632
633	10 CONTINUE
634	IF (i .LT. n) THEN
635	IF (x(i) .LT. x(i-1)) THEN
636	print*, x(i-1), x(i)
637	WRITE(0,*) '>>> ERROR (ADDPNT) <<< x-data must be in ascending order!'
638	STOP
639	ELSE
640	IF (xnew .GT. x(i)) insert = i + 1
641	ENDIF
642	i = i+1
643	GOTO 10
644	ENDIF
645
646	! if <xnew,ynew> needs to be appended at the end, just do so,
647	! otherwise, insert <xnew,ynew> at position INSERT
648
649	IF ( xnew .GT. x(n) ) THEN
650
651	x(n+1) = xnew
652	y(n+1) = ynew
653
654	ELSE
655
656	! shift all existing points one index up
657
658	DO i = n, insert, -1
659	x(i+1) = x(i)
660	y(i+1) = y(i)
661	ENDDO
662
663	! insert new point
664
665	x(insert) = xnew
666	y(insert) = ynew
667
668	ENDIF
669
670	! increase total number of elements in x, y
671
672	n = n+1
673
674	end subroutine addpnt
675
676	!==============================================================================
677
678	subroutine inter2(ng,xg,yg,n,x,y,ierr)
679
680	!-----------------------------------------------------------------------------*
681	!= PURPOSE: =*
682	!= Map input data given on single, discrete points onto a set of target =*
683	!= bins. =*
684	!= The original input data are given on single, discrete points of an =*
685	!= arbitrary grid and are being linearly interpolated onto a specified set =*
686	!= of target bins. In general, this is the case for most of the weighting =*
687	!= functions (action spectra, molecular cross section, and quantum yield =*
688	!= data), which have to be matched onto the specified wavelength intervals. =*
689	!= The average value in each target bin is found by averaging the trapezoi- =*
690	!= dal area underneath the input data curve (constructed by linearly connec-=*
691	!= ting the discrete input values). =*
692	!= Some caution should be used near the endpoints of the grids. If the =*
693	!= input data set does not span the range of the target grid, an error =*
694	!= message is printed and the execution is stopped, as extrapolation of the =*
695	!= data is not permitted. =*
696	!= If the input data does not encompass the target grid, use ADDPNT to =*
697	!= expand the input array. =*
698	!-----------------------------------------------------------------------------*
699	!= PARAMETERS: =*
700	!= NG - INTEGER, number of bins + 1 in the target grid (I)=*
701	!= XG - REAL, target grid (e.g., wavelength grid); bin i is defined (I)=*
702	!= as [XG(i),XG(i+1)] (i = 1..NG-1) =*
703	!= YG - REAL, y-data re-gridded onto XG, YG(i) specifies the value for (O)=*
704	!= bin i (i = 1..NG-1) =*
705	!= N - INTEGER, number of points in input grid (I)=*
706	!= X - REAL, grid on which input data are defined (I)=*
707	!= Y - REAL, input y-data (I)=*
708	!-----------------------------------------------------------------------------*
709
710	IMPLICIT NONE
711
712	! input:
713	INTEGER ng, n
714	REAL x(n), y(n), xg(ng)
715
716	! output:
717	REAL yg(ng)
718
719	! local:
720	REAL area, xgl, xgu
721	REAL darea, slope
722	REAL a1, a2, b1, b2
723	INTEGER ngintv
724	INTEGER i, k, jstart
725	INTEGER ierr
726	!_______________________________________________________________________
727
728	ierr = 0
729
730	! test for correct ordering of data, by increasing value of x
731
732	DO 10, i = 2, n
733	IF (x(i) .LE. x(i-1)) THEN
734	ierr = 1
735	WRITE(,)'data not sorted'
736	WRITE(,) x(i), x(i-1)
737	RETURN
738	ENDIF
739	10 CONTINUE
740
741	DO i = 2, ng
742	IF (xg(i) .LE. xg(i-1)) THEN
743	ierr = 2
744	WRITE(0,*) '>>> ERROR (inter2) <<< xg-grid not sorted!'
745	RETURN
746	ENDIF
747	ENDDO
748
749	! check for xg-values outside the x-range
750
751	IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN
752	WRITE(0,*) '>>> ERROR (inter2) <<< Data do not span grid. '
753	WRITE(0,*) ' Use ADDPNT to expand data and re-run.'
754	STOP
755	ENDIF
756
757	! find the integral of each grid interval and use this to
758	! calculate the average y value for the interval
759	! xgl and xgu are the lower and upper limits of the grid interval
760
761	jstart = 1
762	ngintv = ng - 1
763	DO 50, i = 1,ngintv
764
765	! initialize:
766
767	area = 0.0
768	xgl = xg(i)
769	xgu = xg(i+1)
770
771	! discard data before the first grid interval and after the
772	! last grid interval
773	! for internal grid intervals, start calculating area by interpolating
774	! between the last point which lies in the previous interval and the
775	! first point inside the current interval
776
777	k = jstart
778	IF (k .LE. n-1) THEN
779
780	! if both points are before the first grid, go to the next point
781	30 CONTINUE
782	IF (x(k+1) .LE. xgl) THEN
783	jstart = k - 1
784	k = k+1
785	IF (k .LE. n-1) GO TO 30
786	ENDIF
787
788
789	! if the last point is beyond the end of the grid, complete and go to the next
790	! grid
791	40 CONTINUE
792	IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN
793
794	jstart = k-1
795
796	! compute x-coordinates of increment
797
798	a1 = MAX(x(k),xgl)
799	a2 = MIN(x(k+1),xgu)
800
801	! if points coincide, contribution is zero
802
803	IF (x(k+1).EQ.x(k)) THEN
804	darea = 0.e0
805	ELSE
806	slope = (y(k+1) - y(k))/(x(k+1) - x(k))
807	b1 = y(k) + slope*(a1 - x(k))
808	b2 = y(k) + slope*(a2 - x(k))
809	darea = (a2 - a1)*(b2 + b1)/2.
810	ENDIF
811
812	! find the area under the trapezoid from a1 to a2
813
814	area = area + darea
815
816	! go to next point
817
818	k = k+1
819	GO TO 40
820
821	ENDIF
822	ENDIF
823
824	! calculate the average y after summing the areas in the interval
825
826	yg(i) = area/(xgu - xgl)
827
828	50 CONTINUE
829
830	end subroutine inter2
831
832	!==============================================================================
833
834	subroutine rdxsi(nw,wl,xsi,jparamlinei,xs_tempi,qyi,tdimi,species)
835
836	!-----------------------------------------------------------------------------*
837	!= PURPOSE: =*
838	!= Read molecular absorption cross section. Re-grid data to match =*
839	!= specified wavelength working grid. =*
840	!-----------------------------------------------------------------------------*
841	!= PARAMETERS: =*
842	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
843	!= wavelength grid =*
844	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
845	!= working wavelength grid =*
846	!= XS - REAL, molecular absoprtion cross section (cm^2) at each (O)=*
847	!= specified wavelength =*
848	!-----------------------------------------------------------------------------*
849
850	use datafile_mod, only: datadir
851
852	IMPLICIT NONE
853
854	! input
855
856	integer, intent(in) :: nw ! number of wavelength grid points
857	integer, intent(in) :: tdimi ! temperature dimension
858	real, intent(in) :: wl(nw) ! lower and central wavelength for each interval (nm)
859	character(len = 20) :: species ! species which is photolysed
860	character(len = 300), intent(in) :: jparamlinei ! line of jonline parameters
861
862	! output
863
864	real, intent(inout) :: qyi(nw) ! photodissociation yield
865	real, intent(inout) :: xsi(tdimi,nw) ! absorption cross-section (cm2)
866	real, intent(inout) :: xs_tempi(tdimi) ! absorption cross-section temperature (K)
867
868	! local
869
870	character(len = 50) :: sigfile(tdimi) ! cross section file name
871	character(len = 50) :: qyfile ! quantum yield file name
872	CHARACTER(len = 120) :: fil ! path files
873	integer :: tdimdummy, ifile, itemp, ierr, nheadxs, nheadqy, i, n, naddpnt
874	integer :: kin ! input logical unit
875	real, allocatable :: wlf(:), xsf(:) ! input cross section
876	real, allocatable :: wlqy(:), qyf(:) ! input photodissociation yield
877	real, parameter :: deltax = 1.e-4
878
879	kin = 10 ! input logical unit
880	naddpnt = 4 ! number of adding point (careful: same for xs and qy)
881	nheadxs = 1 ! number of skipping line at the beggining of the cross section files
882	nheadqy = 1 ! number of skipping line at the beggining of the quantum yield files
883	sigfile(:) = ''
884	qyfile = ''
885
886	read(jparamlinei,*) tdimdummy, (xs_tempi(itemp), itemp=1,tdimi), (sigfile(ifile), ifile=1,tdimi), qyfile
887
888	! Check if xs_tempi is sorted from low to high values
889	do i = 1,tdimi-1
890	if (xs_tempi(i)>xs_tempi(i+1)) then
891	print*, 'ERROR: temperature cross section file'
892	print*, ' has to be sorted from low to high values'
893	print*, ' Check reactfile file'
894	stop
895	end if
896	end do
897
898	! Cross section
899	do itemp=1,tdimi
900
901	fil = trim(datadir)//'/cross_sections/'//trim(sigfile(itemp))
902	print*, 'section efficace '//trim(species)//': ', fil
903
904	OPEN(UNIT=kin,FILE=fil,STATUS='old',iostat=ierr)
905
906	if (ierr /= 0) THEN
907	write(,)'Error : cannot open cross_sections file ', trim(sigfile(itemp))
908	write(,)'It should be in :',trim(datadir),'/cross_sections/'
909	write(,)'1) You can change the datadir directory in callphys.def'
910	write(,)' with:'
911	write(,)' datadir=/path/to/the/directory'
912	write(,)'2) You can check if the file is in datadir/cross_sections/'
913	write(,)'3) You can change the input cross_sections file name in'
914	write(,)' chemestry/reactfile file for the specie:',trim(species)
915	stop
916	end if
917
918	! Get number of line in the file
919	DO i = 1, nheadxs
920	READ(kin,*)
921	ENDDO
922	n = 0
923	do
924	read(kin,*,iostat=ierr)
925	if (ierr<0) exit
926	n = n + 1
927	end do
928	rewind(kin)
929	DO i = 1, nheadxs
930	READ(kin,*)
931	ENDDO
932
933	allocate(wlf(n+naddpnt))
934	allocate(xsf(n+naddpnt))
935
936	! Read cross section
937	DO i = 1, n
938	READ(kin,*) wlf(i), xsf(i) !-> xsf in cm2 and wlf in nm
939	ENDDO
940	CLOSE (kin)
941
942	CALL addpnt(wlf,xsf,n+naddpnt,n,wlf(1)*(1.-deltax),0.)
943	CALL addpnt(wlf,xsf,n+naddpnt-1,n, 0.,0.)
944	CALL addpnt(wlf,xsf,n+naddpnt-2,n,wlf(n)*(1.+deltax),0.)
945	CALL addpnt(wlf,xsf,n+naddpnt-3,n, 1.e+38,0.)
946	CALL inter2(nw,wl,xsi(itemp,:),n,wlf,xsf,ierr)
947	IF (ierr .NE. 0) THEN
948	WRITE(,) ierr, fil
949	STOP
950	ENDIF
951
952	deallocate(wlf)
953	deallocate(xsf)
954
955	end do
956
957	! Photodissociation yield
958	fil = trim(datadir)//'/cross_sections/'//trim(qyfile)
959	print*, 'photodissociation yield '//trim(species)//': ', fil
960
961	OPEN(UNIT=kin,FILE=fil,STATUS='old',iostat=ierr)
962
963	if (ierr /= 0) THEN
964	write(,)'Error : cannot open photodissociation yield file ', trim(qyfile)
965	write(,)'It should be in :',trim(datadir),'/cross_sections/'
966	write(,)'1) You can change the datadir directory in callphys.def'
967	write(,)' with:'
968	write(,)' datadir=/path/to/the/directory'
969	write(,)'2) You can check if the file is in datadir/cross_sections/'
970	write(,)'3) You can change the input photodissociation yield file name in'
971	write(,)' chemestry/reactfile file for the specie:',trim(species)
972	stop
973	end if
974
975	! Get number of line in the file
976	DO i = 1, nheadqy
977	READ(kin,*)
978	ENDDO
979	n = 0
980	do
981	read(kin,*,iostat=ierr)
982	if (ierr<0) exit
983	n = n + 1
984	end do
985	rewind(kin)
986	DO i = 1, nheadqy
987	READ(kin,*)
988	ENDDO
989	allocate(wlqy(n+naddpnt))
990	allocate(qyf(n+naddpnt))
991
992	! Read photodissociation yield
993	DO i = 1, n
994	READ(kin,*) wlqy(i), qyf(i) !-> wlqy in nm
995	ENDDO
996	CLOSE (kin)
997
998	CALL addpnt(wlqy,qyf,n+naddpnt,n,wlqy(1)*(1.-deltax),0.)
999	CALL addpnt(wlqy,qyf,n+naddpnt-1,n, 0.,0.)
1000	CALL addpnt(wlqy,qyf,n+naddpnt-2,n,wlqy(n)*(1.+deltax),0.)
1001	CALL addpnt(wlqy,qyf,n+naddpnt-3,n, 1.e+38,0.)
1002	CALL inter2(nw,wl,qyi,n,wlqy,qyf,ierr)
1003	IF (ierr .NE. 0) THEN
1004	WRITE(,) ierr, fil
1005	STOP
1006	ENDIF
1007
1008	end subroutine rdxsi
1009
1010	end module photolysis_mod
1011

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

Download in other formats:

Original Format