Context Navigation

photolysis_online.F @ 2613

Last change on this file since 2613 was 2553, checked in by yjaziri, 3 years ago

Generic GCM:

Chemistry: one input file instead of three for the chemical network
Automatic detection of mono/bi-molecular or quadratic reactions

File size: 42.4 KB

Rev	Line
[2542]	1	!==============================================================================
	2
	3	subroutine photolysis_online(nlayer, alt, press, temp,
	4	$ zmmean, rm,
	5	$ tau, sza, dist_sol, v_phot, e_phot, ig, ngrid,
[2553]	6	$ nreact)
[2542]	7
	8	!***********************************************************************
	9	!
	10	! subject:
	11	! --------
	12	!
	13	! photolysis online
	14	!
	15	! VERSION: Extracted from LMDZ.MARS work of Franck Lefevre (Yassin Jaziri)
	16	! April 2019 - Yassin Jaziri add updates generic input (Yassin Jaziri)
	17	!
	18	!***********************************************************************
	19
	20	use photolysis_mod
	21	use tracer_h
	22	use chimiedata_h, only: indexchim
[2553]	23	use types_asis, only: nb_phot_hv_max, nb_phot_max, jlabel,
	24	$ reactions
[2542]	25
	26	implicit none
	27
	28	! input
	29
	30	integer, intent(in) :: nlayer ! number of atmospheric layers
	31	integer, intent(in) :: ngrid ! number of atmospheric columns
	32
	33	real, intent(in), dimension(nlayer) :: press, temp, zmmean ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
	34	real, intent(in), dimension(nlayer) :: alt ! altitude (km)
	35	real, intent(in), dimension(nlayer,nesp) :: rm ! mixing ratios
	36	real, intent(in) :: tau ! integrated aerosol optical depth at the surface
	37	real, intent(in) :: sza ! solar zenith angle (degrees)
	38	real, intent(in) :: dist_sol ! solar distance (au)
	39	integer, intent(in) :: ig ! grid point index
	40	integer, intent(in) :: nreact ! number of reactions in reactions files
	41
	42	! output
	43
	44	real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot ! photolysis rates (s-1)
	45	real (kind = 8), dimension(nlayer,nb_phot_max) :: e_phot ! photolysis rates by energie (J.mol-1.s-1)
	46
	47	! solar flux at mars
	48
	49	real, dimension(nw) :: fmars ! solar flux (w.m-2.nm-1)
	50	real :: factor
	51
	52	! atmosphere
	53
	54	real, dimension(nw) :: albedo_chim ! Surface albedo calculated on chemistry wavelenght grid
	55	real, dimension(nlayer) :: colinc ! air column increment (molecule.cm-2)
	56	real, dimension(nlayer) :: airlev ! air density at each specified altitude (molec/cc)
	57	real, dimension(nlayer) :: edir, edn, eup ! normalised irradiances
	58	real, dimension(nlayer) :: fdir, fdn, fup ! normalised actinic fluxes
	59	real, dimension(nlayer) :: saflux ! total actinic flux
	60	real, dimension(nlayer,nw) :: dtrl ! rayleigh optical depth
	61	real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth
	62	real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo
	63	real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter
	64	real, dimension(nlayer,nw) :: dtcld ! cloud optical depth
	65	real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo
	66	real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter
	67	real, dimension(nlayer,nw) :: dagas ! total gas optical depth
	68	real, dimension(nlayer,nw,nabs) :: dtgas ! optical depth for each gas
	69
	70	integer, dimension(0:nlayer) :: nid
	71	real, dimension(0:nlayer,nlayer) :: dsdh
	72
[2553]	73	integer :: i, ilay, iw, ialt, iphot, ispe, ij, igas, ireact
[2542]	74	real :: deltaj
	75
	76	!==== air column increments and rayleigh optical depth
	77
	78	call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl,
	79	$ airlev)
	80
	81	!==== set surface albedo
	82
	83	call setalb(nw,wl,ig,ngrid,albedo_chim)
	84
	85	!==== set temperature-dependent cross-sections and optical depths
	86
[2553]	87	iphot = 0
	88	ij = 0
	89	igas = 0
	90	ireact = 0
[2542]	91
	92	do while(iphot<nb_phot_hv_max)
[2553]	93	ij = ij + 1
	94	iphot = iphot + 1
	95	ireact = ireact + 1
[2542]	96
	97	if (tdim(ij).eq.1) then
	98	! Avoid to calculate several times dtgas for a same specie
	99	if (jlabelbis(iphot)) then
	100	ispe = indexchim(trim(jlabel(iphot,2)))
	101	do iw = 1,nw-1
	102	do ilay = 1,nlayer
	103	dtgas(ilay,iw,ij-igas) = colinc(ilay)*rm(ilay,ispe)
	104	$ *xs(1,iw,ij)
	105	end do
	106	end do
	107	else
	108	igas = igas + 1
	109	end if
[2553]	110	do while(reactions(ireact)%rtype/=0)
	111	ireact = ireact + 1
	112	end do
	113	if (reactions(ireact)%three_prod) iphot = iphot + 1
[2542]	114	else
	115	call setsj(nb_phot_hv_max,nlayer,nw,temp,tdim(ij),
	116	$ xs(:tdim(ij),:,ij),
	117	$ xs_temp(:,ij),sj(:,:,iphot))
	118
	119	! Avoid to calculate several times dtgas for a same specie
	120	if (jlabelbis(iphot)) then
	121	ispe = indexchim(trim(jlabel(iphot,2)))
	122	do iw = 1,nw-1
	123	do ilay = 1,nlayer
	124	dtgas(ilay,iw,ij-igas) = colinc(ilay)*rm(ilay,ispe)
	125	$ *sj(ilay,iw,iphot)
	126	end do
	127	end do
	128	else
	129	igas = igas + 1
	130	end if
	131
	132	do iw = 1,nw-1
	133	do ilay = 1,nlayer
	134	sj(ilay,iw,iphot) = sj(ilay,iw,iphot)*qy(iw,ij)
	135	end do
	136	end do
	137
[2553]	138	do while(reactions(ireact)%rtype/=0)
	139	ireact = ireact + 1
	140	end do
	141	if (reactions(ireact)%three_prod) then
[2542]	142	iphot = iphot + 1
	143	do iw = 1,nw-1
	144	do ilay = 1,nlayer
	145	sj(ilay,iw,iphot) = sj(ilay,iw,iphot-1)
	146	end do
	147	end do
	148	end if
	149
	150	endif ! end if tdim .eq. 1
	151
	152	end do ! end while
	153
	154	! total gas optical depth
	155
	156	dagas(:,:) = 0.
	157
	158	do i = 1,nabs
	159	do iw = 1,nw-1
	160	do ilay = 1,nlayer
	161	dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i)
	162	end do
	163	end do
	164	end do
	165
	166	!==== set aerosol properties and optical depth
	167
	168	call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
	169
	170	!==== set cloud properties and optical depth
	171
	172	call setcld(nlayer,nw,dtcld,omcld,gcld)
	173
	174	!==== slant path lengths in spherical geometry
	175
	176	call sphers(nlayer,alt,sza,dsdh,nid)
	177
	178	!==== solar flux at mars
	179
	180	factor = (1./dist_sol)**2.
	181	do iw = 1,nw-1
	182	fmars(iw) = f(iw)*factor
	183	end do
	184
	185	!==== initialise photolysis rates
	186
	187	v_phot(:,1:nb_phot_hv_max) = 0.
	188	e_phot(:,1:nb_phot_hv_max) = 0.
	189
	190	!==== start of wavelength lopp
	191
	192	do iw = 1,nw-1
	193
	194	! monochromatic radiative transfer. outputs are:
	195	! normalized irradiances edir(nlayer), edn(nlayer), eup(nlayer)
	196	! normalized actinic fluxes fdir(nlayer), fdn(nlayer), fup(nlayer)
	197	! where
	198	! dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse
	199
	200	call rtlink(nlayer, nw, iw, albedo_chim(iw),
	201	$ sza, dsdh, nid, dtrl,
	202	$ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
	203	$ edir, edn, eup, fdir, fdn, fup)
	204
	205
	206	! spherical actinic flux
	207
	208	do ilay = 1,nlayer
	209	saflux(ilay) = fmars(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))
	210	end do
	211
	212	! photolysis rate integration
	213	! (0.12/(wc(iw)*1e-9)) E(wc) en J.mol-1
	214
	215	do i = 1,nb_phot_hv_max
	216	do ilay = 1,nlayer
	217	deltaj = saflux(ilay)*sj(ilay,iw,i)
	218	v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw))
	219	if (wc(iw).le.photoheat_lmax) then
	220	e_phot(ilay,i) = e_phot(ilay,i)
	221	$ + deltaj(wu(iw)-wl(iw))(0.12/(wc(iw)*1e-9))
	222	end if
	223	end do
	224	end do
	225
	226	! eliminate small values
	227
	228	where (v_phot(:,1:nb_phot_hv_max) < 1.e-30)
	229	v_phot(:,1:nb_phot_hv_max) = 0.
	230	end where
	231	where (e_phot(:,1:nb_phot_hv_max) < 1.e-30)
	232	e_phot(:,1:nb_phot_hv_max) = 0.
	233	end where
	234
	235	end do ! iw
	236
	237	contains
	238
	239	!==============================================================================
	240
	241	subroutine setair(nlev, nw, wl, wc, press, temp, zmmean,
	242	$ colinc, dtrl, airlev)
	243
	244	-----------------------------------------------------------------------------
	245	= PURPOSE: =
	246	= computes air column increments and rayleigh optical depth =
	247	-----------------------------------------------------------------------------
	248
	249	use comcstfi_mod, only: g, avocado
	250
	251	implicit none
	252
	253	! input:
	254
	255	integer :: nlev, nw
	256
	257	real, dimension(nw) :: wl, wc ! lower and central wavelength grid (nm)
	258	real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1)
	259
	260	! output:
	261
	262	real, dimension(nlev) :: colinc ! air column increments (molecule.cm-2)
	263	real, dimension(nlev) :: airlev ! air density at each specified altitude (molec/cc)
	264	real, dimension(nlev,nw) :: dtrl ! rayleigh optical depth
	265
	266	! local:
	267
	268	real :: dp, nu
	269	real, dimension(nw) :: srayl
	270	integer :: ilev, iw
	271
	272	! compute column increments
	273
	274	do ilev = 1, nlev-1
	275	dp = (press(ilev) - press(ilev+1))*100.
	276	colinc(ilev) = avocado0.1dp/(zmmean(ilev)*g)
	277	end do
	278	colinc(nlev) = avocado0.1press(nlev)100./(zmmean(nlev)g)
	279
	280	do iw = 1, nw - 1
	281
	282	! co2 rayleigh cross-section
	283	! ityaksov et al., chem. phys. lett., 462, 31-34, 2008
	284
	285	! nu = 1./(wc(iw)*1.e-7)
	286	! srayl(iw) = 1.78e-26nu*(4. + 0.625)
	287	! srayl(iw) = srayl(iw)*1.e-20 ! cm2
	288
	289	! calcul Ityaksov et al., 2008 pour N2
	290
	291	nu = 1./(wc(iw)*1.e-7) ! cm-1
	292	srayl(iw) = 1.8e-26nu*(4. + 0.534)
	293	srayl(iw) = srayl(iw)*1.e-20
	294
	295	do ilev = 1, nlev
	296	dtrl(ilev,iw) = colinc(ilev)*srayl(iw) ! cm2
	297	end do
	298	end do
	299
	300	! compute density
	301
	302	do ilev = 1, nlev
	303	airlev(ilev) = press(ilev)/(1.38e-19*temp(ilev))
	304	end do
	305
	306	end subroutine setair
	307
	308	!==============================================================================
	309
	310	subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
	311
	312	!-----------------------------------------------------------------------------*
	313	!= PURPOSE: =*
	314	!= Set aerosol properties for each specified altitude layer. Properties =*
	315	!= may be wavelength dependent. =*
	316	!-----------------------------------------------------------------------------*
	317
	318	implicit none
	319
	320	! input
	321
	322	integer :: nlayer ! number of vertical layers
	323	integer :: nw ! number of wavelength grid points
	324	real, dimension(nlayer) :: alt ! altitude (km)
	325	real :: tau ! integrated aerosol optical depth at the surface
	326
	327	! output
	328
	329	real, dimension(nlayer,nw) :: dtaer ! aerosol optical depth
	330	real, dimension(nlayer,nw) :: omaer ! aerosol single scattering albedo
	331	real, dimension(nlayer,nw) :: gaer ! aerosol asymmetry parameter
	332
	333	! local
	334
	335	integer :: ilay, iw
	336	real, dimension(nlayer) :: aer ! dust extinction
	337	real :: omega, g, scaleh, gamma
	338	real :: dz, tautot, q0
	339
	340	omega = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm
	341	g = 0.88 ! asymmetry factor : mateshvili et al. (2007) at 210 nm
	342	scaleh = 10. ! scale height (km)
	343	gamma = 0.03 ! conrath parameter
	344
	345	dtaer(:,:) = 0.
	346	omaer(:,:) = 0.
	347	gaer(:,:) = 0.
	348
	349	! optical depth profile:
	350
	351	tautot = 0.
	352	do ilay = 1, nlayer-1
	353	dz = alt(ilay+1) - alt(ilay)
	354	tautot = tautot + exp(gamma(1. - exp(alt(ilay)/scaleh)))dz
	355	end do
	356
	357	q0 = tau/tautot
	358	do ilay = 1, nlayer-1
	359	dz = alt(ilay+1) - alt(ilay)
	360	dtaer(ilay,:) = q0exp(gamma(1. - exp(alt(ilay)/scaleh)))*dz
	361	omaer(ilay,:) = omega
	362	gaer(ilay,:) = g
	363	end do
	364
	365	end subroutine setaer
	366
	367	!==============================================================================
	368
	369	subroutine setcld(nlayer,nw,dtcld,omcld,gcld)
	370
	371	!-----------------------------------------------------------------------------*
	372	!= PURPOSE: =*
	373	!= Set cloud properties for each specified altitude layer. Properties =*
	374	!= may be wavelength dependent. =*
	375	!-----------------------------------------------------------------------------*
	376
	377	implicit none
	378
	379	! input
	380
	381	integer :: nlayer ! number of vertical layers
	382	integer :: nw ! number of wavelength grid points
	383
	384	! output
	385
	386	real, dimension(nlayer,nw) :: dtcld ! cloud optical depth
	387	real, dimension(nlayer,nw) :: omcld ! cloud single scattering albedo
	388	real, dimension(nlayer,nw) :: gcld ! cloud asymmetry parameter
	389
	390	! local
	391
	392	integer :: ilay, iw
	393
	394	! dtcld : optical depth
	395	! omcld : single scattering albedo
	396	! gcld : asymmetry factor
	397
	398	do ilay = 1, nlayer - 1
	399	do iw = 1, nw - 1
	400	dtcld(ilay,iw) = 0. ! no clouds for the moment
	401	omcld(ilay,iw) = 0.99
	402	gcld(ilay,iw) = 0.85
	403	end do
	404	end do
	405
	406	end subroutine setcld
	407
	408	!==============================================================================
	409
	410	subroutine sphers(nlev, z, zen, dsdh, nid)
	411
	412	!-----------------------------------------------------------------------------*
	413	!= PURPOSE: =*
	414	!= Calculate slant path over vertical depth ds/dh in spherical geometry. =*
	415	!= Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model =*
	416	!= for computing the radiation field available for photolysis and heating =*
	417	!= at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) =*
	418	!-----------------------------------------------------------------------------*
	419	!= PARAMETERS: =*
	420	!= NZ - INTEGER, number of specified altitude levels in the working (I)=*
	421	!= grid =*
	422	!= Z - REAL, specified altitude working grid (km) (I)=*
	423	!= ZEN - REAL, solar zenith angle (degrees) (I)=*
	424	!= DSDH - REAL, slant path of direct beam through each layer crossed (O)=*
	425	!= when travelling from the top of the atmosphere to layer i; =*
	426	!= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*
	427	!= NID - INTEGER, number of layers crossed by the direct beam when (O)=*
	428	!= travelling from the top of the atmosphere to layer i; =*
	429	!= NID(i), i = 0..NZ-1 =*
	430	!-----------------------------------------------------------------------------*
	431
	432	use comcstfi_mod, only: rad, pi
	433
	434	implicit none
	435
	436	! input
	437
	438	integer, intent(in) :: nlev
	439	real, dimension(nlev), intent(in) :: z
	440	real, intent(in) :: zen
	441
	442	! output
	443
	444	INTEGER nid(0:nlev)
	445	REAL dsdh(0:nlev,nlev)
	446
	447	! more program constants
	448
	449	REAL re, ze(nlev)
	450	REAL dr
	451	real radius ! km
	452
	453	! local
	454
	455	real :: zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
	456	integer :: i, j, k, id, nlay
	457
	458	REAL zd(0:nlev-1)
	459
	460	!-----------------------------------------------------------------------------
	461
	462	radius = rad*1.e-3 ! rad [m] -> radius [km]
	463	dr = pi/180.
	464	zenrad = zen*dr
	465
	466	! number of layers:
	467
	468	nlay = nlev - 1
	469
	470	! include the elevation above sea level to the radius of Mars:
	471
	472	re = radius + z(1)
	473
	474	! correspondingly z changed to the elevation above Mars surface:
	475
	476	DO k = 1, nlev
	477	ze(k) = z(k) - z(1)
	478	END DO
	479
	480	! inverse coordinate of z
	481
	482	zd(0) = ze(nlev)
	483	DO k = 1, nlay
	484	zd(k) = ze(nlev - k)
	485	END DO
	486
	487	! initialise dsdh(i,j), nid(i)
	488
	489	nid(:) = 0.
	490	dsdh(:,:) = 0.
	491
	492	! calculate ds/dh of every layer
	493
	494	do i = 0,nlay
	495	rpsinz = (re + zd(i))*sin(zenrad)
	496
	497	IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
	498	nid(i) = -1
	499	ELSE
	500
	501	! Find index of layer in which the screening height lies
	502
	503	id = i
	504	if (zen > 90.) then
	505	do j = 1,nlay
	506	IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND.
	507	$ (rpsinz .GE. ( zd(j) + re )) ) id = j
	508	end do
	509	end if
	510
	511	do j = 1,id
	512	sm = 1.0
	513	IF (j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0)
	514	$ sm = -1.0
	515
	516	rj = re + zd(j-1)
	517	rjp1 = re + zd(j)
	518
	519	dhj = zd(j-1) - zd(j)
	520
	521	ga = rjrj - rpsinzrpsinz
	522	gb = rjp1rjp1 - rpsinzrpsinz
	523
	524	ga = max(ga, 0.)
	525	gb = max(gb, 0.)
	526
	527	IF (id.GT.i .AND. j.EQ.id) THEN
	528	dsj = sqrt(ga)
	529	ELSE
	530	dsj = sqrt(ga) - sm*sqrt(gb)
	531	END IF
	532	dsdh(i,j) = dsj/dhj
	533	end do
	534	nid(i) = id
	535	end if
	536	end do ! i = 0,nlay
	537
	538	end subroutine sphers
	539
	540	!==============================================================================
	541
	542	SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl,
	543	$ dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
	544	$ edir, edn, eup, fdir, fdn, fup)
	545
	546	implicit none
	547
	548	! input
	549
	550	integer, intent(in) :: nlev, nw, iw ! number of wavelength grid points
	551	REAL ag
	552	REAL zen
	553	REAL dsdh(0:nlev,nlev)
	554	INTEGER nid(0:nlev)
	555
	556	REAL dtrl(nlev,nw)
	557	REAL dagas(nlev,nw)
	558	REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw)
	559	REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw)
	560
	561	! output
	562
	563	REAL edir(nlev), edn(nlev), eup(nlev)
	564	REAL fdir(nlev), fdn(nlev), fup(nlev)
	565
	566	! local:
	567
	568	REAL dt(nlev), om(nlev), g(nlev)
	569	REAL dtabs,dtsct,dscld,dsaer,dacld,daaer
	570	INTEGER i, ii
	571	real, parameter :: largest = 1.e+36
	572
	573	! specific two ps2str
	574
	575	REAL ediri(nlev), edni(nlev), eupi(nlev)
	576	REAL fdiri(nlev), fdni(nlev), fupi(nlev)
	577
	578	logical, save :: delta = .true.
	579
	580	!_______________________________________________________________________
	581
	582	! initialize:
	583
	584	do i = 1, nlev
	585	fdir(i) = 0.
	586	fup(i) = 0.
	587	fdn(i) = 0.
	588	edir(i) = 0.
	589	eup(i) = 0.
	590	edn(i) = 0.
	591	end do
	592
	593	do i = 1, nlev - 1
	594	dscld = dtcld(i,iw)*omcld(i,iw)
	595	dacld = dtcld(i,iw)*(1.-omcld(i,iw))
	596
	597	dsaer = dtaer(i,iw)*omaer(i,iw)
	598	daaer = dtaer(i,iw)*(1.-omaer(i,iw))
	599
	600	dtsct = dtrl(i,iw) + dscld + dsaer
	601	dtabs = dagas(i,iw) + dacld + daaer
	602
	603	dtabs = amax1(dtabs,1./largest)
	604	dtsct = amax1(dtsct,1./largest)
	605
	606	! invert z-coordinate:
	607
	608	ii = nlev - i
	609	dt(ii) = dtsct + dtabs
	610	om(ii) = dtsct/(dtsct + dtabs)
	611	IF(dtsct .EQ. 1./largest) om(ii) = 1./largest
	612	g(ii) = (gcld(i,iw)*dscld +
	613	$ gaer(i,iw)*dsaer)/dtsct
	614	end do
	615
	616	! call rt routine:
	617
	618	call ps2str(nlev, zen, ag, dt, om, g,
	619	$ dsdh, nid, delta,
	620	$ fdiri, fupi, fdni, ediri, eupi, edni)
	621
	622	! output (invert z-coordinate)
	623
	624	do i = 1, nlev
	625	ii = nlev - i + 1
	626	fdir(i) = fdiri(ii)
	627	fup(i) = fupi(ii)
	628	fdn(i) = fdni(ii)
	629	edir(i) = ediri(ii)
	630	eup(i) = eupi(ii)
	631	edn(i) = edni(ii)
	632	end do
	633
	634	end subroutine rtlink
	635
	636	=============================================================================
	637
	638	subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu,
	639	$ dsdh, nid, delta,
	640	$ fdr, fup, fdn, edr, eup, edn)
	641
	642	!-----------------------------------------------------------------------------*
	643	!= PURPOSE: =*
	644	!= Solve two-stream equations for multiple layers. The subroutine is based =*
	645	!= on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
	646	!= It contains 9 two-stream methods to choose from. A pseudo-spherical =*
	647	!= correction has also been added. =*
	648	!-----------------------------------------------------------------------------*
	649	!= PARAMETERS: =*
	650	!= NLEVEL - INTEGER, number of specified altitude levels in the working (I)=*
	651	!= grid =*
	652	!= ZEN - REAL, solar zenith angle (degrees) (I)=*
	653	!= RSFC - REAL, surface albedo at current wavelength (I)=*
	654	!= TAUU - REAL, unscaled optical depth of each layer (I)=*
	655	!= OMU - REAL, unscaled single scattering albedo of each layer (I)=*
	656	!= GU - REAL, unscaled asymmetry parameter of each layer (I)=*
	657	!= DSDH - REAL, slant path of direct beam through each layer crossed (I)=*
	658	!= when travelling from the top of the atmosphere to layer i; =*
	659	!= DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 =*
	660	!= NID - INTEGER, number of layers crossed by the direct beam when (I)=*
	661	!= travelling from the top of the atmosphere to layer i; =*
	662	!= NID(i), i = 0..NZ-1 =*
	663	!= DELTA - LOGICAL, switch to use delta-scaling (I)=*
	664	!= .TRUE. -> apply delta-scaling =*
	665	!= .FALSE.-> do not apply delta-scaling =*
	666	!= FDR - REAL, contribution of the direct component to the total (O)=*
	667	!= actinic flux at each altitude level =*
	668	!= FUP - REAL, contribution of the diffuse upwelling component to (O)=*
	669	!= the total actinic flux at each altitude level =*
	670	!= FDN - REAL, contribution of the diffuse downwelling component to (O)=*
	671	!= the total actinic flux at each altitude level =*
	672	!= EDR - REAL, contribution of the direct component to the total (O)=*
	673	!= spectral irradiance at each altitude level =*
	674	!= EUP - REAL, contribution of the diffuse upwelling component to (O)=*
	675	!= the total spectral irradiance at each altitude level =*
	676	!= EDN - REAL, contribution of the diffuse downwelling component to (O)=*
	677	= the total spectral irradiance at each altitude level =
	678	!-----------------------------------------------------------------------------*
	679
	680	implicit none
	681
	682	! input:
	683
	684	INTEGER nlev
	685	REAL zen, rsfc
	686	REAL tauu(nlev), omu(nlev), gu(nlev)
	687	REAL dsdh(0:nlev,nlev)
	688	INTEGER nid(0:nlev)
	689	LOGICAL delta
	690
	691	! output:
	692
	693	REAL fup(nlev),fdn(nlev),fdr(nlev)
	694	REAL eup(nlev),edn(nlev),edr(nlev)
	695
	696	! local:
	697
	698	REAL tausla(0:nlev), tauc(0:nlev)
	699	REAL mu2(0:nlev), mu, sum
	700
	701	! internal coefficients and matrix
	702
	703	REAL lam(nlev),taun(nlev),bgam(nlev)
	704	REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev)
	705	REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev)
	706	REAL mu1(nlev)
	707	INTEGER row
	708	REAL a(2nlev),b(2nlev),d(2nlev),e(2nlev),y(2*nlev)
	709
	710	! other:
	711
	712	REAL pifs, fdn0
	713	REAL gi(nlev), omi(nlev), tempg
	714	REAL f, g, om
	715	REAL gam1, gam2, gam3, gam4
	716	real, parameter :: largest = 1.e+36
	717	real, parameter :: precis = 1.e-7
	718
	719	! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
	720	! in delta-function, modified quadrature, hemispheric constant,
	721	! Hybrid modified Eddington-delta function metods, p633,Table1.
	722	! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
	723	! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,
	724	! uncomment the following two lines and the appropriate statements further
	725	! down.
	726	! REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
	727	! > BETA1, BETAn, amu1, subd
	728
	729	REAL expon, expon0, expon1, divisr, temp, up, dn
	730	REAL ssfc
	731	INTEGER nlayer, mrows, lev
	732
	733	INTEGER i, j
	734
	735	! Some additional program constants:
	736
	737	real pi, dr
	738	REAL eps
	739	PARAMETER (eps = 1.E-3)
	740	!_______________________________________________________________________
	741
	742	! MU = cosine of solar zenith angle
	743	! RSFC = surface albedo
	744	! TAUU = unscaled optical depth of each layer
	745	! OMU = unscaled single scattering albedo
	746	! GU = unscaled asymmetry factor
	747	! KLEV = max dimension of number of layers in atmosphere
	748	! NLAYER = number of layers in the atmosphere
	749	! NLEVEL = nlayer + 1 = number of levels
	750
	751	! initial conditions: pi*solar flux = 1; diffuse incidence = 0
	752
	753	pifs = 1.
	754	fdn0 = 0.
	755
	756	nlayer = nlev - 1
	757
	758	pi = acos(-1.)
	759	dr = pi/180.
	760	mu = COS(zen*dr)
	761
	762	!************* compute coefficients for each layer:
	763	! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
	764	! EXPON0 = calculation of e when TAU is zero
	765	! EXPON1 = calculation of e when TAU is TAUN
	766	! CUP and CDN = calculation when TAU is zero
	767	! CUPTN and CDNTN = calc. when TAU is TAUN
	768	! DIVISR = prevents division by zero
	769
	770	do j = 0, nlev
	771	tauc(j) = 0.
	772	tausla(j) = 0.
	773	mu2(j) = 1./SQRT(largest)
	774	end do
	775
	776	IF (.NOT. delta) THEN
	777	DO i = 1, nlayer
	778	gi(i) = gu(i)
	779	omi(i) = omu(i)
	780	taun(i) = tauu(i)
	781	END DO
	782	ELSE
	783
	784	! delta-scaling. Have to be done for delta-Eddington approximation,
	785	! delta discrete ordinate, Practical Improved Flux Method, delta function,
	786	! and Hybrid modified Eddington-delta function methods approximations
	787
	788	DO i = 1, nlayer
	789	f = gu(i)*gu(i)
	790	gi(i) = (gu(i) - f)/(1 - f)
	791	omi(i) = (1 - f)omu(i)/(1 - omu(i)f)
	792	taun(i) = (1 - omu(i)f)tauu(i)
	793	END DO
	794	END IF
	795
	796	! calculate slant optical depth at the top of the atmosphere when zen>90.
	797	! in this case, higher altitude of the top layer is recommended.
	798
	799	IF (zen .GT. 90.0) THEN
	800	IF (nid(0) .LT. 0) THEN
	801	tausla(0) = largest
	802	ELSE
	803	sum = 0.0
	804	DO j = 1, nid(0)
	805	sum = sum + 2.taun(j)dsdh(0,j)
	806	END DO
	807	tausla(0) = sum
	808	END IF
	809	END IF
	810
	811	DO 11, i = 1, nlayer
	812	g = gi(i)
	813	om = omi(i)
	814	tauc(i) = tauc(i-1) + taun(i)
	815
	816	! stay away from 1 by precision. For g, also stay away from -1
	817
	818	tempg = AMIN1(abs(g),1. - precis)
	819	g = SIGN(tempg,g)
	820	om = AMIN1(om,1.-precis)
	821
	822	! calculate slant optical depth
	823
	824	IF (nid(i) .LT. 0) THEN
	825	tausla(i) = largest
	826	ELSE
	827	sum = 0.0
	828	DO j = 1, MIN(nid(i),i)
	829	sum = sum + taun(j)*dsdh(i,j)
	830	END DO
	831	DO j = MIN(nid(i),i)+1,nid(i)
	832	sum = sum + 2.taun(j)dsdh(i,j)
	833	END DO
	834	tausla(i) = sum
	835	IF (tausla(i) .EQ. tausla(i-1)) THEN
	836	mu2(i) = SQRT(largest)
	837	ELSE
	838	mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
	839	mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)),
	840	$ mu2(i) )
	841	END IF
	842	END IF
	843
	844	!** the following gamma equations are from pg 16,289, Table 1
	845	!** save mu1 for each approx. for use in converting irradiance to actinic flux
	846
	847	! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
	848
	849	c gam1 = (7. - om(4. + 3.g))/4.
	850	c gam2 = -(1. - om(4. - 3.g))/4.
	851	c gam3 = (2. - 3.gmu)/4.
	852	c gam4 = 1. - gam3
	853	c mu1(i) = 0.5
	854
	855	* quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):
	856
	857	c gam1 = 1.7320508(2. - om(1. + g))/2.
	858	c gam2 = 1.7320508om(1. - g)/2.
	859	c gam3 = (1. - 1.7320508gmu)/2.
	860	c gam4 = 1. - gam3
	861	c mu1(i) = 1./sqrt(3.)
	862
	863	* hemispheric mean (Toon et al., 1089, JGR, 94, 16287):
	864
	865	gam1 = 2. - om*(1. + g)
	866	gam2 = om*(1. - g)
	867	gam3 = (2. - g*mu)/4.
	868	gam4 = 1. - gam3
	869	mu1(i) = 0.5
	870
	871	* PIFM (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
	872	c GAM1 = 0.25(8. - OM(5. + 3.*G))
	873	c GAM2 = 0.75OM(1.-G)
	874	c GAM3 = 0.25(2.-3.G*MU)
	875	c GAM4 = 1. - GAM3
	876	c mu1(i) = 0.5
	877
	878	* delta discrete ordinates (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
	879	c GAM1 = 0.51.7320508(2. - OM*(1. + G))
	880	c GAM2 = 0.51.7320508OM*(1.-G)
	881	c GAM3 = 0.5(1.-1.7320508G*MU)
	882	c GAM4 = 1. - GAM3
	883	c mu1(i) = 1./sqrt(3.)
	884
	885	* Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
	886	* in delta-function, modified quadrature, hemispheric constant,
	887	* Hybrid modified Eddington-delta function metods, p633,Table1.
	888	* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
	889	* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
	890	c YLM0 = 2.
	891	c YLM2 = -3.GMU
	892	c YLM4 = 0.875G3MU(5.MU**2-3.)
	893	c YLM6=-0.171875G5MU(15.-70.MU*2+63.MU**4)
	894	c YLM8=+0.073242G7MU(-35.+315.MU*2-693.MU**4
	895	c +429.MU**6)
	896	c YLM10=-0.008118G9MU(315.-4620.MU*2+18018.MU**4
	897	c -25740.MU*6+12155.MU**8)
	898	c YLM12=0.003685G11MU(-693.+15015.MU*2-90090.MU**4
	899	c +218790.MU*6-230945.MU*8+88179.MU**10)
	900	c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
	901	c YLMS=0.25*YLMS
	902	c BETA0 = YLMS
	903	c
	904	c amu1=1./1.7320508
	905	c YLM0 = 2.
	906	c YLM2 = -3.Gamu1
	907	c YLM4 = 0.875G3amu1(5.amu1**2-3.)
	908	c YLM6=-0.171875G5amu1(15.-70.amu1*2+63.amu1**4)
	909	c YLM8=+0.073242G7amu1(-35.+315.amu1*2-693.amu1**4
	910	c +429.amu1**6)
	911	c YLM10=-0.008118G9amu1(315.-4620.amu1*2+18018.amu1**4
	912	c -25740.amu1*6+12155.amu1**8)
	913	c YLM12=0.003685G11amu1(-693.+15015.amu1*2-90090.amu1**4
	914	c +218790.amu1*6-230945.amu1*8+88179.amu1**10)
	915	c YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
	916	c YLMS=0.25*YLMS
	917	c BETA1 = YLMS
	918	c
	919	c BETAn = 0.25(2. - 1.5G-0.21875G3-0.085938G**5
	920	c -0.045776G**7)
	921
	922
	923	* Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
	924	c subd=4.(1.-GG*(1.-MU))
	925	c GAM1 = (7.-3.GG-OM(4.+3.G)+OMGG(4.BETA0+3.*G))/subd
	926	c GAM2 =-(1.-GG-OM(4.-3.G)-OMGG(4.BETA0+3.G-4.))/subd
	927	c GAM3 = BETA0
	928	c GAM4 = 1. - GAM3
	929	c mu1(i) = (1. - gg(1.- mu) )/(2. - g*g)
	930
	931	*****
	932	* delta function (Meador, and Weaver, 1980, JAS, 37, 630):
	933	c GAM1 = (1. - OM*(1. - beta0))/MU
	934	c GAM2 = OM*BETA0/MU
	935	c GAM3 = BETA0
	936	c GAM4 = 1. - GAM3
	937	c mu1(i) = mu
	938	*****
	939	* modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
	940	c GAM1 = 1.7320508(1. - OM(1. - beta1))
	941	c GAM2 = 1.7320508OMbeta1
	942	c GAM3 = BETA0
	943	c GAM4 = 1. - GAM3
	944	c mu1(i) = 1./sqrt(3.)
	945
	946	* hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
	947	c GAM1 = 2.(1. - OM(1. - betan))
	948	c GAM2 = 2.OMBETAn
	949	c GAM3 = BETA0
	950	c GAM4 = 1. - GAM3
	951	c mu1(i) = 0.5
	952
	953	*****
	954
	955	* lambda = pg 16,290 equation 21
	956	* big gamma = pg 16,290 equation 22
	957	* if gam2 = 0., then bgam = 0.
	958
	959	lam(i) = sqrt(gam1gam1 - gam2gam2)
	960
	961	IF (gam2 .NE. 0.) THEN
	962	bgam(i) = (gam1 - lam(i))/gam2
	963	ELSE
	964	bgam(i) = 0.
	965	END IF
	966
	967	expon = EXP(-lam(i)*taun(i))
	968
	969	* e1 - e4 = pg 16,292 equation 44
	970
	971	e1(i) = 1. + bgam(i)*expon
	972	e2(i) = 1. - bgam(i)*expon
	973	e3(i) = bgam(i) + expon
	974	e4(i) = bgam(i) - expon
	975
	976	* the following sets up for the C equations 23, and 24
	977	* found on page 16,290
	978	* prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
	979	* which is approx equiv to shifting MU by 0.5EPS (MU)**3
	980
	981	expon0 = EXP(-tausla(i-1))
	982	expon1 = EXP(-tausla(i))
	983
	984	divisr = lam(i)lam(i) - 1./(mu2(i)mu2(i))
	985	temp = AMAX1(eps,abs(divisr))
	986	divisr = SIGN(temp,divisr)
	987
	988	up = ompifs((gam1 - 1./mu2(i))gam3 + gam4gam2)/divisr
	989	dn = ompifs((gam1 + 1./mu2(i))gam4 + gam2gam3)/divisr
	990
	991	* cup and cdn are when tau is equal to zero
	992	* cuptn and cdntn are when tau is equal to taun
	993
	994	cup(i) = up*expon0
	995	cdn(i) = dn*expon0
	996	cuptn(i) = up*expon1
	997	cdntn(i) = dn*expon1
	998
	999	11 CONTINUE
	1000
	1001	*************** set up matrix ****
	1002	* ssfc = pg 16,292 equation 37 where pi Fs is one (unity).
	1003
	1004	ssfc = rsfcmuEXP(-tausla(nlayer))*pifs
	1005
	1006	* MROWS = the number of rows in the matrix
	1007
	1008	mrows = 2*nlayer
	1009
	1010	* the following are from pg 16,292 equations 39 - 43.
	1011	* set up first row of matrix:
	1012
	1013	i = 1
	1014	a(1) = 0.
	1015	b(1) = e1(i)
	1016	d(1) = -e2(i)
	1017	e(1) = fdn0 - cdn(i)
	1018
	1019	row=1
	1020
	1021	* set up odd rows 3 thru (MROWS - 1):
	1022
	1023	i = 0
	1024	DO 20, row = 3, mrows - 1, 2
	1025	i = i + 1
	1026	a(row) = e2(i)e3(i) - e4(i)e1(i)
	1027	b(row) = e1(i)e1(i + 1) - e3(i)e3(i + 1)
	1028	d(row) = e3(i)e4(i + 1) - e1(i)e2(i + 1)
	1029	e(row) = e3(i)*(cup(i + 1) - cuptn(i)) +
	1030	$ e1(i)*(cdntn(i) - cdn(i + 1))
	1031	20 CONTINUE
	1032
	1033	* set up even rows 2 thru (MROWS - 2):
	1034
	1035	i = 0
	1036	DO 30, row = 2, mrows - 2, 2
	1037	i = i + 1
	1038	a(row) = e2(i + 1)e1(i) - e3(i)e4(i + 1)
	1039	b(row) = e2(i)e2(i + 1) - e4(i)e4(i + 1)
	1040	d(row) = e1(i + 1)e4(i + 1) - e2(i + 1)e3(i + 1)
	1041	e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) -
	1042	$ (cdn(i + 1) - cdntn(i))*e4(i + 1)
	1043	30 CONTINUE
	1044
	1045	* set up last row of matrix at MROWS:
	1046
	1047	row = mrows
	1048	i = nlayer
	1049
	1050	a(row) = e1(i) - rsfc*e3(i)
	1051	b(row) = e2(i) - rsfc*e4(i)
	1052	d(row) = 0.
	1053	e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)
	1054
	1055	* solve tri-diagonal matrix:
	1056
	1057	CALL tridiag(a, b, d, e, y, mrows)
	1058
	1059	**** unfold solution of matrix, compute output fluxes:
	1060
	1061	row = 1
	1062	lev = 1
	1063	j = 1
	1064
	1065	* the following equations are from pg 16,291 equations 31 & 32
	1066
	1067	fdr(lev) = EXP( -tausla(0) )
	1068	edr(lev) = mu * fdr(lev)
	1069	edn(lev) = fdn0
	1070	eup(lev) = y(row)e3(j) - y(row + 1)e4(j) + cup(j)
	1071	fdn(lev) = edn(lev)/mu1(lev)
	1072	fup(lev) = eup(lev)/mu1(lev)
	1073
	1074	DO 60, lev = 2, nlayer + 1
	1075	fdr(lev) = EXP(-tausla(lev-1))
	1076	edr(lev) = mu *fdr(lev)
	1077	edn(lev) = y(row)e3(j) + y(row + 1)e4(j) + cdntn(j)
	1078	eup(lev) = y(row)e1(j) + y(row + 1)e2(j) + cuptn(j)
	1079	fdn(lev) = edn(lev)/mu1(j)
	1080	fup(lev) = eup(lev)/mu1(j)
	1081
	1082	row = row + 2
	1083	j = j + 1
	1084	60 CONTINUE
	1085
	1086	end subroutine ps2str
	1087
	1088	=============================================================================
	1089
	1090	subroutine tridiag(a,b,c,r,u,n)
	1091
	1092	!_______________________________________________________________________
	1093	! solves tridiagonal system. From Numerical Recipies, p. 40
	1094	!_______________________________________________________________________
	1095
	1096	IMPLICIT NONE
	1097
	1098	! input:
	1099
	1100	INTEGER n
	1101	REAL a, b, c, r
	1102	DIMENSION a(n),b(n),c(n),r(n)
	1103
	1104	! output:
	1105
	1106	REAL u
	1107	DIMENSION u(n)
	1108
	1109	! local:
	1110
	1111	INTEGER j
	1112
	1113	REAL bet, gam
	1114	DIMENSION gam(n)
	1115	!_______________________________________________________________________
	1116
	1117	IF (b(1) .EQ. 0.) STOP 1001
	1118	bet = b(1)
	1119	u(1) = r(1)/bet
	1120	DO 11, j = 2, n
	1121	gam(j) = c(j - 1)/bet
	1122	bet = b(j) - a(j)*gam(j)
	1123	IF (bet .EQ. 0.) STOP 2002
	1124	u(j) = (r(j) - a(j)*u(j - 1))/bet
	1125	11 CONTINUE
	1126	DO 12, j = n - 1, 1, -1
	1127	u(j) = u(j) - gam(j + 1)*u(j + 1)
	1128	12 CONTINUE
	1129	!_______________________________________________________________________
	1130
	1131	end subroutine tridiag
	1132
	1133	!==============================================================================
	1134
	1135	subroutine setalb(nw,wl,ig,ngrid,albedo_chim)
	1136
	1137	!-----------------------------------------------------------------------------*
	1138	!= PURPOSE: =*
	1139	!= Set the albedo of the surface. The albedo is assumed to be Lambertian, =*
	1140	!= i.e., the reflected light is isotropic, and idependt of the direction =*
	1141	!= of incidence of light. Albedo can be chosen to be wavelength dependent. =*
	1142	!-----------------------------------------------------------------------------*
	1143	!= PARAMETERS: =*
	1144	!= NW - INTEGER, number of specified intervals + 1 in working (I)=*
	1145	!= wavelength grid =*
	1146	!= WL - REAL, vector of lower limits of wavelength intervals in (I)=*
	1147	!= working wavelength grid =*
	1148	!= ALBEDO - REAL, surface albedo at each specified wavelength (O)=*
	1149	!-----------------------------------------------------------------------------*
	1150
	1151	use chimiedata_h, only: albedo_snow_chim, albedo_co2_ice_chim
	1152	use slab_ice_h, only: h_alb_ice, alb_ice_min, alb_ice_max
	1153	use tracer_h, only: igcm_h2o_ice, igcm_co2_ice
	1154	use callkeys_mod, only: ok_slab_ocean, co2cond, alb_ocean
	1155	use phys_state_var_mod, only: albedo_bareground,
	1156	& rnat, qsurf, sea_ice,
	1157	& pctsrf_sic, tsurf, capcal
	1158	use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater
	1159
	1160	implicit none
	1161
	1162	! input: (wavelength working grid data)
	1163
	1164	integer, intent(in) :: ngrid ! number of atmospheric columns
	1165	INTEGER nw
	1166	INTEGER ig ! grid point index
	1167	REAL wl(nw)
	1168
	1169	! output:
	1170
	1171	REAL albedo_chim(nw)
	1172
	1173	! local:
	1174
	1175	INTEGER iw
	1176	! REAL alb
	1177	real zfra, alb_ice, twater, hice
	1178	real snowlayer
	1179	parameter (snowlayer=33.0) ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm
	1180
	1181	! 0.015: mean value from clancy et al., icarus, 49-63, 1999.
	1182
	1183	! alb = 0.015
	1184	! alb = albedo_phys
	1185
	1186	! do iw = 1, nw - 1
	1187	! albedo_chim(iw) = alb
	1188	! end do
	1189
	1190	! See hydrol.F90 for taking into account new tendencies
	1191
	1192	if(nint(rnat(ig)).eq.0)then
	1193
	1194	if(ok_slab_ocean) then
	1195
	1196	zfra = MAX(0.0,MIN(1.0,qsurf(ig,igcm_h2o_ice)/45.0)) ! Snow Fraction (Critical height 45kg/m2~15cm)
	1197	alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) ! Ice Albedo
	1198	& *exp(-sea_ice(ig)/h_alb_ice)
	1199	! Albedo final calculation :
	1200	do iw=1,nw - 1
	1201	albedo_chim(iw) = pctsrf_sic(ig)*
	1202	& (albedo_snow_chim(iw)*zfra
	1203	& + alb_ice*(1.0-zfra))
	1204	& + (1.-pctsrf_sic(ig))*alb_ocean
	1205	enddo
	1206
	1207	else !ok_slab_ocean
	1208
	1209
	1210	! calculate oceanic ice height including the latent heat of ice formation
	1211	! hice is the height of oceanic ice with a maximum of maxicethick.
	1212	hice = qsurf(ig,igcm_h2o_ice)/rhowater ! update hice to include recent snowfall
	1213	twater = tsurf(ig) - hiceRLFTTrhowater/capcal(ig)
	1214	! this is temperature water would have if we melted entire ocean ice layer
	1215
	1216	if(twater .lt. T_h2O_ice_liq)then
	1217
	1218	do iw=1,nw - 1
	1219	albedo_chim(iw) = albedo_snow_chim(iw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
	1220	enddo
	1221
	1222	else
	1223
	1224	DO iw=1,nw - 1
	1225	albedo_chim(iw) = alb_ocean
	1226	if(ngrid.eq.1) then
	1227	albedo_chim(iw) = albedo_bareground(ig)
	1228	endif
	1229	ENDDO
	1230
	1231	endif
	1232
	1233	endif!(ok_slab_ocean)
	1234
	1235
	1236	! Continent
	1237	! ---------
	1238	elseif (nint(rnat(ig)).eq.1) then
	1239
	1240	! re-calculate continental albedo
	1241	if (qsurf(ig,igcm_h2o_ice).ge.snowlayer) then
	1242	DO iw=1,nw - 1
	1243	albedo_chim(iw) = albedo_snow_chim(iw)
	1244	ENDDO
	1245	else
	1246	DO iw=1,nw - 1
	1247	albedo_chim(iw) = albedo_bareground(ig)
	1248	& + (albedo_snow_chim(iw)
	1249	& - albedo_bareground(ig))
	1250	& *qsurf(ig,igcm_h2o_ice)/snowlayer
	1251	ENDDO
	1252	endif
	1253
	1254	else
	1255
	1256	print*,'Surface type not recognised in photolysis_online.F!'
	1257	print*,'Exiting...'
	1258	call abort
	1259
	1260	endif
	1261
	1262	! Re-add the albedo effects of CO2 ice if necessary
	1263	! -------------------------------------------------
	1264	if(co2cond)then
	1265
	1266	if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
	1267	DO iw=1,nw - 1
	1268	albedo_chim(iw) = albedo_co2_ice_chim(iw)
	1269	ENDDO
	1270	endif
	1271
	1272	endif ! co2cond
	1273
	1274	end subroutine setalb
	1275
	1276	!==============================================================================
	1277
	1278	subroutine setsj(nd,nlayer,nw,tlay,tdiml,xsl,xs_templ,sjl)
	1279
	1280
	1281	implicit none
	1282
	1283	! input: (wavelength working grid data)
	1284
	1285	integer :: nd ! number of photolysis rates
	1286	integer :: nlayer ! number of vertical layers
	1287	integer :: nw ! number of wavelength grid points
	1288	integer :: tdiml ! number of different temperature cross section file
	1289	real :: xsl(tdiml,nw) ! cross section (cm2) from reading files
	1290	real :: xs_templ(tdiml) ! temperature of the cross section (cm2)
	1291	real :: tlay(nlayer) ! temperature (K)
	1292
	1293	! output:
	1294
	1295	real :: sjl(nlayer,nw) ! output cross section (cm2)
	1296
	1297
	1298	! local:
	1299
	1300	INTEGER ilay,iw,tpos
	1301
	1302	do iw = 1,nw-1
	1303	do ilay = 1,nlayer
	1304	tpos = 1
	1305	do while(tlay(ilay)>xs_templ(tpos) .and. tpos<tdiml)
	1306	tpos = tpos + 1
	1307	end do
	1308	if (tpos.eq.1 .or. tpos.eq.tdiml) then
	1309	sjl(ilay,iw) = xsl(tpos,iw)
	1310	else
	1311	sjl(ilay,iw) = ( (xsl(tpos,iw)-xsl(tpos-1,iw))*
	1312	$ (tlay(ilay)-xs_templ(tpos-1))
	1313	$ /(xs_templ(tpos)-xs_templ(tpos-1))
	1314	$ + xsl(tpos-1,iw) )
	1315	end if
	1316	end do
	1317	end do
	1318
	1319	end subroutine setsj
	1320
	1321	end subroutine photolysis_online

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

Download in other formats:

Original Format