Context Navigation

photolysis_mod.F90 @ 2082

Last change on this file since 2082 was 2044, checked in by flefevre, 7 years ago
Photolyse on-line : desormais par defaut
File size: 57.2 KB

Rev	Line
[2029]	1	module photolysis_mod
	2
	3	implicit none
	4
[2044]	5	! photolysis
	6
	7	logical, parameter :: jonline = .true. ! true: on-line ! false: lookup table
[2030]	8	integer, parameter :: nphot = 13 ! number of photolysis
	9	integer, parameter :: nabs = 10 ! number of absorbing gases
	10
	11	! number of reactions in chemical solver
	12
	13	integer, parameter :: nb_phot_max = nphot + 9 ! photolysis + quenching/heterogeneous
	14	integer, parameter :: nb_reaction_3_max = 6 ! quadratic
	15	integer, parameter :: nb_reaction_4_max = 31 ! bimolecular
	16
[2029]	17	! spectral grid
	18
[2030]	19	integer, parameter :: nw = 162 ! number of spectral intervals (low-res)
[2041]	20	integer, save :: mopt ! high-res/low-res switch
[2029]	21
	22	real, dimension(nw), save :: wl, wc, wu ! lower, center, upper wavelength for each interval
	23
	24	! solar flux
	25
	26	real, dimension(nw), save :: f ! solar flux (w.m-2.nm-1) at 1 au
	27
	28	! cross-sections and quantum yields
	29
	30	real, dimension(nw), save :: xsco2_195, xsco2_295, xsco2_370 ! co2 absorption cross-section at 195-295-370 k (cm2)
	31	real, dimension(nw), save :: yieldco2 ! co2 photodissociation yield
	32	real, dimension(nw), save :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 absorption cross-section at 150-200-250-300 k (cm2)
	33	real, dimension(nw), save :: yieldo2 ! o2 photodissociation yield
	34	real, dimension(nw), save :: xso3_218, xso3_298 ! o3 absorption cross-section at 218-298 k (cm2)
	35	real, dimension(nw), save :: xsh2o ! h2o absorption cross-section (cm2)
	36	real, dimension(nw), save :: xsh2o2 ! h2o2 absorption cross-section (cm2)
	37	real, dimension(nw), save :: xsho2 ! ho2 absorption cross-section (cm2)
	38	real, dimension(nw), save :: xsh2 ! h2 absorption cross-section (cm2)
	39	real, dimension(nw), save :: yieldh2 ! h2 photodissociation yield
	40	real, dimension(nw), save :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2)
	41	real, dimension(nw), save :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k
	42	real, dimension(nw), save :: xsno ! no absorption cross-section (cm2)
	43	real, dimension(nw), save :: yieldno ! no photodissociation yield
	44	real, dimension(nw), save :: xsn2 ! n2 absorption cross-section (cm2)
	45	real, dimension(nw), save :: yieldn2 ! n2 photodissociation yield
	46	real, dimension(nw), save :: albedo ! surface albedo
	47
	48	contains
	49
	50	subroutine init_photolysis
	51
	52	! initialise on-line photolysis
	53
	54	! mopt = 1 high-resolution
	55	! mopt = 2 low-resolution (recommended for gcm use)
	56
	57	mopt = 2
	58
	59	! set wavelength grid
	60
	61	call gridw(nw,wl,wc,wu,mopt)
	62
	63	! read and grid solar flux data
	64
	65	call rdsolarflux(nw,wl,wc,f)
	66
	67	! read and grid o2 cross-sections
	68
	69	call rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2)
	70
	71	! read and grid co2 cross-sections
	72
	73	call rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)
	74
	75	! read and grid o3 cross-sections
	76
	77	call rdxso3(nw,wl,xso3_218,xso3_298)
	78
	79	! read and grid h2o cross-sections
	80
	81	call rdxsh2o(nw,wl,xsh2o)
	82
	83	! read and grid h2o2 cross-sections
	84
	85	call rdxsh2o2(nw,wl,xsh2o2)
	86
	87	! read and grid ho2 cross-sections
	88
	89	call rdxsho2(nw,wl,xsho2)
	90
	91	! read and grid h2 cross-sections
	92
	93	call rdxsh2(nw,wl,wc,xsh2,yieldh2)
	94
	95	! read and grid no cross-sections
	96
	97	call rdxsno(nw,wl,xsno,yieldno)
	98
	99	! read and grid no2 cross-sections
	100
	101	call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248,yldno2_298)
	102
	103	! read and grid n2 cross-sections
	104
	105	call rdxsn2(nw,wl,xsn2,yieldn2)
	106
	107	! set surface albedo
	108
	109	call setalb(nw,wl,albedo)
	110
	111	end subroutine init_photolysis
	112
	113	!==============================================================================
	114
	115	subroutine gridw(nw,wl,wc,wu,mopt)
	116
	117	! Create the wavelength grid for all interpolations and radiative transfer
	118	! calculations. Grid may be irregularly spaced. Wavelengths are in nm.
	119	! No gaps are allowed within the wavelength grid.
	120
	121	implicit none
	122
	123	! input
	124
	125	integer :: nw ! number of wavelength grid points
	126	integer :: mopt ! high-res/low-res switch
	127
	128	! output
	129
	130	real, dimension(nw) :: wl, wc, wu ! lower, center, upper wavelength for each interval
	131
	132	! local
	133
	134	real :: wincr ! wavelength increment
	135	integer :: iw, kw
	136
	137	! mopt = 1 high-resolution mode (3789 intervals)
	138	!
	139	! 0-108 nm : 1.0 nm
	140	! 108-124 nm : 0.1 nm
	141	! 124-175 nm : 0.5 nm
	142	! 175-205 nm : 0.01 nm
	143	! 205-365 nm : 0.5 nm
	144	! 365-850 nm : 5.0 nm
	145	!
	146	! mopt = 2 low-resolution mode
	147	!
	148	! 0-60 nm : 6.0 nm
	149	! 60-80 nm : 2.0 nm
[2030]	150	! 80-85 nm : 5.0 nm
	151	! 85-117 nm : 2.0 nm
	152	! 117-120 nm : 5.0 nm
[2029]	153	! 120-123 nm : 0.2 nm
	154	! 123-163 nm : 5.0 nm
	155	! 163-175 nm : 2.0 nm
	156	! 175-205 nm : 0.5 nm
	157	! 205-245 nm : 5.0 nm
	158	! 245-415 nm : 10.0 nm
	159	! 415-815 nm : 50.0 nm
	160
	161	if (mopt == 1) then ! high-res
	162
	163	! define wavelength intervals of width 1.0 nm from 0 to 108 nm:
	164
	165	kw = 0
	166	wincr = 1.0
	167	do iw = 0, 107
	168	kw = kw + 1
	169	wl(kw) = real(iw)
	170	wu(kw) = wl(kw) + wincr
	171	wc(kw) = (wl(kw) + wu(kw))/2.
	172	end do
	173
	174	! define wavelength intervals of width 0.1 nm from 108 to 124 nm:
	175
	176	wincr = 0.1
	177	do iw = 1080, 1239, 1
	178	kw = kw + 1
	179	wl(kw) = real(iw)/10.
	180	wu(kw) = wl(kw) + wincr
	181	wc(kw) = (wl(kw) + wu(kw))/2.
	182	end do
	183
	184	! define wavelength intervals of width 0.5 nm from 124 to 175 nm:
	185
	186	wincr = 0.5
	187	do iw = 1240, 1745, 5
	188	kw = kw + 1
	189	wl(kw) = real(iw)/10.
	190	wu(kw) = wl(kw) + wincr
	191	wc(kw) = (wl(kw) + wu(kw))/2.
	192	end do
	193
	194	! define wavelength intervals of width 0.01 nm from 175 to 205 nm:
	195
	196	wincr = 0.01
	197	do iw = 17500, 20499, 1
	198	kw = kw + 1
	199	wl(kw) = real(iw)/100.
	200	wu(kw) = wl(kw) + wincr
	201	wc(kw) = (wl(kw) + wu(kw))/2.
	202	end do
	203
	204	! define wavelength intervals of width 0.5 nm from 205 to 365 nm:
	205
	206	wincr = 0.5
	207	do iw = 2050, 3645, 5
	208	kw = kw + 1
	209	wl(kw) = real(iw)/10.
	210	wu(kw) = wl(kw) + wincr
	211	wc(kw) = (wl(kw) + wu(kw))/2.
	212	end do
	213
	214	! define wavelength intervals of width 5.0 nm from 365 to 855 nm:
	215
	216	wincr = 5.0
	217	do iw = 365, 850, 5
	218	kw = kw + 1
	219	wl(kw) = real(iw)
	220	wu(kw) = wl(kw) + wincr
	221	wc(kw) = (wl(kw) + wu(kw))/2.
	222	end do
	223	wl(kw+1) = wu(kw)
	224
[2030]	225	!============================================================
	226
[2029]	227	else if (mopt == 2) then ! low-res
	228
	229	! define wavelength intervals of width 6.0 nm from 0 to 60 nm:
	230
	231	kw = 0
	232	wincr = 6.0
	233	DO iw = 0, 54, 6
	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 2.0 nm from 60 to 80 nm:
	241
	242	wincr = 2.0
	243	DO iw = 60, 78, 2
	244	kw = kw + 1
	245	wl(kw) = real(iw)
	246	wu(kw) = wl(kw) + wincr
	247	wc(kw) = (wl(kw) + wu(kw))/2.
	248	END DO
	249
[2030]	250	! define wavelength intervals of width 5.0 nm from 80 to 85 nm:
[2029]	251
	252	wincr = 5.0
[2041]	253	DO iw = 80, 80
[2029]	254	kw = kw + 1
	255	wl(kw) = real(iw)
	256	wu(kw) = wl(kw) + wincr
	257	wc(kw) = (wl(kw) + wu(kw))/2.
	258	END DO
	259
[2030]	260	! define wavelength intervals of width 2.0 nm from 85 to 117 nm:
	261
	262	wincr = 2.0
	263	DO iw = 85, 115, 2
	264	kw = kw + 1
	265	wl(kw) = real(iw)
	266	wu(kw) = wl(kw) + wincr
	267	wc(kw) = (wl(kw) + wu(kw))/2.
	268	END DO
	269
[2041]	270	! define wavelength intervals of width 3.0 nm from 117 to 120 nm:
[2030]	271
	272	wincr = 3.0
[2041]	273	DO iw = 117, 117
[2030]	274	kw = kw + 1
	275	wl(kw) = real(iw)
	276	wu(kw) = wl(kw) + wincr
	277	wc(kw) = (wl(kw) + wu(kw))/2.
	278	END DO
	279
[2029]	280	! define wavelength intervals of width 0.2 nm from 120 to 123 nm:
	281
	282	wincr = 0.2
	283	DO iw = 1200, 1228, 2
	284	kw = kw + 1
	285	wl(kw) = real(iw)/10.
	286	wu(kw) = wl(kw) + wincr
	287	wc(kw) = (wl(kw) + wu(kw))/2.
[2030]	288	ENDDO
[2029]	289
	290	! define wavelength intervals of width 5.0 nm from 123 to 163 nm:
	291
	292	wincr = 5.0
	293	DO iw = 123, 158, 5
	294	kw = kw + 1
	295	wl(kw) = real(iw)
	296	wu(kw) = wl(kw) + wincr
	297	wc(kw) = (wl(kw) + wu(kw))/2.
[2030]	298	ENDDO
[2029]	299
	300	! define wavelength intervals of width 2.0 nm from 163 to 175 nm:
	301
	302	wincr = 2.0
	303	DO iw = 163, 173, 2
	304	kw = kw + 1
	305	wl(kw) = real(iw)
	306	wu(kw) = wl(kw) + wincr
	307	wc(kw) = (wl(kw) + wu(kw))/2.
[2030]	308	ENDDO
[2029]	309
	310	! define wavelength intervals of width 0.5 nm from 175 to 205 nm:
	311
	312	wincr = 0.5
	313	DO iw = 1750, 2045, 5
	314	kw = kw + 1
	315	wl(kw) = real(iw)/10.
	316	wu(kw) = wl(kw) + wincr
	317	wc(kw) = (wl(kw) + wu(kw))/2.
[2030]	318	ENDDO
[2029]	319
	320	! define wavelength intervals of width 5.0 nm from 205 to 245 nm:
	321
	322	wincr = 5.
	323	DO iw = 205, 240, 5
	324	kw = kw + 1
	325	wl(kw) = real(iw)
	326	wu(kw) = wl(kw) + wincr
	327	wc(kw) = (wl(kw) + wu(kw))/2.
[2030]	328	ENDDO
[2029]	329
	330	! define wavelength intervals of width 10.0 nm from 245 to 415 nm:
	331
	332	wincr = 10.0
	333	DO iw = 245, 405, 10
	334	kw = kw + 1
	335	wl(kw) = real(iw)
	336	wu(kw) = wl(kw) + wincr
	337	wc(kw) = (wl(kw) + wu(kw))/2.
[2030]	338	ENDDO
[2029]	339
	340	! define wavelength intervals of width 50.0 nm from 415 to 815 nm:
	341
	342	wincr = 50.0
	343	DO iw = 415, 815, 50
	344	kw = kw + 1
	345	wl(kw) = real(iw)
	346	wu(kw) = wl(kw) + wincr
	347	wc(kw) = (wl(kw) + wu(kw))/2.
[2030]	348	ENDDO
[2029]	349
	350	wl(kw+1) = wu(kw)
	351
	352	end if ! mopt
[2030]	353
	354	print*, 'number of spectral intervals : ', kw+1
[2029]	355
	356	end subroutine gridw
	357
	358	!==============================================================================
	359
	360	subroutine rdsolarflux(nw,wl,wc,f)
	361
	362	! Read and re-grid solar flux data.
	363
	364	use datafile_mod, only: datadir
	365
	366	implicit none
	367
	368	! input
	369
	370	integer :: nw ! number of wavelength grid points
	371	real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval
	372
	373	! output
	374
	375	real, dimension(nw) :: f ! solar flux (w.m-2.nm-1)
	376
	377	! local
	378
	379	integer, parameter :: kdata = 20000 ! max dimension of input solar flux
	380	integer :: msun ! choice of solar flux
	381	integer :: iw, nhead, ihead, n, i, ierr, kin
	382
	383	real, parameter :: deltax = 1.e-4
	384	real, dimension(kdata) :: x1, y1 ! input solar flux
	385	real, dimension(nw) :: yg1 ! gridded solar flux
	386
	387	character(len=100) :: fil
	388
	389	kin = 10 ! input logical unit
	390
	391	! select desired extra-terrestrial solar irradiance, using msun:
	392
	393	! 18 = atlas3_thuillier_tuv.txt 0-900 nm November 1994
	394	! Thuillier et al., Adv. Space. Res., 34, 256-261, 2004
	395
	396	msun = 18
	397
	398	if (msun == 18) THEN
	399
	400	fil = trim(datadir)//'/solar_fluxes/atlas3_thuillier_tuv.txt'
	401	print*, 'solar flux : ', fil
[2044]	402	open(kin, file=fil, status='old', iostat=ierr)
	403
	404	if (ierr /= 0) THEN
	405	write(,)'cant find solar flux : ', fil
	406	write(,)'It should be in :', trim(datadir),'/solar_fluxes'
	407	write(,)'1) You can change this directory address in '
	408	write(,)' callphys.def with datadir=/path/to/dir'
	409	write(,)'2) If necessary, /solar fluxes (and other datafiles)'
	410	write(,)' can be obtained online on:'
	411	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
	412	stop
	413	end if
	414
[2029]	415	nhead = 9
	416	n = 19193
	417	DO ihead = 1, nhead
	418	READ(kin,*)
	419	ENDDO
	420	DO i = 1, n
	421	READ(kin,*) x1(i), y1(i)
	422	y1(i) = y1(i)*1.e-3 ! mw -> w
	423	ENDDO
	424	CLOSE (kin)
	425	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
	426	CALL addpnt(x1,y1,kdata,n, 0.,0.)
	427	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
	428	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
	429	CALL inter2(nw,wl,yg1,n,x1,y1,ierr)
	430
	431	IF (ierr .NE. 0) THEN
	432	WRITE(,) ierr, fil
	433	STOP
	434	ENDIF
	435
	436	! convert to photon.s-1.nm-1.cm-2
	437	! 5.039e11 = 1.e-41e-9/(hc = 6.62e-342.998e8)
	438
	439	DO iw = 1, nw-1
	440	f(iw) = yg1(iw)wc(iw)5.039e11
	441	! write(25,*) iw, wc(iw), f(iw)
	442	ENDDO
	443
	444	end if
	445
	446	end subroutine rdsolarflux
	447
	448	!==============================================================================
	449
	450	subroutine addpnt ( x, y, ld, n, xnew, ynew )
	451
	452	!-----------------------------------------------------------------------------*
	453	!= PURPOSE: =*
	454	!= Add a point <xnew,ynew> to a set of data pairs <x,y>. x must be in =*
	455	!= ascending order =*
	456	!-----------------------------------------------------------------------------*
	457	!= PARAMETERS: =*
	458	!= X - REAL vector of length LD, x-coordinates (IO)=*
	459	!= Y - REAL vector of length LD, y-values (IO)=*
	460	!= LD - INTEGER, dimension of X, Y exactly as declared in the calling (I)=*
	461	!= program =*
	462	!= N - INTEGER, number of elements in X, Y. On entry, it must be: (IO)=*
	463	!= N < LD. On exit, N is incremented by 1. =*
	464	!= XNEW - REAL, x-coordinate at which point is to be added (I)=*
	465	!= YNEW - REAL, y-value of point to be added (I)=*
	466	!-----------------------------------------------------------------------------*
	467
	468	IMPLICIT NONE
	469
	470	! calling parameters
	471
	472	INTEGER ld, n
	473	REAL x(ld), y(ld)
	474	REAL xnew, ynew
	475	INTEGER ierr
	476
	477	! local variables
	478
	479	INTEGER insert
	480	INTEGER i
	481
	482	!-----------------------------------------------------------------------
	483
	484	! initialize error flag
	485
	486	ierr = 0
	487
	488	! check n<ld to make sure x will hold another point
	489
	490	IF (n .GE. ld) THEN
	491	WRITE(0,*) '>>> ERROR (ADDPNT) <<< Cannot expand array '
	492	WRITE(0,*) ' All elements used.'
	493	STOP
	494	ENDIF
	495
	496	insert = 1
	497	i = 2
	498
	499	! check, whether x is already sorted.
	500	! also, use this loop to find the point at which xnew needs to be inserted
	501	! into vector x, if x is sorted.
	502
	503	10 CONTINUE
	504	IF (i .LT. n) THEN
	505	IF (x(i) .LT. x(i-1)) THEN
	506	print*, x(i-1), x(i)
	507	WRITE(0,*) '>>> ERROR (ADDPNT) <<< x-data must be in ascending order!'
	508	STOP
	509	ELSE
	510	IF (xnew .GT. x(i)) insert = i + 1
	511	ENDIF
	512	i = i+1
	513	GOTO 10
	514	ENDIF
	515
	516	! if <xnew,ynew> needs to be appended at the end, just do so,
	517	! otherwise, insert <xnew,ynew> at position INSERT
	518
	519	IF ( xnew .GT. x(n) ) THEN
	520
	521	x(n+1) = xnew
	522	y(n+1) = ynew
	523
	524	ELSE
	525
	526	! shift all existing points one index up
	527
	528	DO i = n, insert, -1
	529	x(i+1) = x(i)
	530	y(i+1) = y(i)
	531	ENDDO
	532
	533	! insert new point
	534
	535	x(insert) = xnew
	536	y(insert) = ynew
	537
	538	ENDIF
	539
	540	! increase total number of elements in x, y
	541
	542	n = n+1
	543
	544	end subroutine addpnt
	545
	546	!==============================================================================
	547
	548	subroutine inter2(ng,xg,yg,n,x,y,ierr)
	549
	550	!-----------------------------------------------------------------------------*
	551	!= PURPOSE: =*
	552	!= Map input data given on single, discrete points onto a set of target =*
	553	!= bins. =*
	554	!= The original input data are given on single, discrete points of an =*
	555	!= arbitrary grid and are being linearly interpolated onto a specified set =*
	556	!= of target bins. In general, this is the case for most of the weighting =*
	557	!= functions (action spectra, molecular cross section, and quantum yield =*
	558	!= data), which have to be matched onto the specified wavelength intervals. =*
	559	!= The average value in each target bin is found by averaging the trapezoi- =*
	560	!= dal area underneath the input data curve (constructed by linearly connec-=*
	561	!= ting the discrete input values). =*
	562	!= Some caution should be used near the endpoints of the grids. If the =*
	563	!= input data set does not span the range of the target grid, an error =*
	564	!= message is printed and the execution is stopped, as extrapolation of the =*
	565	!= data is not permitted. =*
	566	!= If the input data does not encompass the target grid, use ADDPNT to =*
	567	!= expand the input array. =*
	568	!-----------------------------------------------------------------------------*
	569	!= PARAMETERS: =*
	570	!= NG - INTEGER, number of bins + 1 in the target grid (I)=*
	571	!= XG - REAL, target grid (e.g., wavelength grid); bin i is defined (I)=*
	572	!= as [XG(i),XG(i+1)] (i = 1..NG-1) =*
	573	!= YG - REAL, y-data re-gridded onto XG, YG(i) specifies the value for (O)=*
	574	!= bin i (i = 1..NG-1) =*
	575	!= N - INTEGER, number of points in input grid (I)=*
	576	!= X - REAL, grid on which input data are defined (I)=*
	577	!= Y - REAL, input y-data (I)=*
	578	!-----------------------------------------------------------------------------*
	579
	580	IMPLICIT NONE
	581
	582	! input:
	583	INTEGER ng, n
	584	REAL x(n), y(n), xg(ng)
	585
	586	! output:
	587	REAL yg(ng)
	588
	589	! local:
	590	REAL area, xgl, xgu
	591	REAL darea, slope
	592	REAL a1, a2, b1, b2
	593	INTEGER ngintv
	594	INTEGER i, k, jstart
	595	INTEGER ierr
	596	!_______________________________________________________________________
	597
	598	ierr = 0
	599
	600	! test for correct ordering of data, by increasing value of x
	601
	602	DO 10, i = 2, n
	603	IF (x(i) .LE. x(i-1)) THEN
	604	ierr = 1
	605	WRITE(,)'data not sorted'
	606	WRITE(,) x(i), x(i-1)
	607	RETURN
	608	ENDIF
	609	10 CONTINUE
	610
	611	DO i = 2, ng
	612	IF (xg(i) .LE. xg(i-1)) THEN
	613	ierr = 2
	614	WRITE(0,*) '>>> ERROR (inter2) <<< xg-grid not sorted!'
	615	RETURN
	616	ENDIF
	617	ENDDO
	618
	619	! check for xg-values outside the x-range
	620
	621	IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN
	622	WRITE(0,*) '>>> ERROR (inter2) <<< Data do not span grid. '
	623	WRITE(0,*) ' Use ADDPNT to expand data and re-run.'
	624	STOP
	625	ENDIF
	626
	627	! find the integral of each grid interval and use this to
	628	! calculate the average y value for the interval
	629	! xgl and xgu are the lower and upper limits of the grid interval
	630
	631	jstart = 1
	632	ngintv = ng - 1
	633	DO 50, i = 1,ngintv
	634
	635	! initialize:
	636
	637	area = 0.0
	638	xgl = xg(i)
	639	xgu = xg(i+1)
	640
	641	! discard data before the first grid interval and after the
	642	! last grid interval
	643	! for internal grid intervals, start calculating area by interpolating
	644	! between the last point which lies in the previous interval and the
	645	! first point inside the current interval
	646
	647	k = jstart
	648	IF (k .LE. n-1) THEN
	649
	650	! if both points are before the first grid, go to the next point
	651	30 CONTINUE
	652	IF (x(k+1) .LE. xgl) THEN
	653	jstart = k - 1
	654	k = k+1
	655	IF (k .LE. n-1) GO TO 30
	656	ENDIF
	657
	658
	659	! if the last point is beyond the end of the grid, complete and go to the next
	660	! grid
	661	40 CONTINUE
	662	IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN
	663
	664	jstart = k-1
	665
	666	! compute x-coordinates of increment
	667
	668	a1 = MAX(x(k),xgl)
	669	a2 = MIN(x(k+1),xgu)
	670
	671	! if points coincide, contribution is zero
	672
	673	IF (x(k+1).EQ.x(k)) THEN
	674	darea = 0.e0
	675	ELSE
	676	slope = (y(k+1) - y(k))/(x(k+1) - x(k))
	677	b1 = y(k) + slope*(a1 - x(k))
	678	b2 = y(k) + slope*(a2 - x(k))
	679	darea = (a2 - a1)*(b2 + b1)/2.
	680	ENDIF
	681
	682	! find the area under the trapezoid from a1 to a2
	683
	684	area = area + darea
	685
	686	! go to next point
	687
	688	k = k+1
	689	GO TO 40
	690
	691	ENDIF
	692	ENDIF
	693
	694	! calculate the average y after summing the areas in the interval
	695
	696	yg(i) = area/(xgu - xgl)
	697
	698	50 CONTINUE
	699
	700	end subroutine inter2
	701
	702	!==============================================================================
	703
	704	subroutine rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)
	705
	706	!-----------------------------------------------------------------------------*
	707	!= PURPOSE: =*
	708	!= Read and grid CO2 absorption cross-sections and photodissociation yield =*
	709	!-----------------------------------------------------------------------------*
	710	!= PARAMETERS: =*
	711	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	712	!= wavelength grid =*
	713	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
	714	!= working wavelength grid =*
	715	!= XSCO2 - REAL, molecular absoprtion cross section (cm^2) of CO2 at (O)=*
	716	!= each specified wavelength =*
	717	!-----------------------------------------------------------------------------*
	718
	719	use datafile_mod, only: datadir
	720
	721	implicit none
	722
	723	! input
	724
	725	integer :: nw ! number of wavelength grid points
	726	real, dimension(nw) :: wl ! lower and central wavelength for each interval
	727
	728	! output
	729
	730	real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
	731	real, dimension(nw) :: yieldco2 ! co2 photodissociation yield
	732
	733	! local
	734
	735	integer, parameter :: kdata = 42000
	736	real, parameter :: deltax = 1.e-4
	737	real, dimension(kdata) :: x1, y1, y2, y3, xion, ion
	738	real, dimension(nw) :: yg
	739	real :: xl, xu
	740	integer :: ierr, i, l, n, n1, n2, n3, n4
	741	CHARACTER*100 fil
	742
	743	integer :: kin, kout ! input/ouput logical units
	744
	745	kin = 10
	746	kout = 30
	747
	748	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	749	!
	750	! CO2 absorption cross-sections
	751	!
	752	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	753	!
	754	! 195K: huestis and berkowitz (2010) + Starck et al. (2006)
	755	! + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
	756	!
	757	! 295K: huestis and berkowitz (2010) + Starck et al. (2006)
	758	! + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
	759	!
	760	! 370K: huestis and berkowitz (2010) + Starck et al. (2006)
	761	! + Lewis and Carver (1983) + extrapolation
	762	!
	763	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	764
	765	n1 = 40769
	766	n2 = 41586
	767	n3 = 10110
	768
	769	! 195K:
	770
	771	fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_195k.txt'
	772	print*, 'section efficace CO2 195K: ', fil
	773
	774	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	775	DO i = 1,11
	776	read(kin,*)
	777	END DO
	778
	779	DO i = 1, n1
	780	READ(kin,*) x1(i), y1(i)
	781	END DO
	782	CLOSE (kin)
	783
	784	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
	785	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
	786	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
	787	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
	788	CALL inter2(nw,wl,yg,n1,x1,y1,ierr)
	789	IF (ierr .NE. 0) THEN
	790	WRITE(,) ierr, fil
	791	STOP
	792	ENDIF
	793
	794	DO l = 1, nw-1
	795	xsco2_195(l) = yg(l)
	796	END DO
	797
	798	! 295K:
	799
	800	fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_295k.txt'
	801	print*, 'section efficace CO2 295K: ', fil
	802
	803	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	804	DO i = 1,11
	805	read(kin,*)
	806	END DO
	807
	808	DO i = 1, n2
	809	READ(kin,*) x1(i), y1(i)
	810	END DO
	811	CLOSE (kin)
	812
	813	CALL addpnt(x1,y1,kdata,n2,x1(1)*(1.-deltax),0.)
	814	CALL addpnt(x1,y1,kdata,n2, 0.,0.)
	815	CALL addpnt(x1,y1,kdata,n2,x1(n2)*(1.+deltax),0.)
	816	CALL addpnt(x1,y1,kdata,n2, 1.e+38,0.)
	817	CALL inter2(nw,wl,yg,n2,x1,y1,ierr)
	818	IF (ierr .NE. 0) THEN
	819	WRITE(,) ierr, fil
	820	STOP
	821	ENDIF
	822
	823	DO l = 1, nw-1
	824	xsco2_295(l) = yg(l)
	825	END DO
	826
	827	! 370K:
	828
	829	fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_370k.txt'
	830	print*, 'section efficace CO2 370K: ', fil
	831
	832	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	833	DO i = 1,11
	834	read(kin,*)
	835	END DO
	836
	837	DO i = 1, n3
	838	READ(kin,*) x1(i), y1(i)
	839	END DO
	840	CLOSE (kin)
	841
	842	CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.)
	843	CALL addpnt(x1,y1,kdata,n3, 0.,0.)
	844	CALL addpnt(x1,y1,kdata,n3,x1(n3)*(1.+deltax),0.)
	845	CALL addpnt(x1,y1,kdata,n3, 1.e+38,0.)
	846	CALL inter2(nw,wl,yg,n3,x1,y1,ierr)
	847	IF (ierr .NE. 0) THEN
	848	WRITE(,) ierr, fil
	849	STOP
	850	ENDIF
	851
	852	DO l = 1, nw-1
	853	xsco2_370(l) = yg(l)
	854	END DO
	855
	856	! photodissociation yield:
	857
	858	fil = trim(datadir)//'/cross_sections/efdis_co2-o2_schunkandnagy2000.txt'
	859	print*, 'photodissociation yield CO2: ', fil
	860
	861	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	862
	863	do i = 1,3
	864	read(kin,*)
	865	end do
	866
	867	n4 = 17
	868	do i = 1, n4
	869	read(kin,*) xl, xu, ion(i)
	870	xion(i) = (xl + xu)/2.
	871	ion(i) = max(ion(i), 0.)
	872	end do
	873	close(kin)
	874
	875	CALL addpnt(xion,ion,kdata,n4,xion(1)*(1.-deltax),0.)
	876	CALL addpnt(xion,ion,kdata,n4, 0.,0.)
	877	CALL addpnt(xion,ion,kdata,n4,xion(n4)*(1.+deltax),1.)
	878	CALL addpnt(xion,ion,kdata,n4, 1.e+38,1.)
	879	CALL inter2(nw,wl,yieldco2,n4,xion,ion,ierr)
	880	IF (ierr .NE. 0) THEN
	881	WRITE(,) ierr, fil
	882	STOP
	883	ENDIF
	884
	885	! DO l = 1, nw-1
	886	! write(kout,*) wl(l), xsco2_195(l),
	887	! $ xsco2_295(l),
	888	! $ xsco2_370(l),
	889	! $ yieldco2(l)
	890	! END DO
	891
	892	end subroutine rdxsco2
	893
	894	!==============================================================================
	895
	896	subroutine rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2)
	897
	898	!-----------------------------------------------------------------------------*
	899	!= PURPOSE: =*
	900	!= Read and grid O2 cross-sections and photodissociation yield =*
	901	!-----------------------------------------------------------------------------*
	902	!= PARAMETERS: =*
	903	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	904	!= wavelength grid =*
	905	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
	906	!= working wavelength grid =*
	907	!= XSO2 - REAL, molecular absorption cross section =*
	908	!-----------------------------------------------------------------------------*
	909
	910	use datafile_mod, only: datadir
	911
	912	implicit none
	913
	914	! input
	915
	916	integer :: nw ! number of wavelength grid points
	917	real, dimension(nw) :: wl ! lower and central wavelength for each interval
	918
	919	! output
	920
	921	real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 cross-sections (cm2)
	922	real, dimension(nw) :: yieldo2 ! o2 photodissociation yield
	923
	924	! local
	925
	926	integer, parameter :: kdata = 18000
	927	real, parameter :: deltax = 1.e-4
	928	real, dimension(kdata) :: x1, y1, x2, y2, x3, y3, x4, y4
	929	real, dimension(kdata) :: xion, ion
	930	real :: factor, xl, xu, dummy
	931	integer :: i, ierr, n, n1, n2, n3, n4, nhead
	932	integer :: kin, kout ! input/output logical units
	933	character*100 fil
	934
	935	kin = 10
	936	kout = 30
	937
	938	! read o2 cross section data
	939
	940	nhead = 22
	941	n = 17434
	942
	943	fil = trim(datadir)//'/cross_sections/o2_composite_2018_150K.txt'
	944	print*, 'section efficace O2 150K: ', fil
[2044]	945	open(kin, file=fil, status='old', iostat=ierr)
[2029]	946
[2044]	947	if (ierr /= 0) THEN
	948	write(,)'cant find O2 cross-sections : ', fil
	949	write(,)'It should be in :', trim(datadir),'/cross_sections'
	950	write(,)'1) You can change this directory address in '
	951	write(,)' callphys.def with datadir=/path/to/dir'
	952	write(,)'2) If necessary, /cross_sections (and other datafiles)'
	953	write(,)' can be obtained online on:'
	954	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
	955	stop
	956	end if
	957
[2029]	958	DO i = 1,nhead
	959	read(kin,*)
	960	END DO
	961
	962	n1 = n
	963	DO i = 1, n1
	964	READ(kin,*) x1(i), y1(i)
	965	END DO
	966	CLOSE (kin)
	967
	968	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
	969	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
	970	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
	971	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
	972	CALL inter2(nw,wl,xso2_150,n1,x1,y1,ierr)
	973	IF (ierr .NE. 0) THEN
	974	WRITE(,) ierr, fil
	975	STOP
	976	ENDIF
	977
	978	fil = trim(datadir)//'/cross_sections/o2_composite_2018_200K.txt'
	979	print*, 'section efficace O2 200K: ', fil
	980	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	981
	982	DO i = 1,nhead
	983	read(kin,*)
	984	END DO
	985
	986	n2 = n
	987	DO i = 1, n2
	988	READ(kin,*) x2(i), y2(i)
	989	END DO
	990	CLOSE (kin)
	991
	992	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
	993	CALL addpnt(x2,y2,kdata,n2, 0.,0.)
	994	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
	995	CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.)
	996	CALL inter2(nw,wl,xso2_200,n2,x2,y2,ierr)
	997	IF (ierr .NE. 0) THEN
	998	WRITE(,) ierr, fil
	999	STOP
	1000	ENDIF
	1001
	1002	fil = trim(datadir)//'/cross_sections/o2_composite_2018_250K.txt'
	1003	print*, 'section efficace O2 250K: ', fil
	1004	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1005
	1006	DO i = 1,nhead
	1007	read(kin,*)
	1008	END DO
	1009
	1010	n3 = n
	1011	DO i = 1, n3
	1012	READ(kin,*) x3(i), y3(i)
	1013	END DO
	1014	CLOSE (kin)
	1015
	1016	CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
	1017	CALL addpnt(x3,y3,kdata,n3, 0.,0.)
	1018	CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
	1019	CALL addpnt(x3,y3,kdata,n3, 1.e+38,0.)
	1020	CALL inter2(nw,wl,xso2_250,n3,x3,y3,ierr)
	1021	IF (ierr .NE. 0) THEN
	1022	WRITE(,) ierr, fil
	1023	STOP
	1024	ENDIF
	1025
	1026	fil = trim(datadir)//'/cross_sections/o2_composite_2018_300K.txt'
	1027	print*, 'section efficace O2 300K: ', fil
	1028	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1029
	1030	DO i = 1,nhead
	1031	read(kin,*)
	1032	END DO
	1033
	1034	n4 = n
	1035	DO i = 1, n4
	1036	READ(kin,*) x4(i), y4(i)
	1037	END DO
	1038	CLOSE (kin)
	1039
	1040	CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),0.)
	1041	CALL addpnt(x4,y4,kdata,n4, 0.,0.)
	1042	CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),0.)
	1043	CALL addpnt(x4,y4,kdata,n4, 1.e+38,0.)
	1044	CALL inter2(nw,wl,xso2_300,n4,x4,y4,ierr)
	1045	IF (ierr .NE. 0) THEN
	1046	WRITE(,) ierr, fil
	1047	STOP
	1048	ENDIF
	1049
	1050	! photodissociation yield
	1051
	1052	fil = trim(datadir)//'/cross_sections/efdis_co2-o2_schunkandnagy2000.txt'
	1053	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1054
	1055	do i = 1,11
	1056	read(kin,*)
	1057	end do
	1058
	1059	n = 9
	1060	DO i = 1, n
	1061	READ(kin,*) xl, xu, dummy, ion(i)
	1062	xion(i) = (xl + xu)/2.
	1063	ion(i) = max(ion(i), 0.)
	1064	END DO
	1065	CLOSE (kin)
	1066
	1067	CALL addpnt(xion,ion,kdata,n,xion(1)*(1.-deltax),0.)
	1068	CALL addpnt(xion,ion,kdata,n, 0.,0.)
	1069	CALL addpnt(xion,ion,kdata,n,xion(n)*(1.+deltax),1.)
	1070	CALL addpnt(xion,ion,kdata,n, 1.e+38,1.)
	1071	CALL inter2(nw,wl,yieldo2,n,xion,ion,ierr)
	1072	IF (ierr .NE. 0) THEN
	1073	WRITE(,) ierr, fil
	1074	STOP
	1075	ENDIF
	1076
	1077	end subroutine rdxso2
	1078
	1079	!==============================================================================
	1080
	1081	subroutine rdxso3(nw,wl,xso3_218,xso3_298)
	1082
	1083	!-----------------------------------------------------------------------------*
	1084	!= PURPOSE: =*
	1085	!= Read ozone molecular absorption cross section. Re-grid data to match =*
	1086	!= specified wavelength working grid. =*
	1087	!-----------------------------------------------------------------------------*
	1088	!= PARAMETERS: =*
	1089	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1090	!= wavelength grid =*
	1091	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
	1092	!= working wavelength grid =*
	1093	!= XSO3_218 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=*
	1094	!= each specified wavelength (JPL 2006) 218 K =*
	1095	!= XSO3_298 REAL, molecular absoprtion cross section (cm^2) of O3 at (O)=*
	1096	!= each specified wavelength (JPL 2006) 298 K =*
	1097	!-----------------------------------------------------------------------------*
	1098
	1099	use datafile_mod, only: datadir
	1100
	1101	implicit none
	1102
	1103	! input
	1104
	1105	integer :: nw ! number of wavelength grid points
	1106	real, dimension(nw) :: wl ! lower and central wavelength for each interval
	1107
	1108	! output
	1109
	1110	real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2)
	1111
	1112	! local
	1113
	1114	integer, parameter :: kdata = 200
	1115	real, parameter :: deltax = 1.e-4
	1116	real, dimension(kdata) :: x1, x2, y1, y2
	1117	real, dimension(nw) :: yg
	1118	real :: a1, a2
	1119
	1120	integer :: i, ierr, iw, n, n1, n2
	1121	integer :: kin, kout ! input/output logical units
	1122
	1123	character*100 fil
	1124
[2043]	1125	kin = 10
	1126
[2029]	1127	! JPL 2006 218 K
	1128
	1129	fil = trim(datadir)//'/cross_sections/o3_cross-sections_jpl_2006_218K.txt'
	1130	print*, 'section efficace O3 218K: ', fil
	1131
	1132	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1133	n1 = 167
	1134	DO i = 1, n1
	1135	READ(kin,*) a1, a2, y1(i)
	1136	x1(i) = (a1+a2)/2.
	1137	END DO
	1138	CLOSE (kin)
	1139
	1140	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
	1141	CALL addpnt(x1,y1,kdata,n1, 0.,0.)
	1142	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
	1143	CALL addpnt(x1,y1,kdata,n1, 1.e+38,0.)
	1144	CALL inter2(nw,wl,yg,n1,x1,y1,ierr)
	1145	IF (ierr .NE. 0) THEN
	1146	WRITE(,) ierr, fil
	1147	STOP
	1148	ENDIF
	1149
	1150	DO iw = 1, nw-1
	1151	xso3_218(iw) = yg(iw)
	1152	END DO
	1153
	1154	! JPL 2006 298 K
	1155
	1156	fil = trim(datadir)//'/cross_sections/o3_cross-sections_jpl_2006_298K.txt'
	1157	print*, 'section efficace O3 298K: ', fil
	1158
	1159	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1160	n2 = 167
	1161	DO i = 1, n2
	1162	READ(kin,*) a1, a2, y2(i)
	1163	x2(i) = (a1+a2)/2.
	1164	END DO
	1165	CLOSE (kin)
	1166
	1167	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
	1168	CALL addpnt(x2,y2,kdata,n2, 0.,0.)
	1169	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
	1170	CALL addpnt(x2,y2,kdata,n2, 1.e+38,0.)
	1171	CALL inter2(nw,wl,yg,n2,x2,y2,ierr)
	1172	IF (ierr .NE. 0) THEN
	1173	WRITE(,) ierr, fil
	1174	STOP
	1175	ENDIF
	1176
	1177	DO iw = 1, nw-1
	1178	xso3_298(iw) = yg(iw)
	1179	END DO
	1180
	1181	end subroutine rdxso3
	1182
	1183	!==============================================================================
	1184
	1185	subroutine rdxsh2o(nw, wl, yg)
	1186
	1187	!-----------------------------------------------------------------------------*
	1188	!= PURPOSE: =*
	1189	!= Read H2O molecular absorption cross section. Re-grid data to match =*
	1190	!= specified wavelength working grid. =*
	1191	!-----------------------------------------------------------------------------*
	1192	!= PARAMETERS: =*
	1193	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1194	!= wavelength grid =*
	1195	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
	1196	!= working wavelength grid =*
	1197	!= YG - REAL, molecular absoprtion cross section (cm^2) of H2O at (O)=*
	1198	!= each specified wavelength =*
	1199	!-----------------------------------------------------------------------------*
	1200
	1201	use datafile_mod, only: datadir
	1202
	1203	IMPLICIT NONE
	1204
	1205	! input
	1206
	1207	integer :: nw ! number of wavelength grid points
	1208	real, dimension(nw) :: wl ! lower and central wavelength for each interval
	1209
	1210	! output
	1211
	1212	real, dimension(nw) :: yg ! h2o cross-sections (cm2)
	1213
	1214	! local
	1215
	1216	integer, parameter :: kdata = 500
	1217	real, parameter :: deltax = 1.e-4
	1218	REAL x1(kdata)
	1219	REAL y1(kdata)
	1220	INTEGER ierr
	1221	INTEGER i, n
	1222	CHARACTER*100 fil
	1223	integer :: kin, kout ! input/output logical units
	1224
	1225	kin = 10
	1226
	1227	fil = trim(datadir)//'/cross_sections/h2o_composite_250K.txt'
	1228	print*, 'section efficace H2O: ', fil
	1229
	1230	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1231
	1232	DO i = 1,26
	1233	read(kin,*)
	1234	END DO
	1235
	1236	n = 420
	1237	DO i = 1, n
	1238	READ(kin,*) x1(i), y1(i)
	1239	END DO
	1240	CLOSE (kin)
	1241
	1242	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
	1243	CALL addpnt(x1,y1,kdata,n, 0.,0.)
	1244	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
	1245	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
	1246	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
	1247	IF (ierr .NE. 0) THEN
	1248	WRITE(,) ierr, fil
	1249	STOP
	1250	ENDIF
	1251
	1252	end subroutine rdxsh2o
	1253
	1254	!==============================================================================
	1255
	1256	subroutine rdxsh2o2(nw, wl, xsh2o2)
	1257
	1258	!-----------------------------------------------------------------------------*
	1259	!= PURPOSE: =*
	1260	!= Read and grid H2O2 cross-sections
	1261	!= H2O2 + hv -> 2 OH =*
	1262	!= Cross section: Schuergers and Welge, Z. Naturforsch. 23a (1968) 1508 =*
	1263	!= from 125 to 185 nm, then JPL97 from 190 to 350 nm. =*
	1264	!-----------------------------------------------------------------------------*
	1265	!= PARAMETERS: =*
	1266	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1267	!= wavelength grid =*
	1268	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
	1269	!= working wavelength grid =*
	1270	!-----------------------------------------------------------------------------*
	1271
	1272	use datafile_mod, only: datadir
	1273
	1274	implicit none
	1275
	1276	! input
	1277
	1278	integer :: nw ! number of wavelength grid points
	1279	real, dimension(nw) :: wl ! lower wavelength for each interval
	1280
	1281	! output
	1282
	1283	real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2)
	1284
	1285	! local
	1286
	1287	real, parameter :: deltax = 1.e-4
	1288	integer, parameter :: kdata = 100
	1289	real, dimension(kdata) :: x1, y1
	1290	real, dimension(nw) :: yg
	1291	integer :: i, ierr, iw, n, idum
	1292	integer :: kin, kout ! input/output logical units
	1293	character*100 fil
	1294
	1295	kin = 10
	1296
	1297	! read cross-sections
	1298
	1299	fil = trim(datadir)//'/cross_sections/h2o2_composite.txt'
	1300	print*, 'section efficace H2O2: ', fil
	1301
	1302	OPEN(kin,FILE=fil,STATUS='OLD')
	1303	READ(kin,*) idum,n
	1304	DO i = 1, idum-2
	1305	READ(kin,*)
	1306	ENDDO
	1307	DO i = 1, n
	1308	READ(kin,*) x1(i), y1(i)
	1309	ENDDO
	1310	CLOSE (kin)
	1311
	1312	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
	1313	CALL addpnt(x1,y1,kdata,n, 0.,0.)
	1314	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
	1315	CALL addpnt(x1,y1,kdata,n, 1.e+38,0.)
	1316	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
	1317	IF (ierr .NE. 0) THEN
	1318	WRITE(,) ierr, fil
	1319	STOP
	1320	ENDIF
	1321
	1322	DO iw = 1, nw - 1
	1323	xsh2o2(iw) = yg(iw)
	1324	END DO
	1325
	1326	end subroutine rdxsh2o2
	1327
	1328	!==============================================================================
	1329
	1330	subroutine rdxsho2(nw, wl, yg)
	1331
	1332	!-----------------------------------------------------------------------------*
	1333	!= PURPOSE: =*
	1334	!= Read ho2 cross-sections =*
	1335	!= JPL 2006 recommendation =*
	1336	!-----------------------------------------------------------------------------*
	1337	!= PARAMETERS: =*
	1338	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1339	!= wavelength grid =*
	1340	!-----------------------------------------------------------------------------*
	1341
	1342	use datafile_mod, only: datadir
	1343
	1344	IMPLICIT NONE
	1345
	1346	! input
	1347
	1348	integer :: nw ! number of wavelength grid points
	1349	real, dimension(nw) :: wl ! lower wavelength for each interval
	1350
	1351	! output
	1352
	1353	real, dimension(nw) :: yg ! ho2 cross-sections (cm2)
	1354
	1355	! local
	1356
	1357	real, parameter :: deltax = 1.e-4
	1358	integer, parameter :: kdata = 100
	1359	real, dimension(kdata) :: x1, y1
	1360	integer :: i, n, ierr
	1361	character*100 fil
	1362	integer :: kin, kout ! input/output logical units
	1363
	1364	kin = 10
	1365
	1366	!*** cross sections from Sander et al. [2003]
	1367
	1368	fil = trim(datadir)//'/cross_sections/ho2_jpl2003.txt'
	1369	print*, 'section efficace HO2: ', fil
	1370
	1371	OPEN(kin,FILE=fil,STATUS='OLD')
	1372	READ(kin,*) n
	1373	DO i = 1, n
	1374	READ(kin,*) x1(i), y1(i)
	1375	ENDDO
	1376	CLOSE(kin)
	1377
	1378	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
	1379	CALL addpnt(x1,y1,kdata,n, 0.,0.)
	1380	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
	1381	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
	1382
	1383	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
	1384
	1385	IF (ierr .NE. 0) THEN
	1386	WRITE(,) ierr, fil
	1387	STOP
	1388	ENDIF
	1389
	1390	end subroutine rdxsho2
	1391
	1392	!==============================================================================
	1393
	1394	subroutine rdxsh2(nw, wl, wc, yg, yieldh2)
	1395
	1396	!-----------------------------------------------------------------------------*
	1397	!= PURPOSE: =*
	1398	!= Read h2 cross-sections and photodissociation yield =*
	1399	!-----------------------------------------------------------------------------*
	1400	!= PARAMETERS: =*
	1401	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1402	!= wavelength grid =*
	1403	!-----------------------------------------------------------------------------*
	1404
	1405	use datafile_mod, only: datadir
	1406
	1407	implicit none
	1408
	1409	! input
	1410
	1411	integer :: nw ! number of wavelength grid points
	1412	real, dimension(nw) :: wl, wc ! lower and central wavelength for each interval
	1413
	1414	! output
	1415
	1416	real, dimension(nw) :: yg ! h2 cross-sections (cm2)
	1417	real, dimension(nw) :: yieldh2 ! photodissociation yield
	1418
	1419	! local
	1420
	1421	integer, parameter :: kdata = 1000
	1422	real, parameter :: deltax = 1.e-4
	1423	real, dimension(kdata) :: x1, y1, x2, y2
	1424	real :: xl, xu
	1425	integer :: i, iw, n, ierr
	1426	integer :: kin, kout ! input/output logical units
	1427	character*100 fil
	1428
	1429	kin = 10
	1430
	1431	! h2 cross sections
	1432
	1433	fil = trim(datadir)//'/cross_sections/h2secef.txt'
	1434	print*, 'section efficace H2: ', fil
	1435
	1436	OPEN(kin,FILE=fil,STATUS='OLD')
	1437
	1438	n = 792
	1439	read(kin,*) ! avoid first line with wavelength = 0.
	1440	DO i = 1, n
	1441	READ(kin,*) x1(i), y1(i)
	1442	ENDDO
	1443	CLOSE(kin)
	1444
	1445	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
	1446	CALL addpnt(x1,y1,kdata,n, 0.,0.)
	1447	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
	1448	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
	1449	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
	1450
	1451	IF (ierr .NE. 0) THEN
	1452	WRITE(,) ierr, fil
	1453	STOP
	1454	ENDIF
	1455
	1456	! photodissociation yield
	1457
	1458	fil = trim(datadir)//'/cross_sections/h2_ionef_schunknagy2000.txt'
	1459	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1460
	1461	n = 19
	1462	read(kin,*)
	1463	DO i = 1, n
	1464	READ(kin,*) xl, xu, y2(i)
	1465	x2(i) = (xl + xu)/2.
	1466	y2(i) = max(1. - y2(i),0.)
	1467	END DO
	1468	CLOSE (kin)
	1469
	1470	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
	1471	CALL addpnt(x2,y2,kdata,n, 0.,0.)
	1472	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
	1473	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
	1474	CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr)
	1475	IF (ierr .NE. 0) THEN
	1476	WRITE(,) ierr, fil
	1477	STOP
	1478	ENDIF
	1479
	1480	end subroutine rdxsh2
	1481
	1482	!==============================================================================
	1483
	1484	subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,yldno2_248, yldno2_298)
	1485
	1486	!-----------------------------------------------------------------------------*
	1487	!= PURPOSE: =*
	1488	!= read and grid cross section + quantum yield for NO2 =*
	1489	!= photolysis =*
	1490	!= Jenouvrier et al., 1996 200-238 nm
	1491	!= Vandaele et al., 1998 238-666 nm 220K and 294K
	1492	!= quantum yield from jpl 2006
	1493	!-----------------------------------------------------------------------------*
	1494	!= PARAMETERS: =*
	1495	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1496	!= wavelength grid =*
	1497	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
	1498	!= working wavelength grid =*
	1499	!= SQ - REAL, cross section x quantum yield (cm^2) for each (O)=*
	1500	!= photolysis reaction defined, at each defined wavelength and =*
	1501	!= at each defined altitude level =*
	1502	!-----------------------------------------------------------------------------*
	1503
	1504	use datafile_mod, only: datadir
	1505
	1506	implicit none
	1507
	1508	! input
	1509
	1510	integer :: nw ! number of wavelength grid points
	1511	real, dimension(nw) :: wl ! lower and central wavelength for each interval
	1512
	1513	! output
	1514
	1515	real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2)
	1516	real, dimension(nw) :: yldno2_248, yldno2_298 ! quantum yields at 248-298 k
	1517
	1518	! local
	1519
	1520	integer, parameter :: kdata = 28000
	1521	real, parameter :: deltax = 1.e-4
	1522	real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5
	1523	real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5
	1524	real :: dum, qy
	1525	integer :: i, iw, n, n1, n2, n3, n4, n5, ierr
	1526	character*100 fil
	1527	integer :: kin, kout ! input/output logical units
	1528
	1529	kin = 10
	1530
	1531	!*************** NO2 photodissociation
	1532
	1533	! Jenouvrier 1996 + Vandaele 1998 (JPL 2006)
	1534
	1535	fil = trim(datadir)//'/cross_sections/no2_xs_jenouvrier.txt'
	1536	print*, 'section efficace NO2: ', fil
	1537
	1538	OPEN(UNIT=kin,FILE=fil,status='old')
	1539	DO i = 1, 3
	1540	READ(kin,*)
	1541	ENDDO
	1542	n1 = 10001
	1543	DO i = 1, n1
	1544	READ(kin,*) x1(i), y1(i)
	1545	end do
	1546
	1547	CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax), 0.)
	1548	CALL addpnt(x1,y1,kdata,n1, 0., 0.)
	1549	CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
	1550	CALL addpnt(x1,y1,kdata,n1, 1.e+38, 0.)
	1551	CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
	1552
	1553	fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_294K.txt'
	1554	print*, 'section efficace NO2: ', fil
	1555
	1556	OPEN(UNIT=kin,FILE=fil,status='old')
	1557	DO i = 1, 3
	1558	READ(kin,*)
	1559	ENDDO
	1560	n2 = 27993
	1561	DO i = 1, n2
	1562	READ(kin,*) x2(i), y2(i)
	1563	end do
	1564
	1565	CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax), 0.)
	1566	CALL addpnt(x2,y2,kdata,n2, 0., 0.)
	1567	CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
	1568	CALL addpnt(x2,y2,kdata,n2, 1.e+38, 0.)
	1569	CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
	1570
	1571	fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_220K.txt'
	1572	print*, 'section efficace NO2: ', fil
	1573
	1574	OPEN(UNIT=kin,FILE=fil,status='old')
	1575	DO i = 1, 3
	1576	READ(kin,*)
	1577	ENDDO
	1578	n3 = 27993
	1579	DO i = 1, n3
	1580	READ(kin,*) x3(i), y3(i)
	1581	end do
	1582
	1583	CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.)
	1584	CALL addpnt(x3,y3,kdata,n3, 0., 0.)
	1585	CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
	1586	CALL addpnt(x3,y3,kdata,n3, 1.e+38, 0.)
	1587	CALL inter2(nw,wl,yg3,n3,x3,y3,ierr)
	1588
	1589	do iw = 1, nw - 1
	1590	xsno2(iw) = yg1(iw)
	1591	xsno2_294(iw) = yg2(iw)
	1592	xsno2_220(iw) = yg3(iw)
	1593	end do
	1594
	1595	! photodissociation efficiency from jpl 2006
	1596
	1597	fil = trim(datadir)//'/cross_sections/no2_yield_jpl2006.txt'
	1598	print*, 'quantum yield NO2: ', fil
	1599
	1600	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1601	DO i = 1, 5
	1602	READ(kin,*)
	1603	ENDDO
	1604	n = 25
	1605	n4 = n
	1606	n5 = n
	1607	DO i = 1, n
	1608	READ(kin,*) x4(i), y4(i), y5(i)
	1609	x5(i) = x4(i)
	1610	ENDDO
	1611	CLOSE(kin)
	1612
	1613	CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1))
	1614	CALL addpnt(x4,y4,kdata,n4, 0.,y4(1))
	1615	CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax), 0.)
	1616	CALL addpnt(x4,y4,kdata,n4, 1.e+38, 0.)
	1617	CALL inter2(nw,wl,yg4,n4,x4,y4,ierr)
	1618	IF (ierr .NE. 0) THEN
	1619	WRITE(,) ierr, fil
	1620	STOP
	1621	ENDIF
	1622
	1623	CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1))
	1624	CALL addpnt(x5,y5,kdata,n5, 0.,y5(1))
	1625	CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax), 0.)
	1626	CALL addpnt(x5,y5,kdata,n5, 1.e+38, 0.)
	1627	CALL inter2(nw,wl,yg5,n5,x5,y5,ierr)
	1628	IF (ierr .NE. 0) THEN
	1629	WRITE(,) ierr, fil
	1630	STOP
	1631	ENDIF
	1632
	1633	do iw = 1, nw - 1
	1634	yldno2_298(iw) = yg4(iw)
	1635	yldno2_248(iw) = yg5(iw)
	1636	end do
	1637
	1638	end subroutine rdxsno2
	1639
	1640	!==============================================================================
	1641
	1642	subroutine rdxsno(nw, wl, yg, yieldno)
	1643
	1644	!-----------------------------------------------------------------------------*
	1645	!= PURPOSE: =*
	1646	!= Read NO cross-sections and photodissociation efficiency =*
	1647	!= Lida et al 1986 (provided by Francisco Gonzalez-Galindo) =*
	1648	!-----------------------------------------------------------------------------*
	1649	!= PARAMETERS: =*
	1650	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1651	!= wavelength grid =*
	1652	!-----------------------------------------------------------------------------*
	1653
	1654	use datafile_mod, only: datadir
	1655
	1656	implicit none
	1657
	1658	! input
	1659
	1660	integer :: nw ! number of wavelength grid points
	1661	real, dimension(nw) :: wl ! lower wavelength for each interval
	1662
	1663	! output
	1664
	1665	real, dimension(nw) :: yg ! no cross-sections (cm2)
	1666	real, dimension(nw) :: yieldno ! no photodissociation efficiency
	1667
	1668	! local
	1669
	1670	integer, parameter :: kdata = 110
	1671	real, parameter :: deltax = 1.e-4
	1672	real, dimension(kdata) :: x1, y1, x2, y2
	1673	integer :: i, iw, n, ierr
	1674	character*100 fil
	1675	integer :: kin, kout ! input/output logical units
	1676
	1677	kin = 10
	1678
	1679	! no cross-sections
	1680
	1681	fil = trim(datadir)//'/cross_sections/no_xs_francisco.txt'
	1682	print*, 'section efficace NO: ', fil
	1683	OPEN(kin,FILE=fil,STATUS='OLD')
	1684
	1685	n = 99
	1686	DO i = 1, n
	1687	READ(kin,*) x1(i), y1(i)
	1688	ENDDO
	1689	CLOSE(kin)
	1690
	1691	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
	1692	CALL addpnt(x1,y1,kdata,n, 0.,0.)
	1693	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
	1694	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
	1695
	1696	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
	1697	IF (ierr .NE. 0) THEN
	1698	WRITE(,) ierr, fil
	1699	STOP
	1700	ENDIF
	1701
	1702	! photodissociation yield
	1703
	1704	fil = trim(datadir)//'/cross_sections/noefdis.txt'
	1705	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1706
	1707	n = 33
	1708	DO i = 1, n
	1709	READ(kin,*) x2(n-i+1), y2(n-i+1)
	1710	END DO
	1711	CLOSE (kin)
	1712
	1713	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
	1714	CALL addpnt(x2,y2,kdata,n, 0.,0.)
	1715	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
	1716	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
	1717	CALL inter2(nw,wl,yieldno,n,x2,y2,ierr)
	1718	IF (ierr .NE. 0) THEN
	1719	WRITE(,) ierr, fil
	1720	STOP
	1721	ENDIF
	1722
	1723	end subroutine rdxsno
	1724
	1725	!==============================================================================
	1726
	1727	subroutine rdxsn2(nw, wl, yg, yieldn2)
	1728
	1729	!-----------------------------------------------------------------------------*
	1730	!= PURPOSE: =*
	1731	!= Read n2 cross-sections and photodissociation yield =*
	1732	!-----------------------------------------------------------------------------*
	1733	!= PARAMETERS: =*
	1734	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1735	!= wavelength grid =*
	1736	!-----------------------------------------------------------------------------*
	1737
	1738	use datafile_mod, only: datadir
	1739
	1740	implicit none
	1741
	1742	! input
	1743
	1744	integer :: nw ! number of wavelength grid points
	1745	real, dimension(nw) :: wl ! lower wavelength for each interval
	1746
	1747	! output
	1748
	1749	real, dimension(nw) :: yg ! n2 cross-sections (cm2)
	1750	real, dimension(nw) :: yieldn2 ! n2 photodissociation yield
	1751
	1752	! local
	1753
	1754	integer, parameter :: kdata = 1100
	1755	real, parameter :: deltax = 1.e-4
	1756	real, dimension(kdata) :: x1, y1, x2, y2
	1757	real :: xl, xu
	1758	integer :: i, iw, n, ierr
	1759	integer :: kin, kout ! input/output logical units
	1760	character*100 fil
	1761
	1762	kin = 10
	1763
	1764	! n2 cross sections
	1765
	1766	fil = trim(datadir)//'/cross_sections/n2secef_01nm.txt'
	1767	print*, 'section efficace N2: ', fil
	1768	OPEN(kin,FILE=fil,STATUS='OLD')
	1769
	1770	n = 1020
	1771	DO i = 1, n
	1772	READ(kin,*) x1(i), y1(i)
	1773	ENDDO
	1774	CLOSE(kin)
	1775
	1776	CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
	1777	CALL addpnt(x1,y1,kdata,n, 0.,0.)
	1778	CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
	1779	CALL addpnt(x1,y1,kdata,n, 1E38,0.)
	1780
	1781	CALL inter2(nw,wl,yg,n,x1,y1,ierr)
	1782
	1783	IF (ierr .NE. 0) THEN
	1784	WRITE(,) ierr, fil
	1785	STOP
	1786	ENDIF
	1787
	1788	! photodissociation yield
	1789
	1790	fil = trim(datadir)//'/cross_sections/n2_ionef_schunknagy2000.txt'
	1791	OPEN(UNIT=kin,FILE=fil,STATUS='old')
	1792
	1793	n = 19
	1794	read(kin,*)
	1795	DO i = 1, n
	1796	READ(kin,*) xl, xu, y2(i)
	1797	x2(i) = (xl + xu)/2.
	1798	y2(i) = 1. - y2(i)
	1799	END DO
	1800	CLOSE (kin)
	1801
	1802	CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
	1803	CALL addpnt(x2,y2,kdata,n, 0.,0.)
	1804	CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
	1805	CALL addpnt(x2,y2,kdata,n, 1.e+38,1.)
	1806	CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr)
	1807	IF (ierr .NE. 0) THEN
	1808	WRITE(,) ierr, fil
	1809	STOP
	1810	ENDIF
	1811
	1812	end subroutine rdxsn2
	1813
	1814	!==============================================================================
	1815
	1816	subroutine setalb(nw,wl,albedo)
	1817
	1818	!-----------------------------------------------------------------------------*
	1819	!= PURPOSE: =*
	1820	!= Set the albedo of the surface. The albedo is assumed to be Lambertian, =*
	1821	!= i.e., the reflected light is isotropic, and idependt of the direction =*
	1822	!= of incidence of light. Albedo can be chosen to be wavelength dependent. =*
	1823	!-----------------------------------------------------------------------------*
	1824	!= PARAMETERS: =*
	1825	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1826	!= wavelength grid =*
	1827	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
	1828	!= working wavelength grid =*
	1829	!= ALBEDO - REAL, surface albedo at each specified wavelength (O)=*
	1830	!-----------------------------------------------------------------------------*
	1831
	1832	implicit none
	1833
	1834	! input: (wavelength working grid data)
	1835
	1836	INTEGER nw
	1837	REAL wl(nw)
	1838
	1839	! output:
	1840
	1841	REAL albedo(nw)
	1842
	1843	! local:
	1844
	1845	INTEGER iw
	1846	REAL alb
	1847
	1848	! 0.015: mean value from clancy et al., icarus, 49-63, 1999.
	1849
	1850	alb = 0.015
	1851
	1852	do iw = 1, nw - 1
	1853	albedo(iw) = alb
	1854	end do
	1855
	1856	end subroutine setalb
	1857
	1858	end module photolysis_mod

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

Download in other formats:

Original Format