Context Navigation

photolysis_online.F @ 2883

Last change on this file since 2883 was 2615, checked in by romain.vande, 3 years ago
LMDZ_MARS RV : Open_MP; Put all the "save" variables as "!$OMP THREADPRIVATE" in aeronomars
File size: 51.6 KB

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

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

Download in other formats:

Original Format