Context Navigation

source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_sechiba/albedo.f90 @ 7346

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

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

Download in other formats: