Context Navigation

photolysis_mod.F90 @ 3026

Last change on this file since 3026 was 2614, checked in by romain.vande, 3 years ago
LMDZ_MARS RV : Open_MP; Correction for r2613 correct file is photolysis_mod Reading files in parallel for PHOTOCHEMISTRY parametrisation (photochem = .true.)
File size: 63.1 KB

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

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

Download in other formats:

Original Format