Context Navigation

photolysis_online.F @ 2150

Last change on this file since 2150 was 2041, checked in by flefevre, 7 years ago
Photolyse on-line : Mise a jour de la section efficace Rayleigh (Itiaksov et al., 2008)
File size: 51.0 KB

Rev	Line
[2027]	1	!==============================================================================
	2
[2028]	3	subroutine photolysis_online(nlayer, alt, press, temp,
[2027]	4	$ zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
	5	$ i_h2, i_oh, i_ho2, i_h2o2, i_h2o,
	6	$ i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,
	7	$ tau, sza, dist_sol, v_phot)
	8
[2029]	9	use photolysis_mod
	10
[2027]	11	implicit none
	12
	13	! input
	14
[2029]	15	integer, intent(in) :: nesp ! total number of chemical species
	16	integer, intent(in) :: nlayer
	17	integer, intent(in) :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
	18	$ i_h2, i_oh, i_ho2, i_h2o2, i_h2o,
	19	$ i_n, i_n2d, i_no, i_no2, i_n2
[2027]	20
[2029]	21	real, dimension(nlayer), intent(in) :: press, temp, zmmean ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
	22	real, dimension(nlayer), intent(in) :: alt ! altitude (km)
	23	real, dimension(nlayer,nesp), intent(in) :: rm ! mixing ratios
	24	real, intent(in) :: tau ! integrated aerosol optical depth at the surface
	25	real, intent(in) :: sza ! solar zenith angle (degrees)
	26	real, intent(in) :: dist_sol ! solar distance (au)
[2027]	27
	28	! output
	29
	30	real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot ! photolysis rates (s-1)
	31
[2029]	32	! solar flux at mars
[2027]	33
[2029]	34	real, dimension(nw) :: fmars ! solar flux (w.m-2.nm-1)
[2027]	35	real :: factor
	36
[2029]	37	! cross-sections
[2027]	38
[2028]	39	real, dimension(nlayer,nw,nphot) :: sj ! general cross-section array (cm2)
[2027]	40
	41	! atmosphere
	42
	43	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	44	real, dimension(nlayer,nw) :: dtrl ! rayleigh optical depth
	45	real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth
	46	real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo
	47	real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter
	48	real, dimension(nlayer,nw) :: dtcld ! cloud optical depth
	49	real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo
	50	real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter
	51	real, dimension(nlayer,nw,nabs) :: dtgas ! optical depth for each gas
	52	real, dimension(nlayer,nw) :: dagas ! total gas optical depth
	53	real, dimension(nlayer) :: edir, edn, eup ! normalised irradiances
	54	real, dimension(nlayer) :: fdir, fdn, fup ! normalised actinic fluxes
	55	real, dimension(nlayer) :: saflux ! total actinic flux
	56
	57	integer, dimension(0:nlayer) :: nid
	58	real, dimension(0:nlayer,nlayer) :: dsdh
	59
	60	integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
	61	$ j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2
	62
	63	integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no,
	64	$ a_no2, a_n2
[2028]	65	integer :: i, ilay, iw, ialt
[2027]	66	real :: deltaj
	67
	68	! absorbing gas numbering
	69
	70	a_o2 = 1 ! o2
	71	a_co2 = 2 ! co2
	72	a_o3 = 3 ! o3
	73	a_h2o = 4 ! h2o
	74	a_h2o2 = 5 ! h2o2
	75	a_ho2 = 6 ! ho2
	76	a_h2 = 7 ! h2
	77	a_no = 8 ! no
	78	a_no2 = 9 ! no2
	79	a_n2 = 10 ! no2
	80
[2028]	81	! photodissociation rates numbering.
	82	! photodissociations must be ordered the same way in subroutine "indice"
[2027]	83
	84	j_o2_o = 1 ! o2 + hv -> o + o
	85	j_o2_o1d = 2 ! o2 + hv -> o + o(1d)
	86	j_co2_o = 3 ! co2 + hv -> co + o
	87	j_co2_o1d = 4 ! co2 + hv -> co + o(1d)
	88	j_o3_o1d = 5 ! o3 + hv -> o2 + o(1d)
	89	j_o3_o = 6 ! o3 + hv -> o2 + o
	90	j_h2o = 7 ! h2o + hv -> h + oh
	91	j_h2o2 = 8 ! h2o2 + hv -> oh + oh
	92	j_ho2 = 9 ! ho2 + hv -> oh + o
	93	j_h2 = 10 ! h2 + hv -> h + h
	94	j_no = 11 ! no + hv -> n + o
	95	j_no2 = 12 ! no2 + hv -> no + o
	96	j_n2 = 13 ! n2 + hv -> n + n
	97
	98	!==== air column increments and rayleigh optical depth
	99
	100	call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl)
	101
	102	!==== set temperature-dependent cross-sections and optical depths
	103
	104	! o2:
	105
[2028]	106	call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,
[2027]	107	$ xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d,
	108	$ colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj)
	109
	110	! co2:
	111
[2028]	112	call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295,
[2027]	113	$ xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
	114	$ colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj)
	115
	116	! o3:
	117
[2028]	118	call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298,
[2027]	119	$ j_o3_o, j_o3_o1d,
	120	$ colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj)
	121
	122	! h2o2:
	123
[2028]	124	call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
[2027]	125	$ colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj)
	126
	127	! no2:
	128
[2028]	129	call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
[2027]	130	$ xsno2_294, yldno2_248, yldno2_298, j_no2,
	131	$ colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj)
	132
	133	!==== temperature independent optical depths and cross-sections
	134
	135	! optical depths
	136
	137	do ilay = 1,nlayer
	138	do iw = 1,nw-1
	139	dtgas(ilay,iw,a_h2o) = colinc(ilay)rm(ilay,i_h2o)xsh2o(iw)
	140	dtgas(ilay,iw,a_ho2) = colinc(ilay)rm(ilay,i_ho2)xsho2(iw)
	141	dtgas(ilay,iw,a_h2) = colinc(ilay)rm(ilay,i_h2)xsh2(iw)
	142	dtgas(ilay,iw,a_no) = colinc(ilay)rm(ilay,i_no)xsno(iw)
	143	dtgas(ilay,iw,a_n2) = colinc(ilay)rm(ilay,i_n2)xsn2(iw)
	144	end do
	145	end do
	146
	147	! total gas optical depth
	148
	149	dagas(:,:) = 0.
	150
	151	do ilay = 1,nlayer
	152	do iw = 1,nw-1
	153	do i = 1,nabs
	154	dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i)
	155	end do
	156	end do
	157	end do
	158
	159	! cross-sections
	160
	161	do ilay = 1,nlayer
	162	do iw = 1,nw-1
	163	sj(ilay,iw,j_h2o) = xsh2o(iw) ! h2o
	164	sj(ilay,iw,j_ho2) = xsho2(iw) ! ho2
	165	sj(ilay,iw,j_h2) = xsh2(iw)*yieldh2(iw) ! h2
	166	sj(ilay,iw,j_no) = xsno(iw)*yieldno(iw) ! no
	167	sj(ilay,iw,j_n2) = xsn2(iw)*yieldn2(iw) ! n2
	168	end do
	169	end do
	170
	171	!==== set aerosol properties and optical depth
	172
	173	call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
	174
	175	!==== set cloud properties and optical depth
	176
	177	call setcld(nlayer,nw,dtcld,omcld,gcld)
	178
	179	!==== slant path lengths in spherical geometry
	180
[2029]	181	call sphers(nlayer,alt,sza,dsdh,nid)
[2027]	182
	183	!==== solar flux at mars
	184
	185	factor = (1./dist_sol)**2.
	186	do iw = 1,nw-1
[2029]	187	fmars(iw) = f(iw)*factor
[2027]	188	end do
	189
	190	!==== initialise photolysis rates
	191
[2028]	192	v_phot(:,1:nphot) = 0.
[2027]	193
	194	!==== start of wavelength lopp
	195
	196	do iw = 1,nw-1
	197
	198	! monochromatic radiative transfer. outputs are:
	199	! normalized irradiances edir(nlayer), edn(nlayer), eup(nlayer)
	200	! normalized actinic fluxes fdir(nlayer), fdn(nlayer), fup(nlayer)
	201	! where
	202	! dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse
	203
	204	call rtlink(nlayer, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
	205	$ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
	206	$ edir, edn, eup, fdir, fdn, fup)
	207
	208
	209	! spherical actinic flux
	210
[2028]	211	do ilay = 1,nlayer
[2029]	212	saflux(ilay) = fmars(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))
[2027]	213	end do
	214
	215	! photolysis rate integration
	216
[2028]	217	do i = 1,nphot
	218	do ilay = 1,nlayer
[2027]	219	deltaj = saflux(ilay)*sj(ilay,iw,i)
	220	v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw))
	221	end do
	222	end do
	223
	224	! eliminate small values
	225
[2028]	226	where (v_phot(:,1:nphot) < 1.e-30)
	227	v_phot(:,1:nphot) = 0.
[2027]	228	end where
	229
	230	end do ! iw
	231
[2029]	232	contains
[2027]	233
	234	!==============================================================================
	235
	236	subroutine setair(nlev, nw, wl, wc, press, temp, zmmean,
	237	$ colinc, dtrl)
	238
	239	-----------------------------------------------------------------------------
	240	= PURPOSE: =
	241	= computes air column increments and rayleigh optical depth =
	242	-----------------------------------------------------------------------------
	243
	244	implicit none
	245
	246	! input:
	247
	248	integer :: nlev, nw
	249
	250	real, dimension(nw) :: wl, wc ! lower and central wavelength grid (nm)
	251	real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1)
	252
	253	! output:
	254
	255	real, dimension(nlev) :: colinc ! air column increments (molecule.cm-2)
	256	real, dimension(nlev,nw) :: dtrl ! rayleigh optical depth
	257
	258	! local:
	259
	260	real, parameter :: avo = 6.022e23
	261	real, parameter :: g = 3.72
[2041]	262	real :: dp, nu
[2027]	263	real, dimension(nw) :: srayl
	264	integer :: ilev, iw
	265
[2029]	266	! compute column increments
[2027]	267
	268	do ilev = 1, nlev-1
	269	dp = (press(ilev) - press(ilev+1))*100.
	270	colinc(ilev) = avo0.1dp/(zmmean(ilev)*g)
	271	end do
	272	colinc(nlev) = avo0.1press(nlev)100./(zmmean(nlev)g)
	273
	274	do iw = 1, nw - 1
	275
[2041]	276	! co2 rayleigh cross-section
	277	! ityaksov et al., chem. phys. lett., 462, 31-34, 2008
[2027]	278
[2041]	279	nu = 1./(wc(iw)*1.e-7)
	280	srayl(iw) = 1.78e-26nu*(4. + 0.625)
	281	srayl(iw) = srayl(iw)*1.e-20 ! cm2
[2027]	282
	283	do ilev = 1, nlev
[2041]	284	dtrl(ilev,iw) = colinc(ilev)*srayl(iw) ! cm2
[2027]	285	end do
	286	end do
	287
[2029]	288	end subroutine setair
[2027]	289
	290	!==============================================================================
	291
	292	subroutine setco2(nd, nlayer, nw, wc, tlay, xsco2_195, xsco2_295,
	293	$ xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
	294	$ colinc, rm, dt, sj)
	295
[2029]	296	!-----------------------------------------------------------------------------*
	297	!= PURPOSE: =*
	298	!= Set up the CO2 temperature-dependent cross-sections and optical depth =*
	299	!-----------------------------------------------------------------------------*
[2027]	300
	301	implicit none
	302
	303	! input:
	304
	305	integer :: nd ! number of photolysis rates
	306	integer :: nlayer ! number of vertical layers
	307	integer :: nw ! number of wavelength grid points
	308	integer :: j_co2_o, j_co2_o1d ! photolysis indexes
	309	real, dimension(nw) :: wc ! central wavelength for each interval
	310	real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
	311	real, dimension(nw) :: yieldco2 ! co2 photodissociation yield
	312	real, dimension(nlayer) :: tlay ! temperature (k)
	313	real, dimension(nlayer) :: rm ! co2 mixing ratio
	314	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	315
	316	! output:
	317
	318	real, dimension(nlayer,nw) :: dt ! optical depth
	319	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	320
	321	! local:
	322
	323	integer :: extrapol
	324	integer :: i, l
	325	real :: temp, sco2
	326
	327	! extrapol = 0 no extrapolation below 195 k
	328	! extrapol = 1 extrapolation below 195 k
	329
	330	extrapol = 0
	331
	332	do i = 1, nlayer
	333	if (extrapol == 1) then
	334	temp = tlay(i)
	335	else
	336	temp = max(tlay(i), 195.)
	337	end if
	338	temp = min(temp, 370.)
	339	do l = 1, nw-1
	340	if (temp <= 295.) then
	341	if (xsco2_195(l) /= 0. .and. xsco2_295(l) /= 0.) then
	342	sco2 = alog(xsco2_195(l))
	343	$ + (alog(xsco2_295(l)) - alog(xsco2_195(l)))
	344	$ /(295. - 195.)*(temp - 195.)
	345	sco2 = exp(sco2)
	346	else
	347	sco2 = 0.
	348	end if
	349	else
	350	if (xsco2_295(l) /= 0. .and. xsco2_370(l) /= 0.) then
	351	sco2 = alog(xsco2_295(l))
	352	$ + (alog(xsco2_370(l)) - alog(xsco2_295(l)))
	353	$ /(370. - 295.)*(temp - 295.)
	354	sco2 = exp(sco2)
	355	else
	356	sco2 = 0.
	357	end if
	358	end if
	359
	360	! optical depth
	361
	362	dt(i,l) = colinc(i)rm(i)sco2
	363
	364	! production of o(1d) for wavelengths shorter than 167 nm
	365
	366	if (wc(l) >= 167.) then
	367	sj(i,l,j_co2_o) = sco2*yieldco2(l)
	368	sj(i,l,j_co2_o1d) = 0.
	369	else
	370	sj(i,l,j_co2_o) = 0.
	371	sj(i,l,j_co2_o1d) = sco2*yieldco2(l)
	372	end if
	373	end do
	374	end do
	375
	376	end subroutine setco2
	377
	378	!==============================================================================
	379
	380	subroutine seto2(nd, nlayer, nw, wc, mopt, tlay, xso2_150,
	381	$ xso2_200, xso2_250, xso2_300, yieldo2, j_o2_o,
	382	$ j_o2_o1d, colinc, rm, dt, sj)
	383
[2029]	384	!-----------------------------------------------------------------------------*
	385	!= PURPOSE: =*
	386	!= Set up the O2 temperature-dependent cross-sections and optical depth =*
	387	!-----------------------------------------------------------------------------*
