Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_sechiba/albedo.f90 @ 5242

Last change on this file since 5242 was 2244, checked in by matthew.mcgrath, 10 years ago
DEV: Printing out the noon albedo in the coupled runs
Property svn:executable set to ``*
File size: 105.9 KB

Line
1	! ===============================================================================================================================
2	! MODULE : albedo
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Calculate the surface albedo using a variety of different schemes
10	!!
11	!! \n DESCRIPTION : This module includes the mechanisms for
12	!! 1. :: albedo_calc for computing the albedo with the old method (need to ask
13	!! someone what exactly this is from)...the simplest version
14	!! this does not distinguish between diffuse and direct, but can
15	!! distinguish between NIR and VIS
16	!! 2. :: albedo_two_stream computes the two_stream approach of Pinty et al 2006.
17	!! Requires an effective leaf area index and depends on the solar angle.
18	!! Separates light into NIR and VIS, and diffuse and direct illumination
19	!! 3. :: calculate_surface_albedo_with_snow is a way to calculate the snow albedo
20	!! based on the method of Chalita and Treut (1994)
21	!! This does not distinguish between NIR and VIR light, nor diffuse
22	!! and direct, and should probably only be used with albedo_calc
23	!!
24	!! RECENT CHANGE(S): None
25	!!
26	!! REFERENCES(S) :
27	!! Chalita, S. and H Le Treut (1994), The albedo of temperate
28	!! and boreal forest and the Northern Hemisphere climate: a
29	!! sensitivity experiment using the LMD GCM,
30	!! Climate Dynamics, 10 231-240.
31	!!
32	!! B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski,
33	!! N. Gobron and M. M. Verstraete (2006). Simplifying the
34	!! Interaction of Land Surfaces with Radiation for Relating
35	!! Remote Sensing Products to Climate Models. Journal of
36	!! Geophysical Research. Vol 111, D02116.
37	!!
38	!! SVN :
39	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/perso/matthew.mcgrath/TRUNK/ORCHIDEE/src_sechiba/albedo.f90 $
40	!! $Date: 2012-07-19 15:12:52 +0200 (Thu, 1 Nov 2012) $
41	!! $Revision: 947 $
42	!! \n
43	!_ ================================================================================================================================
44
45	MODULE albedo
46	!
47	USE ioipsl
48	!
49	! modules used :
50	USE constantes
51	USE constantes_soil
52	USE pft_parameters
53	USE structures
54	USE interpol_help
55	USE ioipsl_para
56	USE stomate_laieff
57
58	IMPLICIT NONE
59
60	PRIVATE :: print_debugging_albedo_info,calculate_snow_albedo,&
61	twostream_solver,gammas
62	! public routines :
63	PUBLIC :: albedo_calc, albedo_two_stream, calculate_surface_albedo_with_snow
64
65	!
66	! variables used inside albedo module : declaration and initialisation
67	!
68
69
70	CONTAINS
71
72	!! ==============================================================================================================================
73	!! SUBROUTINE : albedo_calc
74	!!
75	!>\BRIEF This subroutine calculates the albedo without snow.
76	!!
77	!! DESCRIPTION : The albedo is calculated for both the visible and near-infrared
78	!! domain. First the mean albedo of the bare soil is calculated. Two options exist:
79	!! either the soil albedo depends on soil wetness (drysoil_frac variable), or the
80	!! soil albedo is set to a mean soil albedo value. It should be noted that this
81	!! routine essential separates vegetated land into 1) 100Ã¹ plant covered and 2)
82	!! bare soil parts, and calculates the albedo for them seperately.
83	!!
84	!! RECENT CHANGE(S): None
85	!!
86	!! MAIN OUTPUT VARIABLE(S): albedo
87	!!
88	!! REFERENCE(S) : Krinner, G., Viovy, N., Noblet-Ducoudre, N. de, Ogee, J., Polcher,
89	!! J., Friedlingstein, P., Ciais, P., Sitch, S., and Prentice, I. C
90	!! (2005). A dynamic global vegetation model for studies of the
91	!! coupled atmosphere-biosphere system. Global Biogeochem. Cycles, 19, GB1015.
92	!!
93	!! FLOWCHART : None
94	!! \n
95	!_ ================================================================================================================================
96
97	SUBROUTINE albedo_calc (npts, drysoil_frac, veget, soilalb_dry, &
98	soilalb_wet, soilalb_moy, tot_bare_soil, alb_veget, &
99	alb_bare, albedo)
100
101	!! 0. Variable and parameter declaration
102
103	!! 0.1 Input variables
104
105	INTEGER(i_std), INTENT(in) :: npts !! Domain size - Number of land pixels (unitless)
106	REAL(r_std), DIMENSION(npts), INTENT(in) :: drysoil_frac !! Fraction of dry soil (unitless)
107	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
108	!! (m^2 m^{-2})
109	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: soilalb_dry !! Albedo values for the dry bare soil (unitless)
110	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: soilalb_wet !! Albedo values for the wet bare soil (unitless)
111	REAL(r_std), DIMENSION(npts,2), INTENT(in) :: soilalb_moy !! Albedo values for the mean bare soil (unitless)
112	REAL(r_std), DIMENSION (npts), INTENT(in) :: tot_bare_soil !! total evaporating bare soil fraction
113
114	!! 0.2 Output variables
115
116	REAL(r_std), DIMENSION(npts,2), INTENT(out) :: alb_veget !! Mean vegetation albedo for visible and
117	!! near-infrared range (unitless)
118	REAL(r_std), DIMENSION(npts,2), INTENT(out) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
119	REAL(r_std),DIMENSION (npts,2), INTENT (out) :: albedo !! Albedo for visible and near-infrared range
120	!! (unitless)
121
122	!! 0.3 Modified variables
123
124	!! 0.4 Local variables
125
126	REAL(r_std),DIMENSION (nvm,2) :: alb_leaf_tmp !! Variable for albedo values for all PFTs and
127	!! spectral domains (unitless)
128	INTEGER(i_std) :: ks !! Index for visible and near-infraread range
129	INTEGER(i_std) :: jv !! Index for vegetative PFTs
130	!_ ================================================================================================================================
131
132	IF (bavard.GE.2) WRITE(numout,*) 'Entering albedo_calc'
133
134	!! 1. Preliminary calculation
135
136	! Assign values of leaf albedo for visible and near-infrared range
137	! to local variable (constantes_veg.f90)
138	alb_leaf_tmp(:,ivis) = alb_leaf_vis(:)
139	alb_leaf_tmp(:,inir) = alb_leaf_nir(:)
140
141	!! 2. Calculation and assignment of soil albedo
142
143	DO ks = 1, 2! Loop over # of spectra
144
145	! If 'alb_bare_model' is set to TRUE, the soil albedo calculation depends
146	! on soil moisture
147	! If 'alb_bare_model' is set to FALSE, the mean soil albedo is used without
148	! the dependance on soil moisture
149	! see subroutine 'conveg_soilalb'
150	IF ( alb_bare_model ) THEN
151	alb_bare(:,ks) = soilalb_wet(:,ks) + drysoil_frac(:) * &
152	(soilalb_dry(:,ks) - soilalb_wet(:,ks))
153	ELSE
154	alb_bare(:,ks) = soilalb_moy(:,ks)
155	ENDIF
156
157	! Soil albedo is weighed by fraction of bare soil
158	albedo(:,ks) = tot_bare_soil(:) * alb_bare(:,ks)
159
160	!! 3. Calculation of mean albedo of over the grid cell
161
162	! Calculation of mean albedo of over the grid cell and
163	! mean albedo of only vegetative PFTs over the grid cell
164	alb_veget(:,ks) = zero
165
166	DO jv = 2, nvm ! Loop over # of PFTs
167
168	! Mean albedo of grid cell for visible and near-infrared range
169	albedo(:,ks) = albedo(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
170
171	! Mean albedo of vegetation for visible and near-infrared range
172	alb_veget(:,ks) = alb_veget(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
173	ENDDO ! Loop over # of PFTs
174
175	ENDDO ! ks = 1, 2
176
177	IF (bavard.GE.2) WRITE(numout,*) 'Leaving albedo_calc'
178
179	END SUBROUTINE albedo_calc
180
181	!! ==============================================================================================================================\n
182	!! SUBROUTINE : calculate_surface_albedo_with_snow
183	!!
184	!>\BRIEF Calcuating snow albedo and snow cover fraction, which are then used to
185	!! update the gridbox surface albedo following Chalita and Treut (1994). This is the old
186	!! condveg_snow routine, and is only used for calculated the snow albedo with the old
187	!! albedo scheme.
188	!!
189	!! DESCRIPTION : The snow albedo scheme presented below belongs to prognostic albedo
190	!! category, i.e. the snow albedo value at a time step depends on the snow albedo value
191	!! at the previous time step.
192	!!
193	!! First, the following formula (described in Chalita and Treut 1994) is used to describe
194	!! the change in snow albedo with snow age on each PFT and each non-vegetative surfaces,
195	!! i.e. continental ice, lakes, etc.: \n
196	!! \latexonly
197	!! \input{SnowAlbedo.tex}
198	!! \endlatexonly
199	!! \n
200	!! Where snowAge is snow age, tcstSnowa is a critical aging time (tcstSnowa=5 days)
201	!! snowaIni and snowaIni+snowaDec corresponds to albedos measured for aged and
202	!! fresh snow respectively, and their values for each PFT and each non-vegetative surfaces
203	!! is precribed in in constantes_veg.f90.\n
204	!! In order to estimate gridbox snow albedo, snow albedo values for each PFT and
205	!! each non-vegetative surfaces with a grid box are weightedly summed up by their
206	!! respective fractions.\n
207	!! Secondly, the snow cover fraction is computed as:
208	!! \latexonly
209	!! \input{SnowFraction.tex}
210	!! \endlatexonly
211	!! \n
212	!! Where fracSnow is the fraction of snow on total vegetative or total non-vegetative
213	!! surfaces, snow is snow mass (kg/m^2) on total vegetated or total nobio surfaces.\n
214	!! Finally, the surface albedo is then updated as the weighted sum of fracSnow, total
215	!! vegetated fraction, total nobio fraction, gridbox snow albedo, and previous
216	!! time step surface albedo.
217	!!
218	!! RECENT CHANGE(S): None
219	!!
220	!! MAIN OUTPUT VARIABLE(S): :: albedo; surface albedo. :: albedo_snow; snow
221	!! albedo
222	!!
223	!! REFERENCE(S) :
224	!! Chalita, S. and H Le Treut (1994), The albedo of temperate and boreal forest and
225	!! the Northern Hemisphere climate: a sensitivity experiment using the LMD GCM,
226	!! Climate Dynamics, 10 231-240.
227	!!
228	!! FLOWCHART : None
229	!! \n
230	!_ ================================================================================================================================
231
232	SUBROUTINE calculate_surface_albedo_with_snow (npts, veget, veget_max, &
233	frac_nobio, totfrac_nobio, snow, snow_age, &
234	snow_nobio, snow_nobio_age, tot_bare_soil, albedo, &
235	albedo_snow, albedo_glob)
236
237	!! 0. Variable and parameter declaration
238
239	!! 0.1 Input variables
240
241	INTEGER(i_std), INTENT(in) :: npts !! Domain size - Number of land pixels (unitless)
242	REAL(r_std),DIMENSION (npts,nvm), INTENT (in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
243	!! (m^2 m^{-2})
244	REAL(r_std),DIMENSION (npts,nvm), INTENT(in) :: veget_max
245	REAL(r_std),DIMENSION (npts,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
246	!! continental ice, lakes, etc. (unitless)
247	REAL(r_std),DIMENSION (npts), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces, i.e.
248	!! continental ice, lakes, etc. (unitless)
249	REAL(r_std),DIMENSION (npts), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
250	REAL(r_std),DIMENSION (npts,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
251	REAL(r_std),DIMENSION (npts), INTENT(in) :: snow_age !! Snow age (days)
252	REAL(r_std),DIMENSION (npts,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
253	REAL(r_std), DIMENSION (npts), INTENT(in) :: tot_bare_soil !! total evaporating bare soil fraction
254
255	!! 0.2 Output variables
256
257	REAL(r_std),DIMENSION (npts,2), INTENT (inout) :: albedo !! Albedo (unitless ratio)
258	REAL(r_std),DIMENSION (npts), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
259	REAL(r_std),DIMENSION (npts), INTENT (out) :: albedo_glob !! Mean albedo (unitless ratio)
260
261	!! 0.3 Modified variables
262
263	!! 0.4 Local variables
264
265	INTEGER(i_std) :: ji, jv, jb !! indices (unitless)
266	REAL(r_std), DIMENSION(npts) :: frac_snow_veg !! Fraction of snow on vegetation (unitless ratio)
267	REAL(r_std), DIMENSION(npts,nnobio) :: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
268	!! (unitless ratio)
269	REAL(r_std), DIMENSION(npts) :: snowa_veg !! Albedo of snow covered area on vegetation
270	!! (unitless ratio)
271	REAL(r_std), DIMENSION(npts,nnobio) :: snowa_nobio !! Albedo of snow covered area on continental ice,
272	!! lakes, etc. (unitless ratio)
273	REAL(r_std), DIMENSION(npts) :: fraction_veg !! Total vegetation fraction (unitless ratio)
274	REAL(r_std), DIMENSION(npts) :: agefunc_veg !! Age dependency of snow albedo on vegetation
275	!! (unitless)
276	REAL(r_std), DIMENSION(npts,nnobio) :: agefunc_nobio !! Age dependency of snow albedo on ice,
277	!! lakes, .. (unitless)
278	REAL(r_std) :: alb_nobio !! Albedo of continental ice, lakes, etc.
279	!!(unitless ratio)
280	!_ ================================================================================================================================
281
282	IF (bavard.GE.2) WRITE(numout,*) 'Entering calculate_surface_albedo_with_snow '
283
284	!! 1.Reset output variable (::albedo_snow ::albedo_glob)
285
286	DO ji = 1, npts
287	albedo_snow(ji) = zero
288	albedo_glob(ji) = zero
289	ENDDO
290
291	!! 2. Calculate snow albedos on both total vegetated and total nobio surfaces
292
293	! The snow albedo could be either prescribed (in condveg_init.f90) or
294	! calculated following Chalita and Treut (1994).
295	! Check if the precribed value fixed_snow_albedo exists
296	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
297	snowa_veg(:) = fixed_snow_albedo
298	snowa_nobio(:,:) = fixed_snow_albedo
299	ELSE ! calculated following Chalita and Treut (1994)
300
301	!! 2.1 Calculate age dependence
302
303	! On vegetated surfaces
304	DO ji = 1, npts
305	agefunc_veg(ji) = EXP(-snow_age(ji)/tcst_snowa)
306	ENDDO
307
308	! On non-vegtative surfaces
309	DO jv = 1, nnobio ! Loop over # nobio types
310	DO ji = 1, npts
311	agefunc_nobio(ji,jv) = EXP(-snow_nobio_age(ji,jv)/tcst_snowa)
312	ENDDO
313	ENDDO
314
315	!! 2.1 Calculate snow albedo
316
317	! For vegetated surfaces
318	fraction_veg(:) = un - totfrac_nobio(:)
319	snowa_veg(:) = zero
320	IF (control%ok_dgvm) THEN
321	DO ji = 1, npts
322	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
323	snowa_veg(ji) = snowa_veg(ji) + &
324	tot_bare_soil(ji)/fraction_veg(ji) * &
325	( snowa_aged(1)+snowa_dec(1)*agefunc_veg(ji) )
326	ENDIF
327	ENDDO
328
329	DO jv = 2, nvm
330	DO ji = 1, npts
331	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
332	snowa_veg(ji) = snowa_veg(ji) + &
333	veget(ji,jv)/fraction_veg(ji) * &
334	( snowa_aged(jv)+snowa_dec(jv)*agefunc_veg(ji) )
335	ENDIF
336	ENDDO
337	ENDDO
338	ELSE
339	DO jv = 1, nvm
340	DO ji = 1, npts
341	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
342	snowa_veg(ji) = snowa_veg(ji) + &
343	veget_max(ji,jv)/fraction_veg(ji) * &
344	( snowa_aged(jv)+snowa_dec(jv)*agefunc_veg(ji) )
345	ENDIF
346	ENDDO
347	ENDDO
348	ENDIF ! control%ok_dgvm
349	!
350	! snow albedo on other surfaces
351	!
352	DO jv = 1, nnobio
353	DO ji = 1, npts
354	snowa_nobio(ji,jv) = &
355	( snowa_aged(1) + snowa_dec(1) * agefunc_nobio(ji,jv) )
356	ENDDO
357	ENDDO
358	ENDIF ! prescribe or calculate the snow albedo?
359
360	!! 3. Calculate snow cover fraction for both total vegetated and total
361	!! non-vegtative surfaces.\n
362
363	! This has been changed from the trunk version to match what is actually
364	! in the reference...Jan 1, 2013
365	frac_snow_veg(:) = MIN(MAX(snow(:),zero)/&
366	(MAX(snow(:),zero)+snowcri_alb*sn_dens/100.0_r_std),un)
367	DO jv = 1, nnobio
368	frac_snow_nobio(:,jv) = MIN(MAX(snow_nobio(:,jv),zero)/&
369	(MAX(snow_nobio(:,jv),zero)+snowcri_alb*sn_dens/100.0_r_std),un)
370	ENDDO
371
372	!! 4. Update surface albedo
373
374	! Update surface albedo using the weighted sum of previous time step surface
375	! albedo, total vegetated fraction, total nobio fraction, snow cover
376	! fraction (both vegetated and non-vegetative surfaces), and snow albedo
377	! (both vegetated and non-vegetative surfaces). Although both visible and
378	! near-infrared surface albedo are presented, their calculations are the same.
379	DO jb = 1, n_spectralbands
380
381	albedo(:,jb) = ( fraction_veg(:) ) * &
382	( (un-frac_snow_veg(:)) * albedo(:,jb) + &
383	( frac_snow_veg(:) ) * snowa_veg(:) )
384	albedo_snow(:) = albedo_snow(:) + (fraction_veg(:)) * &
385	(frac_snow_veg(:)) * snowa_veg(:)
386	!
387	DO jv = 1, nnobio ! Loop over # nobio surfaces
388	!
389	IF ( jv .EQ. iice ) THEN
390	alb_nobio = alb_ice(jb)
391	ELSE
392	WRITE(numout,*) 'jv=',jv
393	STOP 'DO NOT KNOW ALBEDO OF THIS SURFACE TYPE'
394	ENDIF
395
396	albedo(:,jb) = albedo(:,jb) + &
397	( frac_nobio(:,jv) ) * &
398	( (un-frac_snow_nobio(:,jv)) * alb_nobio + &
399	( frac_snow_nobio(:,jv) ) * snowa_nobio(:,jv) )
400	albedo_snow(:) = albedo_snow(:) + &
401	( frac_nobio(:,jv) ) * ( frac_snow_nobio(:,jv) ) * &
402	snowa_nobio(:,jv)
403	albedo_glob(:) = albedo_glob(:) + albedo(:,jb)
404
405	ENDDO ! jv = 1, nnobio
406
407	END DO ! jb = 1, n_spectralbands
408
409	! this accounts for the fact that the calculations are done for each
410	! spectra
411	DO ji = 1, npts
412	albedo_snow(ji) = (albedo_snow(ji))/REAL(n_spectralbands,r_std)
413	albedo_glob(ji) = (albedo_glob(ji))/REAL(n_spectralbands,r_std)
414	ENDDO
415
416	IF (bavard.GE.2) WRITE(numout,*) 'Leaving calculate_surface_albedo_with_snow '
417
418	END SUBROUTINE calculate_surface_albedo_with_snow
419
420
421
422	!! ==============================================================================================================================
423	!! SUBROUTINE : albedo_two_stream
424	!!
425	!>\BRIEF This subroutine calculates the albedo for two stream radiation transfer model. This routine includes
426	!! the effect of snow on the background albedo, through one of two methods.
427	!!
428	!! DESCRIPTION : The albedo for two stream radiation transfer model is calculated for both the visible and near-infrared
429	!! domain. First the mean albedo of the bare soil is calculated. Two options exist:
430	!! either the soil albedo depends on soil wetness (drysoil_frac variable), or the soil albedo
431	!! is set to a mean soil albedo value.
432	!!
433	!! NOTE: the main output variable, albedo, is an unweighted average of the direct and diffuse albedos
434	!! for each grid point...this is done in this way right now to be consistent with the new scheme,
435	!! but it doesn't have to be combined like that once we have an energy budget which can
436	!! use both types of light directly
437	!!
438	!! RECENT CHANGE(S): None
439	!!
440	!! MAIN OUTPUT VARIABLE(S): ::albedo
441	!!
442	!! REFERENCE(S) : B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski, N. Gobron and M. M. Verstraete (2006).
443	!! Simplifying the Interaction of Land Surfaces with Radiation for Relating Remote Sensing Products to Climate Models.
444	!! Journal of Geophysical Research. Vol 111, D02116.
445	!!
446	!! FLOWCHART : None
447	!! \n
448	!!
449	!!
450	!!
451	!_ ================================================================================================================================
452
453	SUBROUTINE albedo_two_stream(npts, nlevels_loc, &
454	drysoil_frac, lai, &
455	veget_max, sinang, soilalb_dry, soilalb_wet, &
456	frac_nobio, soilalb_moy, alb_bare, albedo, &
457	snow, snow_age, snow_nobio, snow_nobio_age, &
458	albedo_snow, albedo_glob, circ_class_biomass, &
459	circ_class_n, lai_split, z0_veg, Isotrop_Abs_Tot_p,&
460	Isotrop_Tran_Tot_p, laieff_fit, albedo_pft, frac_snow_nobio,&
461	frac_snow_veg)
462
463	!! 0. Variable and parameter declaration
464
465	!! 0.1 Input variables
466
467	INTEGER(i_std), INTENT(in) :: npts !! Domain size - Number of land pixels (unitless)
468	INTEGER(i_std), INTENT(in) :: nlevels_loc !! The number of vertical levels in the canopy (unitless)
469	REAL(r_std), DIMENSION(:), INTENT(in) :: drysoil_frac !! Fraction of dry soil (unitless)
470	REAL(r_std),DIMENSION (:,:), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
471
472	REAL(r_std), DIMENSION(:), INTENT(in) :: sinang !!
473	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! PFT coverage fraction of a PFT (= ind*cn_ind)
474	!! (m^2 m^{-2})
475	REAL(r_std), DIMENSION(:,:), INTENT(in) :: soilalb_dry !! Albedo values for the dry bare soil (unitless)
476	REAL(r_std), DIMENSION(:,:), INTENT(in) :: soilalb_wet !! Albedo values for the wet bare soil (unitless)
477	REAL(r_std), DIMENSION(:,:), INTENT(in) :: soilalb_moy !! Albedo values for the mean bare soil (unitless)
478	REAL(r_std),DIMENSION (:), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
479	REAL(r_std),DIMENSION (:,:), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
480	REAL(r_std),DIMENSION (:), INTENT(in) :: snow_age !! Snow age (days)
481	REAL(r_std),DIMENSION (:,:), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
482	REAL(r_std),DIMENSION (:,:), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
483	!! continental ice, lakes, etc. (unitless)
484	REAL(r_std), DIMENSION(:,:,:,:,:), &
485	INTENT(IN) :: circ_class_biomass !! Stem diameter
486	!! @tex $(m)$ @endtex
487	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: circ_class_n !! Number of trees within each circumference
488	REAL(r_std),DIMENSION(:),INTENT(IN) :: lai_split
489	REAL(r_std), DIMENSION(:), INTENT(in) :: z0_veg !! Roughness height of vegetated part (m)
490	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
491
492
493	!! 0.2 Output variables
494
495	REAL(r_std), DIMENSION(:,:), INTENT(out) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
496	REAL(r_std), DIMENSION (:,:), INTENT (out) :: albedo !! Albedo (two stream radiation transfer model)
497	!! for visible and near-infrared range
498	!! (unitless)
499	REAL(r_std),DIMENSION (:), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
500	REAL(r_std),DIMENSION (:), INTENT (out) :: albedo_glob !! Mean albedo (unitless ratio)
501	REAL(r_std),DIMENSION (:,:,:), INTENT (out) :: Isotrop_Abs_Tot_p !! Absorbed radiation per layer for photosynthesis
502	REAL(r_std),DIMENSION (:,:,:), INTENT (out) :: Isotrop_Tran_Tot_p !! Transmitted radiation per layer for photosynthesis
503	REAL(r_std), DIMENSION (:,:,:), INTENT (out) :: albedo_pft !! Albedo (two stream radiation transfer model)
504	!! for visible and near-infrared range
505	!! for each PFT (unitless)
506	REAL(r_std), DIMENSION(npts,nnobio), INTENT(out) :: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
507	!! (unitless ratio)
508	REAL(r_std), DIMENSION(npts,nvm), INTENT(out) :: frac_snow_veg !! the fraction of ground covered with snow, between 0 and 1
509
510
511
512	!! 0.3 Modified variables
513
514	!! 0.4 Local variables
515
516	REAL(r_std), DIMENSION(nlevels_loc,npts,nvm) :: laieff_collim !! Leaf Area Index Effective for direct light
517	REAL(r_std), DIMENSION(nlevels_loc,npts,nvm) :: laieff_isotrop !! Leaf Area Index Effective converts
518	!! 3D lai into 1D lai for two stream
519	!! radiation transfer model...this is for
520	!! isotropic light and only calculated once per day
521	!! @tex $(m^{2} m^{-2})$ @endtex
522	REAL(r_std),DIMENSION (nlevels_loc) :: laieff_isotrop_pft, laieff_collim_pft !!
523	REAL(r_std) :: cosine_sun_angle !! the cosine of the solar zenith angle
524	INTEGER(i_std) :: ks !! Index for visible and near-infraread range
525	INTEGER(i_std) :: ivm !! Index for vegetative PFTs
526	INTEGER(i_std) :: ipts !! Index for spatial domain
527	INTEGER(i_std) :: ilevel !! Index for canopy levels
528	REAL(r_std) :: leaf_reflectance
529	REAL(r_std) :: leaf_transmittance
530	REAL(r_std) :: br_base_temp,br_base_temp_collim,br_base_temp_isotrop
531	REAL(r_std) :: leaf_psd_temp
532	REAL(r_std) :: leaf_single_scattering_albedo
533	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Alb_Tot !! Collimated (direct) total albedo for a given level
534	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Tran_Uncoll !! collimated total transmission of uncollided light for a given level
535	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Tran_Coll !! collimated total transmission for a given level
536	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Abs_Tot !! collimated total absorption for a given level
537	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Alb_Tot !! isotropic (diffuse) total albedo for a given level
538	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Tran_Uncoll !! isotropic total transmission of uncollided light for a given level
539	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Tran_Coll !! isotropic total transmission for a given level
540	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Abs_Tot !! isotropic total absporbtion for a given level
541	INTEGER :: jv, inno
542	REAL(r_std) :: alb_nobio
543
544	REAL(r_std):: laieff_collim_1, laieff_isotrop_1, Collim_Alb_Tot_1,Collim_Tran_Tot_1,Collim_Abs_Tot_1,&
545	Isotrop_Alb_Tot_1,Isotrop_Tran_Tot_1,Isotrop_Abs_Tot_1,&
546	Collim_Tran_Uncollided_1, Isotrop_Tran_Uncollided_1 !! Values for the one level solution
547
548	REAL(r_std):: converged_albedo
549	REAL(r_std):: isotrop_angle
550	LOGICAL :: lconverged
551	LOGICAL :: lprint !! a flag for printing some debug statements
552	REAL(r_std), DIMENSION(npts,nvm,n_spectralbands) :: snowa_veg !! snow albedo due to vegetated surfaces, between 0 and 1
553	REAL(r_std), DIMENSION(npts,nnobio,n_spectralbands) :: snowa_nobio !! Albedo of snow covered area on continental ice,
554	!! lakes, etc. (unitless ratio)
555
556	!_ ================================================================================================================================
557
558
559	IF (bavard.GE.2) WRITE(numout,*) 'Entering albedo_two_stream'
560
561	!! 1. We will now calculate the background reflectance to be used in the model.
562	!! The parameters read in from the input file do not include the effect of snow,
563	!! so we need to figure out the snow's contribution for each grid space
564
565	!initialize some output variables
566	DO ipts = 1, npts
567	albedo_snow(ipts) = zero
568	ENDDO
569	Isotrop_Abs_Tot_p(:,:,:)=zero
570	Isotrop_Tran_Tot_p(:,:,:)=un
571	albedo_pft(:,:,:)=zero
572
573	! We need to calculate the LAI effective for this solar angle. The other LAI
574	! effective that we calculated was computed at an angle of 60.0 degrees for
575	! the isotropic contribution. Notice that, in many situations, the results
576	! will be exactly the same. One can show that Pgap is proportional to
577	! 1/cos(theta) if the level bottom is below all canopy elements, which
578	! means that the cos(theta) in the calculation of the LAIeff will be
579	! cancelled out
580
581	isotrop_angle=COS(60.0/180.0*pi)
582	DO ipts=1,npts
583	DO ivm=1,nvm
584	DO ilevel=1,nlevels_loc
585	laieff_isotrop(ilevel,ipts,ivm)=calculate_laieff_fit(isotrop_angle,laieff_fit(ipts,ivm,ilevel))
586	laieff_collim(ilevel,ipts,ivm)=calculate_laieff_fit(sinang(ipts),laieff_fit(ipts,ivm,ilevel))
587	! Because we are fitting these values, it's possible that some will
588	! drop below zero. This causes a serious problem here,
589	! so let's make sure this doesn't happen. Unfortunately, by
590	! doing this, we lose a test of things going wrong (an unfitted
591	! negative LAIeff can be the sign of another problem), but we
592	! don't really have a choice.
593	IF(laieff_isotrop(ilevel,ipts,ivm) .LT. min_stomate) &
594	laieff_isotrop(ilevel,ipts,ivm)=zero
595	IF(laieff_collim(ilevel,ipts,ivm) .LT. min_stomate) &
596	laieff_collim(ilevel,ipts,ivm)=zero
597	ENDDO
598	ENDDO
599	ENDDO
600
601	IF(ld_albedo)THEN
602	DO ivm=1,nvm
603	IF(ivm == test_pft)THEN
604	DO ipts = 1, npts
605	IF(ipts == test_grid)THEN
606	WRITE(numout,*) 'Laieff_collim and laieff_isotrop in albedo'
607	WRITE(numout,) 'sinang(ipts): ',ACOS(sinang(ipts))/pi180.0_r_std
608	DO ilevel=nlevels_loc,1,-1
609	WRITE(numout,*) ilevel,laieff_collim(ilevel,ipts,ivm),&
610	laieff_isotrop(ilevel,ipts,ivm)
611	END DO
612	DO ilevel=nlevels_loc,1,-1
613	WRITE(numout,*) 'Fitting parameters: ',laieff_fit(ipts,ivm,ilevel)%a,laieff_fit(ipts,ivm,ilevel)%b,&
614	laieff_fit(ipts,ivm,ilevel)%c,laieff_fit(ipts,ivm,ilevel)%d,laieff_fit(ipts,ivm,ilevel)%e
615	ENDDO
616	ENDIF
617	ENDDO
618	ENDIF
619	ENDDO
620	ENDIF
621
622	CALL calculate_snow_albedo(npts, sinang, snow, snow_nobio, &
623	snow_age, snow_nobio_age, frac_nobio, albedo_snow, &
624	snowa_veg, frac_snow_veg, snowa_nobio, frac_snow_nobio, &
625	veget_max, z0_veg)
626
627	!!$ !+++++ CHECK +++++++
628	!!$ DO ipts=1,npts
629	!!$ DO ivm=1,nvm
630	!!$ IF(frac_snow_veg(ipts,ivm) .LT. 0.0)THEN
631	!!$ WRITE(numout,*) 'SNOWFRAC ipts ',ipts,frac_snow_veg(ipts,ivm)
632	!!$ CALL ipslerr_p (3,'albedo', 'snow frac error','','')
633	!!$ ENDIF
634	!!$ ENDDO
635	!!$ DO inno=1,nnobio
636	!!$ IF(frac_snow_nobio(ipts,inno) .LT. 0.0)THEN
637	!!$ WRITE(numout,*) 'SNOWFRAC ipts ',ipts,frac_snow_nobio(ipts,inno)
638	!!$ CALL ipslerr_p (3,'albedo', 'nobio snow frac error','','')
639	!!$ ENDIF
640	!!$ ENDDO
641	!!$ ENDDO
642	!!$
643	!!$ !++++++++++++++++++
644
645	! initialize the albedo in every square for both spectra by calculating the
646	! bare soil albedo
647
648	DO ks = 1, n_spectralbands
649	DO ipts = 1, npts
650	! If 'alb_bare_model' is set to TRUE, the soil albedo calculation
651	! depends on soil moisture
652	! If 'alb_bare_model' is set to FALSE, the mean soil albedo is used
653	! without the dependance on soil moisture
654	! see subroutine 'conveg_soilalb'
655	IF ( alb_bare_model ) THEN
656	alb_bare(ipts,ks) = soilalb_wet(ipts,ks) + drysoil_frac(ipts) * &
657	(soilalb_dry(ipts,ks) - soilalb_wet(ipts,ks))
658	ELSE
659	alb_bare(ipts,ks) = soilalb_moy(ipts,ks)
660	ENDIF
661
662	! we need to take into account any snow that has fallen on the bare soil
663	! since there may be some bare soil, initialize the albedo with this fraction
664	albedo(ipts,ks) = veget_max(ipts,1) * (alb_bare(ipts,ks)*&
665	(un-frac_snow_veg(ipts,1))+frac_snow_veg(ipts,1)*snowa_veg(ipts,1,ks))
666	ENDDO ! ipts = 1, npts
667	ENDDO ! ks = 1, n_spectralbands
668
669	!! 3. Calculation of albedo of the whole grid cell, taking into account all
670	!! PFTs (except bare soil) and spectral bands
671
672	DO ipts = 1, npts ! loop over the grid squares
673
674	! it appears that sinang is the sin of the angle from the horizon
675	! this is equal to the cos of the zenith angle, which is what we need
676	cosine_sun_angle = sinang(ipts)
677
678	! For the coupled model, the angles can be negative.
679	! Offline, they are always between zero and 90 degrees.
680	IF(cosine_sun_angle .LE. min_sechiba)THEN
681	! it's night, so no sunlight...don't need to calculate albedo. Set it equal to
682	! zero, because the snow albedo has already been calculated and it looks funny on
683	! the coupled run maps otherwise.
684	albedo(ipts,:)=zero
685	CYCLE
686	ENDIF
687
688
689	! are there any non-bio PFTs here that we need to take into account?
690	DO ks = 1, n_spectralbands ! Loop over # of spectra
691	DO jv = 1, nnobio
692	! now update the albedo
693	IF ( jv .EQ. iice ) THEN
694	alb_nobio = alb_ice(ks)
695	ELSE
696	WRITE(numout,*) 'jv=',jv
697	STOP 'DO NOT KNOW ALBEDO OF THIS SURFACE TYPE'
698	ENDIF
699
700	! this takes into account both snow covered and non-snow covered non
701	! bio regios in all grid squares
702	albedo(ipts,ks) = albedo(ipts,ks) + &
703	( frac_nobio(ipts,jv) ) * &
704	( (un-frac_snow_nobio(ipts,jv)) * alb_nobio + &
705	( frac_snow_nobio(ipts,jv) ) * snowa_nobio(ipts,jv,ks) )
706	albedo_snow(ipts) = albedo_snow(ipts) + &
707	( frac_nobio(ipts,jv) ) * ( frac_snow_nobio(ipts,jv) ) * &
708	snowa_nobio(ipts,jv,ks)/REAL(n_spectralbands,r_std)
709
710	ENDDO ! jv = 1, nnobio
711	ENDDO ! ks = 1, n_spectralbands
712
713	DO ivm = 2, nvm ! Loop over # of PFTs
714
715	! for this grid square and this vegetation type, we have a set LAI
716	! effective
717	laieff_collim_pft(1:nlevels_loc) = laieff_collim(1:nlevels_loc,ipts,ivm)
718	laieff_isotrop_pft(1:nlevels_loc) = laieff_isotrop(1:nlevels_loc,ipts,ivm)
719
720	IF(veget_max(ipts,ivm) == zero)THEN
721	! this vegetation type is not present, so no reason to do the
722	! calculation
723	CYCLE
724	ENDIF
725
726	DO ks = 1, n_spectralbands ! Loop over # of spectra
727
728	! set single scattering albedo and preferred scattering direction
729	leaf_single_scattering_albedo = leaf_ssa(ivm,ks)
730	leaf_psd_temp = leaf_psd(ivm,ks)
731
732	! calculate the rleaf and tleaf from wl=rl+tl and dl=rl/tl
733
734	leaf_transmittance = leaf_single_scattering_albedo/ &
735	(leaf_psd_temp+1)
736	leaf_reflectance = leaf_psd_temp * &
737	leaf_transmittance
738
739	! We need to take into account the effect of snow on the background
740	! reflectance here
741
742	! from the snow above I can calculate the true background reflectance
743	! for this PFT
744	br_base_temp= (un-frac_snow_veg(ipts,ivm)) * bgd_reflectance(ivm,ks)&
745	+ ( frac_snow_veg(ipts,ivm) ) * snowa_veg(ipts,ivm,ks)
746
747	! At this step, we assume that there is no difference between the
748	! direct and the diffuse background reflectance, which is true for
749	! the background but not the canopy albedo
750	br_base_temp_collim=br_base_temp
751	br_base_temp_isotrop=br_base_temp
752
753	! This prints out some debugging information
754	IF(ld_albedo .AND. ivm == test_pft .AND. ipts == test_grid)THEN
755	DO ilevel=nlevels_loc,1,-1
756	laieff_collim_1=laieff_collim_pft(ilevel)
757	laieff_isotrop_1=laieff_isotrop_pft(ilevel)
758	call twostream_solver(leaf_reflectance,leaf_transmittance,&
759	br_base_temp_collim,br_base_temp_isotrop,&
760	cosine_sun_angle,Collim_Alb_Tot_1,Collim_Tran_Tot_1,&
761	Collim_Abs_Tot_1,Isotrop_Alb_Tot_1,Isotrop_Tran_Tot_1,&
762	Isotrop_Abs_Tot_1,laieff_collim_1, laieff_isotrop_1, &
763	Collim_Tran_Uncollided_1, Isotrop_Tran_Uncollided_1)
764
765
766	CALL print_debugging_albedo_info(ilevel,cosine_sun_angle,&
767	leaf_reflectance,leaf_transmittance,&
768	laieff_collim_1,laieff_isotrop_1,&
769	br_base_temp_collim,br_base_temp_isotrop,Collim_Alb_Tot_1,&
770	Collim_Tran_Tot_1,Collim_Abs_Tot_1,&
771	Isotrop_Alb_Tot_1,Isotrop_Tran_Tot_1,Isotrop_Abs_Tot_1)
772
773	ENDDO ! loop over levels
774	ENDIF ! debugging
775
776	! Calculate the one level approach to use as a template. Notice that
777	! this is not actually needed for the multilevel approach, since there
778	! are no parameters that we can tune to get the multilevel result
779	! closer to the one level result.
780
781	! First, we need to recombine all the levels into a one level LAI
782	laieff_collim_1=SUM(laieff_collim_pft(1:nlevels_loc))
783	laieff_isotrop_1=SUM(laieff_isotrop_pft(1:nlevels_loc))
784	call twostream_solver(leaf_reflectance, leaf_transmittance, &
785	br_base_temp_collim, br_base_temp_isotrop, &
786	cosine_sun_angle, Collim_Alb_Tot_1, Collim_Tran_Tot_1, &
787	Collim_Abs_Tot_1, &
788	Isotrop_Alb_Tot_1, Isotrop_Tran_Tot_1, Isotrop_Abs_Tot_1, &
789	laieff_collim_1, &
790	laieff_isotrop_1, Collim_Tran_Uncollided_1, &
791	Isotrop_Tran_Uncollided_1)
792
793	! Print out more debugging information
794	IF(ld_albedo .AND. ipts == test_grid .AND. ivm == test_pft) &
795	CALL print_debugging_albedo_info(1, cosine_sun_angle, &
796	leaf_reflectance, leaf_transmittance, laieff_collim_1, &
797	laieff_isotrop_1, br_base_temp_collim, br_base_temp_isotrop, &
798	Collim_Alb_Tot_1, Collim_Tran_Tot_1, Collim_Abs_Tot_1, &
799	Isotrop_Alb_Tot_1, Isotrop_Tran_Tot_1, Isotrop_Abs_Tot_1)
800
801	! If we only have one canopy level, there is no reason to try to fit
802	! the nlevel case, since we just calculated the answer above
803	IF(nlevels_loc == 1)THEN
804	! notice that this albedo assumes that diffuse and direct albedo
805	! are of equal importance
806	converged_albedo=(Collim_Alb_Tot_1+Isotrop_Alb_Tot_1)/2.0_r_std
807	albedo(ipts,ks) = albedo(ipts,ks) + veget_max(ipts,ivm)*converged_albedo
808	albedo_pft(ipts,ivm,ks)=veget_max(ipts,ivm)*converged_albedo
809	Isotrop_Abs_Tot_p(ipts,ivm,:) = Isotrop_Abs_Tot_1
810	Isotrop_Tran_Tot_p(ipts,ivm,:) = Isotrop_Tran_Tot_1
811
812	CYCLE
813	ENDIF
814
815
816	! now solve the multilevel scheme
817	IF(ld_albedo .AND. test_pft == ivm .AND. test_grid == ipts)THEN
818	lprint=.TRUE.
819	ELSE
820	lprint=.FALSE.
821	ENDIF
822
823	CALL multilevel_albedo(nlevels_loc, cosine_sun_angle, leaf_single_scattering_albedo, &
824	leaf_psd_temp, br_base_temp_collim, br_base_temp_isotrop, &
825	laieff_collim_pft, laieff_isotrop_pft, lconverged, &
826	Collim_Alb_Tot, Collim_Tran_Coll, Collim_Abs_Tot, Isotrop_Alb_Tot, &
827	Isotrop_Tran_Coll, Isotrop_Abs_Tot, Collim_Tran_Uncoll, Isotrop_Tran_Uncoll, lprint)
828
829	! This is a lot of information printed out for the paper
830	! on the multilevel albedo scheme. Probably not useful
831	! for anyone else, but I'm keeping it in until the paper is published.
832	! This is why I'm protecting it with an IF statement. I do not use
833	! ld_albedo since it really is too much information for normal usage.
834	IF(.FALSE.)THEN
835
836	WRITE(6,'(A,2F14.10)') 'Background reflectance ',&
837	br_base_temp_collim,br_base_temp_isotrop
838	WRITE(6,'(A,F14.10)') 'Single scattering albedo (wl): ',&
839	leaf_single_scattering_albedo
840	WRITE(6,'(A,F14.10)') 'Preferred scattering direction (dl): ',&
841	leaf_psd_temp
842
843	WRITE(6,) 'Collimated light ',ACOS(cosine_sun_angle)/pi180.0
844	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' 1 & ',laieff_collim_1,&
845	' & ',Collim_Alb_Tot_1,' & ',&
846	Collim_Tran_Tot_1-Collim_Tran_Uncollided_1,' & ',&
847	Collim_Tran_Uncollided_1,' & ',&
848	Collim_Alb_Tot_1+Collim_Tran_Tot_1,' \\ '
849	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' 2 & ',laieff_collim_1,&
850	' & ',Collim_Alb_Tot(nlevels_loc),' & ',&
851	Collim_Tran_Coll(1),' & ',Collim_Tran_Uncoll(1),' & ',&
852	Collim_Alb_Tot(nlevels_loc)+Collim_Tran_Coll(1)+&
853	Collim_Tran_Uncoll(1),' \\ '
854	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' Diff & ',laieff_collim_1,&
855	' & ',ABS(Collim_Alb_Tot_1-Collim_Alb_Tot(nlevels_loc)),' & ',&
856	ABS((Collim_Tran_Tot_1-Collim_Tran_Uncollided_1)-&
857	Collim_Tran_Coll(1)),' & ',&
858	ABS(Collim_Tran_Uncollided_1-Collim_Tran_UnColl(1)),' & ',&
859	ABS((Collim_Alb_Tot_1+Collim_Tran_Tot_1)-&
860	(Collim_Alb_Tot(nlevels_loc)+Collim_Tran_Coll(1)+&
861	Collim_Tran_Uncoll(1))),' \\ '
862
863	WRITE(6,*) 'Isotropic light '
864	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' 1 & ',&
865	laieff_isotrop_1,' & ',Isotrop_Alb_Tot_1,' & ',&
866	Isotrop_Tran_Tot_1-Isotrop_Tran_Uncollided_1,' & ',&
867	Isotrop_Tran_Uncollided_1,' & ',&
868	Isotrop_Alb_Tot_1+Isotrop_Tran_Tot_1,' \\ '
869	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' 2 & ',&
870	laieff_isotrop_1,' & ',Isotrop_Alb_Tot(nlevels_loc),' & ',&
871	Isotrop_Tran_Coll(1),' & ',Isotrop_Tran_Uncoll(1),' & ',&
872	Isotrop_Alb_Tot(nlevels_loc)+Isotrop_Tran_Coll(1)+&
873	Isotrop_Tran_Uncoll(1),' \\ '
874	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' Diff & ',&
875	laieff_isotrop_1,&
876	' & ',ABS(Isotrop_Alb_Tot_1-Isotrop_Alb_Tot(nlevels_loc)),' & ',&
877	ABS((Isotrop_Tran_Tot_1-Isotrop_Tran_Uncollided_1)-&
878	Isotrop_Tran_Coll(1)),' & ',&
879	ABS(Isotrop_Tran_Uncollided_1-Isotrop_Tran_UnColl(1)),' & ',&
880	ABS((Isotrop_Alb_Tot_1+Isotrop_Tran_Tot_1)-&
881	(Isotrop_Alb_Tot(nlevels_loc)+Isotrop_Tran_Coll(1)+&
882	Isotrop_Tran_Uncoll(1))),' \\ '
883
884	WRITE(6,1010) 'COLLIMTRAN ',laieff_collim_1, &
885	laieff_collim_pft(nlevels_loc),1,&
886	leaf_single_scattering_albedo,br_base_temp_collim,&
887	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
888	Collim_Tran_Tot_1,(Collim_Tran_Coll(1)+Collim_Tran_Uncoll(1))
889	WRITE(6,1010) 'COLLIMALB ',laieff_collim_1, &
890	laieff_collim_pft(nlevels_loc),2,&
891	leaf_single_scattering_albedo,br_base_temp_collim,&
892	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
893	Collim_Alb_Tot_1,Collim_Alb_Tot(nlevels_loc)
894	WRITE(6,1010) 'COLLIMABSCAN ',laieff_collim_1, &
895	laieff_collim_pft(nlevels_loc),3,&
896	leaf_single_scattering_albedo,br_base_temp_collim,&
897	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
898	Collim_Abs_Tot_1,SUM(Collim_Abs_Tot(:))
899	WRITE(6,1010) 'COLLIMABSSOIL ',laieff_collim_1, &
900	laieff_collim_pft(nlevels_loc),4,&
901	leaf_single_scattering_albedo,br_base_temp_collim,&
902	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
903	un-Collim_Alb_Tot_1-Collim_Abs_Tot_1,un-&
904	Collim_Alb_Tot(nlevels_loc)-SUM(Collim_Abs_Tot(:))
905
906	WRITE(6,1010) 'ISOTROPTRAN ',laieff_isotrop_1, &
907	laieff_isotrop_pft(nlevels_loc),1,&
908	leaf_single_scattering_albedo,br_base_temp_isotrop,&
909	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
910	Isotrop_Tran_Tot_1,(Isotrop_Tran_Coll(1)+Isotrop_Tran_Uncoll(1))
911	WRITE(6,1010) 'ISOTROPALB ',laieff_isotrop_1, &
912	laieff_isotrop_pft(nlevels_loc),2,&
913	leaf_single_scattering_albedo,br_base_temp_isotrop,&
914	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
915	Isotrop_Alb_Tot_1,Isotrop_Alb_Tot(nlevels_loc)
916	WRITE(6,1010) 'ISOTROPABSCAN ',laieff_isotrop_1, &
917	laieff_isotrop_pft(nlevels_loc),3,&
918	leaf_single_scattering_albedo,br_base_temp_isotrop,&
919	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
920	Isotrop_Abs_Tot_1,SUM(Isotrop_Abs_Tot(:))
921	WRITE(6,1010) 'ISOTROPABSSOIL ',laieff_isotrop_1, &
922	laieff_isotrop_pft(nlevels_loc),4,&
923	leaf_single_scattering_albedo,br_base_temp_isotrop,&
924	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
925	un-Isotrop_Alb_Tot_1-Isotrop_Abs_Tot_1,un-&
926	Isotrop_Alb_Tot(nlevels_loc)-SUM(Isotrop_Abs_Tot(:))
927
928	1010 FORMAT(A,2(F12.6,2X),I4,2X,6(F14.8,2X))
929
930	ENDIF ! debugging
931
932	! For our total albedo, we are just taking a simple mean of the diffuse
933	! and the direct light, which means we lose some information. This
934	! might be changed in the future. We also need to weight this by the
935	! percentage of nonbio land, but this appears to be included in veget_max
936
937	converged_albedo =(Collim_Alb_Tot(nlevels_loc)+Isotrop_Alb_Tot(nlevels_loc))/2.0_r_std
938	albedo(ipts,ks) = albedo(ipts,ks) + veget_max(ipts,ivm)*converged_albedo
939	albedo_pft(ipts,ivm,ks)=veget_max(ipts,ivm)*converged_albedo
940
941	! Save the absorbed radiation for photosynthesis, we need only the
942	! visible range and use the isotropic part, this can be change later
943	! to distinguish between sunlit and shaded leaves.
944	! Be careful here
945	IF (ks == 1) THEN
946	! Test if Collim_Abs_Tot has negative values
947	! If so, set it to min_sechiba
948	DO ilevel=1,nlevels_loc
949	IF (Isotrop_Abs_Tot(ilevel) .LT. zero) THEN
950	Isotrop_Abs_Tot(ilevel) = min_sechiba
951	ENDIF
952	ENDDO
953	!
954	Isotrop_Abs_Tot_p(ipts,ivm,:) = Isotrop_Abs_Tot(:)
955	Isotrop_Tran_Tot_p(ipts,ivm,:) = Isotrop_Tran_Coll(:) + Isotrop_Tran_Uncoll(:)
956
957	! Notice that this light is actually cumulative, not per
958	! level! This was needed for debugging purposes and running
959	! tests. However, the photosynthesis routines need the light
960	! transmitted per level, i.e. if there is zero LAI
961	! in a level, the transmitted light will be one.
962	DO ilevel=1,nlevels_loc-1
963	IF(Isotrop_Tran_Tot_p(ipts,ivm,ilevel+1) .GT. alb_threshold)THEN
964	Isotrop_Tran_Tot_p(ipts,ivm,ilevel)=&
965	Isotrop_Tran_Tot_p(ipts,ivm,ilevel)/Isotrop_Tran_Tot_p(ipts,ivm,ilevel+1)
966	ELSE
967	! Here, we really don't know anything about how much light is transmitted
968	! in this layer, but there is no light reaching it from above so we
969	! can safely assume no photosynthesis takes place. This is equivalent
970	! to assuming it has no LAI, which means the transmission will be unity.
971	Isotrop_Tran_Tot_p(ipts,ivm,ilevel)=un
972	ENDIF
973	ENDDO
974
975	! Debugging information
976	IF (ld_albedo .AND. ivm == test_pft .AND. ipts == test_grid) THEN
977	WRITE(numout,*) 'Absorbed light in albedo.f90: ', ipts, ivm
978	DO ilevel=nlevels_loc,1,-1
979	WRITE(numout,'(I5,10F20.10)') ilevel, Isotrop_Abs_Tot_p(ipts,ivm,ilevel),&
980	Isotrop_Tran_Tot_p(ipts,ivm,ilevel),&
981	laieff_isotrop(ilevel,ipts,ivm),laieff_collim(ilevel,ipts,ivm)
982	ENDDO
983	ENDIF
984
985	ENDIF ! IF ks==1
986
987	ENDDO ! ks = 1, n_spectralbands
988
989	ENDDO ! ivm = 2, nvm
990
991	ENDDO ! ipts = 1, npts
992
993	! now we need to average the albedo over all our spectra, so we can pass it to
994	! methods which do not distinguish between different spectral bands
995	albedo_glob(:) = zero
996	DO ks = 1, n_spectralbands ! Loop over # of spectra
997	albedo_glob(:) = albedo_glob(:) + albedo(:,ks)
998	ENDDO
999	albedo_glob(:)=albedo_glob(:)/REAL(n_spectralbands,r_std)
1000
1001	! WRITE(,) '281013 in albedo_two_stream, Isotrop_Abs_Tot_p is: ', Isotrop_Abs_Tot_p
1002
1003	IF (bavard.GE.2) WRITE(numout,*) 'Leaving albedo_two_stream'
1004
1005	END SUBROUTINE albedo_two_stream
1006
1007	!! ==============================================================================================================================
1008	!! SUBROUTINE : twostream_solver
1009	!!
1010	!>\BRIEF : Computes the two-stream albedo solution for a level with given single scatterer properties
1011	!! and a defined background albedo
1012	!!
1013	!! DESCRIPTION : This solution is the two-stream solution of Pinty et al (2006) for a vegetation level above
1014	!! an isotropically reflecting background. It breaks down the problem into three parts, solving
1015	!! each part for the case of diffuse (isotropic) and direct (collimated) light. The three parts are,
1016	!! 1) The term due to light that does not interact at all with the canopy (Black Canopy)
1017	!! 2) Light that does not interact at all with the background (Black Background)
1018	!! 3) Light that bounces between the background and the canopy
1019	!! This routine (and the routines it uses) were received directly from Bernard Pinty, and only some
1020	!! minor modifcations were made, in addition to more documentation
1021	!!
1022	!!
1023	!!
1024	!! RECENT CHANGE(S): None
1025	!!
1026	!! MAIN OUTPUT VARIABLE(S): Collim_Alb_Tot,Collim_Tran_Tot,
1027	!! Collim_Abs_Tot,Isotrop_Alb_Tot,Isotrop_Tran_Tot,Isotrop_Abs_To
1028	!!
1029	!! REFERENCE(S) : B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski, N. Gobron and M. M. Verstraete (2006).
1030	!! Simplifying the Interaction of Land Surfaces with Radiation for Relating Remote Sensing Products to Climate Models.
1031	!! Journal of Geophysical Research. Vol 111, D02116.
1032	!!
1033	!! FLOWCHART : None
1034	!! \n
1035	!_ ================================================================================================================================
1036
1037
1038	SUBROUTINE twostream_solver(leaf_reflectance, leaf_transmittance, background_reflectance_collim, background_reflectance_isotrop, &
1039	cosine_sun_angle, Collim_Alb_Tot, Collim_Tran_Tot, Collim_Abs_Tot, &
1040	Isotrop_Alb_Tot, Isotrop_Tran_Tot, Isotrop_Abs_Tot, laieff_collim, &
1041	laieff_isotrop, Collim_Tran_Uncollided, Isotrop_Tran_Uncollided)
1042
1043
1044	!! 0. Variables and parameter declaration
1045	!! 0.1 Input variables
1046	REAL(r_std), INTENT(IN) :: leaf_reflectance !! effective leaf reflectance, between 0 and 1
1047	!! @tex $()$ @endtex
1048	REAL(r_std), INTENT(IN) :: leaf_transmittance !! effective leaf transmittance, between 0 and 1
1049	!! @tex $()$ @endtex
1050	REAL(r_std), INTENT(IN) :: background_reflectance_collim !! background reflectance for direct radiation,
1051	!! between 0 and 1
1052	REAL(r_std), INTENT(IN) :: background_reflectance_isotrop !! background reflectance for diffuse radiation,
1053	!! between 0 and 1
1054	REAL(r_std), INTENT(IN) :: cosine_sun_angle !! cosine of the solar zenith angle, between 0 and 1
1055	!! @tex $()$ @endtex
1056	REAL(r_std), INTENT(IN) :: laieff_collim !! Effective Leaf Area Index, computed at current sun angle
1057	REAL(r_std), INTENT(IN) :: laieff_isotrop !! Effective Leaf Area Index, computed at 60 degrees
1058	!! @tex $(m^{2} m^{-2})$ @endtex
1059
1060	!! 0.2 Output variables
1061	!! notice that these variables (absorption + transmission + reflection) do not necessarily add up to 1 due to multiple scattering
1062	REAL(r_std), INTENT(OUT) :: Collim_Alb_Tot !! collimated total albedo from this level, between 0 and 1
1063	!! @tex $()$ @endtex
1064	REAL(r_std), INTENT(OUT) :: Collim_Tran_Tot !! collimated total transmission through this level, between 0 and 1
1065	!! @tex $()$ @endtex
1066	REAL(r_std), INTENT(OUT) :: Collim_Abs_Tot !! collimated total absorption by this level, between 0 and 1
1067	!! @tex $()$ @endtex
1068	REAL(r_std), INTENT(OUT) :: Isotrop_Alb_Tot !! isotropic total albedo from this level, between 0 and 1
1069	!! @tex $()$ @endtex
1070	REAL(r_std), INTENT(OUT) :: Isotrop_Tran_Tot !! isotropic total transmission through this level, between 0 and 1
1071	!! @tex $()$ @endtex
1072	REAL(r_std), INTENT(OUT) :: Isotrop_Abs_Tot !! isotropic total absorption by this level, between 0 and 1
1073	!! @tex $()$ @endtex
1074	REAL(r_std), INTENT(OUT) :: Collim_Tran_Uncollided !! collimated uncollied transmission through this level, between 0 and 1
1075	!! @tex $()$ @endtex
1076	REAL(r_std), INTENT(OUT) :: Isotrop_Tran_Uncollided !! isotropic uncollied transmission through this level, between 0 and 1
1077	!! @tex $()$ @endtex
1078
1079	!! 0.3 Modified variables
1080
1081	!! 0.4 Local variables
1082	LOGICAL :: ok
1083	REAL(r_std), DIMENSION(4) :: gammaCoeffs
1084	REAL(r_std), DIMENSION(4) :: gammaCoeffs_star
1085	REAL(r_std) :: tauprimetilde
1086	REAL(r_std) :: tauprimestar
1087	REAL(r_std) :: sun_zenith_angle_radians
1088	REAL(r_std), PARAMETER :: isotropic_cosine_constant=0.5/0.705
1089
1090	!calculated fluxes
1091
1092	REAL(r_std) :: Collim_Alb_BB
1093	REAL(r_std) :: Collim_Tran_BB
1094	REAL(r_std) :: Collim_Abs_BB
1095	REAL(r_std) :: Isotrop_Alb_BB
1096	REAL(r_std) :: Isotrop_Tran_BB
1097	REAL(r_std) :: Isotrop_Abs_BB
1098	REAL(r_std) :: Collim_Tran_BC
1099	REAL(r_std) :: Isotrop_Tran_BC
1100	REAL(r_std) :: Collim_Tran_TotalOneWay
1101	REAL(r_std) :: Isotrop_Tran_TotalOneWay
1102	REAL(r_std) :: Below_reinject_rad_collim,Below_reinject_rad_isotrop
1103
1104	!_ ================================================================================================================================
1105
1106	IF (bavard.GE.2) WRITE(numout,*) 'Entering twostream_solver'
1107
1108	! convert angular values
1109
1110	sun_zenith_angle_radians = acos(cosine_sun_angle)
1111
1112	! calculate the 4 gamma coefficients both in isotropic and collimated illumination
1113
1114	call gammas(leaf_reflectance, leaf_transmittance, cosine_sun_angle, gammaCoeffs)
1115	call gammas(leaf_reflectance,leaf_transmittance,isotropic_cosine_constant,&
1116	gammaCoeffs_star)
1117
1118	! estimate the effective value of the optical thickness
1119
1120	tauprimetilde = 0.5_r_std * laieff_collim
1121	!tauprimetilde = 1.
1122
1123	! Should become zenit angle dependent (=60Â°) calculated by the Monte Carlo
1124	! photon shooting routine
1125
1126	tauprimestar = 0.5_r_std * laieff_isotrop
1127	!tauprimestar = 1.
1128
1129	! +++++++++++++ BLACK BACKGROUND ++++++++++++++++++++++
1130	! * Apply the black-background 2stream solution
1131	! * These equations are written for the part of the incoming radiation
1132	! * that never hits the background but does interact with the vegetation
1133	! * canopy.
1134	! *
1135	! * Note : the same routine dhrT1() is used for both the isotropic and
1136	! * collimated illumination conditions but the calling arguments differ.
1137	! * (especially the solar angle used).
1138	! */
1139
1140	! /* 1) collimated source */
1141	ok=dhrT1(leaf_reflectance,leaf_transmittance,&
1142	gammaCoeffs(1),gammaCoeffs(2),gammaCoeffs(3),gammaCoeffs(4),&
1143	sun_zenith_angle_radians, tauprimetilde,&
1144	Collim_Alb_BB,Collim_Tran_BB,Collim_Abs_BB)
1145	! /* 2) isotropic source */
1146	ok=dhrT1(leaf_reflectance,leaf_transmittance,&
1147	gammaCoeffs_star(1),gammaCoeffs_star(2),gammaCoeffs_star(3),&
1148	gammaCoeffs_star(4),&
1149	acos(isotropic_cosine_constant),tauprimestar,&
1150	Isotrop_Alb_BB,Isotrop_Tran_BB,Isotrop_Abs_BB)
1151
1152
1153	! +++++++++++++ BLACK CANOPY ++++++++++++++++++++++
1154	! * Apply the black canopy solution.
1155	! * These equations hold for the part of the incoming radiation
1156	! * that do not interact with the vegetation, travelling through
1157	! * its gaps.
1158	! */
1159
1160	! /* 1) collimated source */
1161	IF (cosine_sun_angle .NE. 0) THEN
1162	Collim_Tran_BC = exp( - tauprimetilde/cosine_sun_angle)
1163	Collim_Tran_Uncollided = Collim_Tran_BC
1164	ENDIF
1165	! /* 2) isotropic source */
1166	Isotrop_Tran_BC = TBarreUncoll_exact(tauprimestar)
1167	Isotrop_Tran_Uncollided = Isotrop_Tran_BC
1168
1169	! /* Total one-way transmissions:
1170	! * The vegetation canopy is crossed (one way) by the uncollided radiation
1171	! * (black canopy) and the collided one (black background). */
1172
1173	! /* 1) collimated source */
1174	Collim_Tran_TotalOneWay = Collim_Tran_BC + Collim_Tran_BB
1175
1176	! /* 2) isotropic source */
1177	Isotrop_Tran_TotalOneWay = Isotrop_Tran_BC + Isotrop_Tran_BB
1178
1179	! * The Below_reinject_rad describes the process of reflecting toward the
1180	! * background the upward travelling radiation (re-emitted from below the
1181	! * canopy). It appears in the coupling equations as the limit of the series:
1182	! * 1 + rgrbv + (rgrbv)^2 + (rg*rbv)^3 + ...
1183	! * where rg is the background_reflectance and rbv is Isotrop_Alb_BB
1184	! * (with Isotrop describing the Lambertian reflectance of the background).
1185	! */
1186	! this might be involved in Eq. 27,
1187	Below_reinject_rad_collim = un / (un - background_reflectance_collim*Isotrop_Alb_BB)
1188	Below_reinject_rad_isotrop = un / (un - background_reflectance_isotrop*Isotrop_Alb_BB)
1189
1190	!/* TOTAL ALBEDO */
1191	! /* 1) collimated source */
1192	Collim_Alb_Tot = Collim_Alb_BB + &
1193	background_reflectance_collim * Collim_Tran_TotalOneWay * &
1194	Isotrop_Tran_TotalOneWay * Below_reinject_rad_collim
1195	! /* 2) isotropic source */
1196	Isotrop_Alb_Tot = Isotrop_Alb_BB + &
1197	background_reflectance_isotrop * Isotrop_Tran_TotalOneWay * &
1198	Isotrop_Tran_TotalOneWay * Below_reinject_rad_isotrop ! this seems to be
1199	! exactly Eq. 33
1200
1201	!/* TOTAL TRANSMITION TO THE BACKGROUND LEVEL */
1202	! /* 1) collimated source */
1203	Collim_Tran_Tot = Collim_Tran_TotalOneWay * Below_reinject_rad_collim ;
1204
1205	! /* 2) isotropic source */
1206	Isotrop_Tran_Tot = Isotrop_Tran_TotalOneWay * Below_reinject_rad_isotrop ;
1207
1208	!/* TOTAL ABSORPTION BY THE VEGETATION LEVEL */
1209	! /* 1) collimated source */
1210	Collim_Abs_Tot = un - (Collim_Tran_Tot + Collim_Alb_Tot) + &
1211	background_reflectance_collim * Collim_Tran_Tot;
1212	! /* 2) isotropic source */
1213	Isotrop_Abs_Tot = un - (Isotrop_Tran_Tot + Isotrop_Alb_Tot) + &
1214	background_reflectance_isotrop * Isotrop_Tran_Tot;
1215
1216	IF (bavard.GE.2) WRITE(numout,*) 'Exiting twostream_solver'
1217
1218	END SUBROUTINE twostream_solver
1219
1220
1221	!! ============================================================================\n
1222	!! FUNCTION : dhrT1
1223	!!
1224	!>\BREF
1225	!!
1226	!! DESCRIPTION : Taken directly from the Fortan code provided by Pinty et al for their \n
1227	!! two-stream model. The only changes are the REAL types. This is used
1228	!! to compute the black background absorption, transmission, and reflectance
1229	!! in Pintys scheme and it's based on the older model of Meador and Weaver...
1230	!! Pinty 2006 also includes a short discussion of it in Appendix A
1231	!!
1232	!! RECENT CHANGE(S): None\n
1233	!!
1234	!! RETURN VALUE : dhrT1
1235	!!
1236	!! REFERENCE(S) : Meador and Weaver, 'Two-stream approximations to radiative transfer in
1237	!! planetary atmosphere: a unified description of existing methods and a new improvement'
1238	!! J. Atmospheric Sciences, VOL 37, p. 630--643, 1980.
1239	!!
1240	!! FLOWCHART : None
1241	!_ =============================================================================
1242
1243	FUNCTION dhrT1(rl,tl,gamma1,gamma2,gamma3,gamma4,tta,tau,AlbBS,Tdif,AbsVgt)
1244
1245
1246	!! 0. Variables and parameter declaration
1247	!! 0.1 Input variables
1248	REAL(r_std), INTENT(IN) :: rl ! the effective reflectance of a single scatterer, between 0 and 1
1249	!! @tex $()$ @endtex
1250	REAL(r_std), INTENT(IN) :: tl ! the effective transmittance of a single scatterer, between 0 and 1
1251	!! @tex $()$ @endtex
1252	REAL(r_std), INTENT(IN) :: gamma1 ! the gamma coefficients from Meador and Weaver
1253	REAL(r_std), INTENT(IN) :: gamma2
1254	REAL(r_std), INTENT(IN) :: gamma3
1255	REAL(r_std), INTENT(IN) :: gamma4
1256	REAL(r_std), INTENT(IN) :: tta ! the solar zenith angle, between 0 and pi/2
1257	!! @tex $(radians)$ @endtex
1258
1259	REAL(r_std), INTENT(IN) :: tau ! Effective LAI * G(theta)
1260	!! @tex $()$ @endtex
1261
1262	!! 0.2 Output variables
1263	REAL(r_std), INTENT(OUT) :: AlbBS ! the albedo of level, between 0 and 1
1264	!! @tex $()$ @endtex
1265	REAL(r_std), INTENT(OUT) :: Tdif ! the transmitted light (all diffuse) from this light source, between 0 and 1
1266	!! @tex $()$ @endtex
1267	REAL(r_std), INTENT(OUT) :: AbsVgt ! the light absorbed by the level, between 0 and 1
1268	!! @tex $()$ @endtex
1269	LOGICAL :: dhrT1 ! the output flag, which always appears to be true
1270
1271	!! 0.3 Modified variables
1272	!! 0.4 Local variables
1273
1274	REAL(r_std) :: alpha1,alpha2,ksquare,k
1275	REAL(r_std) :: first_term,secnd_term1,secnd_term2,secnd_term3
1276	REAL(r_std) :: expktau,Tdir
1277	REAL(r_std) :: mu0,w0
1278
1279	!_ ================================================================================================================================
1280
1281	mu0=cos(tta)
1282	w0=rl+tl
1283
1284
1285	Tdir = exp(-tau/mu0) ! direct transmission
1286
1287	! There is a difference between conservative and non-conservative scattering
1288	! conditions */
1289	IF (w0 .ne. 1.0 .AND. w0 .ne. 0.0) THEN
1290	!NON_CONSERVATIVE SCATTERING
1291
1292	! some additional parameters, taken from Eq. 16--18 of Meador and Weaver, J. Atmospheric Sciences,
1293	! Vol 37, p. 630--643 (1980)
1294	alpha1 = gamma1gamma4 + gamma2gamma3
1295	alpha2 = gamma1gamma3 + gamma2gamma4
1296	ksquare = gamma1gamma1 - gamma2gamma2
1297	k = sqrt(ksquare)
1298
1299	expktau = exp(k*tau)
1300
1301	!Black Soil Albedo...Eq. 14 in Meador and Weaver
1302	first_term = ((un-ksquaremu0mu0)((k+gamma1)expktau + (k-gamma1)/expktau))
1303	IF (first_term .eq. 0.0) THEN
1304	!we will be dividing by zero : cannot continue.
1305	dhrT1 = .false.
1306	ELSE
1307	first_term = un/first_term
1308	secnd_term1 = (un - kmu0)(alpha2 + kgamma3)expktau
1309	secnd_term2 = (un + kmu0)(alpha2 - k*gamma3)/expktau
1310	secnd_term3 = 2.0_r_std * k * (gamma3 - alpha2mu0)Tdir
1311	AlbBS = (w0 * first_term * (secnd_term1 - secnd_term2 - secnd_term3))
1312
1313	!Transmission...Eq. 15 in Meador and Weaver, for diffuse light?
1314	IF (ksquare .eq. 0.0) THEN
1315	first_term = un
1316	ENDIF
1317	secnd_term1 = (un+kmu0)(alpha1+kgamma4)expktau
1318	secnd_term2 = (un-kmu0)(alpha1-k*gamma4)/expktau
1319	secnd_term3 = 2.0_r_std * k * (gamma4 + alpha1*mu0)
1320	Tdif = - w0first_term(Tdir*(secnd_term1 - secnd_term2) - secnd_term3)
1321
1322	! Absorption by vegetation...whatever is not transmitted or reflected
1323	! must be absorbed
1324	AbsVgt = (un- (Tdif+Tdir) - AlbBS)
1325	ENDIF ! first_term .eq. 0.0
1326	ELSE IF (w0 .eq. 0.) THEN
1327	!BLACK CANOPY
1328	AlbBS = zero
1329	Tdif = zero
1330	AbsVgt = un - Tdir
1331	ELSE
1332	!CONSERATIVE SCATTERING...Eq. 24 in Meador and Weaver
1333	AlbBS = (un/(un + gamma1tau))(gamma1tau + (gamma3-gamma1mu0)*&
1334	(un-exp(-tau/mu0)));
1335	Tdif = un - AlbBS - Tdir;
1336	AbsVgt = zero;
1337	ENDIF ! w0 .ne. 1.0 .AND. w0 .ne. 0.0
1338
1339	! not sure what the purpose of this flag is, as it will always be set to true here
1340	dhrT1 = .true.
1341
1342	END FUNCTION dhrT1
1343
1344	!! ============================================================================\n
1345	!! FUNCTION : TBarreUncoll_exact
1346	!!
1347	!>\BRIEF Computes the transmission of the diffuse black canopy light
1348	!!
1349	!! DESCRIPTION : Taken directly from the Fortan code provided by Pinty et al for their \n
1350	!! two-stream model. The only changes are the REAL types. This appears
1351	!! to be solving Eq. 16 in Pinty 2006, which is the transmission which
1352	!! does not collide with the canopy
1353	!!
1354	!! RECENT CHANGE(S): None\n
1355	!!
1356	!! RETURN VALUE : TBarreUncoll_exact
1357	!!
1358	!! REFERENCE(S) : None
1359	!!
1360	!! FLOWCHART : None
1361	!_ =============================================================================
1362
1363	FUNCTION TBarreUncoll_exact(tau)
1364
1365
1366	!! 0. Variables and parameter declaration
1367
1368	!! 0.1 Input variables
1369	REAL(r_std), INTENT(IN) :: tau ! Effective LAI * G(theta)
1370	!! @tex $()$ @endtex
1371
1372	!! 0.2 Output variables
1373	REAL(r_std) :: TBarreUncoll_exact ! the isotropic light transmission uncollided
1374	! with the canopy, between 0 and 1
1375	!! @tex $()$ @endtex
1376
1377
1378	!! 0.3 Modified variables
1379
1380	!! 0.4 Local variables
1381	! INTEGER :: j,ind
1382	! REAL(r_std) :: iGammaloc
1383	! INTEGER :: order
1384
1385	!_ ================================================================================================================================
1386
1387	!+++++ CHECK +++++
1388
1389	!!$ iGammaloc = zero
1390	!!$ order=20
1391	!!$
1392	!!$ ! in the case where tau is equal to zero, this crashes in the first loop where
1393	!!$ ! iGammaloc is also zero...quick fix, and might even be accurate
1394	!!$
1395	!!$ IF(tau .GT. 1e-10_r_std)THEN
1396	!!$
1397	!!$ DO j=0,order-1
1398	!!$ ind = order - j
1399	!!$ iGammaloc = ind / (un + ind/(tau+iGammaloc))
1400	!!$ END DO
1401	!!$ iGammaloc=un/(tau + iGammaloc)
1402	!!$ ENDIF
1403	!!$
1404	!!$ TBarreUncoll_exact = exp(-tau)(un - tau + tautau*iGammaloc)
1405
1406	! this is a change suggested by Bernard to improve the matching between the one
1407	! level and multilevel case
1408
1409	TBarreUncoll_exact = exp(-tau2.0_r_std0.705_r_std)
1410	!++++++++++
1411
1412	END FUNCTION TBarreUncoll_exact
1413
1414	!! ==============================================================================================================================
1415	!! SUBROUTINE : gammas
1416	!!
1417	!>\BRIEF Computes a set of gamma coefficients for use in the black background equations
1418	!!
1419	!! DESCRIPTION : Taken directly from the Fortan code provided by Pinty et al for their \n
1420	!! two-stream model. The only changes are the REAL types. This calculates
1421	!! the gamma coefficients based on Appendix A in the paper. This seems to
1422	!! make the assumption of the Ross's G function and a spherical leaf
1423	!! angle distribution.
1424	!!
1425	!! RECENT CHANGE(S): None\n
1426	!!
1427	!! MAIN OUTPUT VARIABLE(S): gammaCoeff
1428	!!
1429	!! REFERENCE(S) : B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski, N. Gobron and M. M. Verstraete (2006).
1430	!! Simplifying the Interaction of Land Surfaces with Radiation for Relating Remote Sensing Products to Climate Models.
1431	!! Journal of Geophysical Research. Vol 111, D02116.
1432	!!
1433	!! FLOWCHART : None
1434	!_ ================================================================================================================================
1435
1436	SUBROUTINE gammas(rl, tl, mu0, gammaCoeff)
1437
1438
1439	!! 0. Variables and parameter declaration
1440	!! 0.1 Input variables
1441	REAL(r_std), INTENT(IN) :: rl ! the effective reflectance of a single scatterer, between 0 and 1
1442	!! @tex $()$ @endtex
1443	REAL(r_std), INTENT(IN) :: tl ! the effective transmittance of a single scatterer, between 0 and 1
1444	!! @tex $()$ @endtex
1445	REAL(r_std), INTENT(IN) :: mu0 ! the cosine of the solar zenith angle, between 0 and 1
1446	!! @tex $()$ @endtex
1447
1448	!! 0.2 Output variables
1449	REAL(r_std), DIMENSION(1:4), INTENT(OUT) :: gammaCoeff ! the set of gamma coefficients in the above reference
1450	!! @tex $()$ @endtex
1451
1452	!! 0.3 Modified variables
1453
1454	!! 0.4 Local variables
1455	REAL(r_std) :: wd,w0,w0half,wdsixth
1456
1457	!_ ================================================================================================================================
1458
1459
1460	w0 = rl+tl
1461	wd = rl-tl
1462	w0half = w0*0.5_r_std
1463	wdsixth = wd/6.0_r_std
1464
1465	gammaCoeff(1)=2._r_std*(un - w0half + wdsixth)
1466	gammaCoeff(2)=2._r_std*(w0half + wdsixth)
1467	gammaCoeff(3)=2._r_std(w0half0.5_r_std + mu0*wdsixth)/w0
1468	gammaCoeff(4)=un-gammaCoeff(3)
1469
1470	END SUBROUTINE gammas
1471
1472
1473	!! ==============================================================================================================================
1474	!! SUBROUTINE : calculate_snow_albedo
1475	!!
1476	!>\BRIEF Computes some of the information needed to calculate the effect of the snow albedo
1477	!! on the background reflectance in the two stream model. Right now, this can be done
1478	!! with the old snow albedo scheme from Krinner et al 2005 as well as the snow albedo
1479	!! scheme from Community Land Model 3 (CLM3). The calculation of
1480	!! snow cover fraction is taken from Yang et al. 1997
1481	!!
1482	!! DESCRIPTION : In order to compute the background albedo in Pinty's two stream model, we have to take
1483	!! into account any snow that has fallen through the canopy and landed on the ground. In particular,
1484	!! we need the amount of ground covered by the snow and the albedo of this snow. Both of these
1485	!! quantities are calculated here, using a function that changes the albedo of the snow based on
1486	!! its age. This is done in two ways, but the model of CLM3 is better suited to be used with Pinty's
1487	!! albedo scheme, as it divides the radiation into VIS and NIR light. The calculation of snow cover
1488	!! fraction depends on roughness lenght.
1489	!!
1490	!! RECENT CHANGE(S): None\n
1491	!!
1492	!! MAIN OUTPUT VARIABLE(S): ::snowa_veg, ::frac_snow_veg, ::albedo_snow
1493	!!
1494	!! REFERENCE(S) : "Technical description of the community land model CLM)", K. W. Oleson, Y. Dai, G. Bonan, M. Bosilovich, et al.,
1495	!! NCAR Technical Note, May 2004.
1496	!!
1497	!! Yang Z, Dickinson R, Robock A, Vinnikov K (1997) Validation of the snow submodel of the
1498	!! biosphere-atmosphere transfer scheme with Russian
1499	!! snow cover and meteorological observational data. Journal of Climate, 10, 353â373.
1500	!!
1501	!! FLOWCHART : None
1502	!_ ================================================================================================================================
1503
1504	SUBROUTINE calculate_snow_albedo(npts, sinang, snow, snow_nobio, &
1505	snow_age, snow_nobio_age, frac_nobio, albedo_snow, &
1506	snowa_veg, frac_snow_veg, snowa_nobio, frac_snow_nobio, &
1507	veget_max, z0_veg)
1508
1509	!! 0. Variables and parameter declaration
1510	!! 0.1 Input variables
1511	INTEGER,INTENT(IN) :: npts !! Domain size - Number of land pixels (unitless)
1512	REAL(r_std), DIMENSION(npts), INTENT(in) :: sinang !! cosine of the solar zenith angle, between 0 and 1
1513	!! @tex $()$ @endtex
1514	REAL(r_std),DIMENSION (npts), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
1515	REAL(r_std),DIMENSION (npts,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
1516	REAL(r_std),DIMENSION (npts), INTENT(in) :: snow_age !! Snow age (days)
1517	REAL(r_std),DIMENSION (npts,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
1518	REAL(r_std),DIMENSION (npts,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
1519	!! continental ice, lakes, etc. (unitless)
1520	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! PFT coverage fraction of a PFT (= ind*cn_ind)
1521	!! (m^2 m^{-2})
1522	REAL(r_std), DIMENSION(npts), INTENT(in) :: z0_veg !! Roughness height of vegetated part (m)
1523
1524
1525	!! 0.2 Output variables
1526	REAL(r_std), DIMENSION(npts),INTENT(OUT) :: albedo_snow !! the averaged snow albedo for a given grid cell
1527	REAL(r_std), DIMENSION(npts,nvm,n_spectralbands),INTENT(OUT) :: snowa_veg !! the albedo of snow...this is seperated into
1528	!! PFT types, even though the calculation is identical
1529	!! for all PFTs right now
1530	REAL(r_std), DIMENSION(npts,nvm),INTENT(OUT) :: frac_snow_veg !! the fraction of the surface for each PFT covered by snow
1531	REAL(r_std), DIMENSION(npts,nnobio,n_spectralbands),INTENT(OUT) :: snowa_nobio !! the albedo of snow on nonbiological land types
1532	REAL(r_std), DIMENSION(npts,nnobio),INTENT(OUT) :: frac_snow_nobio !! the fraction of nonbiological land types covered by snow
1533
1534
1535	!! 0.3 Modified variables
1536
1537	!! 0.4 Local variables
1538	REAL(r_std) :: agefunc,agefunc_nobio,snow_alb_direct, snow_alb_diffuse,f_mu
1539	REAL(r_std),DIMENSION(n_spectralbands) :: c_albedo,alb_snow_0 ! parameter values, taken directly from CLM3
1540
1541	INTEGER :: ipts,ks,ivm,jv ! indices
1542
1543	!_ ================================================================================================================================
1544
1545
1546	! initialize the output
1547	albedo_snow=val_exp
1548	snowa_veg=val_exp
1549	frac_snow_veg=val_exp
1550	snowa_nobio=val_exp
1551	frac_snow_nobio=val_exp
1552
1553	! set some values for the new snow albedo
1554	alb_snow_0(ivis)=0.95_r_std
1555	alb_snow_0(inir)=0.65_r_std
1556	c_albedo(ivis)=0.2_r_std
1557	c_albedo(inir)=0.5_r_std
1558
1559
1560	DO ipts = 1, npts ! loop over all the grid squares
1561
1562
1563	! calculate the snow cover fraction...right now this doesn't depend on PFT
1564	! snow fraction calculated following Yang et al 1997, eq. 32
1565	frac_snow_veg(ipts,:) = TANH((snow(ipts)/sn_dens)/ &
1566	MAX((2.5_r_std*z0_veg(ipts)), &
1567	min_sechiba))
1568
1569	DO ivm=1,nvm ! loop over all the PFTs
1570
1571
1572	IF(veget_max(ipts,ivm) == zero)THEN
1573	! this vegetation type is not present, so no reason to do the
1574	! calculation
1575	CYCLE
1576	ENDIF
1577
1578
1579	DO ks=1,n_spectralbands ! loop over spectra
1580
1581	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. &
1582	EPSILON(undef_sechiba)) THEN
1583	snowa_veg(ipts,ivm,ks) = fixed_snow_albedo
1584	ELSE
1585
1586	IF(control%do_new_snow_albedo)THEN
1587	! first we calculate diffuse albedo
1588	! NOTE: The difference between these two methods has not been
1589	! tested. The only constraint is that it must be dimensionless
1590	agefunc = un - EXP(-snow_age(ipts)/tcst_snowa) ! already in ORCHIDEE
1591	!agefunc = zero
1592
1593	! using the constant from CLM3 conveted into days**-1
1594	!agefunc=un-un/(un+snow_age(ipts)*0.864_r_std)
1595	!agefunc=un/(un+snow_age(ipts)*0.864_r_std)
1596	snow_alb_diffuse=(un-c_albedo(ks)agefunc)alb_snow_0(ks)
1597
1598	! now the direct albedo
1599	IF(sinang(ipts) .GT. 0.5_r_std)THEN
1600	f_mu=zero
1601	ELSE
1602	! substituting b=2.0 into the equation...recommended by CLM3
1603	f_mu=(3.0_r_std)/(un+sinang(ipts))-2.0_r_std
1604	ENDIF
1605
1606	snow_alb_direct=snow_alb_diffuse+0.4_r_stdf_mu&
1607	(un-snow_alb_diffuse)
1608
1609	! assume that the total snow albedo is a mix of
1610	! 50% direct and 50% diffuse
1611	snowa_veg(ipts,ivm,ks)=(snow_alb_direct+snow_alb_diffuse)&
1612	/2.0_r_std
1613
1614	ELSE
1615	! calculate the age of the snow for this grid square
1616	agefunc = EXP(-snow_age(ipts)/tcst_snowa)
1617
1618	! now the albedo...use only the bare soil snow decay parameters,
1619	! since the true albedo is calculated
1620	! by the PFT-dependent two-stream solver
1621	snowa_veg(ipts,ivm,ks) = snowa_aged(1)+snowa_dec(1)*agefunc
1622
1623	ENDIF ! control%do_new_snow_albedo
1624	ENDIF ! prescribe or calculate albedo
1625
1626
1627
1628	! Notice that the fraction of nobio land is included in veget_max.
1629	! We have to average by the number of spectras, too, for this array,
1630	! because it doesn't distinguish between PFTs or spectra
1631	albedo_snow(ipts) = albedo_snow(ipts) + &
1632	frac_snow_veg(ipts,ivm) * veget_max(ipts,ivm) * &
1633	snowa_veg(ipts,ivm,ks)/REAL(n_spectralbands,r_std)
1634
1635	ENDDO ! ks=1,n_spectralbands
1636
1637	ENDDO ! ivm=1,nvm
1638
1639	! there is no distinction between NIR and VIS for non bio surfaces
1640	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
1641	snowa_nobio(:,:,:) = fixed_snow_albedo
1642	ELSE
1643
1644	DO jv = 1, nnobio
1645	! calculate the snow age
1646	agefunc_nobio = EXP(-snow_nobio_age(ipts,jv)/tcst_snowa)
1647
1648	! and the albedo due to snow of this age on this nonbio surface
1649	snowa_nobio(ipts,jv,:) = ( snowa_aged(1) + snowa_dec(1) * agefunc_nobio )
1650	ENDDO
1651	ENDIF ! prescribe or calculate
1652
1653	DO jv = 1, nnobio
1654	! calculate the fraction of the surface covered by snow
1655	! snow fraction calculated following Yang et al 1997, eq. 32
1656	frac_snow_nobio(ipts,jv) = TANH((snow_nobio(ipts,jv)/sn_dens)/ &
1657	MAX((2.5_r_std*z0_ice), &
1658	min_sechiba))
1659
1660	! For some reason, occasionally snow_nobio is negative. This is driver
1661	! data that we cannot control, but it should not influence our snowcover, so
1662	! set it equal to zero in that case.
1663	IF(snow_nobio(ipts,jv) .LT. zero) frac_snow_nobio(ipts,jv)=zero
1664
1665	albedo_snow(ipts) = albedo_snow(ipts) + &
1666	( frac_nobio(ipts,jv) ) * ( frac_snow_nobio(ipts,jv) ) * &
1667	snowa_nobio(ipts,jv,1) ! not sensitive to the spectral band.
1668	! Don't need to average here since we are outside the loop over the bands
1669
1670	ENDDO ! jv = 1, nnobio
1671
1672	ENDDO ! ipts = 1, npts
1673
1674
1675
1676
1677	END SUBROUTINE CALCULATE_SNOW_ALBEDO
1678
1679	!! ==============================================================================================================================
1680	!! SUBROUTINE : optimize_albedo_values
1681	!!
1682	!>\BRIEF Follow the radiation scattered through the canopy in the
1683	!! case of multiple levels.
1684	!!
1685	!! DESCRIPTION : Right now, we know that using Pintys method in the case of a
1686	!! single canopy level will result in different top of the canopy
1687	!! albedos than iwhat is found if we use multiple canopy levels.
1688	!! We trust that the single level case gives the 'true' values.
1689	!!
1690	!! This algorithm follows light as it scatters through multiple
1691	!! levels in the canopy. At each level, the scattering
1692	!! solution is found by solving Pinty's two stream model. This
1693	!! means that light which enters a level can either be transmitted
1694	!! without colliding with the vegetation, transmitted after
1695	!! collision with the vegetation, reflecting off the vegetation,
1696	!! or being absorbed. We follow the fluxes until they get
1697	!! really small.
1698	!!
1699	!! RECENT CHANGE(S): None\n
1700	!!
1701	!! MAIN OUTPUT VARIABLE(S): ::lconverged, ::Collim_Alb_Tot, ::Collim_Tran_Tot,
1702	!! ::Collim_Abs_Tot, ::Isotrop_Alb_Tot, ::Isotrop_Tran_Tot, ::Isotrop_Abs_Tot
1703	!!
1704	!! REFERENCE(S) : B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski, N. Gobron
1705	!! and M. M. Verstraete (2006). Simplifying the Interaction of Land
1706	!! Surfaces with Radiation for Relating Remote Sensing Products to
1707	!! Climate Models. Journal of Geophysical Research. Vol 111, D02116.
1708	!!
1709	!! FLOWCHART : None
1710	!_ ================================================================================================================================
1711	SUBROUTINE multilevel_albedo(nlevels_loc, cosine_sun_angle, leaf_single_scattering_albedo_start,&
1712	leaf_psd_start, br_base_temp_collim, br_base_temp_isotrop, &
1713	laieff_collim_pft, laieff_isotrop_pft, lconverged, &
1714	Collim_Alb_Coll, Collim_Tran_Coll, Collim_Abs_Tot, Isotrop_Alb_Coll, &
1715	Isotrop_Tran_Coll, Isotrop_Abs_Tot, Collim_Tran_Uncoll, Isotrop_Tran_Uncoll, lprint)
1716
1717	!! 0. Variables and parameter declaration
1718	!! 0.1 Input variables
1719	INTEGER,INTENT(IN) :: nlevels_loc
1720	REAL(r_std), INTENT(IN) :: cosine_sun_angle !! cosine of the solar zenith angle, between 0 and 1
1721	!! @tex $()$ @endtex
1722	REAL(r_std), INTENT(IN) :: leaf_single_scattering_albedo_start !! cosine of the solar zenith angle, between 0 and 1
1723	REAL(r_std), INTENT(IN) :: leaf_psd_start !! cosine of the solar zenith angle, between 0 and 1
1724	REAL(r_std), INTENT(IN) :: br_base_temp_collim !! cosine of the solar zenith angle, between 0 and 1
1725	REAL(r_std), INTENT(IN) :: br_base_temp_isotrop !! cosine of the solar zenith angle, between 0 and 1
1726	REAL(r_std),DIMENSION(:),INTENT(IN) :: laieff_collim_pft !! Effective lai for a single pixel and pft to be
1727	REAL(r_std),DIMENSION(:),INTENT(IN) :: laieff_isotrop_pft !! Effective lai for a single pixel and pft to be
1728	LOGICAL,INTENT(IN) :: lprint !! A flag to print some debug statements
1729	!! used in the two stream approach...every level
1730
1731	!! 0.2 Output variables
1732	LOGICAL,INTENT(OUT) :: lconverged !! did the optimization converge?
1733	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Collim_Alb_Coll !! collimated total albedo for the converged case
1734	!! unitless, between 0 and 1
1735	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Collim_Tran_Coll !! collimated total transmission
1736	!! unitless, between 0 and 1
1737	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Collim_Abs_Tot !! collimated total absorption
1738	!! unitless, between 0 and 1
1739	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Isotrop_Alb_Coll !! isotropic total albedo
1740	!! unitless, between 0 and 1
1741	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Isotrop_Tran_Coll !! isotropic total transmission
1742	!! unitless, between 0 and 1
1743	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Isotrop_Abs_Tot !! isotropic total absorption
1744	!! unitless, between 0 and 1
1745	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Collim_Tran_Uncoll !! collimated uncollided transmission
1746	!! unitless, between 0 and 1
1747	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Isotrop_Tran_Uncoll !! isotropic uncollided transmission
1748	!! unitless, between 0 and 1
1749
1750	!! 0.3 Modified variables
1751
1752	!! 0.4 Local variables
1753	! use an extra level here...this is basically to store the transmission from the sun, which is normalized to 1.0
1754	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Abs_Coll_Unscaled
1755	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Alb_Coll_Unscaled
1756	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Tran_Coll_Unscaled
1757	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Tran_UnColl_Unscaled
1758	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Abs_Coll_Unscaled
1759	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Alb_Coll_Unscaled
1760	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Tran_Coll_Unscaled
1761	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Tran_UnColl_Unscaled
1762	REAL(r_std),DIMENSION(nlevels_loc) :: Collim_Tran_Tot
1763	REAL(r_std),DIMENSION(nlevels_loc) :: Isotrop_Tran_Tot
1764
1765	REAL(r_std),DIMENSION(0:nlevels_loc+1) :: &
1766	collim_down_cs,isotrop_down_cs,isotrop_up_cs
1767
1768	INTEGER :: istep,ilevel
1769
1770	REAL(r_std) :: leaf_reflectance
1771	REAL(r_std) :: leaf_transmittance
1772
1773	LOGICAL :: lexit
1774	LOGICAL :: ok
1775	REAL(r_std), DIMENSION(4) :: gammaCoeffs
1776	REAL(r_std), DIMENSION(4) :: gammaCoeffs_star
1777	REAL(r_std) :: sun_zenith_angle_radians
1778	REAL(r_std), PARAMETER :: isotropic_cosine_constant=0.5/0.705
1779
1780	INTEGER,SAVE :: icalls=0
1781	!_ ================================================================================================================================
1782
1783
1784	icalls=icalls+1
1785
1786	! convet the sun angle
1787	sun_zenith_angle_radians = acos(cosine_sun_angle)
1788
1789
1790
1791	! initialize some values that are used in the optimization loop
1792	istep=0
1793	lexit=.FALSE.
1794	lconverged=.FALSE.
1795
1796	! calculate the albedo at our starting point
1797
1798	leaf_transmittance = leaf_single_scattering_albedo_start/ &
1799	( leaf_psd_start+un)
1800	leaf_reflectance = leaf_psd_start * &
1801	leaf_transmittance
1802
1803	! some debugging stuff
1804	! DO ilevel=1,nlevels_loc
1805
1806	! CALL print_debugging_albedo_info(ilevel,cosine_sun_angle,leaf_reflectance,&
1807	! leaf_transmittance,&
1808	! laieff_collim_pft(ilevel),laieff_isotrop_pft(ilevel), &
1809	! reflectance_collim(ilevel-1),reflectance_isotrop(ilevel-1),&
1810	! Collim_Alb_Tot_temp(ilevel),Collim_Tran_Tot_temp(ilevel),&
1811	! Collim_Abs_Tot_temp(ilevel),&
1812	! Isotrop_Alb_Tot_temp(ilevel),Isotrop_Tran_Tot_temp(ilevel),&
1813	! Isotrop_Abs_Tot_temp(ilevel))
1814	! ENDDO
1815
1816
1817	! calculate the gamma coefficients used in the case of the black background
1818	call gammas(leaf_reflectance, leaf_transmittance, cosine_sun_angle, gammaCoeffs)
1819	call gammas(leaf_reflectance,leaf_transmittance,isotropic_cosine_constant,&
1820	gammaCoeffs_star)
1821
1822
1823
1824	!************************ step one **** !
1825	! compute all the unscaled quantities
1826	DO ilevel=nlevels_loc,1,-1
1827
1828	Collim_Tran_UnColl_Unscaled(ilevel)=exp( - 0.5_r_std*&
1829	laieff_collim_pft(ilevel)/cosine_sun_angle)
1830
1831	Isotrop_Tran_UnColl_Unscaled(ilevel)= &
1832	TBarreUncoll_exact(0.5_r_std*laieff_isotrop_pft(ilevel))
1833
1834
1835	! /* 1) collimated source */
1836	ok=dhrT1(leaf_reflectance,leaf_transmittance,&
1837	gammaCoeffs(1),gammaCoeffs(2),gammaCoeffs(3),gammaCoeffs(4),&
1838	sun_zenith_angle_radians, 0.5_r_std*laieff_collim_pft(ilevel),&
1839	Collim_Alb_Coll_Unscaled(ilevel),Collim_Tran_Coll_Unscaled(ilevel),&
1840	Collim_Abs_Coll_Unscaled(ilevel))
1841
1842	! /* 2) isotropic source */
1843	ok=dhrT1(leaf_reflectance,leaf_transmittance,&
1844	gammaCoeffs_star(1),gammaCoeffs_star(2),gammaCoeffs_star(3),gammaCoeffs_star(4),&
1845	acos(isotropic_cosine_constant),0.5_r_std*laieff_isotrop_pft(ilevel),&
1846	Isotrop_Alb_Coll_Unscaled(ilevel),Isotrop_Tran_Coll_Unscaled(ilevel),&
1847	Isotrop_Abs_Coll_Unscaled(ilevel))
1848
1849	! Collim_Tran_OneWay_Unscaled(ilevel)=Collim_Tran_Coll_Unscaled(ilevel)+&
1850	! Collim_Tran_UnColl_Unscaled(ilevel)
1851	! Isotrop_Tran_OneWay_Unscaled(ilevel)=Isotrop_Tran_Coll_Unscaled(ilevel)+&
1852	! Isotrop_Tran_UnColl_Unscaled(ilevel)
1853
1854	! WRITE(6,'(A,I5,3F10.6)') 'unscaled ',ilevel,Collim_Tran_UnColl_Unscaled(ilevel),&
1855	! Collim_Tran_Coll_Unscaled(ilevel),&
1856	! Collim_Alb_Coll_Unscaled(ilevel)
1857	ENDDO ! ilevel=nlevels_loc,1,-1
1858
1859
1860
1861
1862	! Following the fate of the light at every step
1863
1864	! The downwelling array indicates the quantity of light flowing into this level
1865	! from above therefore, to start our system for collimated light, we give a
1866	! unit of light coming into the top level from a collimated source.
1867	! The upwelling array is light entering this level from below.
1868	! cs is the current step, while ns is the next step for the iteration.
1869	! Level 0 is the background. Light can enter this level from above but not below
1870	! Level nlevels_loc+1 is the atomosphere. Light can enter this level from below
1871	! but not above.
1872	collim_down_cs(:)=zero
1873	collim_down_cs(nlevels_loc)=un
1874	isotrop_down_cs(:)=zero
1875	isotrop_up_cs(:)=zero
1876
1877	CALL propagate_fluxes(nlevels_loc, collim_down_cs, isotrop_down_cs, isotrop_up_cs, &
1878	Collim_Tran_UnColl_Unscaled, Collim_Tran_Coll_Unscaled, &
1879	Collim_Alb_Coll_Unscaled, Isotrop_Tran_UnColl_Unscaled, &
1880	Isotrop_Tran_Coll_Unscaled, Isotrop_Alb_Coll_Unscaled, &
1881	br_base_temp_collim, br_base_temp_isotrop, .FALSE., &
1882	Collim_Tran_Uncoll, Collim_Tran_Coll, Collim_Alb_Coll, lconverged, lprint)
1883
1884	! Now for the isotropic light
1885	collim_down_cs(:)=zero
1886	isotrop_down_cs(:)=zero
1887	isotrop_down_cs(nlevels_loc)=un
1888	isotrop_up_cs(:)=zero
1889
1890	! The colliminated light is never used here, but it's passed to keep things
1891	! clean between the two cases. Might have to redo this if we are having
1892	! peformance issues.
1893	CALL propagate_fluxes(nlevels_loc, collim_down_cs, isotrop_down_cs, isotrop_up_cs, &
1894	Collim_Tran_UnColl_Unscaled, Collim_Tran_Coll_Unscaled, &
1895	Collim_Alb_Coll_Unscaled, Isotrop_Tran_UnColl_Unscaled, &
1896	Isotrop_Tran_Coll_Unscaled, Isotrop_Alb_Coll_Unscaled, &
1897	br_base_temp_collim, br_base_temp_isotrop, .TRUE., &
1898	Isotrop_Tran_Uncoll, Isotrop_Tran_Coll, Isotrop_Alb_Coll, lconverged, lprint)
1899
1900	! Calculate the absorption profile
1901	Collim_Tran_Tot(:)=Collim_Tran_Coll(:)+Collim_Tran_Uncoll(:)
1902	Isotrop_Tran_Tot(:)=Isotrop_Tran_Coll(:)+Isotrop_Tran_Uncoll(:)
1903
1904	! bottom level
1905	Collim_Abs_Tot(1)=Collim_Tran_Tot(1+1)+Collim_Tran_Tot(1)*br_base_temp_collim&
1906	-Collim_Tran_Tot(1)-Collim_Alb_Coll(1)
1907	Isotrop_Abs_Tot(1)=Isotrop_Tran_Tot(1+1)+Isotrop_Tran_Tot(1)*br_base_temp_isotrop&
1908	-Isotrop_Tran_Tot(1)-Isotrop_Alb_Coll(1)
1909
1910	! all middle levels
1911	DO ilevel=2,nlevels_loc-1
1912	Collim_Abs_Tot(ilevel)=Collim_Tran_Tot(ilevel+1)+Collim_Alb_Coll(ilevel-1)&
1913	-Collim_Tran_Tot(ilevel)-Collim_Alb_Coll(ilevel)
1914	Isotrop_Abs_Tot(ilevel)=Isotrop_Tran_Tot(ilevel+1)+Isotrop_Alb_Coll(ilevel-1)&
1915	-Isotrop_Tran_Tot(ilevel)-Isotrop_Alb_Coll(ilevel)
1916	ENDDO
1917
1918	! top level
1919	Collim_Abs_Tot(nlevels_loc)=un+Collim_Alb_Coll(nlevels_loc-1)&
1920	-Collim_Tran_Tot(nlevels_loc)-Collim_Alb_Coll(nlevels_loc)
1921	Isotrop_Abs_Tot(nlevels_loc)=un+Isotrop_Alb_Coll(nlevels_loc-1)&
1922	-Isotrop_Tran_Tot(nlevels_loc)-Isotrop_Alb_Coll(nlevels_loc)
1923
1924
1925
1926	END SUBROUTINE multilevel_albedo
1927
1928	!! ==============================================================================================================================
1929	!! SUBROUTINE : print_debugging_albedo_info
1930	!!
1931	!>\BRIEF Prints out some albedo information in a nice format.
1932	!! Should only be used for debugging, never for production runs.
1933	!!
1934	!! DESCRIPTION :
1935	!!
1936	!! RECENT CHANGE(S): None\n
1937	!!
1938	!! MAIN OUTPUT VARIABLE(S): None.
1939	!!
1940	!! REFERENCE(S) :
1941	!!
1942	!! FLOWCHART : None
1943	!_ ================================================================================================================================
1944	SUBROUTINE print_debugging_albedo_info(ilevel, cosine_sun_angle, &
1945	leaf_reflectance, leaf_transmittance, laieff_collim_temp, laieff_isotrop_temp, &
1946	br_base_temp_collim, br_base_temp_isotrop, Collim_Alb_Tot, &
1947	Collim_Tran_Tot, Collim_Abs_Tot, Isotrop_Alb_Tot, Isotrop_Tran_Tot, Isotrop_Abs_Tot)
1948
1949	!! 0. Variables and parameter declaration
1950	!! 0.1 Input variables
1951	INTEGER,INTENT(IN) :: ilevel
1952	REAL(r_std), INTENT(IN) :: cosine_sun_angle !! cosine of the solar zenith angle, between 0 and 1
1953	!! @tex $()$ @endtex
1954	REAL(r_std), INTENT(IN) :: laieff_collim_temp !! cosine of the solar zenith angle, between 0 and 1
1955	REAL(r_std), INTENT(IN) :: laieff_isotrop_temp !! cosine of the solar zenith angle, between 0 and 1
1956	REAL(r_std), INTENT(IN) :: leaf_reflectance !! cosine of the solar zenith angle, between 0 and 1
1957	REAL(r_std), INTENT(IN) :: leaf_transmittance !! cosine of the solar zenith angle, between 0 and 1
1958	REAL(r_std), INTENT(IN) :: br_base_temp_collim !! cosine of the solar zenith angle, between 0 and 1
1959	REAL(r_std), INTENT(IN) :: br_base_temp_isotrop !! cosine of the solar zenith angle, between 0 and 1
1960	REAL(r_std),INTENT(IN) :: Collim_Alb_Tot
1961	REAL(r_std),INTENT(IN) :: Collim_Tran_Tot
1962	REAL(r_std),INTENT(IN) :: Collim_Abs_Tot
1963	REAL(r_std),INTENT(IN) :: Isotrop_Alb_Tot
1964	REAL(r_std),INTENT(IN) :: Isotrop_Tran_Tot
1965	REAL(r_std),INTENT(IN) :: Isotrop_Abs_Tot
1966
1967	!! 0.2 Output variables
1968
1969	!! 0.3 Modified variables
1970
1971	!! 0.4 Local variables
1972
1973	!_ ================================================================================================================================
1974
1975	1223 FORMAT(I5,I3,7(F14.7))
1976
1977	WRITE (6,'(8(11A))') ' Level ',' Angle ',' laieff_c ',' laieff_i ',&
1978	' rleaf ',' tleaf ',' rbgd_c ',' rbgd_i '
1979
1980	WRITE(6,FMT='(I6,5X,7(F11.6))') &
1981	ilevel,180/pi*ACOS(cosine_sun_angle),&
1982	laieff_collim_temp,laieff_isotrop_temp,&
1983	leaf_reflectance,leaf_transmittance,br_base_temp_collim,&
1984	br_base_temp_isotrop
1985
1986
1987	WRITE (6,'(10X,6(11A))') ' Rtot(sun) ',' Ttot(sun) ',&
1988	' Atot(sun) ',' Rtot(iso) ',' Ttot(iso) ',' Atot(iso) '
1989	WRITE(6,FMT='(10X,6(F11.6))') &
1990	Collim_Alb_Tot,Collim_Tran_Tot,Collim_Abs_Tot,&
1991	Isotrop_Alb_Tot,Isotrop_Tran_Tot,Isotrop_Abs_Tot
1992
1993
1994	END SUBROUTINE print_debugging_albedo_info
1995
1996
1997	!! ==============================================================================================================================
1998	!! SUBROUTINE : propagate_fluxes
1999	!!
2000	!>\BRIEF Propogates the radiation fluxes through each level of the
2001	!! canopy.
2002	!!
2003	!! DESCRIPTION : This is an algorithm to follow radition fluxes through the
2004	!! canopy. At each level, the path probabilities are determined by the
2005	!! raditional transfer scheme of Pinty et al (2006). Notice that
2006	!! the fluxes given by this routine are all cumulative fluxes, not per
2007	!! level.
2008	!!
2009	!! RECENT CHANGE(S): None\n
2010	!!
2011	!! MAIN OUTPUT VARIABLE(S): ::Tran_Uncoll_Tot, ::Tran_Coll_Tot, ::Alb_Coll_Tot
2012	!!
2013	!! REFERENCE(S) :
2014	!!
2015	!! FLOWCHART : None
2016	!_ ================================================================================================================================
2017	SUBROUTINE propagate_fluxes(nlevels_loc, collim_down_cs, isotrop_down_cs, isotrop_up_cs, &
2018	Collim_Tran_UnColl_Unscaled, Collim_Tran_Coll_Unscaled, &
2019	Collim_Alb_Coll_Unscaled, Isotrop_Tran_UnColl_Unscaled, &
2020	Isotrop_Tran_Coll_Unscaled, Isotrop_Alb_Coll_Unscaled, br_base_temp_collim, &
2021	br_base_temp_isotrop, lisotrop, Tran_Uncoll_Tot, Tran_Coll_Tot, &
2022	Alb_Coll_Tot, lconverged, lprint)
2023
2024	!! 0. Variables and parameter declaration
2025	!! 0.1 Input variables
2026	INTEGER,INTENT(IN) :: nlevels_loc
2027	REAL(r_std), DIMENSION(nlevels_loc),INTENT(IN) :: Collim_Tran_UnColl_Unscaled
2028	REAL(r_std), DIMENSION(nlevels_loc),INTENT(IN) :: Collim_Tran_Coll_Unscaled
2029	REAL(r_std), DIMENSION(nlevels_loc),INTENT(IN) :: Collim_Alb_Coll_Unscaled
2030	REAL(r_std), DIMENSION(nlevels_loc),INTENT(IN) :: Isotrop_Tran_UnColl_Unscaled
2031	REAL(r_std), DIMENSION(nlevels_loc),INTENT(IN) :: Isotrop_Tran_Coll_Unscaled
2032	REAL(r_std), DIMENSION(nlevels_loc),INTENT(IN) :: Isotrop_Alb_Coll_Unscaled
2033	REAL(r_std), INTENT(IN) :: br_base_temp_collim !! cosine of the solar zenith angle, between 0 and 1
2034	REAL(r_std), INTENT(IN) :: br_base_temp_isotrop !! cosine of the solar zenith angle, between 0 and 1
2035	LOGICAL, INTENT(IN) :: lisotrop !! are we dealing with an isotropic source? only needed for correct partitioning
2036	!! of the collided and uncollided transmitted light...the total is not affected
2037	LOGICAL, INTENT(IN) :: lprint !! a flag to print
2038
2039	!! 0.2 Output variables
2040	REAL(r_std), DIMENSION(nlevels_loc),INTENT(OUT) :: Tran_Uncoll_Tot
2041	REAL(r_std), DIMENSION(nlevels_loc),INTENT(OUT) :: Tran_Coll_Tot
2042	REAL(r_std), DIMENSION(nlevels_loc),INTENT(OUT) :: Alb_Coll_Tot
2043	LOGICAL,INTENT(OUT) :: lconverged
2044
2045	!! 0.3 Modified variables
2046	REAL(r_std), DIMENSION(0:nlevels_loc+1),INTENT(INOUT) :: collim_down_cs
2047	REAL(r_std), DIMENSION(0:nlevels_loc+1),INTENT(INOUT) :: isotrop_down_cs
2048	REAL(r_std), DIMENSION(0:nlevels_loc+1),INTENT(INOUT) :: isotrop_up_cs
2049
2050	!! 0.4 Local variables
2051	INTEGER :: istep,ilevel
2052	REAL(r_std),DIMENSION(0:nlevels_loc+1) :: collim_down_ns,isotrop_down_ns,&
2053	isotrop_up_ns,Tran_Uncoll_Tot_temp
2054
2055	!_ ================================================================================================================================
2056
2057	Tran_Uncoll_Tot(:)=zero
2058	Tran_Coll_Tot(:)=zero
2059	Alb_Coll_Tot(:)=zero
2060
2061	Tran_Uncoll_Tot_temp(:)=un
2062
2063	istep=0
2064	DO
2065
2066	istep=istep+1
2067
2068	lconverged=.TRUE.
2069
2070	! Zero out the counters for the next step
2071	collim_down_ns(:)=zero
2072	isotrop_down_ns(:)=zero
2073	isotrop_up_ns(:)=zero
2074
2075
2076	! Now we need to loop over all levels and see what light is entering the
2077	! level, and how it will propagate in the next step
2078	DO ilevel=1,nlevels_loc
2079
2080	! For collimated downwelling light into the level, it can be scattered
2081	! up, down, or pass through uncollided
2082	IF(collim_down_cs(ilevel) .GT. zero)THEN
2083	collim_down_ns(ilevel-1)=collim_down_ns(ilevel-1)+&
2084	collim_down_cs(ilevel)*Collim_Tran_UnColl_Unscaled(ilevel)
2085
2086	! This statement checks to see if this light has been previously
2087	! scattered or not. This term is only present in level nlevels_loc
2088	! for the first step, level nlevels_loc-1 for the second step,
2089	! level nlevels_loc-2 for the third step, etc., and it only happens
2090	! in the case of a collimated light source.
2091	IF(ilevel == nlevels_loc-istep+1 .AND. .NOT. lisotrop) THEN
2092	Tran_Uncoll_Tot_temp(ilevel)=Tran_Uncoll_Tot_temp(ilevel+1)*&
2093	Collim_Tran_UnColl_Unscaled(ilevel)
2094	ENDIF
2095
2096	isotrop_down_ns(ilevel-1)=isotrop_down_ns(ilevel-1)+&
2097	collim_down_cs(ilevel)*Collim_Tran_Coll_Unscaled(ilevel)
2098	isotrop_up_ns(ilevel+1)=isotrop_up_ns(ilevel+1)+&
2099	collim_down_cs(ilevel)*Collim_Alb_Coll_Unscaled(ilevel)
2100	ENDIF ! collim_down_cs(ilevel) .GT. zero
2101
2102	! For isotropic downwelling light, it can also be scattered up, down,
2103	! or pass through uncollided
2104	IF(isotrop_down_cs(ilevel) .GT. zero)THEN
2105	isotrop_down_ns(ilevel-1)=isotrop_down_ns(ilevel-1)+&
2106	isotrop_down_cs(ilevel)*Isotrop_Tran_UnColl_Unscaled(ilevel)
2107
2108	! This is the same check as above, but this time the light has
2109	! an isotropic source and not a collimated source
2110	IF(ilevel == nlevels_loc-istep+1 .AND. lisotrop) THEN
2111	Tran_Uncoll_Tot_temp(ilevel)=Tran_Uncoll_Tot_temp(ilevel+1)&
2112	*Isotrop_Tran_UnColl_Unscaled(ilevel)
2113	ENDIF
2114
2115	isotrop_down_ns(ilevel-1)=isotrop_down_ns(ilevel-1)+&
2116	isotrop_down_cs(ilevel)*Isotrop_Tran_Coll_Unscaled(ilevel)
2117	isotrop_up_ns(ilevel+1)=isotrop_up_ns(ilevel+1)+&
2118	isotrop_down_cs(ilevel)*Isotrop_Alb_Coll_Unscaled(ilevel)
2119	ENDIF ! isotrop_down_cs(ilevel) .GT. zero
2120
2121	! Isotropic upwelling light can pass through upwards collided or
2122	! uncollided with vegetation in this level, or it can be reflected downwards
2123	IF(isotrop_up_cs(ilevel) .GT. zero)THEN
2124	isotrop_up_ns(ilevel+1)=isotrop_up_ns(ilevel+1)+isotrop_up_cs(ilevel)&
2125	*Isotrop_Tran_UnColl_Unscaled(ilevel)
2126	isotrop_up_ns(ilevel+1)=isotrop_up_ns(ilevel+1)+isotrop_up_cs(ilevel)&
2127	*Isotrop_Tran_Coll_Unscaled(ilevel)
2128	isotrop_down_ns(ilevel-1)=isotrop_down_ns(ilevel-1)+&
2129	isotrop_up_cs(ilevel)*Isotrop_Alb_Coll_Unscaled(ilevel)
2130	ENDIF
2131
2132	! there can be no collimated upwards light, since upwards light must
2133	! always have reflected off something
2134
2135
2136	ENDDO ! ilevel=1,nlevels_loc
2137
2138	! The background is a bit special. there is no transmission, but there is
2139	! reflection, which leads to a source term to the bottom level from below.
2140	isotrop_up_ns(1)=isotrop_up_ns(1)+collim_down_cs(0)*br_base_temp_collim
2141	isotrop_up_ns(1)=isotrop_up_ns(1)+isotrop_down_cs(0)*br_base_temp_isotrop
2142
2143
2144	! Now we add the light generated here to our cumulative counters to track
2145	! the total amount that leaves the canopy, either through being absorbed
2146	! by the background or being reflected back into the atmosphere.
2147	! We keep track of the uncoll above, so here we track the total light that
2148	! is transmitted to the soil, and then at the end we taken the difference
2149	! to get the collided radation.
2150	Tran_Coll_Tot(1:nlevels_loc)=Tran_Coll_Tot(1:nlevels_loc)+&
2151	isotrop_down_ns(0:nlevels_loc-1)+collim_down_ns(0:nlevels_loc-1)
2152	Alb_Coll_Tot(1:nlevels_loc)=Alb_Coll_Tot(1:nlevels_loc)+isotrop_up_ns(2:nlevels_loc+1)
2153
2154	! now we update the light we are currently tracking
2155	collim_down_cs(:)=collim_down_ns(:)
2156	isotrop_down_cs(:)=isotrop_down_ns(:)
2157	isotrop_up_cs(:)=isotrop_up_ns(:)
2158
2159
2160
2161
2162	! check for convegence...if all the values of light are currently below a
2163	! threshold, we're probably good assuming the threshold is low enough.
2164	IF(MAXVAL(collim_down_cs,1) .GT. alb_threshold) lconverged=.FALSE.
2165	IF(MAXVAL(isotrop_down_cs,1) .GT. alb_threshold) lconverged=.FALSE.
2166	IF(MAXVAL(isotrop_up_cs,1) .GT. alb_threshold) lconverged=.FALSE.
2167
2168	IF(lconverged)THEN
2169	IF(lprint) WRITE(numout,*) 'Converged after this many steps: ',istep
2170	EXIT
2171	ENDIF
2172
2173	! This number here could also be externalized.
2174	IF(istep .GE. 1000)THEN
2175	WRITE(numout,) '********************************************************'
2176	WRITE(numout,'(A,I6,F14.10)') 'Albedo not converging!',istep,alb_threshold
2177	WRITE(numout,'(A)') ' collim_down_cs ' // &
2178	'isotrop_down_cs isotrop_up_cs'
2179	DO ilevel=0,nlevels_loc+1
2180	WRITE(numout,'(I4,3F14.10)') ilevel, &
2181	collim_down_cs(ilevel),isotrop_down_cs(ilevel),isotrop_up_cs(ilevel)
2182	ENDDO
2183	WRITE(numout,*) 'You should increase either the number' // &
2184	'of steps or the alb_threshold.'
2185	WRITE(numout,) '********************************************************'
2186	EXIT
2187	ENDIF ! istep .GE. 1000
2188	ENDDO ! convergence loop
2189
2190
2191	! now separate the collided from the uncollided light. Notice that
2192	! this is not really needed for any purposes other than debugging, as the
2193	! important quantity is the total amount of light striking the ground.
2194	Tran_UnColl_Tot(1:nlevels_loc)=Tran_UnColl_Tot_temp(1:nlevels_loc)
2195	Tran_Coll_Tot(1:nlevels_loc)=Tran_Coll_Tot(1:nlevels_loc)-Tran_UnColl_Tot(1:nlevels_loc)
2196
2197	! Some debugging information
2198	IF(lprint)THEN
2199	WRITE(numout,'(A)') ' Tran_Uncoll_Tot Tran_Coll_Tot Alb_Coll_Tot'
2200	DO ilevel=nlevels_loc,1,-1
2201	WRITE(numout,'(I4,3F10.6)') ilevel, &
2202	Tran_Uncoll_Tot(ilevel),Tran_Coll_Tot(ilevel),Alb_Coll_Tot(ilevel)
2203	ENDDO
2204	ENDIF
2205
2206	END SUBROUTINE propagate_fluxes
2207
2208
2209	END MODULE albedo

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

Download in other formats: