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

Rev	Line
[2542]	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
[3309]	34	real, allocatable, save :: fstar1AU(:) ! stellar flux (photon.s-1.nm-1.cm-2) at 1 au
[2542]	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
[2553]	47	subroutine init_photolysis(nlayer,nreact)
[2542]	48
[2553]	49	use types_asis, only: reactions
[2542]	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
[2553]	56	integer :: iphot,tdim_max,i3prod,ij,iw,ilay,ireact
[2542]	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
[3309]	71	allocate(fstar1AU(nw))
[2542]	72
	73	! read and grid solar flux data
	74
[3309]	75	call rdsolarflux(nw,wl,wc,fstar1AU)
[2542]	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
[2553]	98	allocate(tdim(nb_hv_max))
[2542]	99	tdim(:) = 1
	100	tdim_max = 1
[2553]	101	do iphot=1,nb_hv_max
[2542]	102	read(jparamline(iphot),*) tdim(iphot)
	103	if (tdim(iphot) > tdim_max) tdim_max = tdim(iphot)
	104	end do
	105
	106	! allocation
[2553]	107	allocate(qy(nw,nb_hv_max))
	108	allocate(xs(tdim_max,nw,nb_hv_max))
[2542]	109	allocate(sj(nlayer,nw,nb_phot_hv_max))
[2553]	110	allocate(xs_temp(tdim_max,nb_hv_max))
[2542]	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
[2553]	120	ireact = 0
[2542]	121
[2553]	122	! read and grid cross-sections
	123	do iphot=1,nb_hv_max
	124	ireact = ireact + 1
[2542]	125	call rdxsi(nw,wl,xs(:tdim(iphot),:,iphot),jparamline(iphot),xs_temp(:tdim(iphot),iphot),qy(:,iphot),tdim(iphot),jlabel(iphot+i3prod,2))
[2553]	126	do while(reactions(ireact)%rtype/=0)
	127	ireact = ireact + 1
	128	end do
	129	if (reactions(ireact)%three_prod) i3prod = i3prod + 1
[2542]	130	end do
	131
	132	! init sj for no temperature dependent cross sections (tdim.eq.1)
[2553]	133	iphot = 0
	134	ij = 0
	135	ireact = 0
[2542]	136	do while(iphot<nb_phot_hv_max)
[2553]	137	ij = ij + 1
	138	iphot = iphot + 1
	139	ireact = ireact + 1
[2542]	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
[2553]	147	do while(reactions(ireact)%rtype/=0)
	148	ireact = ireact + 1
	149	end do
	150	if (reactions(ireact)%three_prod) then
[2542]	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
[3309]	476	subroutine rdsolarflux(nw,wl,wc,fstar1AU)
[2542]	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
[3309]	492	real, dimension(nw) :: fstar1AU ! stellar flux (photon.s-1.nm-1.cm-2)
[2542]	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:"
[3370]	511	stellarflux = "Sun_Claire2012_0.0Ga.txt" ! default
[2542]	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
[3370]	517	fil = trim(datadir)//'/stellar_spectra/photochem_stellar_spectra/'//stellarflux
[2542]	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)
[3370]	523	write(,)'It should be in :',trim(datadir),'/stellar_spectra/photochem_stellar_spectra/'
[2542]	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
[3309]	570	! fstar1AU need to be in photon.s-1.nm-1.cm-2
[2542]	571	DO iw = 1, nw-1
	572	!f(iw) = yg1(iw)wc(iw)5.039e11 ! If yg1 in w.m-2.nm-1
[3309]	573	fstar1AU(iw) = yg1(iw)
[2542]	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'
[2553]	893	print*, ' Check reactfile file'
[2542]	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'
[2553]	914	write(,)' chemestry/reactfile file for the specie:',trim(species)
[2542]	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'
[2553]	971	write(,)' chemestry/reactfile file for the specie:',trim(species)
[2542]	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