[2027]	388
	389	implicit none
	390
	391	! input:
	392
	393	integer :: nd ! number of photolysis rates
	394	integer :: nlayer ! number of vertical layers
	395	integer :: nw ! number of wavelength grid points
	396	integer :: mopt ! high-res/low-res switch
	397	integer :: j_o2_o, j_o2_o1d ! photolysis indexes
	398	real, dimension(nw) :: wc ! central wavelength for each interval
	399	real, dimension(nw) :: xso2_150, xso2_200, xso2_250, ! o2 cross-sections (cm2)
	400	$ xso2_300
	401	real, dimension(nw) :: yieldo2 ! o2 photodissociation yield
	402	real, dimension(nlayer) :: tlay ! temperature (k)
	403	real, dimension(nlayer) :: rm ! o2 mixing ratio
	404	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	405
	406	! output:
	407
	408	real, dimension(nlayer,nw) :: dt ! optical depth
	409	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	410
	411	! local:
	412
	413	integer :: ilev, iw
	414	real :: temp
	415	real :: xso2, factor
	416
	417	! correction by factor if low-resolution in schumann-runge bands
	418
	419	if (mopt == 1) then
	420	factor = 1.
	421	else if (mopt == 2) then
	422	factor = 0.8
	423	end if
	424
	425	! calculate temperature dependance
	426
	427	do ilev = 1,nlayer
	428	temp = max(tlay(ilev),150.)
	429	temp = min(temp, 300.)
	430	do iw = 1, nw-1
	431	if (tlay(ilev) > 250.) then
	432	xso2 = xso2_250(iw) + (xso2_300(iw) - xso2_250(iw))
	433	$ /(300. - 250.)*(temp - 250.)
	434	else if (tlay(ilev) > 200.) then
	435	xso2 = xso2_200(iw) + (xso2_250(iw) - xso2_200(iw))
	436	$ /(250. - 200.)*(temp - 200.)
	437	else
	438	xso2 = xso2_150(iw) + (xso2_200(iw) - xso2_150(iw))
	439	$ /(200. - 150.)*(temp - 150.)
	440	end if
	441
	442	if (wc(iw) > 180. .and. wc(iw) < 200.) then
	443	xso2 = xso2*factor
	444	end if
	445
	446	! optical depth
	447
	448	dt(ilev,iw) = colinc(ilev)rm(ilev)xso2
	449
	450	! production of o(1d) for wavelengths shorter than 175 nm
	451
	452	if (wc(iw) >= 175.) then
	453	sj(ilev,iw,j_o2_o) = xso2*yieldo2(iw)
	454	sj(ilev,iw,j_o2_o1d) = 0.
	455	else
	456	sj(ilev,iw,j_o2_o) = 0.
	457	sj(ilev,iw,j_o2_o1d) = xso2*yieldo2(iw)
	458	end if
	459
	460	end do
	461	end do
	462
	463	end subroutine seto2
	464
	465	!==============================================================================
	466
	467	subroutine seto3(nd, nlayer, nw, wc, tlay, xso3_218, xso3_298,
	468	$ j_o3_o, j_o3_o1d,
	469	$ colinc, rm, dt, sj)
	470
[2029]	471	!-----------------------------------------------------------------------------*
	472	!= PURPOSE: =*
	473	!= Set up the O3 temperature dependent cross-sections and optical depth =*
	474	!-----------------------------------------------------------------------------*
[2027]	475
	476	implicit none
	477
	478	! input:
	479
	480	integer :: nd ! number of photolysis rates
	481	integer :: nlayer ! number of vertical layers
	482	integer :: nw ! number of wavelength grid points
	483	integer :: j_o3_o, j_o3_o1d ! photolysis indexes
	484	real, dimension(nw) :: wc ! central wavelength for each interval
	485	real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2)
	486	real, dimension(nlayer) :: tlay ! temperature (k)
	487	real, dimension(nlayer) :: rm ! o3 mixing ratio
	488	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	489
	490	! output:
	491
	492	real, dimension(nlayer,nw) :: dt ! optical depth
	493	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	494
	495	! local:
	496	!
	497	integer :: ilev, iw
	498	real :: temp
	499	real, dimension(nw) :: xso3(nw)
	500	real, dimension(nw) :: qy1d ! quantum yield for o(1d) production
	501	real :: q1, q2, a1, a2, a3
	502
	503	do ilev = 1, nlayer
	504	temp = max(tlay(ilev), 218.)
	505	temp = min(temp,298.)
	506	do iw = 1, nw-1
	507	xso3(iw) = xso3_218(iw) + (xso3_298(iw) - xso3_218(iw))
	508	$ /(298. - 218.) *(temp - 218.)
	509
	510	! optical depth
	511
	512	dt(ilev,iw) = colinc(ilev)rm(ilev)xso3(iw)
	513
	514	end do
	515
	516	! calculate quantum yield for o(1d) production (jpl 2006)
	517
	518	temp = max(tlay(ilev),200.)
	519	temp = min(temp,320.)
	520	do iw = 1, nw-1
	521	if (wc(iw) <= 306.) then
	522	qy1d(iw) = 0.90
	523	else if (wc(iw) > 306. .and. wc(iw) < 328.) then
	524	q1 = 1.
	525	q2 = exp(-825.518/(0.695*temp))
	526	a1 = (304.225 - wc(iw))/5.576
	527	a2 = (314.957 - wc(iw))/6.601
	528	a3 = (310.737 - wc(iw))/2.187
	529	qy1d(iw) = (q1/(q1 + q2))0.8036exp(-(a1a1a1*a1))
	530	$ + (q2/(q1 + q2))8.9061(temp/300.)**2.
	531	$ exp(-(a2a2))
	532	$ + 0.1192(temp/300.)1.5exp(-(a3*a3))
	533	$ + 0.0765
	534	else if (wc(iw) >= 328. .and. wc(iw) <= 340.) then
	535	qy1d(iw) = 0.08
	536	else
	537	qy1d(iw) = 0.
	538	endif
	539	end do
	540	do iw = 1, nw-1
	541	sj(ilev,iw,j_o3_o) = xso3(iw)*(1. - qy1d(iw))
	542	sj(ilev,iw,j_o3_o1d) = xso3(iw)*qy1d(iw)
	543	end do
	544	end do
	545
	546	end subroutine seto3
	547
	548	!==============================================================================
	549
	550	subroutine seth2o2(nd, nlayer, nw, wc, tlay, xsh2o2, j_h2o2,
	551	$ colinc, rm, dt, sj)
	552
[2029]	553	!-----------------------------------------------------------------------------*
	554	!= PURPOSE: =*
	555	!= Set up the h2o2 temperature dependent cross-sections and optical depth =*
	556	!-----------------------------------------------------------------------------*
[2027]	557
	558	implicit none
	559
	560	! input:
	561
	562	integer :: nd ! number of photolysis rates
	563	integer :: nlayer ! number of vertical layers
	564	integer :: nw ! number of wavelength grid points
	565	integer :: j_h2o2 ! photolysis index
	566	real, dimension(nw) :: wc ! central wavelength for each interval
	567	real, dimension(nw) :: xsh2o2 ! h2o2 cross-sections (cm2)
	568	real, dimension(nlayer) :: tlay ! temperature (k)
	569	real, dimension(nlayer) :: rm ! h2o2 mixing ratio
	570	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	571
	572	! output:
	573
	574	real, dimension(nlayer,nw) :: dt ! optical depth
	575	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	576
	577	! local:
	578
	579	integer :: ilev, iw
	580	real :: a0, a1, a2, a3, a4, a5, a6, a7
	581	real :: b0, b1, b2, b3, b4
	582	real :: lambda, suma, sumb, chi, temp, xs
	583
	584	A0 = 6.4761E+04
	585	A1 = -9.2170972E+02
	586	A2 = 4.535649
	587	A3 = -4.4589016E-03
	588	A4 = -4.035101E-05
	589	A5 = 1.6878206E-07
	590	A6 = -2.652014E-10
	591	A7 = 1.5534675E-13
	592
	593	B0 = 6.8123E+03
	594	B1 = -5.1351E+01
	595	B2 = 1.1522E-01
	596	B3 = -3.0493E-05
	597	B4 = -1.0924E-07
	598
	599	! temperature dependance: jpl 2006
	600
	601	do ilev = 1,nlayer
	602	temp = min(max(tlay(ilev),200.),400.)
	603	chi = 1./(1. + exp(-1265./temp))
	604	do iw = 1, nw-1
	605	if ((wc(iw) >= 260.) .and. (wc(iw) < 350.)) then
	606	lambda = wc(iw)
	607	sumA = ((((((A7lambda + A6)lambda + A5)*lambda +
	608	$ A4)lambda +A3)lambda + A2)*lambda +
	609	$ A1)*lambda + A0
	610	sumB = (((B4lambda + B3)lambda + B2)*lambda +
	611	$ B1)*lambda + B0
	612	xs = (chisumA + (1. - chi)sumB)*1.e-21
	613	sj(ilev,iw,j_h2o2) = xs
	614	else
	615	sj(ilev,iw,j_h2o2) = xsh2o2(iw)
	616	end if
	617
	618	! optical depth
	619
	620	dt(ilev,iw) = colinc(ilev)rm(ilev)sj(ilev,iw,j_h2o2)
	621	end do
	622	end do
	623
	624	end subroutine seth2o2
	625
	626	!==============================================================================
	627
	628	subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220,
	629	$ xsno2_294, yldno2_248, yldno2_298, j_no2,
	630	$ colinc, rm, dt, sj)
	631
[2029]	632	!-----------------------------------------------------------------------------*
	633	!= PURPOSE: =*
	634	!= Set up the no2 temperature-dependent cross-sections and optical depth =*
	635	!-----------------------------------------------------------------------------*
[2027]	636
	637	implicit none
	638
	639	! input:
	640
	641	integer :: nd ! number of photolysis rates
	642	integer :: nlayer ! number of vertical layers
	643	integer :: nw ! number of wavelength grid points
	644	integer :: j_no2 ! photolysis index
	645	real, dimension(nw) :: wc ! central wavelength for each interval
	646	real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2)
	647	real, dimension(nw) :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k
	648	real, dimension(nlayer) :: tlay ! temperature (k)
	649	real, dimension(nlayer) :: rm ! no2 mixing ratio
	650	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	651
	652	! output:
	653
	654	real, dimension(nlayer,nw) :: dt ! optical depth
	655	real, dimension(nlayer,nw,nd) :: sj ! cross-section array (cm2)
	656
	657	! local:
	658
	659	integer :: ilev, iw
	660	real :: temp, qy
	661
	662	! temperature dependance: jpl 2006
	663
	664	do ilev = 1,nlayer
	665	temp = max(220.,min(tlay(ilev),294.))
	666	do iw = 1, nw - 1
	667	if (wc(iw) < 238.) then
	668	sj(ilev,iw,j_no2) = xsno2(iw)
	669	else
	670	sj(ilev,iw,j_no2) = xsno2_220(iw)
	671	$ + (xsno2_294(iw) - xsno2_220(iw))
	672	$ /(294. - 220.)*(temp - 220.)
	673	end if
	674
	675	! optical depth
	676
	677	dt(ilev,iw) = colinc(ilev)rm(ilev)sj(ilev,iw,j_no2)
	678	end do
	679	end do
	680
	681	! quantum yield: jpl 2006
	682
	683	do ilev = 1,nlayer
	684	temp = max(248.,min(tlay(ilev),298.))
	685	do iw = 1, nw - 1
	686	qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw))
	687	$ /(298. - 248.)*(temp - 248.)
	688	sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy
	689	end do
	690	end do
	691
	692	end subroutine setno2
	693
	694	!==============================================================================
	695
	696	subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
	697
[2029]	698	!-----------------------------------------------------------------------------*
	699	!= PURPOSE: =*
	700	!= Set aerosol properties for each specified altitude layer. Properties =*
	701	!= may be wavelength dependent. =*
	702	!-----------------------------------------------------------------------------*
[2027]	703
	704	implicit none
	705
	706	! input
	707
	708	integer :: nlayer ! number of vertical layers
	709	integer :: nw ! number of wavelength grid points
	710	real, dimension(nlayer) :: alt ! altitude (km)
	711	real :: tau ! integrated aerosol optical depth at the surface
	712
	713	! output
	714
	715	real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth
	716	real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo
	717	real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter
	718
	719	! local
	720
	721	integer :: ilay, iw
	722	real, dimension(nlayer) :: aer ! dust extinction
	723	real :: omega, g, scaleh, gamma
	724	real :: dz, tautot, q0
	725
	726	omega = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm
	727	g = 0.88 ! asymmetry factor : mateshvili et al. (2007) at 210 nm
	728	scaleh = 10. ! scale height (km)
	729	gamma = 0.03 ! conrath parameter
	730
	731	dtaer(:,:) = 0.
	732	omaer(:,:) = 0.
	733	gaer(:,:) = 0.
	734
	735	! optical depth profile:
	736
	737	tautot = 0.
	738	do ilay = 1, nlayer-1
	739	dz = alt(ilay+1) - alt(ilay)
	740	tautot = tautot + exp(gamma(1. - exp(alt(ilay)/scaleh)))dz
	741	end do
	742
	743	q0 = tau/tautot
	744	do ilay = 1, nlayer-1
	745	dz = alt(ilay+1) - alt(ilay)
	746	dtaer(ilay,:) = q0exp(gamma(1. - exp(alt(ilay)/scaleh)))*dz
	747	omaer(ilay,:) = omega
	748	gaer(ilay,:) = g
	749	end do
	750
	751	end subroutine setaer
	752
	753	!==============================================================================
	754
	755	subroutine setcld(nlayer,nw,dtcld,omcld,gcld)
	756
[2029]	757	!-----------------------------------------------------------------------------*
	758	!= PURPOSE: =*
	759	!= Set cloud properties for each specified altitude layer. Properties =*
	760	!= may be wavelength dependent. =*
	761	!-----------------------------------------------------------------------------*
[2027]	762
	763	implicit none
	764
	765	! input
	766
	767	integer :: nlayer ! number of vertical layers
	768	integer :: nw ! number of wavelength grid points
	769
	770	! output
	771
	772	real, dimension(nlayer,nw) :: dtcld ! cloud optical depth
	773	real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo
	774	real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter
	775
	776	! local
	777
	778	integer :: ilay, iw
	779
[2029]	780	! dtcld : optical depth
	781	! omcld : single scattering albedo
	782	! gcld : asymmetry factor
[2027]	783
	784	do ilay = 1, nlayer - 1
	785	do iw = 1, nw - 1
	786	dtcld(ilay,iw) = 0. ! no clouds for the moment
	787	omcld(ilay,iw) = 0.99
	788	gcld(ilay,iw) = 0.85
	789	end do
	790	end do
	791
[2029]	792	end subroutine setcld
[2027]	793
	794	!==============================================================================
	795
	796	subroutine sphers(nlev, z, zen, dsdh, nid)
	797
[2029]	798	!-----------------------------------------------------------------------------*
	799	!= PURPOSE: =*
	800	!= Calculate slant path over vertical depth ds/dh in spherical geometry. =*
	801	!= Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model =*
	802	!= for computing the radiation field available for photolysis and heating =*
	803	!= at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) =*
	804	!-----------------------------------------------------------------------------*
	805	!= PARAMETERS: =*
	806	!= NZ - INTEGER, number of specified altitude levels in the working (I)=*
	807	!= grid =*
	808	!= Z - REAL, specified altitude working grid (km) (I)=*
	809	!= ZEN - REAL, solar zenith angle (degrees) (I)=*
	810	!= DSDH - REAL, slant path of direct beam through each layer crossed (O)=*
	811	!= when travelling from the top of the atmosphere to layer i; =*
	812	!= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*
	813	!= NID - INTEGER, number of layers crossed by the direct beam when (O)=*
	814	!= travelling from the top of the atmosphere to layer i; =*
	815	!= NID(i), i = 0..NZ-1 =*
	816	!-----------------------------------------------------------------------------*
	817
[2027]	818	implicit none
	819
[2029]	820	! input
[2027]	821
[2029]	822	integer, intent(in) :: nlev
	823	real, dimension(nlev), intent(in) :: z
	824	real, intent(in) :: zen
	825
	826	! output
	827
[2027]	828	INTEGER nid(0:nlev)
	829	REAL dsdh(0:nlev,nlev)
	830
[2029]	831	! more program constants
	832
[2027]	833	REAL re, ze(nlev)
	834	REAL dr
	835	real radius
	836	parameter (radius = 3393.)
	837
[2029]	838	! local
[2027]	839
[2029]	840	real :: pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
	841	integer :: i, j, k, id, nlay
[2027]	842
	843	REAL zd(0:nlev-1)
	844
[2029]	845	!-----------------------------------------------------------------------------
[2027]	846
	847	pi = acos(-1.0)
	848	dr = pi/180.
	849	zenrad = zen*dr
	850
[2029]	851	! number of layers:
	852
[2027]	853	nlay = nlev - 1
	854
[2029]	855	! include the elevation above sea level to the radius of Mars:
	856
[2027]	857	re = radius + z(1)
[2029]	858
	859	! correspondingly z changed to the elevation above Mars surface:
	860
[2027]	861	DO k = 1, nlev
	862	ze(k) = z(k) - z(1)
	863	END DO
	864
[2029]	865	! inverse coordinate of z
	866
[2027]	867	zd(0) = ze(nlev)
	868	DO k = 1, nlay
	869	zd(k) = ze(nlev - k)
	870	END DO
	871
[2029]	872	! initialise dsdh(i,j), nid(i)
[2027]	873
[2029]	874	nid(:) = 0.
	875	dsdh(:,:) = 0.
[2027]	876
[2029]	877	! calculate ds/dh of every layer
	878
	879	do i = 0,nlay
	880	rpsinz = (re + zd(i))*sin(zenrad)
[2027]	881
	882	IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
	883	nid(i) = -1
	884	ELSE
	885
[2029]	886	! Find index of layer in which the screening height lies
	887
[2027]	888	id = i
[2029]	889	if (zen > 90.) then
	890	do j = 1,nlay
[2027]	891	IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND.
	892	$ (rpsinz .GE. ( zd(j) + re )) ) id = j
[2029]	893	end do
	894	end if
[2027]	895
[2029]	896	do j = 1,id
	897	sm = 1.0
	898	IF (j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0)
	899	$ sm = -1.0
[2027]	900
[2029]	901	rj = re + zd(j-1)
	902	rjp1 = re + zd(j)
[2027]	903
[2029]	904	dhj = zd(j-1) - zd(j)
[2027]	905
[2029]	906	ga = rjrj - rpsinzrpsinz
	907	gb = rjp1rjp1 - rpsinzrpsinz
[2027]	908
[2030]	909	ga = max(ga, 0.)
	910	gb = max(gb, 0.)
[2027]	911
[2029]	912	IF (id.GT.i .AND. j.EQ.id) THEN
	913	dsj = sqrt(ga)
	914	ELSE
	915	dsj = sqrt(ga) - sm*sqrt(gb)
	916	END IF
	917	dsdh(i,j) = dsj/dhj
	918	end do
	919	nid(i) = id
	920	end if
	921	end do ! i = 0,nlay
[2027]	922
[2029]	923	end subroutine sphers
[2027]	924
	925	!==============================================================================
	926
	927	SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl,
	928	$ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
	929	$ edir, edn, eup, fdir, fdn, fup)
	930
	931	implicit none
	932
[2029]	933	! input
[2027]	934
[2029]	935	integer, intent(in) :: nlev, nw, iw ! number of wavelength grid points
[2027]	936	REAL ag
	937	REAL zen
	938	REAL dsdh(0:nlev,nlev)
	939	INTEGER nid(0:nlev)
	940
	941	REAL dtrl(nlev,nw)
	942	REAL dagas(nlev,nw)
	943	REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw)
	944	REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw)
	945
[2029]	946	! output
[2027]	947
	948	REAL edir(nlev), edn(nlev), eup(nlev)
	949	REAL fdir(nlev), fdn(nlev), fup(nlev)
	950
[2029]	951	! local:
[2027]	952
	953	REAL dt(nlev), om(nlev), g(nlev)
	954	REAL dtabs,dtsct,dscld,dsaer,dacld,daaer
	955	INTEGER i, ii
	956	real, parameter :: largest = 1.e+36
	957
[2029]	958	! specific two ps2str
[2027]	959
	960	REAL ediri(nlev), edni(nlev), eupi(nlev)
	961	REAL fdiri(nlev), fdni(nlev), fupi(nlev)
	962
	963	logical, save :: delta = .true.
	964
[2029]	965	!_______________________________________________________________________
[2027]	966
[2029]	967	! initialize:
[2027]	968
	969	do i = 1, nlev
	970	fdir(i) = 0.
	971	fup(i) = 0.
	972	fdn(i) = 0.
	973	edir(i) = 0.
	974	eup(i) = 0.
	975	edn(i) = 0.
	976	end do
[2029]	977
[2027]	978	do i = 1, nlev - 1
	979	dscld = dtcld(i,iw)*omcld(i,iw)
	980	dacld = dtcld(i,iw)*(1.-omcld(i,iw))
[2029]	981
[2027]	982	dsaer = dtaer(i,iw)*omaer(i,iw)
	983	daaer = dtaer(i,iw)*(1.-omaer(i,iw))
[2029]	984
[2027]	985	dtsct = dtrl(i,iw) + dscld + dsaer
	986	dtabs = dagas(i,iw) + dacld + daaer
[2029]	987
[2027]	988	dtabs = amax1(dtabs,1./largest)
	989	dtsct = amax1(dtsct,1./largest)
[2029]	990
	991	! invert z-coordinate:
	992
[2027]	993	ii = nlev - i
	994	dt(ii) = dtsct + dtabs
	995	om(ii) = dtsct/(dtsct + dtabs)
	996	IF(dtsct .EQ. 1./largest) om(ii) = 1./largest
	997	g(ii) = (gcld(i,iw)*dscld +
	998	$ gaer(i,iw)*dsaer)/dtsct
	999	end do
[2029]	1000
	1001	! call rt routine:
	1002
[2027]	1003	call ps2str(nlev, zen, ag, dt, om, g,
	1004	$ dsdh, nid, delta,
	1005	$ fdiri, fupi, fdni, ediri, eupi, edni)
[2029]	1006
	1007	! output (invert z-coordinate)
	1008
[2027]	1009	do i = 1, nlev
	1010	ii = nlev - i + 1
	1011	fdir(i) = fdiri(ii)
	1012	fup(i) = fupi(ii)
	1013	fdn(i) = fdni(ii)
	1014	edir(i) = ediri(ii)
	1015	eup(i) = eupi(ii)
	1016	edn(i) = edni(ii)
	1017	end do
	1018
[2029]	1019	end subroutine rtlink
	1020
[2027]	1021	=============================================================================
	1022
[2029]	1023	subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu,
[2027]	1024	$ dsdh, nid, delta,
	1025	$ fdr, fup, fdn, edr, eup, edn)
	1026
[2029]	1027	!-----------------------------------------------------------------------------*
	1028	!= PURPOSE: =*
	1029	!= Solve two-stream equations for multiple layers. The subroutine is based =*
	1030	!= on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
	1031	!= It contains 9 two-stream methods to choose from. A pseudo-spherical =*
	1032	!= correction has also been added. =*
	1033	!-----------------------------------------------------------------------------*
	1034	!= PARAMETERS: =*
	1035	!= NLEVEL - INTEGER, number of specified altitude levels in the working (I)=*
	1036	!= grid =*
	1037	!= ZEN - REAL, solar zenith angle (degrees) (I)=*
	1038	!= RSFC - REAL, surface albedo at current wavelength (I)=*
	1039	!= TAUU - REAL, unscaled optical depth of each layer (I)=*
	1040	!= OMU - REAL, unscaled single scattering albedo of each layer (I)=*
	1041	!= GU - REAL, unscaled asymmetry parameter of each layer (I)=*
	1042	!= DSDH - REAL, slant path of direct beam through each layer crossed (I)=*
	1043	!= when travelling from the top of the atmosphere to layer i; =*
	1044	!= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*
	1045	!= NID - INTEGER, number of layers crossed by the direct beam when (I)=*
	1046	!= travelling from the top of the atmosphere to layer i; =*
	1047	!= NID(i), i = 0..NZ-1 =*
	1048	!= DELTA - LOGICAL, switch to use delta-scaling (I)=*
	1049	!= .TRUE. -> apply delta-scaling =*
	1050	!= .FALSE.-> do not apply delta-scaling =*
	1051	!= FDR - REAL, contribution of the direct component to the total (O)=*
	1052	!= actinic flux at each altitude level =*
	1053	!= FUP - REAL, contribution of the diffuse upwelling component to (O)=*
	1054	!= the total actinic flux at each altitude level =*
	1055	!= FDN - REAL, contribution of the diffuse downwelling component to (O)=*
	1056	!= the total actinic flux at each altitude level =*
	1057	!= EDR - REAL, contribution of the direct component to the total (O)=*
	1058	!= spectral irradiance at each altitude level =*
	1059	!= EUP - REAL, contribution of the diffuse upwelling component to (O)=*
	1060	!= the total spectral irradiance at each altitude level =*
	1061	!= EDN - REAL, contribution of the diffuse downwelling component to (O)=*
[2027]	1062	= the total spectral irradiance at each altitude level =
[2029]	1063	!-----------------------------------------------------------------------------*
	1064
[2027]	1065	implicit none
[2029]	1066
	1067	! input:
	1068
[2027]	1069	INTEGER nlev
	1070	REAL zen, rsfc
	1071	REAL tauu(nlev), omu(nlev), gu(nlev)
	1072	REAL dsdh(0:nlev,nlev)
	1073	INTEGER nid(0:nlev)
	1074	LOGICAL delta
	1075
[2029]	1076	! output:
	1077
[2027]	1078	REAL fup(nlev),fdn(nlev),fdr(nlev)
	1079	REAL eup(nlev),edn(nlev),edr(nlev)
	1080
[2029]	1081	! local:
	1082
[2027]	1083	REAL tausla(0:nlev), tauc(0:nlev)
	1084	REAL mu2(0:nlev), mu, sum
	1085
[2029]	1086	! internal coefficients and matrix
	1087
[2027]	1088	REAL lam(nlev),taun(nlev),bgam(nlev)
	1089	REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev)
	1090	REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev)
	1091	REAL mu1(nlev)
	1092	INTEGER row
	1093	REAL a(2nlev),b(2nlev),d(2nlev),e(2nlev),y(2*nlev)
	1094
[2029]	1095	! other:
	1096
[2027]	1097	REAL pifs, fdn0
	1098	REAL gi(nlev), omi(nlev), tempg
	1099	REAL f, g, om
	1100	REAL gam1, gam2, gam3, gam4
	1101	real, parameter :: largest = 1.e+36
	1102	real, parameter :: precis = 1.e-7
	1103
[2029]	1104	! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
	1105	! in delta-function, modified quadrature, hemispheric constant,
	1106	! Hybrid modified Eddington-delta function metods, p633,Table1.
	1107	! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
	1108	! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,
	1109	! uncomment the following two lines and the appropriate statements further
	1110	! down.
	1111	! REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
	1112	! > BETA1, BETAn, amu1, subd
[2027]	1113
	1114	REAL expon, expon0, expon1, divisr, temp, up, dn
	1115	REAL ssfc
	1116	INTEGER nlayer, mrows, lev
	1117
	1118	INTEGER i, j
	1119
[2029]	1120	! Some additional program constants:
[2027]	1121
	1122	real pi, dr
	1123	REAL eps
	1124	PARAMETER (eps = 1.E-3)
[2029]	1125	!_______________________________________________________________________
[2027]	1126
[2029]	1127	! MU = cosine of solar zenith angle
	1128	! RSFC = surface albedo
	1129	! TAUU = unscaled optical depth of each layer
	1130	! OMU = unscaled single scattering albedo
	1131	! GU = unscaled asymmetry factor
	1132	! KLEV = max dimension of number of layers in atmosphere
	1133	! NLAYER = number of layers in the atmosphere
	1134	! NLEVEL = nlayer + 1 = number of levels
[2027]	1135
[2029]	1136	! initial conditions: pi*solar flux = 1; diffuse incidence = 0
[2027]	1137
	1138	pifs = 1.
	1139	fdn0 = 0.
	1140
	1141	nlayer = nlev - 1
	1142
	1143	pi = acos(-1.)
	1144	dr = pi/180.
	1145	mu = COS(zen*dr)
	1146
[2029]	1147	!************* compute coefficients for each layer:
	1148	! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
	1149	! EXPON0 = calculation of e when TAU is zero
	1150	! EXPON1 = calculation of e when TAU is TAUN
	1151	! CUP and CDN = calculation when TAU is zero
	1152	! CUPTN and CDNTN = calc. when TAU is TAUN
	1153	! DIVISR = prevents division by zero
[2027]	1154
	1155	do j = 0, nlev
	1156	tauc(j) = 0.
	1157	tausla(j) = 0.
	1158	mu2(j) = 1./SQRT(largest)
	1159	end do
[2029]	1160
[2027]	1161	IF (.NOT. delta) THEN
	1162	DO i = 1, nlayer
	1163	gi(i) = gu(i)
	1164	omi(i) = omu(i)
	1165	taun(i) = tauu(i)
	1166	END DO
	1167	ELSE
	1168
[2029]	1169	! delta-scaling. Have to be done for delta-Eddington approximation,
	1170	! delta discrete ordinate, Practical Improved Flux Method, delta function,
	1171	! and Hybrid modified Eddington-delta function methods approximations
[2027]	1172
	1173	DO i = 1, nlayer
	1174	f = gu(i)*gu(i)
	1175	gi(i) = (gu(i) - f)/(1 - f)
	1176	omi(i) = (1 - f)omu(i)/(1 - omu(i)f)
	1177	taun(i) = (1 - omu(i)f)tauu(i)
	1178	END DO
	1179	END IF
[2029]	1180
	1181	! calculate slant optical depth at the top of the atmosphere when zen>90.
	1182	! in this case, higher altitude of the top layer is recommended.
	1183
[2027]	1184	IF (zen .GT. 90.0) THEN
	1185	IF (nid(0) .LT. 0) THEN
	1186	tausla(0) = largest
	1187	ELSE
	1188	sum = 0.0
	1189	DO j = 1, nid(0)
	1190	sum = sum + 2.taun(j)dsdh(0,j)
	1191	END DO
	1192	tausla(0) = sum
	1193	END IF
	1194	END IF
[2029]	1195
[2027]	1196	DO 11, i = 1, nlayer
	1197	g = gi(i)
	1198	om = omi(i)
	1199	tauc(i) = tauc(i-1) + taun(i)
	1200
[2029]	1201	! stay away from 1 by precision. For g, also stay away from -1
[2027]	1202
	1203	tempg = AMIN1(abs(g),1. - precis)
	1204	g = SIGN(tempg,g)
	1205	om = AMIN1(om,1.-precis)
	1206
[2029]	1207	! calculate slant optical depth
	1208
[2027]	1209	IF (nid(i) .LT. 0) THEN
	1210	tausla(i) = largest
	1211	ELSE
	1212	sum = 0.0
	1213	DO j = 1, MIN(nid(i),i)
	1214	sum = sum + taun(j)*dsdh(i,j)
	1215	END DO
	1216	DO j = MIN(nid(i),i)+1,nid(i)
	1217	sum = sum + 2.taun(j)dsdh(i,j)
	1218	END DO
	1219	tausla(i) = sum
	1220	IF (tausla(i) .EQ. tausla(i-1)) THEN
	1221	mu2(i) = SQRT(largest)
	1222	ELSE
	1223	mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
	1224	mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)),
	1225	$ mu2(i) )
	1226	END IF
	1227	END IF
	1228
[2029]	1229	!** the following gamma equations are from pg 16,289, Table 1
	1230	!** save mu1 for each approx. for use in converting irradiance to actinic flux
[2027]	1231
[2029]	1232	! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
	1233
[2027]	1234	c gam1 = (7. - om(4. + 3.g))/4.
	1235	c gam2 = -(1. - om(4. - 3.g))/4.
	1236	c gam3 = (2. - 3.gmu)/4.
	1237	c gam4 = 1. - gam3
	1238	c mu1(i) = 0.5
	1239
	1240	* quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):
	1241
	1242	c gam1 = 1.7320508(2. - om(1. + g))/2.
	1243	c gam2 = 1.7320508om(1. - g)/2.
	1244	c gam3 = (1. - 1.7320508gmu)/2.
	1245	c gam4 = 1. - gam3
	1246	c mu1(i) = 1./sqrt(3.)
	1247
	1248	* hemispheric mean (Toon et al., 1089, JGR, 94, 16287):
	1249
	1250	gam1 = 2. - om*(1. + g)
	1251	gam2 = om*(1. - g)
	1252	gam3 = (2. - g*mu)/4.
	1253	gam4 = 1. - gam3
	1254	mu1(i) = 0.5
	1255
	1256	* PIFM (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
	1257	c GAM1 = 0.25(8. - OM(5. + 3.*G))
	1258	c GAM2 = 0.75OM(1.-G)
	1259	c GAM3 = 0.25(2.-3.G*MU)
	1260	c GAM4 = 1. - GAM3
	1261	c mu1(i) = 0.5
	1262
	1263	* delta discrete ordinates (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
	1264	c GAM1 = 0.51.7320508(2. - OM*(1. + G))
	1265	c GAM2 = 0.51.7320508OM*(1.-G)
	1266	c GAM3 = 0.5(1.-1.7320508G*MU)
	1267	c GAM4 = 1. - GAM3
	1268	c mu1(i) = 1./sqrt(3.)
	1269
	1270	* Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
	1271	* in delta-function, modified quadrature, hemispheric constant,
	1272	* Hybrid modified Eddington-delta function metods, p633,Table1.
	1273	* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
	1274	* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
	1275	c YLM0 = 2.
	1276	c YLM2 = -3.GMU
	1277	c YLM4 = 0.875G3MU(5.MU**2-3.)
	1278	c YLM6=-0.171875G5MU(15.-70.MU*2+63.MU**4)
	1279	c YLM8=+0.073242G7MU(-35.+315.MU*2-693.MU**4
	1280	c +429.MU**6)
	1281	c YLM10=-0.008118G9MU(315.-4620.MU*2+18018.MU**4
	1282	c -25740.MU*6+12155.MU**8)
	1283	c YLM12=0.003685G11MU(-693.+15015.MU*2-90090.MU**4
	1284	c +218790.MU*6-230945.MU*8+88179.MU**10)
	1285	c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
	1286	c YLMS=0.25*YLMS
	1287	c BETA0 = YLMS
	1288	c
	1289	c amu1=1./1.7320508
	1290	c YLM0 = 2.
	1291	c YLM2 = -3.Gamu1
	1292	c YLM4 = 0.875G3amu1(5.amu1**2-3.)
	1293	c YLM6=-0.171875G5amu1(15.-70.amu1*2+63.amu1**4)
	1294	c YLM8=+0.073242G7amu1(-35.+315.amu1*2-693.amu1**4
	1295	c +429.amu1**6)
	1296	c YLM10=-0.008118G9amu1(315.-4620.amu1*2+18018.amu1**4
	1297	c -25740.amu1*6+12155.amu1**8)
	1298	c YLM12=0.003685G11amu1(-693.+15015.amu1*2-90090.amu1**4
	1299	c +218790.amu1*6-230945.amu1*8+88179.amu1**10)
	1300	c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
	1301	c YLMS=0.25*YLMS
	1302	c BETA1 = YLMS
	1303	c
	1304	c BETAn = 0.25(2. - 1.5G-0.21875G3-0.085938G**5
	1305	c -0.045776G**7)
	1306
	1307
	1308	* Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
	1309	c subd=4.(1.-GG*(1.-MU))
	1310	c GAM1 = (7.-3.GG-OM(4.+3.G)+OMGG(4.BETA0+3.*G))/subd
	1311	c GAM2 =-(1.-GG-OM(4.-3.G)-OMGG(4.BETA0+3.G-4.))/subd
	1312	c GAM3 = BETA0
	1313	c GAM4 = 1. - GAM3
	1314	c mu1(i) = (1. - gg(1.- mu) )/(2. - g*g)
	1315
	1316	*****
	1317	* delta function (Meador, and Weaver, 1980, JAS, 37, 630):
	1318	c GAM1 = (1. - OM*(1. - beta0))/MU
	1319	c GAM2 = OM*BETA0/MU
	1320	c GAM3 = BETA0
	1321	c GAM4 = 1. - GAM3
	1322	c mu1(i) = mu
	1323	*****
	1324	* modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
	1325	c GAM1 = 1.7320508(1. - OM(1. - beta1))
	1326	c GAM2 = 1.7320508OMbeta1
	1327	c GAM3 = BETA0
	1328	c GAM4 = 1. - GAM3
	1329	c mu1(i) = 1./sqrt(3.)
	1330
	1331	* hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
	1332	c GAM1 = 2.(1. - OM(1. - betan))
	1333	c GAM2 = 2.OMBETAn
	1334	c GAM3 = BETA0
	1335	c GAM4 = 1. - GAM3
	1336	c mu1(i) = 0.5
	1337
	1338	*****
	1339
	1340	* lambda = pg 16,290 equation 21
	1341	* big gamma = pg 16,290 equation 22
	1342	* if gam2 = 0., then bgam = 0.
	1343
	1344	lam(i) = sqrt(gam1gam1 - gam2gam2)
	1345
	1346	IF (gam2 .NE. 0.) THEN
	1347	bgam(i) = (gam1 - lam(i))/gam2
	1348	ELSE
	1349	bgam(i) = 0.
	1350	END IF
	1351
	1352	expon = EXP(-lam(i)*taun(i))
	1353
	1354	* e1 - e4 = pg 16,292 equation 44
	1355
	1356	e1(i) = 1. + bgam(i)*expon
	1357	e2(i) = 1. - bgam(i)*expon
	1358	e3(i) = bgam(i) + expon
	1359	e4(i) = bgam(i) - expon
	1360
	1361	* the following sets up for the C equations 23, and 24
	1362	* found on page 16,290
	1363	* prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
	1364	* which is approx equiv to shifting MU by 0.5EPS (MU)**3
	1365
	1366	expon0 = EXP(-tausla(i-1))
	1367	expon1 = EXP(-tausla(i))
	1368
	1369	divisr = lam(i)lam(i) - 1./(mu2(i)mu2(i))
	1370	temp = AMAX1(eps,abs(divisr))
	1371	divisr = SIGN(temp,divisr)
	1372
	1373	up = ompifs((gam1 - 1./mu2(i))gam3 + gam4gam2)/divisr
	1374	dn = ompifs((gam1 + 1./mu2(i))gam4 + gam2gam3)/divisr
	1375
	1376	* cup and cdn are when tau is equal to zero
	1377	* cuptn and cdntn are when tau is equal to taun
	1378
	1379	cup(i) = up*expon0
	1380	cdn(i) = dn*expon0
	1381	cuptn(i) = up*expon1
	1382	cdntn(i) = dn*expon1
	1383
	1384	11 CONTINUE
	1385
	1386	*************** set up matrix ****
	1387	* ssfc = pg 16,292 equation 37 where pi Fs is one (unity).
	1388
	1389	ssfc = rsfcmuEXP(-tausla(nlayer))*pifs
	1390
	1391	* MROWS = the number of rows in the matrix
	1392
	1393	mrows = 2*nlayer
	1394
	1395	* the following are from pg 16,292 equations 39 - 43.
	1396	* set up first row of matrix:
	1397
	1398	i = 1
	1399	a(1) = 0.
	1400	b(1) = e1(i)
	1401	d(1) = -e2(i)
	1402	e(1) = fdn0 - cdn(i)
	1403
	1404	row=1
	1405
	1406	* set up odd rows 3 thru (MROWS - 1):
	1407
	1408	i = 0
	1409	DO 20, row = 3, mrows - 1, 2
	1410	i = i + 1
	1411	a(row) = e2(i)e3(i) - e4(i)e1(i)
	1412	b(row) = e1(i)e1(i + 1) - e3(i)e3(i + 1)
	1413	d(row) = e3(i)e4(i + 1) - e1(i)e2(i + 1)
	1414	e(row) = e3(i)*(cup(i + 1) - cuptn(i)) +
	1415	$ e1(i)*(cdntn(i) - cdn(i + 1))
	1416	20 CONTINUE
	1417
	1418	* set up even rows 2 thru (MROWS - 2):
	1419
	1420	i = 0
	1421	DO 30, row = 2, mrows - 2, 2
	1422	i = i + 1
	1423	a(row) = e2(i + 1)e1(i) - e3(i)e4(i + 1)
	1424	b(row) = e2(i)e2(i + 1) - e4(i)e4(i + 1)
	1425	d(row) = e1(i + 1)e4(i + 1) - e2(i + 1)e3(i + 1)
	1426	e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) -
	1427	$ (cdn(i + 1) - cdntn(i))*e4(i + 1)
	1428	30 CONTINUE
	1429
	1430	* set up last row of matrix at MROWS:
	1431
	1432	row = mrows
	1433	i = nlayer
	1434
	1435	a(row) = e1(i) - rsfc*e3(i)
	1436	b(row) = e2(i) - rsfc*e4(i)
	1437	d(row) = 0.
	1438	e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)
	1439
	1440	* solve tri-diagonal matrix:
	1441
	1442	CALL tridiag(a, b, d, e, y, mrows)
	1443
	1444	**** unfold solution of matrix, compute output fluxes:
	1445
	1446	row = 1
	1447	lev = 1
	1448	j = 1
	1449
	1450	* the following equations are from pg 16,291 equations 31 & 32
	1451
	1452	fdr(lev) = EXP( -tausla(0) )
	1453	edr(lev) = mu * fdr(lev)
	1454	edn(lev) = fdn0
	1455	eup(lev) = y(row)e3(j) - y(row + 1)e4(j) + cup(j)
	1456	fdn(lev) = edn(lev)/mu1(lev)
	1457	fup(lev) = eup(lev)/mu1(lev)
	1458
	1459	DO 60, lev = 2, nlayer + 1
	1460	fdr(lev) = EXP(-tausla(lev-1))
	1461	edr(lev) = mu *fdr(lev)
	1462	edn(lev) = y(row)e3(j) + y(row + 1)e4(j) + cdntn(j)
	1463	eup(lev) = y(row)e1(j) + y(row + 1)e2(j) + cuptn(j)
	1464	fdn(lev) = edn(lev)/mu1(j)
	1465	fup(lev) = eup(lev)/mu1(j)
	1466
	1467	row = row + 2
	1468	j = j + 1
	1469	60 CONTINUE
	1470
[2029]	1471	end subroutine ps2str
	1472
[2027]	1473	=============================================================================
	1474
[2029]	1475	subroutine tridiag(a,b,c,r,u,n)
[2027]	1476
[2029]	1477	!_______________________________________________________________________
	1478	! solves tridiagonal system. From Numerical Recipies, p. 40
	1479	!_______________________________________________________________________
[2027]	1480
	1481	IMPLICIT NONE
	1482
[2029]	1483	! input:
	1484
[2027]	1485	INTEGER n
	1486	REAL a, b, c, r
	1487	DIMENSION a(n),b(n),c(n),r(n)
	1488
[2029]	1489	! output:
	1490
[2027]	1491	REAL u
	1492	DIMENSION u(n)
	1493
[2029]	1494	! local:
	1495
[2027]	1496	INTEGER j
	1497
	1498	REAL bet, gam
	1499	DIMENSION gam(n)
[2029]	1500	!_______________________________________________________________________
[2027]	1501
	1502	IF (b(1) .EQ. 0.) STOP 1001
	1503	bet = b(1)
	1504	u(1) = r(1)/bet
	1505	DO 11, j = 2, n
	1506	gam(j) = c(j - 1)/bet
	1507	bet = b(j) - a(j)*gam(j)
	1508	IF (bet .EQ. 0.) STOP 2002
	1509	u(j) = (r(j) - a(j)*u(j - 1))/bet
	1510	11 CONTINUE
	1511	DO 12, j = n - 1, 1, -1
	1512	u(j) = u(j) - gam(j + 1)*u(j + 1)
	1513	12 CONTINUE
[2029]	1514	!_______________________________________________________________________
[2027]	1515
[2029]	1516	end subroutine tridiag
	1517
	1518	end subroutine photolysis_online

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

Download in other formats:

Original Format