Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_sechiba/albedo_surface.f90 @ 8119

Last change on this file since 8119 was 7647, checked in by xiaoni.wang, 2 years ago
Modified globgrd.f90 in order that SAFRAN forcing files can be used in regional simulations. First offline test works on Obelix. More test would be better.
Property svn:executable set to ``*
File size: 147.0 KB

Line
1	! ===============================================================================================================================
2	! MODULE : albedo_surface
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_surface_main 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_surface
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	USE xios_orchidee
50	USE grid
51
52	IMPLICIT NONE
53
54	PRIVATE
55	PUBLIC :: albedo_surface_main, albedo_surface_initialize, albedo_surface_finalize, &
56	albedo_surface_clear, albedo_surface_xios_initialize
57
58	INTEGER, SAVE :: printlev_loc !! Local level of text output for current module
59	!$OMP THREADPRIVATE(printlev_loc)
60	LOGICAL, SAVE :: firstcall_albedo_surface=.TRUE. !! Initialization phase of this module
61	!$OMP THREADPRIVATE(firstcall_albedo_surface)
62	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: soilalb_dry !! Albedo values for the dry bare soil (unitless)
63	!$OMP THREADPRIVATE(soilalb_dry)
64	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: soilalb_wet !! Albedo values for the wet bare soil (unitless)
65	!$OMP THREADPRIVATE(soilalb_wet)
66	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: soilalb_moy !! Albedo values for the mean bare soil (unitless)
67	!$OMP THREADPRIVATE(soilalb_moy)
68	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bckgrnd_alb !! Albedo values for the background bare soil (unitless)
69	!$OMP THREADPRIVATE(bckgrnd_alb)
70
71	CONTAINS
72
73	!! =============================================================================================================================
74	!! SUBROUTINE: albedo_surface_xios_initialize
75	!!
76	!>\BRIEF Initialize xios dependant defintion before closing context defintion
77	!!
78	!! DESCRIPTION: Initialize xios dependant defintion needed for the interpolations done in albedo.
79	!! Reading is deactivated if the sechiba restart file exists because the variable
80	!! should be in the restart file already.
81	!! This subruting is called before closing context with xios_orchidee_close_definition in intersurf
82	!! via the subroutine sechiba_xios_initialize.
83	!!
84	!! \n
85	!_ ==============================================================================================================================
86
87	SUBROUTINE albedo_surface_xios_initialize
88
89	CHARACTER(LEN=255) :: filename !! Filename read from run.def
90	CHARACTER(LEN=255) :: name !! Filename without suffix .nc
91	LOGICAL :: lerr !! Flag to dectect error
92
93
94	! Read the file name for the albedo input file from run.def
95	filename = 'alb_bg.nc'
96	CALL getin_p('ALB_BG_FILE',filename)
97
98	! Remove suffix .nc from filename
99	name = filename(1:LEN_TRIM(FILENAME)-3)
100
101	! Define the file name in XIOS
102	CALL xios_orchidee_set_file_attr("albedo_file",name=name)
103
104	! Define default values for albedo
105	lerr=xios_orchidee_setvar('albbg_vis_default',0.129)
106	lerr=xios_orchidee_setvar('albbg_nir_default',0.247)
107
108	! Check if the albedo file will be read by XIOS, by IOIPSL or not at all
109	IF (xios_interpolation .AND. alb_bg_modis .AND. restname_in=='NONE') THEN
110	! The albedo file will be read using XIOS
111	IF (printlev>=2) WRITE(numout,*) 'Reading of albedo file will be done later using XIOS. The filename is ', filename
112
113	ELSE
114	IF (.NOT. alb_bg_modis) THEN
115	IF (printlev>=2) WRITE (numout,*) 'No reading of albedo will be done because alb_bg_modis=FALSE'
116	ELSE IF (restname_in=='NONE') THEN
117	IF (printlev>=2) WRITE (numout,*) 'The albedo file will be read later by IOIPSL'
118	ELSE
119	IF (printlev>=2) WRITE (numout,*) 'The albedo file will not be read because the restart file exists.'
120	END IF
121
122	! The albedo file will not be read by XIOS. Now deactivate albedo for XIOS.
123	CALL xios_orchidee_set_file_attr("albedo_file",enabled=.FALSE.)
124	CALL xios_orchidee_set_field_attr("bg_alb_vis_interp",enabled=.FALSE.)
125	CALL xios_orchidee_set_field_attr("bg_alb_nir_interp",enabled=.FALSE.)
126	END IF
127
128	END SUBROUTINE albedo_surface_xios_initialize
129
130	!! =============================================================================================================================
131	!! SUBROUTINE : albedo_surface_initialize
132	!!
133	!>\BRIEF Initialize albedo module
134	!!
135	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
136	!!
137	!! MAIN OUTPUT VARIABLE(S)
138	!!
139	!! REFERENCE(S) : None
140	!!
141	!! \n
142	!_ ==============================================================================================================================
143	!+++CHECK+++
144	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
145	! when running a larger domain. Needs to be corrected when implementing
146	! a global use of the multi-layer energy budget.
147	SUBROUTINE albedo_surface_initialize (kjit, kjpindex, rest_id, &
148	lalo, neighbours, resolution, contfrac, &
149	drysoil_frac, veget_max, coszang, frac_nobio, &
150	snow, snow_age, snow_nobio, &
151	snow_nobio_age, frac_snow_veg, frac_snow_nobio,&
152	z0m, laieff_fit, &
153	albedo, &
154	Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
155	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
156	laieff_isotrop, veget)
157	!+++++++++++
158
159	!! 0. Variable and parameter declaration
160	!! 0.1 Input variables
161	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
162	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
163	INTEGER(i_std), INTENT(in) :: rest_id !! _Restart_ file identifier
164	REAL(r_std), DIMENSION(kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
165	INTEGER(i_std), DIMENSION(kjpindex,NbNeighb), INTENT(in) :: neighbours !! neighoring grid points if land
166	REAL(r_std), DIMENSION(kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
167	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: contfrac !! Fraction of land in each grid box.
168
169	REAL(r_std), DIMENSION(:), INTENT(in) :: drysoil_frac !! Fraction of dry soil (unitless)
170	REAL(r_std), DIMENSION(:), INTENT(in) :: coszang !! The cosine of the solar zenith angle (unitless)
171	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! PFT coverage fraction of a PFT (= ind*cn_ind) (m^2 m^{-2})
172	REAL(r_std),DIMENSION (:,:), INTENT(in) :: veget !! Fraction of vegetation types
173	REAL(r_std), DIMENSION(:), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
174	REAL(r_std), DIMENSION(:,:), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
175	REAL(r_std), DIMENSION(:), INTENT(in) :: snow_age !! Snow age (days)
176	REAL(r_std), DIMENSION(:,:), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
177	REAL(r_std), DIMENSION(:,:), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces
178	REAL(r_std), DIMENSION(:), INTENT(in) :: z0m !! Roughness height for momentum of vegetated part (m)
179	TYPE(laieff_type),DIMENSION(:,:,:), INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
180	REAL(r_std), DIMENSION(:,:), INTENT(in) :: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc. (unitless ratio)
181	REAL(r_std), DIMENSION(:), INTENT(in) :: frac_snow_veg !! The fraction of ground covered with snow, between 0 and 1
182
183
184	!! 0.2 Output variables
185
186	REAL(r_std), DIMENSION(kjpindex,n_spectralbands), INTENT(out) :: albedo !! Albedo (two stream radiation transfer model)
187	!! for visible and near-infrared range
188	!! (unitless)
189	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(out) :: Isotrop_Abs_Tot_p !! Absorbed radiation per layer for photosynthesis
190	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(out) :: Isotrop_Tran_Tot_p !! Transmitted radiation per layer for photosynthesis
191	!+++CHECK+++
192	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
193	! when running a larger domain. Needs to be corrected when implementing
194	! a global use of the multi-layer energy budget.
195	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT(out) :: Light_Abs_Tot_mean !! total absorption for a given level * changed for enerbil integration
196	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT(out) :: Light_Alb_Tot_mean !! total albedo for a given level * changed for enerbil integration
197	!+++++++++++
198	REAL(r_std), DIMENSION(kjpindex,nlevels_tot,nvm), INTENT(out) :: laieff_isotrop !! Leaf Area Index Effective converts
199	!! 3D lai into 1D lai for two stream
200	!! radiation transfer model...this is for
201	!! isotropic light and only calculated once per day
202	!! @tex $(m^{2} m^{-2})$ @endtex
203
204
205
206	!! 0.4 Local variables
207	INTEGER :: ier, ji, ilevel
208	LOGICAL :: init_needed !! Variable indicating true if one or several restart
209	!! variables were missing
210	INTEGER(i_std), DIMENSION(kjpindex) :: index_dummy !! Dummy array
211
212	!_ ================================================================================================================================
213
214	! Get local printlev variable
215	printlev_loc=get_printlev('albedo')
216
217	IF (alb_bg_modis) THEN
218	! Allocate background soil albedo
219	ALLOCATE (bckgrnd_alb(kjpindex,2),stat=ier)
220	IF (ier /= 0) CALL ipslerr_p(3,'albedo_surface_initialize','Pb in allocation for bckgrnd_alb','','')
221
222	! Read background albedo from restart file
223	CALL ioconf_setatt_p('UNITS', '-')
224	CALL ioconf_setatt_p('LONG_NAME','Background soil albedo for visible and near-infrared range')
225	CALL restget_p (rest_id, 'bckgrnd_alb', nbp_glo, 2, 1, kjit, .TRUE., bckgrnd_alb, "gather", nbp_glo, index_g)
226
227	! Initialize by interpolating from file if the variable was not in restart file
228
229	IF ( ALL(bckgrnd_alb(:,:) == val_exp) ) THEN
230	CALL read_background_albedo(kjpindex, lalo, neighbours, resolution, contfrac)
231	END IF
232	CALL xios_orchidee_send_field("bckgrnd_alb",bckgrnd_alb)
233
234	ELSE
235	! Dry soil albedo
236	ALLOCATE (soilalb_dry(kjpindex,2),stat=ier)
237	IF (ier /= 0) CALL ipslerr_p(3,'albedo_surface_initialize','Pb in allocation for soilalb_dry','','')
238
239	! Wet soil albedo
240	ALLOCATE (soilalb_wet(kjpindex,2),stat=ier)
241	IF (ier /= 0) CALL ipslerr_p(3,'albedo_surface_initialize','Pb in allocation for soilalb_wet','','')
242
243	! Mean soil albedo
244	ALLOCATE (soilalb_moy(kjpindex,2),stat=ier)
245	IF (ier /= 0) CALL ipslerr_p(3,'albedo_surface_initialize','Pb in allocation for soilalb_moy','','')
246
247
248	!! Read variables from restart file or initialize them
249	! dry soil albedo
250	CALL ioconf_setatt_p('UNITS', '-')
251	CALL ioconf_setatt_p('LONG_NAME','Dry bare soil albedo')
252	CALL restget_p (rest_id,'soilalbedo_dry' , nbp_glo, 2, 1, kjit, .TRUE., soilalb_dry, "gather", nbp_glo, index_g)
253
254	! wet soil albedo
255	CALL ioconf_setatt_p('UNITS', '-')
256	CALL ioconf_setatt_p('LONG_NAME','Wet bare soil albedo')
257	CALL restget_p (rest_id, 'soilalbedo_wet', nbp_glo, 2, 1, kjit, .TRUE., soilalb_wet, "gather", nbp_glo, index_g)
258
259	! mean soil aledo
260	CALL ioconf_setatt_p('UNITS', '-')
261	CALL ioconf_setatt_p('LONG_NAME','Mean bare soil albedo')
262	CALL restget_p (rest_id, 'soilalbedo_moy', nbp_glo, 2, 1, kjit, .TRUE., soilalb_moy, "gather", nbp_glo, index_g)
263
264	! Initialize the variables if not found in restart file
265	IF ( ALL(soilalb_wet(:,:) == val_exp) .OR. &
266	ALL(soilalb_dry(:,:) == val_exp) .OR. &
267	ALL(soilalb_moy(:,:) == val_exp)) THEN
268	! One or more of the variables were not in the restart file.
269	! Call routine albedo_surface_soilalb to calculate them.
270	CALL albedo_surface_soilalb(kjpindex, lalo, neighbours, resolution, contfrac)
271	WRITE(numout,*) '---> val_exp ', val_exp
272	WRITE(numout,*) '---> ALBEDO_wet VIS:', MINVAL(soilalb_wet(:,ivis)), MAXVAL(soilalb_wet(:,ivis))
273	WRITE(numout,*) '---> ALBEDO_wet NIR:', MINVAL(soilalb_wet(:,inir)), MAXVAL(soilalb_wet(:,inir))
274	WRITE(numout,*) '---> ALBEDO_dry VIS:', MINVAL(soilalb_dry(:,ivis)), MAXVAL(soilalb_dry(:,ivis))
275	WRITE(numout,*) '---> ALBEDO_dry NIR:', MINVAL(soilalb_dry(:,inir)), MAXVAL(soilalb_dry(:,inir))
276	WRITE(numout,*) '---> ALBEDO_moy VIS:', MINVAL(soilalb_moy(:,ivis)), MAXVAL(soilalb_moy(:,ivis))
277	WRITE(numout,*) '---> ALBEDO_moy NIR:', MINVAL(soilalb_moy(:,inir)), MAXVAL(soilalb_moy(:,inir))
278	ENDIF
279
280	END IF
281
282	!! Read variables from restart file
283	! If a variable is not found in the restart file, the flag init_needed is set to true.
284	init_needed=.FALSE.
285	! albedo : DIM(kjpindex, n_spectralbands)
286	CALL ioconf_setatt_p('UNITS', '-')
287	CALL ioconf_setatt_p('LONG_NAME','albedo')
288	CALL restget_p (rest_id,'albedo' , nbp_glo, n_spectralbands, 1, kjit, .TRUE., &
289	albedo, "gather", nbp_glo, index_g)
290	IF (ALL(albedo(:,:) == val_exp)) init_needed=.TRUE.
291
292	! Isotrop_Abs_Tot_p : DIM(kjpindex,nvm,nlevels_tot)
293	CALL ioconf_setatt_p('UNITS', '')
294	CALL ioconf_setatt_p('LONG_NAME','Absorbed radiation per layer for photosynthesis')
295	CALL restget_p (rest_id, 'Isotrop_Abs_Tot_p', nbp_glo, nvm, nlevels_tot, kjit, .TRUE., &
296	Isotrop_Abs_Tot_p, "gather", nbp_glo, index_g)
297	IF ( ALL(Isotrop_Abs_Tot_p(:,:,:) == val_exp ) ) init_needed=.TRUE.
298
299	! Isotrop_Tran_Tot_p : DIM(kjpindex,nvm,nlevels_tot)
300	CALL ioconf_setatt_p('UNITS', '')
301	CALL ioconf_setatt_p('LONG_NAME','Transmitted radiation per layer for photosynthesis')
302	CALL restget_p (rest_id, 'Isotrop_Tran_Tot_p', nbp_glo, nvm, nlevels_tot, kjit, .TRUE., &
303	Isotrop_Tran_Tot_p, "gather", nbp_glo, index_g)
304	IF ( ALL(Isotrop_Tran_Tot_p(:,:,:) == val_exp ) ) init_needed=.TRUE.
305
306	!+++CHECK+++
307	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
308	! when running a larger domain. Needs to be corrected when implementing
309	! a global use of the multi-layer energy budget.
310	!!$ ! Light_Abs_Tot_mean, Light_Alb_Tot_mean : DIM(nlevels_tot) : Needed in mleb_main
311	!!$ IF (is_root_prc) THEN
312	!!$ CALL ioconf_setatt_p('UNITS', '')
313	!!$ CALL ioconf_setatt_p('LONG_NAME','')
314	!!$ CALL restget (rest_id, 'Light_Abs_Tot_mean', nlevels_tot, 1, 1, kjit, .TRUE., Light_Abs_Tot_mean)
315	!!$
316	!!$ CALL ioconf_setatt_p('UNITS', '')
317	!!$ CALL ioconf_setatt_p('LONG_NAME','')
318	!!$ CALL restget (rest_id, 'Light_Alb_Tot_mean', nlevels_tot, 1, 1, kjit, .TRUE., Light_Alb_Tot_mean)
319	!!$ END IF
320	!!$ CALL bcast(Light_Abs_Tot_mean)
321	!!$ CALL bcast(Light_Alb_Tot_mean)
322	!!$ IF ( ALL(Light_Abs_Tot_mean(:) == val_exp ) ) init_needed=.TRUE.
323	!!$ IF ( ALL(Light_Alb_Tot_mean(:) == val_exp ) ) init_needed=.TRUE.
324
325	! laieff_isotrop : DIM(kjpindex, nlevels_tot,nvm) : Needed in mleb_main called before condveg_main
326	!+++++++++++
327	CALL ioconf_setatt_p('UNITS', '')
328	CALL ioconf_setatt_p('LONG_NAME','')
329	CALL restget_p (rest_id, 'laieff_isotrop', nbp_glo, nlevels_tot, nvm, kjit, .TRUE., &
330	laieff_isotrop, "gather", nbp_glo, index_g)
331	IF ( ALL(laieff_isotrop(:,:,:) == val_exp ) ) init_needed=.TRUE.
332
333
334	IF (init_needed) THEN
335	! One or more of the variables were not in the restart file.
336	! Call routine albedo_surface_main to calculate them.
337	! In the initialization phase, the arguments for hist_id, hist_id2 and index will never be used.
338	! Therefor send dummy variables.
339	index_dummy(:)=0
340	!+++CHECK+++
341	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
342	! when running a larger domain. Needs to be corrected when implementing
343	! a global use of the multi-layer energy budget.
344	CALL albedo_surface_main( &
345	kjit, kjpindex, 0, 0, &
346	index_dummy, drysoil_frac, veget_max, coszang, &
347	frac_nobio, frac_snow_veg, frac_snow_nobio, &
348	snow, snow_age, snow_nobio, snow_nobio_age, &
349	z0m, laieff_fit, &
350	albedo, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
351	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
352	laieff_isotrop, veget)
353	!+++++++++++
354	END IF
355
356
357	! Initalization finished
358	firstcall_albedo_surface=.FALSE.
359
360	END SUBROUTINE albedo_surface_initialize
361
362	!! ==============================================================================================================================
363	!! SUBROUTINE : albedo_surface_main
364	!!
365	!>\BRIEF This subroutine calculates the albedo for two stream radiation transfer model. This routine includes
366	!! the effect of snow on the background albedo, through one of two methods.
367	!!
368	!! DESCRIPTION : The albedo for two stream radiation transfer model is calculated for both the visible and near-infrared
369	!! domain. First the mean albedo of the bare soil is calculated. Two options exist:
370	!! either the soil albedo depends on soil wetness (drysoil_frac variable), or the soil albedo
371	!! is set to a mean soil albedo value.
372	!!
373	!! NOTE: the main output variable, albedo, is an unweighted average of the direct and diffuse albedos
374	!! for each grid point...this is done in this way right now to be consistent with the new scheme,
375	!! but it doesn't have to be combined like that once we have an energy budget which can
376	!! use both types of light directly
377	!!
378	!! RECENT CHANGE(S): None
379	!!
380	!! MAIN OUTPUT VARIABLE(S): ::albedo
381	!!
382	!! REFERENCE(S) : B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski, N. Gobron and M. M. Verstraete (2006).
383	!! Simplifying the Interaction of Land Surfaces with Radiation for Relating Remote Sensing Products to Climate Models.
384	!! Journal of Geophysical Research. Vol 111, D02116.
385	!!
386	!! FLOWCHART : None
387	!! \n
388	!!
389	!!
390	!!
391	!_ ================================================================================================================================
392	!+++CHECK+++
393	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
394	! when running a larger domain. Needs to be corrected when implementing
395	! a global use of the multi-layer energy budget.
396	SUBROUTINE albedo_surface_main( &
397	kjit, kjpindex, hist_id, hist2_id, &
398	index, drysoil_frac, veget_max, coszang, &
399	frac_nobio, frac_snow_veg, frac_snow_nobio, &
400	snow, snow_age, snow_nobio, snow_nobio_age, &
401	z0m, laieff_fit, &
402	albedo, Light_Abs_Tot, Light_Tran_Tot, &
403	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
404	laieff_isotrop, veget)
405	!+++++++++++
406
407	!! 0. Variable and parameter declaration
408
409	!! 0.1 Input variables
410	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
411	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
412	INTEGER(i_std), INTENT (in) :: hist_id !! _History_ file identifier
413	INTEGER(i_std), INTENT (in) :: hist2_id !! _History_ file 2 identifier
414	INTEGER(i_std), DIMENSION(kjpindex), INTENT (in) :: index !! Indeces of the points on the map
415	REAL(r_std), DIMENSION(:), INTENT(in) :: drysoil_frac !! Fraction of dry soil (unitless)
416	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! PFT coverage fraction of a PFT (= ind*cn_ind)
417	!! (m^2 m^{-2})
418	REAL(r_std),DIMENSION (:,:), INTENT(in) :: veget !! Fraction of vegetation types
419	REAL(r_std), DIMENSION(:), INTENT(in) :: coszang !! The cosine of the solar zenith angle (unitless)
420	REAL(r_std),DIMENSION (:,:), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
421	!! continental ice, lakes, etc. (unitless)
422	REAL(r_std),DIMENSION(:), INTENT(in) :: frac_snow_veg !! The fraction of ground covered with snow, between 0 and 1
423	REAL(r_std), DIMENSION(:,:), INTENT(in) :: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
424	!! (unitless ratio)
425	REAL(r_std),DIMENSION (:), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
426	REAL(r_std),DIMENSION (:,:), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
427	REAL(r_std),DIMENSION (:), INTENT(in) :: snow_age !! Snow age (days)
428	REAL(r_std),DIMENSION (:,:), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
429	REAL(r_std), DIMENSION(:), INTENT(in) :: z0m !! Roughness height for momentum of vegetated part (m)
430	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
431
432	!! 0.2 Output variables
433
434	REAL(r_std), DIMENSION (kjpindex,n_spectralbands), &
435	INTENT (out) :: albedo !! Albedo (two stream radiation transfer model)
436	!! for visible and near-infrared range
437	!! (unitless)
438	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT(out) :: Light_Abs_Tot !! Absorbed radiation per layer, averaged between light
439	!! from a collimated (direct) source and that from an
440	!! isotropic (diffuse) source. Expressed as the fraction
441	!! of overall light hitting the canopy.
442	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT(out) :: Light_Tran_Tot !! Transmitted radiation per layer, averaged between light
443	!! from a collimated (direct) source and that from an
444	!! isotropic (diffuse) source. Expressed as the fraction
445	!! of overall light hitting the canopy.
446	!+++CHECK+++
447	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
448	! when running a larger domain. Needs to be corrected when implementing
449	! a global use of the multi-layer energy budget.
450	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT (out) :: Light_Abs_Tot_mean !! total absorption for a given level * changed for enerbil integration
451	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT (out) :: Light_Alb_Tot_mean !! total albedo for a given level * changed for enerbil integration
452	!+++++++++++
453	REAL(r_std), DIMENSION(kjpindex,nlevels_tot,nvm), INTENT(out) :: laieff_isotrop !! Leaf Area Index Effective converts
454	!! 3D lai into 1D lai for two stream
455	!! radiation transfer model...this is for
456	!! isotropic light and only calculated once per day
457	!! @tex $(m^{2} m^{-2})$ @endtex
458
459
460
461	!! 0.3 Modified variables
462
463
464	!! 0.4 Local variables
465
466	REAL(r_std),DIMENSION (nlevels_tot) :: laieff_isotrop_pft, laieff_collim_pft !!
467	REAL(r_std) :: cosine_sun_angle !! the cosine of the solar zenith angle
468	INTEGER(i_std) :: ks !! Index for visible and near-infraread range
469	INTEGER(i_std) :: ivm !! Index for vegetative PFTs
470	INTEGER(i_std) :: ipts !! Index for spatial domain
471	INTEGER(i_std) :: ilevel !! Index for canopy levels
472	REAL(r_std) :: leaf_reflectance
473	REAL(r_std) :: leaf_transmittance
474	REAL(r_std) :: br_base_temp,br_base_temp_collim,br_base_temp_isotrop
475	REAL(r_std) :: leaf_psd_temp
476	REAL(r_std) :: leaf_single_scattering_albedo
477	REAL(r_std), DIMENSION(kjpindex,n_spectralbands) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
478	REAL(r_std), DIMENSION(kjpindex,n_spectralbands) :: albedo_snow !! Snow albedo (unitless ratio)
479	REAL(r_std), DIMENSION(kjpindex) :: albedo_glob !! Mean albedo (unitless ratio)
480	REAL(r_std), DIMENSION(kjpindex,nvm, n_spectralbands) :: albedo_pft !! Albedo (two stream radiation transfer model)
481	!! for visible and near-infrared range for each PFT (unitless)
482	REAL(r_std), DIMENSION(kjpindex,nlevels_tot,nvm) :: laieff_collim !! Leaf Area Index Effective for direct light
483	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Tran_Uncoll !! collimated total transmission of uncollided light for a given level
484	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Tran_Coll !! collimated total transmission for a given level
485	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Abs_Tot !! collimated total absorption for a given level * changed for enerbil integration
486	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Alb_Tot !! Collimated (direct) total albedo for a given level * changed for enerbil integration
487
488	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Alb_Tot !! isotropic (diffuse) total albedo for a given level
489	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Tran_Uncoll !! isotropic total transmission of uncollided light for a given level
490	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Tran_Coll !! isotropic total transmission for a given level
491	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Abs_Tot !! isotropic total absporbtion for a given level
492	INTEGER :: jv, inno
493	REAL(r_std) :: alb_nobio
494	REAL(r_std):: laieff_collim_1, laieff_isotrop_1, Collim_Alb_Tot_1,Collim_Tran_Tot_1,Collim_Abs_Tot_1,&
495	Isotrop_Alb_Tot_1,Isotrop_Tran_Tot_1,Isotrop_Abs_Tot_1,&
496	Collim_Tran_Uncollided_1, Isotrop_Tran_Uncollided_1 !! Values for the one level solution
497
498	REAL(r_std):: converged_albedo
499	REAL(r_std):: isotrop_angle
500	LOGICAL :: lconverged
501	LOGICAL :: lprint !! a flag for printing some debug statements
502	REAL(r_std), DIMENSION(kjpindex,nvm,n_spectralbands) :: snowa_veg !! snow albedo due to vegetated surfaces, between 0 and 1
503	REAL(r_std), DIMENSION(kjpindex,nnobio,n_spectralbands) :: snowa_nobio !! Albedo of snow covered area on continental ice,
504	!! lakes, etc. (unitless ratio)
505
506	REAL(r_std), DIMENSION(kjpindex,n_spectralbands) :: albedo_diag !! Array used to filter the values that will be send to xios
507	REAL(r_std), DIMENSION(kjpindex,n_spectralbands) :: alb_bare_diag !! Array used to filter the values that will be send to xios
508	REAL(r_std), DIMENSION(kjpindex,n_spectralbands) :: albedo_snow_diag !! Array used to filter the values that will be send to xios
509
510	!_ ================================================================================================================================
511
512	IF (printlev_loc>=3) WRITE(numout,*) 'Entering albedo_surface_main'
513
514	IF (printlev_loc>=3) WRITE(numout,*) 'nlevels_tot',nlevels_tot
515	IF (printlev_loc>=4) THEN
516	WRITE(numout,*) 'laieff_isotrop start albedo for test_pft',test_pft,'is', &
517	laieff_isotrop(test_grid,:,test_pft)
518
519	WRITE(numout,*)'albedo.f90, coszang',coszang(:)
520	ENDIF
521
522	!! 1. We will now calculate the background reflectance to be used in the model.
523	! The parameters read in from the input file do not include the effect
524	! of snow, so we need to figure out the snow's contribution for each
525	! grid space
526	! Initialize some output variables
527	albedo_pft(:,:,:)=zero
528	Light_Abs_Tot(:,:,:)=zero
529	Light_Tran_Tot(:,:,:)=un
530	!+++CHECK+++
531	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
532	! when running a larger domain. Needs to be corrected when implementing
533	! a global use of the multi-layer energy budget.
534	!!$ Light_Abs_Tot_mean(:) = zero
535	!!$ Light_Alb_Tot_mean(:) = zero
536	!+++++++++++
537
538	! We need to calculate the LAI effective for this solar angle. The other LAI
539	! effective that we calculated was computed at an angle of 60.0 degrees for
540	! the isotropic contribution. Notice that, in many situations, the results
541	! will be exactly the same. One can show that Pgap is proportional to
542	! 1/cos(theta) if the level bottom is below all canopy elements, which
543	! means that the cos(theta) in the calculation of the LAIeff will be
544	! cancelled out
545	isotrop_angle=COS(60.0/180.0*pi)
546	DO ipts=1,kjpindex
547	DO ivm=1,nvm
548	DO ilevel=1,nlevels_tot
549	laieff_isotrop(ipts,ilevel,ivm) = &
550	calculate_laieff_fit(isotrop_angle,laieff_fit(ipts,ivm,ilevel))
551	laieff_collim(ipts,ilevel,ivm) = &
552	calculate_laieff_fit(coszang(ipts),laieff_fit(ipts,ivm,ilevel))
553	! Because we are fitting these values, it's possible that some
554	! will drop below zero. This causes a serious problem here,
555	! so let's make sure this doesn't happen. Unfortunately, by
556	! doing this, we lose a test of things going wrong (an unfitted
557	! negative LAIeff can be the sign of another problem), but we
558	! don't really have a choice.
559	IF(laieff_isotrop(ipts,ilevel,ivm) .LT. min_stomate) &
560	laieff_isotrop(ipts,ilevel,ivm)=zero
561	IF(laieff_collim(ipts,ilevel,ivm) .LT. min_stomate) &
562	laieff_collim(ipts,ilevel,ivm)=zero
563	ENDDO
564	ENDDO
565	ENDDO
566
567	! Debug
568	IF(printlev_loc>=4)THEN
569	DO ivm=1,nvm
570	IF(ivm == test_pft)THEN
571	DO ipts = 1, kjpindex
572	IF(ipts == test_grid)THEN
573	WRITE(numout,*) 'Laieff_collim and laieff_isotrop in albedo'
574	WRITE(numout,) 'ACOS(coszang(ipts))/pi180: ',&
575	ACOS(coszang(ipts))/pi*180.0_r_std
576	WRITE(numout,*) 'coszang(ipst)', coszang
577	DO ilevel=nlevels_tot,1,-1
578	WRITE(numout,*) ilevel,laieff_collim(ipts,ilevel,ivm),&
579	laieff_isotrop(ipts,ilevel,ivm)
580	END DO
581	DO ilevel=nlevels_tot,1,-1
582	WRITE(numout,*) 'Fitting parameters: ',&
583	laieff_fit(ipts,ivm,ilevel)%a,&
584	laieff_fit(ipts,ivm,ilevel)%b,&
585	laieff_fit(ipts,ivm,ilevel)%c,&
586	laieff_fit(ipts,ivm,ilevel)%d,&
587	laieff_fit(ipts,ivm,ilevel)%e
588	ENDDO
589	ENDIF
590	ENDDO
591	ENDIF
592	ENDDO
593	ENDIF
594	!-
595
596	! Calculate the snow albedo
597	CALL calculate_snow_albedo(kjpindex, coszang, snow, snow_nobio, &
598	snow_age, snow_nobio_age, frac_nobio, albedo_snow, &
599	snowa_veg, frac_snow_veg, snowa_nobio, frac_snow_nobio, &
600	veget_max, z0m, veget)
601
602	! Check for possible problems in the calcualtion of the fraction
603	! of the grid cell covered by snow. Fractions should be positive
604	! numbers
605	DO ipts=1,kjpindex
606	IF(frac_snow_veg(ipts) .LT. 0.0)THEN
607	WRITE(numout,*) 'SNOWFRAC ipts ',ipts,frac_snow_veg(ipts)
608	CALL ipslerr_p (3,'albedo', 'snow frac error','','')
609	END IF
610
611	DO inno=1,nnobio
612	IF(frac_snow_nobio(ipts,inno) .LT. 0.0)THEN
613	WRITE(numout,*) 'SNOWFRAC ipts inno',ipts,inno,frac_snow_nobio(ipts,inno)
614	CALL ipslerr_p (3,'albedo', 'nobio snow frac error','','')
615	ENDIF
616	ENDDO
617	ENDDO
618
619	! Initialize the albedo in every square for both spectra by calculating the
620	! bare soil albedo
621	DO ipts = 1, kjpindex
622	DO ks = 1, n_spectralbands
623
624	! If 'alb_bare_model' is set to TRUE, the soil albedo calculation
625	! depends on soil moisture
626	! If 'alb_bare_model' is set to FALSE, the mean soil albedo is used
627	! without the dependance on soil moisture
628	! see subroutine 'conveg_soilalb'
629	! Note that these two methods are considered old. The albedo
630	! background map has been done more recently, as is now
631	! used as the default in simulations with TAG 2.2. The
632	! lines below ensure that the albedos over the Sahara are identical
633	! for Tags 2.2, 3.0, and 4.0 when alb_bg_modis=.TRUE.
634	IF ( alb_bg_modis ) THEN
635	alb_bare(ipts,ks) = bckgrnd_alb(ipts,ks)
636	ELSEIF ( alb_bare_model ) THEN
637	alb_bare(ipts,ks) = soilalb_wet(ipts,ks) + drysoil_frac(ipts) * &
638	(soilalb_dry(ipts,ks) - soilalb_wet(ipts,ks))
639	ELSE
640	alb_bare(ipts,ks) = soilalb_moy(ipts,ks)
641	ENDIF
642
643
644	! we need to take into account any snow that has fallen on the bare
645	! soil since there may be some bare soil, initialize the albedo with
646	! this fraction
647	albedo(ipts,ks) = veget_max(ipts,1) * (alb_bare(ipts,ks)*&
648	(un-frac_snow_veg(ipts)) + &
649	frac_snow_veg(ipts)*snowa_veg(ipts,1,ks))
650
651	ENDDO ! ks = 1, n_spectralbands
652
653	ENDDO ! ipts = 1, kjpindex
654
655	!! 2. Calculation of albedo of the whole grid cell, taking into account all
656	!! PFTs (except bare soil) and spectral bands
657	DO ipts = 1, kjpindex ! loop over the grid squares
658
659	!+++CHECK+++
660	! observations from the MERGE, sinang is calculated as cos of the zenith
661	! angle in solar.f90 (i.e. it is the cosine of the zenith angle).
662	! Then maybe also the cosine_sun_angle is a redundant variable.
663	cosine_sun_angle = coszang(ipts)
664	!+++++++++++
665
666	! Debug
667	IF(printlev_loc >= 4) THEN
668	WRITE(numout,*)'cosine_sun_angle',cosine_sun_angle
669	ENDIF
670
671	! For the coupled model, the angles can be negative.
672	! Offline, they are always between zero and 90 degrees.
673	IF(cosine_sun_angle .LE. min_sechiba)THEN
674
675	! It's night, so no sunlight...don't need to calculate albedo.
676	! Set it equal to zero, because the snow albedo has already been
677	! calculated and it looks funny on the coupled run maps otherwise.
678	albedo(ipts,:)=zero
679	CYCLE
680
681	ENDIF
682
683	! Are there any non-bio PFTs here that we need to take into account?
684	DO ks = 1, n_spectralbands ! Loop over # of spectra
685	DO jv = 1, nnobio
686	! now update the albedo
687	IF ( jv .EQ. iice ) THEN
688	alb_nobio = alb_ice(ks)
689	ELSE
690	WRITE(numout,*) 'jv=',jv
691	CALL ipslerr_p (3,'Albedo.f90', &
692	'DO NOT KNOW ALBEDO OF THIS SURFACE TYPE','','')
693	ENDIF
694
695	! this takes into account both snow covered and non-snow
696	! covered non-bio regios in all grid squares
697	albedo(ipts,ks) = albedo(ipts,ks) + &
698	( frac_nobio(ipts,jv) ) * &
699	( (un-frac_snow_nobio(ipts,jv)) * alb_nobio + &
700	( frac_snow_nobio(ipts,jv) ) * snowa_nobio(ipts,jv,ks) )
701
702	!+++CHECK+++
703	! Already calculated in calculate_snow_albedo. SL thinks
704	! this should not be repeated here.
705	!!$ albedo_snow(ipts,ks) = albedo_snow(ipts,ks) + &
706	!!$ ( frac_nobio(ipts,jv) ) * ( frac_snow_nobio(ipts,jv) ) * &
707	!!$ snowa_nobio(ipts,jv,ks)
708	!+++++++++++
709
710	ENDDO ! jv = 1, nnobio
711	ENDDO ! ks = 1, n_spectralbands
712
713	! Calculate the albedo of the vegetated fractions of the pixel
714	DO ivm = 2, nvm ! Loop over # of PFTs
715
716	! for this grid square and this vegetation type, we have a set LAI
717	! effective
718	laieff_collim_pft(1:nlevels_tot) = &
719	laieff_collim(ipts,1:nlevels_tot,ivm)
720	laieff_isotrop_pft(1:nlevels_tot) = &
721	laieff_isotrop(ipts,1:nlevels_tot,ivm)
722
723	IF(veget_max(ipts,ivm) == zero)THEN
724	! this vegetation type is not present, so no reason to do the
725	! calculation
726	CYCLE
727	ENDIF
728
729	DO ks = 1, n_spectralbands ! Loop over # of spectra
730
731	! set single scattering albedo and preferred scattering direction
732	leaf_single_scattering_albedo = leaf_ssa(ivm,ks)
733	leaf_psd_temp = leaf_psd(ivm,ks)
734
735	! calculate the rleaf and tleaf from wl=rl+tl and dl=rl/tl
736	leaf_transmittance = leaf_single_scattering_albedo/ &
737	(leaf_psd_temp+1)
738	leaf_reflectance = leaf_psd_temp * &
739	leaf_transmittance
740
741	! We need to take into account the effect of snow on the
742	! background reflectance here. From the snow above I can
743	! calculate the true background reflectance for this PFT.
744	! The flag alb_bg_modis is a bit confusing because in
745	! both cases the data source is the MODIS albedo product
746	IF (alb_bg_modis) THEN
747
748	! If the flag is set to TRUE the model will read a
749	! spatially explicit map of the background albedo for
750	! each PIXEL. The background albedo of the pixel is
751	! then used to calculate the background albedo of the
752	! PFT accounting for snow fraction and snow albedo. Note
753	! that both snow fraction and snow albedo vary over time.
754	br_base_temp = (un-frac_snow_veg(ipts)) * &
755	bckgrnd_alb(ipts,ks) + &
756	( frac_snow_veg(ipts) ) * snowa_veg(ipts,ivm,ks)
757
758	ELSE
759
760	! If the flag is set to FALSE the model uses the original
761	! parameterization at the PFT level. Each PFT has a fixed
762	! background albedo. The background albedo is recalculated
763	! taking into account the snowfraction and the snow albedo.
764	! Note that snow fraction and snoe albedo both vary over
765	! time.
766	br_base_temp= (un-frac_snow_veg(ipts)) * &
767	bgd_reflectance(ivm,ks) + &
768	( frac_snow_veg(ipts) ) * snowa_veg(ipts,ivm,ks)
769
770	ENDIF
771
772	! At this step, we assume that there is no difference between the
773	! direct and the diffuse background reflectance, which is true for
774	! the background but not the canopy albedo
775	br_base_temp_collim=br_base_temp
776	br_base_temp_isotrop=br_base_temp
777
778	! Debug - Set flag for extra output
779	IF(printlev_loc>=4 .AND. test_pft == ivm .AND. &
780	test_grid == ipts)THEN
781	lprint=.TRUE.
782	ELSE
783	lprint=.FALSE.
784	ENDIF
785	!-
786
787	! Now solve the multilevel scheme
788	CALL multilevel_albedo(cosine_sun_angle, &
789	leaf_single_scattering_albedo, &
790	leaf_psd_temp, br_base_temp_collim, br_base_temp_isotrop, &
791	laieff_collim_pft, laieff_isotrop_pft, lconverged, &
792	Collim_Alb_Tot, Collim_Tran_Coll, Collim_Abs_Tot, &
793	Isotrop_Alb_Tot, &
794	Isotrop_Tran_Coll, Isotrop_Abs_Tot, Collim_Tran_Uncoll, &
795	Isotrop_Tran_Uncoll, lprint)
796
797	!+++CHECK+++
798	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
799	! when running a larger domain. Needs to be corrected when implementing
800	! a global use of the multi-layer energy budget.
801	!!$ ! Here we weight the light in a consistent way between diffuse and
802	!!$ ! direct sources
803	!!$ Light_Abs_Tot_mean(:) = Light_Abs_Tot_mean(:) + &
804	!!$ direct_light_weightCollim_Abs_Tot(:)+ (un-direct_light_weight)Isotrop_Abs_Tot
805	!!$ Light_Alb_Tot_mean(:) = Light_Alb_Tot_mean(:) + &
806	!!$ direct_light_weightCollim_Alb_Tot(:)+ (un-direct_light_weight)Isotrop_Alb_Tot
807	!+++++++++++
808
809	! This is a lot of information printed out for the paper
810	! on the multilevel albedo scheme (McGrath et al 2016, GMDD).
811	! Probably not useful for anyone else, but I'm keeping it in
812	! until the paper is published. This is why I'm protecting it
813	! with an IF statement. I do not use printlev_loc>=4 since it
814	! really is too much information for normal usage.
815	IF(.FALSE.)THEN
816
817	WRITE(6,'(A,2F14.10)') 'Background reflectance ',&
818	br_base_temp_collim,br_base_temp_isotrop
819	WRITE(6,'(A,F14.10)') 'Single scattering albedo (wl): ',&
820	leaf_single_scattering_albedo
821	WRITE(6,'(A,F14.10)') 'Preferred scattering direction (dl): ',&
822	leaf_psd_temp
823
824	WRITE(6,) 'Collimated light ',ACOS(cosine_sun_angle)/pi180.0
825	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' 1 & ',laieff_collim_1,&
826	' & ',Collim_Alb_Tot_1,' & ',&
827	Collim_Tran_Tot_1-Collim_Tran_Uncollided_1,' & ',&
828	Collim_Tran_Uncollided_1,' & ',&
829	Collim_Alb_Tot_1+Collim_Tran_Tot_1,' \\ '
830	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' 2 & ',laieff_collim_1,&
831	' & ',Collim_Alb_Tot(nlevels_tot),' & ',&
832	Collim_Tran_Coll(1),' & ',Collim_Tran_Uncoll(1),' & ',&
833	Collim_Alb_Tot(nlevels_tot)+Collim_Tran_Coll(1)+&
834	Collim_Tran_Uncoll(1),' \\ '
835	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' Diff & ',laieff_collim_1,&
836	' & ',ABS(Collim_Alb_Tot_1-Collim_Alb_Tot(nlevels_tot)),' & ',&
837	ABS((Collim_Tran_Tot_1-Collim_Tran_Uncollided_1)-&
838	Collim_Tran_Coll(1)),' & ',&
839	ABS(Collim_Tran_Uncollided_1-Collim_Tran_UnColl(1)),' & ',&
840	ABS((Collim_Alb_Tot_1+Collim_Tran_Tot_1)-&
841	(Collim_Alb_Tot(nlevels_tot)+Collim_Tran_Coll(1)+&
842	Collim_Tran_Uncoll(1))),' \\ '
843
844	WRITE(6,*) 'Isotropic light '
845	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' 1 & ',&
846	laieff_isotrop_1,' & ',Isotrop_Alb_Tot_1,' & ',&
847	Isotrop_Tran_Tot_1-Isotrop_Tran_Uncollided_1,' & ',&
848	Isotrop_Tran_Uncollided_1,' & ',&
849	Isotrop_Alb_Tot_1+Isotrop_Tran_Tot_1,' \\ '
850	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' 2 & ',&
851	laieff_isotrop_1,' & ',Isotrop_Alb_Tot(nlevels_tot),' & ',&
852	Isotrop_Tran_Coll(1),' & ',Isotrop_Tran_Uncoll(1),' & ',&
853	Isotrop_Alb_Tot(nlevels_tot)+Isotrop_Tran_Coll(1)+&
854	Isotrop_Tran_Uncoll(1),' \\ '
855	WRITE(6,'(A,F4.1,A,4(F10.6,A))') ' Diff & ',&
856	laieff_isotrop_1,&
857	' & ',ABS(Isotrop_Alb_Tot_1-Isotrop_Alb_Tot(nlevels_tot)),' & ',&
858	ABS((Isotrop_Tran_Tot_1-Isotrop_Tran_Uncollided_1)-&
859	Isotrop_Tran_Coll(1)),' & ',&
860	ABS(Isotrop_Tran_Uncollided_1-Isotrop_Tran_UnColl(1)),' & ',&
861	ABS((Isotrop_Alb_Tot_1+Isotrop_Tran_Tot_1)-&
862	(Isotrop_Alb_Tot(nlevels_tot)+Isotrop_Tran_Coll(1)+&
863	Isotrop_Tran_Uncoll(1))),' \\ '
864
865	WRITE(6,1010) 'COLLIMTRAN ',laieff_collim_1, &
866	laieff_collim_pft(nlevels_tot),1,&
867	leaf_single_scattering_albedo,br_base_temp_collim,&
868	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
869	Collim_Tran_Tot_1,(Collim_Tran_Coll(1)+Collim_Tran_Uncoll(1))
870	WRITE(6,1010) 'COLLIMALB ',laieff_collim_1, &
871	laieff_collim_pft(nlevels_tot),2,&
872	leaf_single_scattering_albedo,br_base_temp_collim,&
873	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
874	Collim_Alb_Tot_1,Collim_Alb_Tot(nlevels_tot)
875	WRITE(6,1010) 'COLLIMABSCAN ',laieff_collim_1, &
876	laieff_collim_pft(nlevels_tot),3,&
877	leaf_single_scattering_albedo,br_base_temp_collim,&
878	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
879	Collim_Abs_Tot_1,SUM(Collim_Abs_Tot(:))
880	WRITE(6,1010) 'COLLIMABSSOIL ',laieff_collim_1, &
881	laieff_collim_pft(nlevels_tot),4,&
882	leaf_single_scattering_albedo,br_base_temp_collim,&
883	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
884	un-Collim_Alb_Tot_1-Collim_Abs_Tot_1,un-&
885	Collim_Alb_Tot(nlevels_tot)-SUM(Collim_Abs_Tot(:))
886
887	WRITE(6,1010) 'ISOTROPTRAN ',laieff_isotrop_1, &
888	laieff_isotrop_pft(nlevels_tot),1,&
889	leaf_single_scattering_albedo,br_base_temp_isotrop,&
890	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
891	Isotrop_Tran_Tot_1,(Isotrop_Tran_Coll(1)+Isotrop_Tran_Uncoll(1))
892	WRITE(6,1010) 'ISOTROPALB ',laieff_isotrop_1, &
893	laieff_isotrop_pft(nlevels_tot),2,&
894	leaf_single_scattering_albedo,br_base_temp_isotrop,&
895	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
896	Isotrop_Alb_Tot_1,Isotrop_Alb_Tot(nlevels_tot)
897	WRITE(6,1010) 'ISOTROPABSCAN ',laieff_isotrop_1, &
898	laieff_isotrop_pft(nlevels_tot),3,&
899	leaf_single_scattering_albedo,br_base_temp_isotrop,&
900	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
901	Isotrop_Abs_Tot_1,SUM(Isotrop_Abs_Tot(:))
902	WRITE(6,1010) 'ISOTROPABSSOIL ',laieff_isotrop_1, &
903	laieff_isotrop_pft(nlevels_tot),4,&
904	leaf_single_scattering_albedo,br_base_temp_isotrop,&
905	leaf_psd_temp,ACOS(cosine_sun_angle)/pi*180.0,&
906	un-Isotrop_Alb_Tot_1-Isotrop_Abs_Tot_1,un-&
907	Isotrop_Alb_Tot(nlevels_tot)-SUM(Isotrop_Abs_Tot(:))
908
909	1010 FORMAT(A,2(F12.6,2X),I4,2X,6(F14.8,2X))
910
911	ENDIF ! debugging
912
913	! For our total albedo, we take a weighted average of the
914	! diffuse and the direct light, which means we lose some
915	! information. This might be changed in the future. We also need
916	! to weight this by the percentage of nonbio land, but this
917	! appears to be included in veget_max
918	converged_albedo = direct_light_weight*Collim_Alb_Tot(nlevels_tot) + &
919	(un-direct_light_weight)*Isotrop_Alb_Tot(nlevels_tot)
920	albedo(ipts,ks) = albedo(ipts,ks) + &
921	veget_max(ipts,ivm)*converged_albedo
922	albedo_pft(ipts,ivm,ks)=veget_max(ipts,ivm)*converged_albedo
923
924	! Save the absorbed radiation for photosynthesis. We need only the
925	! visible range. We weight the isotropic and collimated parts
926	! according to an external parameter, but 0.5 is a sensible value
927	! if we don't actually know the fraction of light coming from direct
928	! and diffuse sources. This makes it consistent with the albedo
929	! calculated above. This can be changed later to distinguish between
930	! sunlit and shaded leaves.
931	! Notice that sunlight leaves will only get light from a Collimated
932	! source, whereas shaded leaves will get light from both direct
933	! sunlight (that has bounced off other leaves...Collim in this
934	! subroutine) and sunilght that has already reflected off of aerosols
935	! and clouds (Isotrop in this subroutine). Therefore, it is not a
936	! simple matter of taking Isotrop_Abs_Tot or Collim_Abs_Tot, but
937	! the Uncollided fraction of Collim_Abs_Tot (for the sunlit leaves)
938	! and the Collided fraction of Collum_Abs_Tot and the collided
939	! and uncollided fraction of Isotrop_Abs_Tot (for shaded leaves).
940	IF (ks == 1) THEN
941	! Test if Collim_Abs_Tot has negative values
942	! If so, set it to min_sechiba
943	DO ilevel=1,nlevels_tot
944	IF (Isotrop_Abs_Tot(ilevel) .LT. zero) THEN
945	Isotrop_Abs_Tot(ilevel) = min_sechiba
946	ENDIF
947	IF (Collim_Abs_Tot(ilevel) .LT. zero) THEN
948	Collim_Abs_Tot(ilevel) = min_sechiba
949	ENDIF
950	ENDDO
951	!
952	Light_Abs_Tot(ipts,ivm,:) = direct_light_weight*Collim_Abs_Tot(:) + &
953	(un-direct_light_weight)*Isotrop_Abs_Tot(:)
954	Light_Tran_Tot(ipts,ivm,:) = direct_light_weight*(Isotrop_Tran_Coll(:) + &
955	Isotrop_Tran_Uncoll(:)) + &
956	(un-direct_light_weight)*(Isotrop_Tran_Coll(:) + &
957	Isotrop_Tran_Uncoll(:))
958
959	! Notice that this light is actually cumulative, not per
960	! level! This was needed for debugging purposes and running
961	! tests. However, the photosynthesis routines need the light
962	! transmitted per level, i.e. if there is zero LAI
963	! in a level, the transmitted light will be one.
964	DO ilevel=1,nlevels_tot-1
965	IF(Light_Tran_Tot(ipts,ivm,ilevel+1) .GT. alb_threshold)THEN
966	Light_Tran_Tot(ipts,ivm,ilevel)=&
967	Light_Tran_Tot(ipts,ivm,ilevel) / &
968	Light_Tran_Tot(ipts,ivm,ilevel+1)
969	ELSE
970
971	! Here, we really don't know anything about how much
972	! light is transmitted in this layer, but there is no
973	! light reaching it from above so we can safely assume
974	! no photosynthesis takes place. This is equivalent
975	! to assuming it has no LAI, which means the
976	! transmission will be unity.
977	Light_Tran_Tot(ipts,ivm,ilevel)=un
978
979	ENDIF
980	ENDDO
981
982	! Debug
983	! Notice that the sum of the transmissions may not equal one. This
984	! is due to multiple scattering (light can be reflected up from a
985	! lower layer, and then back down).
986	IF (printlev_loc>=4 .AND. ivm == test_pft .AND. &
987	ipts == test_grid) THEN
988	WRITE(numout,*) 'Light values passed out of albedo_surface_main: ', ipts, ivm
989	WRITE(numout,'(7X,4(A15,3X))') ' Light_Abs_Tot',' Light_Tran_Tot',&
990	' laieff_isotrop',' laieff_collim'
991	DO ilevel=nlevels_tot,1,-1
992	WRITE(numout,'(I4,3X,4(F15.8,3X))') ilevel, &
993	Light_Abs_Tot(ipts,ivm,ilevel),&
994	Light_Tran_Tot(ipts,ivm,ilevel),&
995	laieff_isotrop(ipts,ilevel,ivm),&
996	laieff_collim(ipts,ilevel,ivm)
997	ENDDO
998
999	ENDIF
1000	!-
1001
1002	ENDIF ! IF ks==1
1003
1004	ENDDO ! ks = 1, n_spectralbands
1005
1006	ENDDO ! ivm = 2, nvm
1007
1008	ENDDO ! ipts = 1, kjpindex
1009
1010	! now we need to average the albedo over all our spectra, so we can
1011	! pass it to modules which do not distinguish between different
1012	! spectral bands
1013	albedo_glob(:) = zero
1014	DO ks = 1, n_spectralbands ! Loop over # of spectra
1015	albedo_glob(:) = albedo_glob(:) + albedo(:,ks)
1016	ENDDO
1017	albedo_glob(:)=albedo_glob(:)/REAL(n_spectralbands,r_std)
1018
1019	!+++CHECK+++
1020	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1021	! when running a larger domain. Needs to be corrected when implementing
1022	! a global use of the multi-layer energy budget.
1023	!!$ Light_Abs_Tot_mean(:) = Light_Abs_Tot_mean(:) / n_spectralbands
1024	!!$ Light_Alb_Tot_mean(:) = Light_Alb_Tot_mean(:) / n_spectralbands
1025	!+++++++++++
1026
1027	! Error checking
1028	IF (abs(albedo(1,1)) .gt. 1.0d0) THEN
1029
1030	CALL print_debugging_albedo_info(1, cosine_sun_angle, &
1031	leaf_reflectance, leaf_transmittance, laieff_collim_1, &
1032	laieff_isotrop_1, br_base_temp_collim, br_base_temp_isotrop, &
1033	Collim_Alb_Tot_1, Collim_Tran_Tot_1, Collim_Abs_Tot_1, &
1034	Isotrop_Alb_Tot_1, Isotrop_Tran_Tot_1, Isotrop_Abs_Tot_1)
1035	STOP
1036
1037	END IF
1038	!-
1039
1040	! Debug
1041	IF (printlev_loc>=4) THEN
1042	WRITE(numout,*) 'laieff_isotrop end albedo for test_pft',test_pft,'is', &
1043	laieff_isotrop(test_grid,:,test_pft)
1044	ENDIF
1045	!-
1046
1047	!! 5. Output diagnostics
1048	!! Return if current call is done from the initialization phase
1049	IF (firstcall_albedo_surface) RETURN
1050
1051	! Add XIOS default value where no sun
1052	DO ipts=1,kjpindex
1053	IF(coszang(ipts) .LE. min_sechiba)THEN
1054	albedo_diag(ipts,:) = xios_default_val
1055	alb_bare_diag(ipts,:) = xios_default_val
1056	albedo_snow_diag(ipts,:) = xios_default_val
1057	ELSE
1058	albedo_diag(ipts,:) = albedo(ipts,:)
1059	alb_bare_diag(ipts,:) = alb_bare(ipts,:)
1060	albedo_snow_diag(ipts,:) = albedo_snow(ipts,:)
1061	END IF
1062	END DO
1063
1064	IF (.NOT. impaze) THEN
1065	CALL xios_orchidee_send_field("soilalb_vis",alb_bare_diag(:,1))
1066	CALL xios_orchidee_send_field("soilalb_nir",alb_bare_diag(:,2))
1067	END IF
1068	CALL xios_orchidee_send_field("albedo_vis",albedo_diag(:,1))
1069	CALL xios_orchidee_send_field("albedo_nir",albedo_diag(:,2))
1070	CALL xios_orchidee_send_field("albedo_snow_vis",(albedo_snow_diag(:,1)))
1071	CALL xios_orchidee_send_field("albedo_snow_nir",(albedo_snow_diag(:,2)))
1072
1073	IF ( almaoutput ) THEN
1074	CALL histwrite_p(hist_id, 'Albedo', kjit, (albedo(:,1) + albedo(:,2))/2, kjpindex, index)
1075	CALL histwrite_p(hist_id, 'SAlbedo', kjit, (albedo_snow(:,1)+albedo_snow(:,2))/2, kjpindex, index)
1076	IF ( hist2_id > 0 ) THEN
1077	CALL histwrite_p(hist2_id, 'Albedo', kjit, (albedo(:,1) + albedo(:,2))/2, kjpindex, index)
1078	CALL histwrite_p(hist2_id, 'SAlbedo', kjit, (albedo_snow(:,1)+albedo_snow(:,2))/2, kjpindex, index)
1079	ENDIF
1080	ELSE
1081	IF (.NOT. impaze) THEN
1082	CALL histwrite_p(hist_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
1083	CALL histwrite_p(hist_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
1084	IF ( hist2_id > 0 ) THEN
1085	CALL histwrite_p(hist2_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
1086	CALL histwrite_p(hist2_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
1087	ENDIF
1088	END IF
1089	ENDIF
1090
1091	IF (printlev_loc>=3) WRITE(numout,*) 'Leaving albedo_surface_main'
1092
1093	END SUBROUTINE albedo_surface_main
1094
1095
1096	!! =============================================================================================================================
1097	!! SUBROUTINE : albedo_surface_finalize
1098	!!
1099	!>\BRIEF Write to restart file
1100	!!
1101	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in albedo
1102	!! to restart file
1103	!!
1104	!! RECENT CHANGE(S) : None
1105	!!
1106	!! REFERENCE(S) : None
1107	!!
1108	!! FLOWCHART : None
1109	!! \n
1110	!_ ==============================================================================================================================
1111	!+++CHECK+++
1112	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1113	! when running a larger domain. Needs to be corrected when implementing
1114	! a global use of the multi-layer energy budget.
1115	SUBROUTINE albedo_surface_finalize (kjit, kjpindex, rest_id, &
1116	albedo, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
1117	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
1118	laieff_isotrop)
1119	!+++++++++++
1120
1121	!! 0. Variable and parameter declaration
1122	!! 0.1 Input variables
1123	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
1124	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1125	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
1126	REAL(r_std), DIMENSION(kjpindex,n_spectralbands), INTENT(in) :: albedo !! Albedo (two stream radiation transfer model)
1127	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) :: Isotrop_Abs_Tot_p !! Absorbed radiation per layer for photosynthesis
1128	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot), INTENT(in) :: Isotrop_Tran_Tot_p !! Transmitted radiation per layer for photosynthesis
1129	!+++CHECK+++
1130	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1131	! when running a larger domain. Needs to be corrected when implementing
1132	! a global use of the multi-layer energy budget.
1133	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT(in) :: Light_Abs_Tot_mean !! total absorption for a given level * changed for enerbil integration
1134	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT(in) :: Light_Alb_Tot_mean !! total albedo for a given level * changed for enerbil integration
1135	!+++++++++++
1136	REAL(r_std), DIMENSION(kjpindex,nlevels_tot,nvm), INTENT(in) :: laieff_isotrop !! Leaf Area Index Effective converts
1137	!! 3D lai into 1D lai for two stream
1138
1139	!_ ================================================================================================================================
1140
1141	IF ( alb_bg_modis ) THEN
1142	CALL restput_p (rest_id, 'bckgrnd_alb', nbp_glo, 2, 1, kjit, bckgrnd_alb, 'scatter', nbp_glo, index_g)
1143	ELSE
1144	CALL restput_p (rest_id, 'soilalbedo_dry', nbp_glo, 2, 1, kjit, soilalb_dry, 'scatter', nbp_glo, index_g)
1145	!-
1146	CALL restput_p (rest_id, 'soilalbedo_wet', nbp_glo, 2, 1, kjit, soilalb_wet, 'scatter', nbp_glo, index_g)
1147	!-
1148	CALL restput_p (rest_id, 'soilalbedo_moy', nbp_glo, 2, 1, kjit, soilalb_moy, 'scatter', nbp_glo, index_g)
1149	ENDIF
1150
1151	CALL restput_p (rest_id, 'albedo' , nbp_glo, n_spectralbands, 1, kjit, albedo, "scatter", nbp_glo, index_g)
1152	!-
1153	CALL restput_p (rest_id, 'Isotrop_Abs_Tot_p', nbp_glo, nvm, nlevels_tot, kjit, Isotrop_Abs_Tot_p, "scatter", nbp_glo, index_g)
1154	!-
1155	CALL restput_p (rest_id, 'Isotrop_Tran_Tot_p', nbp_glo, nvm, nlevels_tot, kjit, Isotrop_Tran_Tot_p, "scatter", nbp_glo, index_g)
1156
1157	!+++CHECK+++
1158	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
1159	! when running a larger domain. Needs to be corrected when implementing
1160	! a global use of the multi-layer energy budget.
1161	!!$ IF (is_root_prc) THEN
1162	!!$ CALL restput (rest_id, 'Light_Abs_Tot_mean', nlevels_tot, 1, 1, kjit, Light_Abs_Tot_mean)
1163	!!$ CALL restput (rest_id, 'Light_Alb_Tot_mean', nlevels_tot, 1, 1, kjit, Light_Alb_Tot_mean)
1164	!!$ END IF
1165	!+++++++++++
1166
1167	CALL restput_p (rest_id, 'laieff_isotrop', nbp_glo, nlevels_tot, nvm, kjit, laieff_isotrop, "scatter", nbp_glo, index_g)
1168
1169
1170	END SUBROUTINE albedo_surface_finalize
1171
1172
1173	!! ==============================================================================================================================
1174	!! SUBROUTINE : albedo_surface_clear
1175	!!
1176	!>\BRIEF Deallocate albedo variables
1177	!!
1178	!! DESCRIPTION : None
1179	!!
1180	!! RECENT CHANGE(S): None
1181	!!
1182	!! MAIN OUTPUT VARIABLE(S): None
1183	!!
1184	!! REFERENCES : None
1185	!!
1186	!! FLOWCHART : None
1187	!! \n
1188	!_ ================================================================================================================================
1189
1190	SUBROUTINE albedo_surface_clear ()
1191
1192	IF (ALLOCATED (soilalb_dry)) DEALLOCATE (soilalb_dry)
1193	IF (ALLOCATED(soilalb_wet)) DEALLOCATE (soilalb_wet)
1194	IF (ALLOCATED(soilalb_moy)) DEALLOCATE (soilalb_moy)
1195	IF (ALLOCATED(bckgrnd_alb)) DEALLOCATE (bckgrnd_alb)
1196
1197	END SUBROUTINE albedo_surface_clear
1198
1199
1200	!! ==============================================================================================================================
1201	!! SUBROUTINE : twostream_solver
1202	!!
1203	!>\BRIEF : Computes the two-stream albedo solution for a level with given single scatterer properties
1204	!! and a defined background albedo
1205	!!
1206	!! DESCRIPTION : This solution is the two-stream solution of Pinty et al (2006) for a vegetation level above
1207	!! an isotropically reflecting background. It breaks down the problem into three parts, solving
1208	!! each part for the case of diffuse (isotropic) and direct (collimated) light. The three parts are,
1209	!! 1) The term due to light that does not interact at all with the canopy (Black Canopy)
1210	!! 2) Light that does not interact at all with the background (Black Background)
1211	!! 3) Light that bounces between the background and the canopy
1212	!! This routine (and the routines it uses) were received directly from Bernard Pinty, and only some
1213	!! minor modifcations were made, in addition to more documentation
1214	!!
1215	!! NOTE: the single layer solution is no longer used. The code is kept here in case it is decided that
1216	!! the multi layer solution is two expensive if only a single layer is used for the energy budget.
1217	!! In that case a solution needs to be implemeted to calculate the light distribution within
1218	!! the canopy.
1219	!!
1220	!! RECENT CHANGE(S): None
1221	!!
1222	!! MAIN OUTPUT VARIABLE(S): Collim_Alb_Tot,Collim_Tran_Tot,
1223	!! Collim_Abs_Tot,Isotrop_Alb_Tot,Isotrop_Tran_Tot,Isotrop_Abs_To
1224	!!
1225	!! REFERENCE(S) : B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski, N. Gobron and M. M. Verstraete (2006).
1226	!! Simplifying the Interaction of Land Surfaces with Radiation for Relating Remote Sensing Products to Climate Models.
1227	!! Journal of Geophysical Research. Vol 111, D02116.
1228	!!
1229	!! FLOWCHART : None
1230	!! \n
1231	!_ ================================================================================================================================
1232
1233	SUBROUTINE twostream_solver(leaf_reflectance, leaf_transmittance, background_reflectance_collim, background_reflectance_isotrop, &
1234	cosine_sun_angle, Collim_Alb_Tot, Collim_Tran_Tot, Collim_Abs_Tot, &
1235	Isotrop_Alb_Tot, Isotrop_Tran_Tot, Isotrop_Abs_Tot, laieff_collim, &
1236	laieff_isotrop, Collim_Tran_Uncollided, Isotrop_Tran_Uncollided)
1237
1238
1239	!! 0. Variables and parameter declaration
1240	!! 0.1 Input variables
1241	REAL(r_std), INTENT(IN) :: leaf_reflectance !! effective leaf reflectance, between 0 and 1
1242	!! @tex $()$ @endtex
1243	REAL(r_std), INTENT(IN) :: leaf_transmittance !! effective leaf transmittance, between 0 and 1
1244	!! @tex $()$ @endtex
1245	REAL(r_std), INTENT(IN) :: background_reflectance_collim !! background reflectance for direct radiation,
1246	!! between 0 and 1
1247	REAL(r_std), INTENT(IN) :: background_reflectance_isotrop !! background reflectance for diffuse radiation,
1248	!! between 0 and 1
1249	REAL(r_std), INTENT(IN) :: cosine_sun_angle !! cosine of the solar zenith angle, between 0 and 1
1250	!! @tex $()$ @endtex
1251	REAL(r_std), INTENT(IN) :: laieff_collim !! Effective Leaf Area Index, computed at current sun angle
1252	REAL(r_std), INTENT(IN) :: laieff_isotrop !! Effective Leaf Area Index, computed at 60 degrees
1253	!! @tex $(m^{2} m^{-2})$ @endtex
1254
1255	!! 0.2 Output variables
1256	!! notice that these variables (absorption + transmission + reflection) do not necessarily add up to 1 due to multiple scattering
1257	REAL(r_std), INTENT(OUT) :: Collim_Alb_Tot !! collimated total albedo from this level, between 0 and 1
1258	!! @tex $()$ @endtex
1259	REAL(r_std), INTENT(OUT) :: Collim_Tran_Tot !! collimated total transmission through this level, between 0 and 1
1260	!! @tex $()$ @endtex
1261	REAL(r_std), INTENT(OUT) :: Collim_Abs_Tot !! collimated total absorption by this level, between 0 and 1
1262	!! @tex $()$ @endtex
1263	REAL(r_std), INTENT(OUT) :: Isotrop_Alb_Tot !! isotropic total albedo from this level, between 0 and 1
1264	!! @tex $()$ @endtex
1265	REAL(r_std), INTENT(OUT) :: Isotrop_Tran_Tot !! isotropic total transmission through this level, between 0 and 1
1266	!! @tex $()$ @endtex
1267	REAL(r_std), INTENT(OUT) :: Isotrop_Abs_Tot !! isotropic total absorption by this level, between 0 and 1
1268	!! @tex $()$ @endtex
1269	REAL(r_std), INTENT(OUT) :: Collim_Tran_Uncollided !! collimated uncollied transmission through this level, between 0 and 1
1270	!! @tex $()$ @endtex
1271	REAL(r_std), INTENT(OUT) :: Isotrop_Tran_Uncollided !! isotropic uncollied transmission through this level, between 0 and 1
1272	!! @tex $()$ @endtex
1273
1274	!! 0.3 Modified variables
1275
1276	!! 0.4 Local variables
1277	LOGICAL :: ok
1278	REAL(r_std), DIMENSION(4) :: gammaCoeffs
1279	REAL(r_std), DIMENSION(4) :: gammaCoeffs_star
1280	REAL(r_std) :: tauprimetilde
1281	REAL(r_std) :: tauprimestar
1282	REAL(r_std) :: sun_zenith_angle_radians
1283	REAL(r_std), PARAMETER :: isotropic_cosine_constant=0.5/0.705
1284
1285	!calculated fluxes
1286
1287	REAL(r_std) :: Collim_Alb_BB
1288	REAL(r_std) :: Collim_Tran_BB
1289	REAL(r_std) :: Collim_Abs_BB
1290	REAL(r_std) :: Isotrop_Alb_BB
1291	REAL(r_std) :: Isotrop_Tran_BB
1292	REAL(r_std) :: Isotrop_Abs_BB
1293	REAL(r_std) :: Collim_Tran_BC
1294	REAL(r_std) :: Isotrop_Tran_BC
1295	REAL(r_std) :: Collim_Tran_TotalOneWay
1296	REAL(r_std) :: Isotrop_Tran_TotalOneWay
1297	REAL(r_std) :: Below_reinject_rad_collim,Below_reinject_rad_isotrop
1298
1299	!_ ================================================================================================================================
1300
1301	IF (printlev_loc>=3) WRITE(numout,*) 'Entering twostream_solver'
1302
1303	! convert angular values
1304
1305	sun_zenith_angle_radians = acos(cosine_sun_angle)
1306
1307	! calculate the 4 gamma coefficients both in isotropic and collimated illumination
1308
1309	call gammas(leaf_reflectance, leaf_transmittance, cosine_sun_angle, gammaCoeffs)
1310	call gammas(leaf_reflectance,leaf_transmittance,isotropic_cosine_constant,&
1311	gammaCoeffs_star)
1312
1313	! estimate the effective value of the optical thickness
1314
1315	tauprimetilde = 0.5_r_std * laieff_collim
1316	!tauprimetilde = 1.
1317
1318	! Should become zenit angle dependent (=60Â°) calculated by the Monte Carlo
1319	! photon shooting routine
1320
1321	tauprimestar = 0.5_r_std * laieff_isotrop
1322	!tauprimestar = 1.
1323
1324	! +++++++++++++ BLACK BACKGROUND ++++++++++++++++++++++
1325	! * Apply the black-background 2stream solution
1326	! * These equations are written for the part of the incoming radiation
1327	! * that never hits the background but does interact with the vegetation
1328	! * canopy.
1329	! *
1330	! * Note : the same routine dhrT1() is used for both the isotropic and
1331	! * collimated illumination conditions but the calling arguments differ.
1332	! * (especially the solar angle used).
1333	! */
1334
1335	! /* 1) collimated source */
1336	ok=dhrT1(leaf_reflectance,leaf_transmittance,&
1337	gammaCoeffs(1),gammaCoeffs(2),gammaCoeffs(3),gammaCoeffs(4),&
1338	sun_zenith_angle_radians, tauprimetilde,&
1339	Collim_Alb_BB,Collim_Tran_BB,Collim_Abs_BB)
1340	! /* 2) isotropic source */
1341	ok=dhrT1(leaf_reflectance,leaf_transmittance,&
1342	gammaCoeffs_star(1),gammaCoeffs_star(2),gammaCoeffs_star(3),&
1343	gammaCoeffs_star(4),&
1344	acos(isotropic_cosine_constant),tauprimestar,&
1345	Isotrop_Alb_BB,Isotrop_Tran_BB,Isotrop_Abs_BB)
1346
1347
1348	! +++++++++++++ BLACK CANOPY ++++++++++++++++++++++
1349	! * Apply the black canopy solution.
1350	! * These equations hold for the part of the incoming radiation
1351	! * that do not interact with the vegetation, travelling through
1352	! * its gaps.
1353	! */
1354
1355	! /* 1) collimated source */
1356	IF (cosine_sun_angle .NE. 0) THEN
1357	Collim_Tran_BC = exp( - tauprimetilde/cosine_sun_angle)
1358	Collim_Tran_Uncollided = Collim_Tran_BC
1359	ENDIF
1360
1361
1362	! /* 2) isotropic source */
1363	Isotrop_Tran_BC = TBarreUncoll_exact(tauprimestar)
1364	Isotrop_Tran_Uncollided = Isotrop_Tran_BC
1365
1366
1367	! /* Total one-way transmissions:
1368	! * The vegetation canopy is crossed (one way) by the uncollided radiation
1369	! * (black canopy) and the collided one (black background). */
1370
1371	! /* 1) collimated source */
1372	Collim_Tran_TotalOneWay = Collim_Tran_BC + Collim_Tran_BB
1373
1374	! /* 2) isotropic source */
1375	Isotrop_Tran_TotalOneWay = Isotrop_Tran_BC + Isotrop_Tran_BB
1376
1377	! * The Below_reinject_rad describes the process of reflecting toward the
1378	! * background the upward travelling radiation (re-emitted from below the
1379	! * canopy). It appears in the coupling equations as the limit of the series:
1380	! * 1 + rgrbv + (rgrbv)^2 + (rg*rbv)^3 + ...
1381	! * where rg is the background_reflectance and rbv is Isotrop_Alb_BB
1382	! * (with Isotrop describing the Lambertian reflectance of the background).
1383	! */
1384	! this might be involved in Eq. 27,
1385	Below_reinject_rad_collim = un / (un - background_reflectance_collim*Isotrop_Alb_BB)
1386	Below_reinject_rad_isotrop = un / (un - background_reflectance_isotrop*Isotrop_Alb_BB)
1387
1388	!/* TOTAL ALBEDO */
1389	! /* 1) collimated source */
1390	Collim_Alb_Tot = Collim_Alb_BB + &
1391	background_reflectance_collim * Collim_Tran_TotalOneWay * &
1392	Isotrop_Tran_TotalOneWay * Below_reinject_rad_collim
1393	! /* 2) isotropic source */
1394	Isotrop_Alb_Tot = Isotrop_Alb_BB + &
1395	background_reflectance_isotrop * Isotrop_Tran_TotalOneWay * &
1396	Isotrop_Tran_TotalOneWay * Below_reinject_rad_isotrop ! this seems to be
1397	! exactly Eq. 33
1398
1399
1400	!/* TOTAL TRANSMITION TO THE BACKGROUND LEVEL */
1401	! /* 1) collimated source */
1402	Collim_Tran_Tot = Collim_Tran_TotalOneWay * Below_reinject_rad_collim ;
1403
1404	! /* 2) isotropic source */
1405	Isotrop_Tran_Tot = Isotrop_Tran_TotalOneWay * Below_reinject_rad_isotrop ;
1406
1407	!/* TOTAL ABSORPTION BY THE VEGETATION LEVEL */
1408	! /* 1) collimated source */
1409	Collim_Abs_Tot = un - (Collim_Tran_Tot + Collim_Alb_Tot) + &
1410	background_reflectance_collim * Collim_Tran_Tot;
1411	! /* 2) isotropic source */
1412	Isotrop_Abs_Tot = un - (Isotrop_Tran_Tot + Isotrop_Alb_Tot) + &
1413	background_reflectance_isotrop * Isotrop_Tran_Tot;
1414
1415	!+++ TEMP +++
1416	! Combine the collimated & isotrop together as Collimate for energybudget
1417	Collim_Alb_Tot = Collim_Alb_Tot + Isotrop_Alb_Tot
1418	Collim_Abs_Tot = Collim_Abs_Tot + Isotrop_Abs_Tot
1419	!+++ TEMP +++
1420
1421
1422	IF (printlev_loc>=3) WRITE(numout,*) 'Exiting twostream_solver'
1423
1424	END SUBROUTINE twostream_solver
1425
1426
1427	!! ============================================================================\n
1428	!! FUNCTION : dhrT1
1429	!!
1430	!>\BREF
1431	!!
1432	!! DESCRIPTION : Taken directly from the Fortan code provided by Pinty et al for their \n
1433	!! two-stream model. The only changes are the REAL types. This is used
1434	!! to compute the black background absorption, transmission, and reflectance
1435	!! in Pintys scheme and it's based on the older model of Meador and Weaver...
1436	!! Pinty 2006 also includes a short discussion of it in Appendix A
1437	!!
1438	!! RECENT CHANGE(S): None\n
1439	!!
1440	!! RETURN VALUE : dhrT1
1441	!!
1442	!! REFERENCE(S) : Meador and Weaver, 'Two-stream approximations to radiative transfer in
1443	!! planetary atmosphere: a unified description of existing methods and a new improvement'
1444	!! J. Atmospheric Sciences, VOL 37, p. 630--643, 1980.
1445	!!
1446	!! FLOWCHART : None
1447	!_ =============================================================================
1448
1449	FUNCTION dhrT1(rl,tl,gamma1,gamma2,gamma3,gamma4,tta,tau,AlbBS,Tdif,AbsVgt)
1450
1451
1452	!! 0. Variables and parameter declaration
1453	!! 0.1 Input variables
1454	REAL(r_std), INTENT(IN) :: rl ! the effective reflectance of a single scatterer, between 0 and 1
1455	!! @tex $()$ @endtex
1456	REAL(r_std), INTENT(IN) :: tl ! the effective transmittance of a single scatterer, between 0 and 1
1457	!! @tex $()$ @endtex
1458	REAL(r_std), INTENT(IN) :: gamma1 ! the gamma coefficients from Meador and Weaver
1459	REAL(r_std), INTENT(IN) :: gamma2
1460	REAL(r_std), INTENT(IN) :: gamma3
1461	REAL(r_std), INTENT(IN) :: gamma4
1462	REAL(r_std), INTENT(IN) :: tta ! the solar zenith angle, between 0 and pi/2
1463	!! @tex $(radians)$ @endtex
1464
1465	REAL(r_std), INTENT(IN) :: tau ! Effective LAI * G(theta)
1466	!! @tex $()$ @endtex
1467
1468	!! 0.2 Output variables
1469	REAL(r_std), INTENT(OUT) :: AlbBS ! the albedo of level, between 0 and 1
1470	!! @tex $()$ @endtex
1471	REAL(r_std), INTENT(OUT) :: Tdif ! the transmitted light (all diffuse) from this light source, between 0 and 1
1472	!! @tex $()$ @endtex
1473	REAL(r_std), INTENT(OUT) :: AbsVgt ! the light absorbed by the level, between 0 and 1
1474	!! @tex $()$ @endtex
1475	LOGICAL :: dhrT1 ! the output flag, which always appears to be true
1476
1477	!! 0.3 Modified variables
1478	!! 0.4 Local variables
1479
1480	REAL(r_std) :: alpha1,alpha2,ksquare,k
1481	REAL(r_std) :: first_term,secnd_term1,secnd_term2,secnd_term3
1482	REAL(r_std) :: expktau,Tdir
1483	REAL(r_std) :: mu0,w0
1484
1485	!_ ================================================================================================================================
1486
1487	mu0=cos(tta)
1488	w0=rl+tl
1489
1490
1491	Tdir = exp(-tau/mu0) ! direct transmission
1492
1493	! There is a difference between conservative and non-conservative scattering
1494	! conditions */
1495	IF (w0 .ne. 1.0 .AND. w0 .ne. 0.0) THEN
1496	!NON_CONSERVATIVE SCATTERING
1497
1498	! some additional parameters, taken from Eq. 16--18 of Meador and Weaver, J. Atmospheric Sciences,
1499	! Vol 37, p. 630--643 (1980)
1500
1501	alpha1 = gamma1gamma4 + gamma2gamma3
1502	alpha2 = gamma1gamma3 + gamma2gamma4
1503	ksquare = gamma1gamma1 - gamma2gamma2
1504	k = sqrt(ksquare)
1505
1506	expktau = exp(k*tau)
1507
1508	!Black Soil Albedo...Eq. 14 in Meador and Weaver
1509	first_term = ((un-ksquaremu0mu0)((k+gamma1)expktau + (k-gamma1)/expktau))
1510	IF (first_term .eq. 0.0) THEN
1511	!we will be dividing by zero : cannot continue.
1512	dhrT1 = .false.
1513	ELSE
1514	first_term = un/first_term
1515	secnd_term1 = (un - kmu0)(alpha2 + kgamma3)expktau
1516	secnd_term2 = (un + kmu0)(alpha2 - k*gamma3)/expktau
1517	secnd_term3 = 2.0_r_std * k * (gamma3 - alpha2mu0)Tdir
1518
1519	AlbBS = (w0 * first_term * (secnd_term1 - secnd_term2 - secnd_term3))
1520
1521	!Transmission...Eq. 15 in Meador and Weaver, for diffuse light?
1522	IF (ksquare .eq. 0.0) THEN
1523	first_term = un
1524	ENDIF
1525	secnd_term1 = (un+kmu0)(alpha1+kgamma4)expktau
1526	secnd_term2 = (un-kmu0)(alpha1-k*gamma4)/expktau
1527	secnd_term3 = 2.0_r_std * k * (gamma4 + alpha1*mu0)
1528	Tdif = - w0first_term(Tdir*(secnd_term1 - secnd_term2) - secnd_term3)
1529
1530	! Absorption by vegetation...whatever is not transmitted or reflected
1531	! must be absorbed
1532	AbsVgt = (un- (Tdif+Tdir) - AlbBS)
1533	ENDIF ! first_term .eq. 0.0
1534	ELSE IF (w0 .eq. 0.) THEN
1535	!BLACK CANOPY
1536	AlbBS = zero
1537	Tdif = zero
1538	AbsVgt = un - Tdir
1539	ELSE
1540	!CONSERATIVE SCATTERING...Eq. 24 in Meador and Weaver
1541	AlbBS = (un/(un + gamma1tau))(gamma1tau + (gamma3-gamma1mu0)*&
1542	(un-exp(-tau/mu0)));
1543	Tdif = un - AlbBS - Tdir;
1544	AbsVgt = zero;
1545	ENDIF ! w0 .ne. 1.0 .AND. w0 .ne. 0.0
1546
1547	! not sure what the purpose of this flag is, as it will always be set to true here
1548	dhrT1 = .true.
1549
1550	END FUNCTION dhrT1
1551
1552	!! ============================================================================\n
1553	!! FUNCTION : TBarreUncoll_exact
1554	!!
1555	!>\BRIEF Computes the transmission of the diffuse black canopy light
1556	!!
1557	!! DESCRIPTION : Taken directly from the Fortan code provided by Pinty et al for their \n
1558	!! two-stream model. The only changes are the REAL types. This appears
1559	!! to be solving Eq. 16 in Pinty 2006, which is the transmission which
1560	!! does not collide with the canopy
1561	!!
1562	!! RECENT CHANGE(S): None\n
1563	!!
1564	!! RETURN VALUE : TBarreUncoll_exact
1565	!!
1566	!! REFERENCE(S) : None
1567	!!
1568	!! FLOWCHART : None
1569	!_ =============================================================================
1570
1571	FUNCTION TBarreUncoll_exact(tau)
1572
1573
1574	!! 0. Variables and parameter declaration
1575
1576	!! 0.1 Input variables
1577	REAL(r_std), INTENT(IN) :: tau ! Effective LAI * G(theta)
1578	!! @tex $()$ @endtex
1579
1580	!! 0.2 Output variables
1581	REAL(r_std) :: TBarreUncoll_exact ! the isotropic light transmission uncollided
1582	! with the canopy, between 0 and 1
1583	!! @tex $()$ @endtex
1584
1585
1586	!! 0.3 Modified variables
1587
1588	!! 0.4 Local variables
1589	! INTEGER :: j,ind
1590	! REAL(r_std) :: iGammaloc
1591	! INTEGER :: order
1592
1593	!_ ================================================================================================================================
1594
1595	!+++++ CHECK +++++
1596
1597	!!$ iGammaloc = zero
1598	!!$ order=20
1599	!!$
1600	!!$ ! in the case where tau is equal to zero, this crashes in the first loop where
1601	!!$ ! iGammaloc is also zero...quick fix, and might even be accurate
1602	!!$
1603	!!$ IF(tau .GT. 1e-10_r_std)THEN
1604	!!$
1605	!!$ DO j=0,order-1
1606	!!$ ind = order - j
1607	!!$ iGammaloc = ind / (un + ind/(tau+iGammaloc))
1608	!!$ END DO
1609	!!$ iGammaloc=un/(tau + iGammaloc)
1610	!!$ ENDIF
1611	!!$
1612	!!$ TBarreUncoll_exact = exp(-tau)(un - tau + tautau*iGammaloc)
1613
1614	! this is a change suggested by Bernard to improve the matching between the one
1615	! level and multilevel case
1616
1617	TBarreUncoll_exact = exp(-tau2.0_r_std0.705_r_std)
1618	!++++++++++
1619
1620	END FUNCTION TBarreUncoll_exact
1621
1622	!! ==============================================================================================================================
1623	!! SUBROUTINE : gammas
1624	!!
1625	!>\BRIEF Computes a set of gamma coefficients for use in the black background equations
1626	!!
1627	!! DESCRIPTION : Taken directly from the Fortan code provided by Pinty et al for their \n
1628	!! two-stream model. The only changes are the REAL types. This calculates
1629	!! the gamma coefficients based on Appendix A in the paper. This seems to
1630	!! make the assumption of the Ross's G function and a spherical leaf
1631	!! angle distribution.
1632	!!
1633	!! RECENT CHANGE(S): None\n
1634	!!
1635	!! MAIN OUTPUT VARIABLE(S): gammaCoeff
1636	!!
1637	!! REFERENCE(S) : B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski, N. Gobron and M. M. Verstraete (2006).
1638	!! Simplifying the Interaction of Land Surfaces with Radiation for Relating Remote Sensing Products to Climate Models.
1639	!! Journal of Geophysical Research. Vol 111, D02116.
1640	!!
1641	!! FLOWCHART : None
1642	!_ ================================================================================================================================
1643
1644	SUBROUTINE gammas(rl, tl, mu0, gammaCoeff)
1645
1646
1647	!! 0. Variables and parameter declaration
1648	!! 0.1 Input variables
1649	REAL(r_std), INTENT(IN) :: rl ! the effective reflectance of a single scatterer, between 0 and 1
1650	!! @tex $()$ @endtex
1651	REAL(r_std), INTENT(IN) :: tl ! the effective transmittance of a single scatterer, between 0 and 1
1652	!! @tex $()$ @endtex
1653	REAL(r_std), INTENT(IN) :: mu0 ! the cosine of the solar zenith angle, between 0 and 1
1654	!! @tex $()$ @endtex
1655
1656	!! 0.2 Output variables
1657	REAL(r_std), DIMENSION(1:4), INTENT(OUT) :: gammaCoeff ! the set of gamma coefficients in the above reference
1658	!! @tex $()$ @endtex
1659
1660	!! 0.3 Modified variables
1661
1662	!! 0.4 Local variables
1663	REAL(r_std) :: wd,w0,w0half,wdsixth
1664
1665	!_ ================================================================================================================================
1666
1667
1668	w0 = rl+tl
1669	wd = rl-tl
1670	w0half = w0*0.5_r_std
1671	wdsixth = wd/6.0_r_std
1672
1673	gammaCoeff(1)=2._r_std*(un - w0half + wdsixth)
1674	gammaCoeff(2)=2._r_std*(w0half + wdsixth)
1675	gammaCoeff(3)=2._r_std(w0half0.5_r_std + mu0*wdsixth)/w0
1676	gammaCoeff(4)=un-gammaCoeff(3)
1677
1678	END SUBROUTINE gammas
1679
1680
1681	!! ==============================================================================================================================
1682	!! SUBROUTINE : calculate_snow_albedo
1683	!!
1684	!>\BRIEF Computes some of the information needed to calculate the effect of the snow albedo
1685	!! on the background reflectance in the two stream model. Right now, this can be done
1686	!! with the old snow albedo scheme from Krinner et al 2005 as well as the snow albedo
1687	!! scheme from Community Land Model 3 (CLM3). The calculation of
1688	!! snow cover fraction is taken from Yang et al. 1997
1689	!!
1690	!! DESCRIPTION : In order to compute the background albedo in Pinty's two stream model, we have to take
1691	!! into account any snow that has fallen through the canopy and landed on the ground. In particular,
1692	!! we need the amount of ground covered by the snow and the albedo of this snow. Both of these
1693	!! quantities are calculated here, using a function that changes the albedo of the snow based on
1694	!! its age. This is done in two ways, but the model of CLM3 is better suited to be used with Pinty's
1695	!! albedo scheme, as it divides the radiation into VIS and NIR light. The calculation of snow cover
1696	!! fraction depends on roughness lenght.
1697	!!
1698	!! RECENT CHANGE(S): None\n
1699	!!
1700	!! MAIN OUTPUT VARIABLE(S): ::snowa_veg, ::frac_snow_veg, ::albedo_snow
1701	!!
1702	!! REFERENCE(S) : "Technical description of the community land model CLM)", K. W. Oleson, Y. Dai, G. Bonan, M. Bosilovich, et al.,
1703	!! NCAR Technical Note, May 2004.
1704	!!
1705	!! Yang Z, Dickinson R, Robock A, Vinnikov K (1997) Validation of the snow submodel of the
1706	!! biosphere-atmosphere transfer scheme with Russian
1707	!! snow cover and meteorological observational data. Journal of Climate, 10, 353â373.
1708	!!
1709	!! FLOWCHART : None
1710	!_ ================================================================================================================================
1711
1712	SUBROUTINE calculate_snow_albedo(kjpindex, coszang, snow, snow_nobio, &
1713	snow_age, snow_nobio_age, frac_nobio, albedo_snow, &
1714	snowa_veg, frac_snow_veg, snowa_nobio, frac_snow_nobio, &
1715	veget_max, z0m, veget)
1716
1717	!! 0. Variables and parameter declaration
1718	!! 0.1 Input variables
1719	INTEGER,INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1720	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: coszang !! cosine of the solar zenith angle, between 0 and 1
1721	!! @tex $()$ @endtex
1722	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
1723	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
1724	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age (days)
1725	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
1726	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
1727	!! continental ice, lakes, etc. (unitless)
1728	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT coverage fraction of a PFT (= ind*cn_ind)
1729	!! (m^2 m^{-2})
1730	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
1731	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: z0m !! Roughness height for momentum of vegetated part (m)
1732	REAL(r_std), DIMENSION(:),INTENT(in) :: frac_snow_veg !! The fraction of the surface for each PFT covered by snow
1733	REAL(r_std), DIMENSION(:,:),INTENT(in) :: frac_snow_nobio !! The fraction of nonbiological land types covered by snow
1734
1735	!! 0.2 Output variables
1736	REAL(r_std), DIMENSION (kjpindex,n_spectralbands), &
1737	INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
1738	REAL(r_std), DIMENSION(kjpindex,nvm,n_spectralbands),&
1739	INTENT(out) :: snowa_veg !! the albedo of snow...this is seperated into
1740	!! PFT types, even though the calculation is identical
1741	!! for all PFTs right now
1742	REAL(r_std), DIMENSION(kjpindex,nnobio,n_spectralbands),&
1743	INTENT(out) :: snowa_nobio !! the albedo of snow on nonbiological land types
1744
1745
1746	!! 0.3 Modified variables
1747
1748	!! 0.4 Local variables
1749	REAL(r_std), DIMENSION(kjpindex) :: agefunc_veg
1750	REAL(r_std), DIMENSION(kjpindex,nnobio) :: agefunc_nobio
1751	REAL(r_std) :: snow_alb_direct, snow_alb_diffuse,f_mu
1752	REAL(r_std), DIMENSION(kjpindex) :: tot_bare_soil_notree !! Total bare soil fraction without accounting for Trees
1753	REAL(r_std), DIMENSION(kjpindex) :: fraction_veg !! fraction of the pixel with vegetation (bare soil is considered vegetation)
1754
1755	INTEGER :: ipts,ks,ivm,jv ! indices
1756
1757	REAL(r_std), DIMENSION (nvm,2) :: snowa_aged_tmp !! spectral domains (unitless)
1758	REAL(r_std), DIMENSION (nvm,2) :: snowa_dec_tmp
1759
1760	!_ ================================================================================================================================
1761
1762	! initialize the output
1763	albedo_snow(:,:) = zero
1764	snowa_veg(:,:,:) = val_exp
1765	snowa_nobio(:,:,:) = val_exp
1766	fraction_veg(:) = un - SUM(frac_nobio(:,:),2)
1767	tot_bare_soil_notree(:) = zero
1768	snowa_aged_tmp(:,ivis) = snowa_aged_vis(:)
1769	snowa_aged_tmp(:,inir) = snowa_aged_nir(:)
1770	snowa_dec_tmp(:,ivis) = snowa_dec_vis(:)
1771	snowa_dec_tmp(:,inir) = snowa_dec_nir(:)
1772
1773	DO ipts = 1, kjpindex ! loop over all the grid squares
1774
1775	IF (SUM(veget_max(ipts,:)).LT.min_stomate) THEN
1776	! The pixel is inconsistent. It is present in climate forcing
1777	! but it is absent in the PFT map. If we do not want the model
1778	! to crash or generate NaNs. Some variables need to be initialized
1779	agefunc_nobio(ipts,:) = un
1780	END IF
1781
1782	DO ivm=1,nvm
1783	! Use tot_bare_soil_notree to quantify the total bare soil
1784	! fraction outside the forest PFTs
1785	IF ( fraction_veg(ipts) .GT. min_sechiba .AND. &
1786	.NOT. is_tree(ivm) ) THEN
1787	tot_bare_soil_notree(ipts) = tot_bare_soil_notree(ipts) + &
1788	veget_max(ipts,ivm) - veget(ipts,ivm)
1789	ENDIF
1790	END DO
1791
1792	DO ivm=1,nvm
1793
1794	IF(veget_max(ipts,ivm) == zero)THEN
1795	! No vegetation or bare soil present, so no reason to do the
1796	! calculation
1797	CYCLE
1798	ENDIF
1799
1800	DO ks=1,n_spectralbands ! loop over spectra
1801
1802	! The snow albedo could be either prescribed or
1803	! calculated following Chalita and Treut (1994)(ok_snow_albedo_clm3 = n)
1804	! or calculated following CLM3 (ok_snow_albedo_clm3 = y). NOTE: The
1805	! difference between these two methods has not been tested.
1806
1807	! Check if the precribed value fixed_snow_albedo exists.
1808	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
1809
1810	snowa_veg(ipts,ivm,ks) = fixed_snow_albedo
1811	snowa_nobio(ipts,ivm,ks) = fixed_snow_albedo
1812
1813	ELSE
1814
1815	! Calculate age dependence
1816	! On vegetated surfaces
1817	agefunc_veg(ipts) = EXP(-snow_age(ipts)/tcst_snowa)
1818
1819	! On non-vegtative surfaces
1820	DO jv = 1, nnobio ! Loop over # nobio types
1821	agefunc_nobio(ipts,jv) = EXP(-snow_nobio_age(ipts,jv)/tcst_snowa)
1822	ENDDO
1823
1824	IF(ok_snow_albedo_clm3)THEN
1825
1826	! Albedo of the snow under the vegetation. These
1827	! values are used to calculate the background albedo which
1828	! in turn is used in the two-stream solver to calculate
1829	! the radiative transfer through the canopy.
1830	! First we calculate diffuse albedo using the constant from
1831	! CLM3 conveted into days**-1
1832	!agefunc=un-un/(un+snow_age(ipts)*0.864_r_std)
1833	!agefunc=un/(un+snow_age(ipts)*0.864_r_std)
1834	snow_alb_diffuse = (un-c_albedo(ks) * &
1835	(un-agefunc_veg(ipts)))*alb_snow_0(ks)
1836
1837	! now the direct albedo
1838	IF(coszang(ipts) .GT. 0.5_r_std)THEN
1839	f_mu=zero
1840	ELSE
1841	! substituting b=2.0 into the equation...recommended by CLM3
1842	f_mu=(3.0_r_std)/(un+coszang(ipts))-2.0_r_std
1843	ENDIF
1844	snow_alb_direct=snow_alb_diffuse+0.4_r_stdf_mu&
1845	(un-snow_alb_diffuse)
1846
1847	! assume that the total snow albedo is a mix of
1848	! 50% direct and 50% diffuse
1849	snowa_veg(ipts,ivm,ks) = (snow_alb_direct+snow_alb_diffuse) / 2.0_r_std
1850
1851	ELSE
1852
1853	!+++CHECK+++
1854	! Two simulations of north America, one with CLM3 and one with
1855	! the other albedo scheme result in very different snowalbedo
1856	! estimates. The CLM3 scheme results in OK estimates. The
1857	! original scheme contains bugs. Only very few pixels (even
1858	! over Greenland have non-zero value for snow albedo.
1859	CALL ipslerr_p(3,'calculate_snow_albedo','This code compiles and runs',&
1860	'but the results have never been checked after the merge',&
1861	'remove this stop and check the results')
1862
1863	! Snow albedo following ORCHIDEE 3.0. Albedo of the snow
1864	! under the vegetation. These values are used to calculate
1865	! the background albedo which in turn is used in the
1866	! two-stream solver to calculate the radiative transfer
1867	! through the canopy
1868	snowa_veg(ipts,ivm,ks) = zero
1869
1870	! The new rational for the snow albedo and its dependency to
1871	! snow age is (see ticket 223):
1872	! i) Short vegetation (grass, crop): we assume that when there
1873	! are no leaves (LAI=0), the snow aging is not impacted by
1874	! the vegetation. The difference between veget and vegetmax
1875	! (bare soil fraction within the PFT) should be considered in
1876	! the same way as the bare soil fraction (PFT1).
1877	! ii) High vegetation (trees): In this case we consider that
1878	! the vegetation impacts the snow aging even when LAI=0
1879	! because of the trunks, branches. So these PFTs should be
1880	! treated differently without any consideration of the bare
1881	! soil fraction within the PFT.
1882	IF ( ivm .EQ. ibare_sechiba) THEN
1883
1884	IF( fraction_veg(ipts) .GT. min_sechiba ) THEN
1885
1886	! Bare soil use tot_bare_soil_notree
1887	snowa_veg(ipts,ivm,ks) = tot_bare_soil_notree(ipts) / &
1888	fraction_veg(ipts) * &
1889	( snowa_aged_tmp(1,ks)+snowa_dec_tmp(1,ks) * &
1890	agefunc_veg(ipts) )
1891
1892	ENDIF
1893
1894	ELSE
1895
1896	! Vegetated PFTs
1897	IF ( is_tree(ivm) .AND. fraction_veg(ipts) .GT. min_sechiba ) THEN
1898	! For forested PFTs use veget_max
1899	snowa_veg(ipts,ivm,ks) = veget_max(ipts,ivm) / &
1900	fraction_veg(ipts) * &
1901	( snowa_aged_tmp(ivm,ks)+snowa_dec_tmp(ivm,ks) * &
1902	agefunc_veg(ipts) )
1903
1904	ELSEIF (fraction_veg(ipts) .GT. min_sechiba ) THEN
1905
1906	! For grasses and crops use veget
1907	snowa_veg(ipts,ivm,ks) = veget(ipts,ivm) / &
1908	fraction_veg(ipts) * &
1909	( snowa_aged_tmp(ivm,ks)+snowa_dec_tmp(ivm,ks) * &
1910	agefunc_veg(ipts) )
1911	ENDIF
1912
1913	ENDIF
1914	!++++++++++
1915
1916	ENDIF ! control ok_snow_albedo_clm3
1917
1918	ENDIF ! prescribe or calculate albedo
1919
1920	! Albedo of the snow under the vegetation. Notice that the fraction
1921	! of nobio land is accounted for in veget_max. This value is only
1922	! used as an output variable
1923	albedo_snow(ipts,ks) = albedo_snow(ipts,ks) + &
1924	frac_snow_veg(ipts) * veget_max(ipts,ivm) * &
1925	snowa_veg(ipts,ivm,ks)
1926
1927	ENDDO ! ks=1,n_spectralbands
1928
1929	ENDDO ! ivm=1,nvm
1930
1931	! NIR and VIS albedo for the snow on bio surfaces
1932	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
1933
1934	snowa_nobio(ipts,:,:) = fixed_snow_albedo
1935
1936	ELSE
1937
1938	DO ks = 1, n_spectralbands
1939
1940	DO jv = 1, nnobio
1941
1942	! The albedo due to snow of this age on this nonbio surface
1943	! This value is used in albedo_main to calculate
1944	! the albedo of the entire pixel (= vegetation + nobio + bare soil)
1945	snowa_nobio(ipts,jv,ks) = ( snowa_aged_tmp(1,ks) + &
1946	snowa_dec_tmp(1,ks) * agefunc_nobio(ipts,jv) )
1947
1948	ENDDO
1949
1950	ENDDO
1951
1952	ENDIF ! prescribe or calculate
1953
1954	! NIR and VIS albedo for the snow on non-biological surfaces
1955	DO ks = 1, n_spectralbands
1956
1957	DO jv = 1, nnobio
1958
1959	! This value is only used as an output variable
1960	albedo_snow(ipts,ks) = albedo_snow(ipts,ks) + &
1961	( frac_nobio(ipts,jv) ) * ( frac_snow_nobio(ipts,jv) ) * &
1962	snowa_nobio(ipts,jv,ks)
1963
1964	ENDDO ! jv = 1, nnobio
1965
1966	ENDDO ! spectralbands
1967
1968	ENDDO ! ipts = 1, kjpindex
1969
1970	END SUBROUTINE calculate_snow_albedo
1971
1972	!! ==============================================================================================================================
1973	!! SUBROUTINE : optimize_albedo_values
1974	!!
1975	!>\BRIEF Follow the radiation scattered through the canopy in the
1976	!! case of multiple levels.
1977	!!
1978	!! DESCRIPTION : Right now, we know that using Pintys method in the case of a
1979	!! single canopy level will result in different top of the canopy
1980	!! albedos than iwhat is found if we use multiple canopy levels.
1981	!! We trust that the single level case gives the 'true' values.
1982	!!
1983	!! This algorithm follows light as it scatters through multiple
1984	!! levels in the canopy. At each level, the scattering
1985	!! solution is found by solving Pinty's two stream model. This
1986	!! means that light which enters a level can either be transmitted
1987	!! without colliding with the vegetation, transmitted after
1988	!! collision with the vegetation, reflecting off the vegetation,
1989	!! or being absorbed. We follow the fluxes until they get
1990	!! really small.
1991	!!
1992	!! RECENT CHANGE(S): None\n
1993	!!
1994	!! MAIN OUTPUT VARIABLE(S): ::lconverged, ::Collim_Alb_Tot, ::Collim_Tran_Tot,
1995	!! ::Collim_Abs_Tot, ::Isotrop_Alb_Tot, ::Isotrop_Tran_Tot, ::Isotrop_Abs_Tot
1996	!!
1997	!! REFERENCE(S) : B. Pinty, T. Lavergne, R.E. Dickinson, J-L. Widlowski, N. Gobron
1998	!! and M. M. Verstraete (2006). Simplifying the Interaction of Land
1999	!! Surfaces with Radiation for Relating Remote Sensing Products to
2000	!! Climate Models. Journal of Geophysical Research. Vol 111, D02116.
2001	!!
2002	!! FLOWCHART : None
2003	!_ ================================================================================================================================
2004	SUBROUTINE multilevel_albedo(cosine_sun_angle, leaf_single_scattering_albedo_start,&
2005	leaf_psd_start, br_base_temp_collim, br_base_temp_isotrop, &
2006	laieff_collim_pft, laieff_isotrop_pft, lconverged, &
2007	Collim_Alb_Coll, Collim_Tran_Coll, Collim_Abs_Tot, Isotrop_Alb_Coll, &
2008	Isotrop_Tran_Coll, Isotrop_Abs_Tot, Collim_Tran_Uncoll, Isotrop_Tran_Uncoll, lprint)
2009
2010	!! 0. Variables and parameter declaration
2011	!! 0.1 Input variables
2012	REAL(r_std), INTENT(IN) :: cosine_sun_angle !! cosine of the solar zenith angle, between 0 and 1
2013	!! @tex $()$ @endtex
2014	REAL(r_std), INTENT(IN) :: leaf_single_scattering_albedo_start !! cosine of the solar zenith angle, between 0 and 1
2015	REAL(r_std), INTENT(IN) :: leaf_psd_start !! cosine of the solar zenith angle, between 0 and 1
2016	REAL(r_std), INTENT(IN) :: br_base_temp_collim !! cosine of the solar zenith angle, between 0 and 1
2017	REAL(r_std), INTENT(IN) :: br_base_temp_isotrop !! cosine of the solar zenith angle, between 0 and 1
2018	REAL(r_std),DIMENSION(:),INTENT(IN) :: laieff_collim_pft !! Effective lai for a single pixel and pft to be
2019	REAL(r_std),DIMENSION(:),INTENT(IN) :: laieff_isotrop_pft !! Effective lai for a single pixel and pft to be
2020	LOGICAL,INTENT(IN) :: lprint !! A flag to print some debug statements
2021	!! used in the two stream approach...every level
2022
2023	!! 0.2 Output variables
2024	LOGICAL,INTENT(OUT) :: lconverged !! did the optimization converge?
2025	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Collim_Alb_Coll !! collimated total albedo for the converged case
2026	!! unitless, between 0 and 1
2027	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Collim_Tran_Coll !! collimated total transmission
2028	!! unitless, between 0 and 1
2029	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Collim_Abs_Tot !! collimated total absorption
2030	!! unitless, between 0 and 1
2031	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Isotrop_Alb_Coll !! isotropic total albedo
2032	!! unitless, between 0 and 1
2033	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Isotrop_Tran_Coll !! isotropic total transmission
2034	!! unitless, between 0 and 1
2035	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Isotrop_Abs_Tot !! isotropic total absorption
2036	!! unitless, between 0 and 1
2037	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Collim_Tran_Uncoll !! collimated uncollided transmission
2038	!! unitless, between 0 and 1
2039	REAL(r_std),DIMENSION(:),INTENT(OUT) :: Isotrop_Tran_Uncoll !! isotropic uncollided transmission
2040	!! unitless, between 0 and 1
2041
2042	!! 0.3 Modified variables
2043
2044	!! 0.4 Local variables
2045	! use an extra level here...this is basically to store the transmission from the sun, which is normalized to 1.0
2046	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Abs_Coll_Unscaled
2047	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Alb_Coll_Unscaled
2048	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Tran_Coll_Unscaled
2049	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Tran_UnColl_Unscaled
2050	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Abs_Coll_Unscaled
2051	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Alb_Coll_Unscaled
2052	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Tran_Coll_Unscaled
2053	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Tran_UnColl_Unscaled
2054	REAL(r_std),DIMENSION(nlevels_tot) :: Collim_Tran_Tot
2055	REAL(r_std),DIMENSION(nlevels_tot) :: Isotrop_Tran_Tot
2056
2057	REAL(r_std),DIMENSION(0:nlevels_tot+1) :: &
2058	collim_down_cs,isotrop_down_cs,isotrop_up_cs
2059
2060	INTEGER :: istep,ilevel
2061
2062	REAL(r_std) :: leaf_reflectance
2063	REAL(r_std) :: leaf_transmittance
2064
2065	LOGICAL :: lexit
2066	LOGICAL :: ok
2067	REAL(r_std), DIMENSION(4) :: gammaCoeffs
2068	REAL(r_std), DIMENSION(4) :: gammaCoeffs_star
2069	REAL(r_std) :: sun_zenith_angle_radians
2070	REAL(r_std), PARAMETER :: isotropic_cosine_constant=0.5/0.705
2071
2072	!_ ================================================================================================================================
2073
2074
2075	! convet the sun angle
2076	sun_zenith_angle_radians = acos(cosine_sun_angle)
2077
2078
2079
2080	! initialize some values that are used in the optimization loop
2081	istep=0
2082	lexit=.FALSE.
2083	lconverged=.FALSE.
2084
2085	! calculate the albedo at our starting point
2086
2087	leaf_transmittance = leaf_single_scattering_albedo_start/ &
2088	( leaf_psd_start+un)
2089	leaf_reflectance = leaf_psd_start * &
2090	leaf_transmittance
2091
2092	! some debugging stuff
2093	! DO ilevel=1,nlevels_tot
2094	!
2095	! CALL print_debugging_albedo_info(ilevel,cosine_sun_angle,leaf_reflectance,&
2096	! leaf_transmittance,&
2097	! laieff_collim_pft(ilevel),laieff_isotrop_pft(ilevel), &
2098	! reflectance_collim(ilevel-1),reflectance_isotrop(ilevel-1),&
2099	! Collim_Alb_Tot_temp(ilevel),Collim_Tran_Tot_temp(ilevel),&
2100	! Collim_Abs_Tot_temp(ilevel),&
2101	! Isotrop_Alb_Tot_temp(ilevel),Isotrop_Tran_Tot_temp(ilevel),&
2102	! Isotrop_Abs_Tot_temp(ilevel))
2103	! ENDDO
2104
2105	! calculate the gamma coefficients used in the case of the black background
2106	call gammas(leaf_reflectance, leaf_transmittance, cosine_sun_angle, gammaCoeffs)
2107	call gammas(leaf_reflectance,leaf_transmittance,isotropic_cosine_constant,&
2108	gammaCoeffs_star)
2109
2110
2111
2112	!************************ step one **** !
2113	! compute all the unscaled quantities
2114	DO ilevel=nlevels_tot,1,-1
2115
2116	Collim_Tran_UnColl_Unscaled(ilevel)=exp( - 0.5_r_std*&
2117	laieff_collim_pft(ilevel)/cosine_sun_angle)
2118
2119	Isotrop_Tran_UnColl_Unscaled(ilevel)= &
2120	TBarreUncoll_exact(0.5_r_std*laieff_isotrop_pft(ilevel))
2121
2122
2123	! /* 1) collimated source */
2124	ok=dhrT1(leaf_reflectance,leaf_transmittance,&
2125	gammaCoeffs(1),gammaCoeffs(2),gammaCoeffs(3),gammaCoeffs(4),&
2126	sun_zenith_angle_radians, 0.5_r_std*laieff_collim_pft(ilevel),&
2127	Collim_Alb_Coll_Unscaled(ilevel),Collim_Tran_Coll_Unscaled(ilevel),&
2128	Collim_Abs_Coll_Unscaled(ilevel))
2129
2130	! /* 2) isotropic source */
2131	ok=dhrT1(leaf_reflectance,leaf_transmittance,&
2132	gammaCoeffs_star(1),gammaCoeffs_star(2),gammaCoeffs_star(3),gammaCoeffs_star(4),&
2133	acos(isotropic_cosine_constant),0.5_r_std*laieff_isotrop_pft(ilevel),&
2134	Isotrop_Alb_Coll_Unscaled(ilevel),Isotrop_Tran_Coll_Unscaled(ilevel),&
2135	Isotrop_Abs_Coll_Unscaled(ilevel))
2136	ENDDO ! ilevel=nlevels_tot,1,-1
2137
2138
2139	! Following the fate of the light at every step
2140
2141	! The downwelling array indicates the quantity of light flowing into this level
2142	! from above therefore, to start our system for collimated light, we give a
2143	! unit of light coming into the top level from a collimated source.
2144	! The upwelling array is light entering this level from below.
2145	! cs is the current step, while ns is the next step for the iteration.
2146	! Level 0 is the background. Light can enter this level from above but not below
2147	! Level nlevels_tot+1 is the atomosphere. Light can enter this level from below
2148	! but not above.
2149	collim_down_cs(:)=zero
2150	collim_down_cs(nlevels_tot)=un
2151	isotrop_down_cs(:)=zero
2152	isotrop_up_cs(:)=zero
2153
2154	IF(lprint)THEN
2155	WRITE(numout,*) "Solving fluxes for a collimated light source (i.e., the sun) in the canopy."
2156	ENDIF
2157	CALL propagate_fluxes(collim_down_cs, isotrop_down_cs, isotrop_up_cs, &
2158	Collim_Tran_UnColl_Unscaled, Collim_Tran_Coll_Unscaled, &
2159	Collim_Alb_Coll_Unscaled, Isotrop_Tran_UnColl_Unscaled, &
2160	Isotrop_Tran_Coll_Unscaled, Isotrop_Alb_Coll_Unscaled, &
2161	br_base_temp_collim, br_base_temp_isotrop, .FALSE., &
2162	Collim_Tran_Uncoll, Collim_Tran_Coll, Collim_Alb_Coll, lconverged, lprint)
2163
2164	! Now for the isotropic light
2165	collim_down_cs(:)=zero
2166	isotrop_down_cs(:)=zero
2167	isotrop_down_cs(nlevels_tot)=un
2168	isotrop_up_cs(:)=zero
2169
2170	IF(lprint)THEN
2171	WRITE(numout,*) "Solving fluxes for diffuse light sources (e.g., clouds, aerosols) in the canopy."
2172	ENDIF
2173	! The colliminated light is never used here in propagate_fluxes for the
2174	! isotropic light source, but it's passed to keep things
2175	! clean between the two cases.
2176	CALL propagate_fluxes(collim_down_cs, isotrop_down_cs, isotrop_up_cs, &
2177	Collim_Tran_UnColl_Unscaled, Collim_Tran_Coll_Unscaled, &
2178	Collim_Alb_Coll_Unscaled, Isotrop_Tran_UnColl_Unscaled, &
2179	Isotrop_Tran_Coll_Unscaled, Isotrop_Alb_Coll_Unscaled, &
2180	br_base_temp_collim, br_base_temp_isotrop, .TRUE., &
2181	Isotrop_Tran_Uncoll, Isotrop_Tran_Coll, Isotrop_Alb_Coll, lconverged, lprint)
2182
2183	! Calculate the absorption profile
2184	Collim_Tran_Tot(:)=Collim_Tran_Coll(:)+Collim_Tran_Uncoll(:)
2185	Isotrop_Tran_Tot(:)=Isotrop_Tran_Coll(:)+Isotrop_Tran_Uncoll(:)
2186
2187	! bottom level
2188	Collim_Abs_Tot(1)=Collim_Tran_Tot(1+1)+Collim_Tran_Tot(1)*br_base_temp_collim&
2189	-Collim_Tran_Tot(1)-Collim_Alb_Coll(1)
2190	Isotrop_Abs_Tot(1)=Isotrop_Tran_Tot(1+1)+Isotrop_Tran_Tot(1)*br_base_temp_isotrop&
2191	-Isotrop_Tran_Tot(1)-Isotrop_Alb_Coll(1)
2192
2193	! all middle levels
2194	DO ilevel=2,nlevels_tot-1
2195	Collim_Abs_Tot(ilevel)=Collim_Tran_Tot(ilevel+1)+Collim_Alb_Coll(ilevel-1)&
2196	-Collim_Tran_Tot(ilevel)-Collim_Alb_Coll(ilevel)
2197	Isotrop_Abs_Tot(ilevel)=Isotrop_Tran_Tot(ilevel+1)+Isotrop_Alb_Coll(ilevel-1)&
2198	-Isotrop_Tran_Tot(ilevel)-Isotrop_Alb_Coll(ilevel)
2199	ENDDO
2200
2201	! top level
2202	Collim_Abs_Tot(nlevels_tot)=un+Collim_Alb_Coll(nlevels_tot-1)&
2203	-Collim_Tran_Tot(nlevels_tot)-Collim_Alb_Coll(nlevels_tot)
2204	Isotrop_Abs_Tot(nlevels_tot)=un+Isotrop_Alb_Coll(nlevels_tot-1)&
2205	-Isotrop_Tran_Tot(nlevels_tot)-Isotrop_Alb_Coll(nlevels_tot)
2206
2207	END SUBROUTINE multilevel_albedo
2208
2209	!! ==============================================================================================================================
2210	!! SUBROUTINE : print_debugging_albedo_info
2211	!!
2212	!>\BRIEF Prints out some albedo information in a nice format.
2213	!! Should only be used for debugging, never for production runs.
2214	!!
2215	!! DESCRIPTION :
2216	!!
2217	!! RECENT CHANGE(S): None\n
2218	!!
2219	!! MAIN OUTPUT VARIABLE(S): None.
2220	!!
2221	!! REFERENCE(S) :
2222	!!
2223	!! FLOWCHART : None
2224	!_ ================================================================================================================================
2225	SUBROUTINE print_debugging_albedo_info(ilevel, cosine_sun_angle, &
2226	leaf_reflectance, leaf_transmittance, laieff_collim_temp, laieff_isotrop_temp, &
2227	br_base_temp_collim, br_base_temp_isotrop, Collim_Alb_Tot, &
2228	Collim_Tran_Tot, Collim_Abs_Tot, Isotrop_Alb_Tot, Isotrop_Tran_Tot, Isotrop_Abs_Tot)
2229
2230	!! 0. Variables and parameter declaration
2231	!! 0.1 Input variables
2232	INTEGER,INTENT(IN) :: ilevel
2233	REAL(r_std), INTENT(IN) :: cosine_sun_angle !! cosine of the solar zenith angle, between 0 and 1
2234	!! @tex $()$ @endtex
2235	REAL(r_std), INTENT(IN) :: laieff_collim_temp !! cosine of the solar zenith angle, between 0 and 1
2236	REAL(r_std), INTENT(IN) :: laieff_isotrop_temp !! cosine of the solar zenith angle, between 0 and 1
2237	REAL(r_std), INTENT(IN) :: leaf_reflectance !! cosine of the solar zenith angle, between 0 and 1
2238	REAL(r_std), INTENT(IN) :: leaf_transmittance !! cosine of the solar zenith angle, between 0 and 1
2239	REAL(r_std), INTENT(IN) :: br_base_temp_collim !! cosine of the solar zenith angle, between 0 and 1
2240	REAL(r_std), INTENT(IN) :: br_base_temp_isotrop !! cosine of the solar zenith angle, between 0 and 1
2241	REAL(r_std),INTENT(IN) :: Collim_Alb_Tot
2242	REAL(r_std),INTENT(IN) :: Collim_Tran_Tot
2243	REAL(r_std),INTENT(IN) :: Collim_Abs_Tot
2244	REAL(r_std),INTENT(IN) :: Isotrop_Alb_Tot
2245	REAL(r_std),INTENT(IN) :: Isotrop_Tran_Tot
2246	REAL(r_std),INTENT(IN) :: Isotrop_Abs_Tot
2247
2248	!! 0.2 Output variables
2249
2250	!! 0.3 Modified variables
2251
2252	!! 0.4 Local variables
2253
2254	!_ ================================================================================================================================
2255
2256
2257	WRITE (numout,'(8(11A))') ' Level ',' Angle ',' laieff_c ',' laieff_i ',&
2258	' rleaf ',' tleaf ',' rbgd_c ',' rbgd_i '
2259
2260	WRITE(numout,FMT='(I6,5X,7(F11.6))') &
2261	ilevel,180/pi*ACOS(cosine_sun_angle),&
2262	laieff_collim_temp,laieff_isotrop_temp,&
2263	leaf_reflectance,leaf_transmittance,br_base_temp_collim,&
2264	br_base_temp_isotrop
2265
2266
2267	WRITE (numout,'(10X,6(11A))') ' Rtot(sun) ',' Ttot(sun) ',&
2268	' Atot(sun) ',' Rtot(iso) ',' Ttot(iso) ',' Atot(iso) '
2269	WRITE(numout,FMT='(10X,6(F11.6))') &
2270	Collim_Alb_Tot,Collim_Tran_Tot,Collim_Abs_Tot,&
2271	Isotrop_Alb_Tot,Isotrop_Tran_Tot,Isotrop_Abs_Tot
2272
2273
2274	END SUBROUTINE print_debugging_albedo_info
2275
2276
2277	!! ==============================================================================================================================
2278	!! SUBROUTINE : propagate_fluxes
2279	!!
2280	!>\BRIEF Propogates the radiation fluxes through each level of the
2281	!! canopy.
2282	!!
2283	!! DESCRIPTION : This is an algorithm to follow radition fluxes through the
2284	!! canopy. At each level, the path probabilities are determined by the
2285	!! raditional transfer scheme of Pinty et al (2006). Notice that
2286	!! the fluxes given by this routine are all cumulative fluxes, not per
2287	!! level.
2288	!!
2289	!! RECENT CHANGE(S): None\n
2290	!!
2291	!! MAIN OUTPUT VARIABLE(S): ::Tran_Uncoll_Tot, ::Tran_Coll_Tot, ::Alb_Coll_Tot
2292	!!
2293	!! REFERENCE(S) :
2294	!!
2295	!! FLOWCHART : None
2296	!_ ================================================================================================================================
2297	SUBROUTINE propagate_fluxes(collim_down_cs, isotrop_down_cs, isotrop_up_cs, &
2298	Collim_Tran_UnColl_Unscaled, Collim_Tran_Coll_Unscaled, &
2299	Collim_Alb_Coll_Unscaled, Isotrop_Tran_UnColl_Unscaled, &
2300	Isotrop_Tran_Coll_Unscaled, Isotrop_Alb_Coll_Unscaled, br_base_temp_collim, &
2301	br_base_temp_isotrop, lisotrop, Tran_Uncoll_Tot, Tran_Coll_Tot, &
2302	Alb_Coll_Tot, lconverged, lprint)
2303
2304	!! 0. Variables and parameter declaration
2305	!! 0.1 Input variables
2306	REAL(r_std), DIMENSION(nlevels_tot),INTENT(IN) :: Collim_Tran_UnColl_Unscaled
2307	REAL(r_std), DIMENSION(nlevels_tot),INTENT(IN) :: Collim_Tran_Coll_Unscaled
2308	REAL(r_std), DIMENSION(nlevels_tot),INTENT(IN) :: Collim_Alb_Coll_Unscaled
2309	REAL(r_std), DIMENSION(nlevels_tot),INTENT(IN) :: Isotrop_Tran_UnColl_Unscaled
2310	REAL(r_std), DIMENSION(nlevels_tot),INTENT(IN) :: Isotrop_Tran_Coll_Unscaled
2311	REAL(r_std), DIMENSION(nlevels_tot),INTENT(IN) :: Isotrop_Alb_Coll_Unscaled
2312	REAL(r_std), INTENT(IN) :: br_base_temp_collim !! cosine of the solar zenith angle, between 0 and 1
2313	REAL(r_std), INTENT(IN) :: br_base_temp_isotrop !! cosine of the solar zenith angle, between 0 and 1
2314	LOGICAL, INTENT(IN) :: lisotrop !! are we dealing with an isotropic source? only needed for correct partitioning
2315	!! of the collided and uncollided transmitted light...the total is not affected
2316	LOGICAL, INTENT(IN) :: lprint !! a flag to print
2317
2318	!! 0.2 Output variables
2319	REAL(r_std), DIMENSION(nlevels_tot),INTENT(OUT) :: Tran_Uncoll_Tot
2320	REAL(r_std), DIMENSION(nlevels_tot),INTENT(OUT) :: Tran_Coll_Tot
2321	REAL(r_std), DIMENSION(nlevels_tot),INTENT(OUT) :: Alb_Coll_Tot
2322	LOGICAL,INTENT(OUT) :: lconverged
2323
2324	!! 0.3 Modified variables
2325	REAL(r_std), DIMENSION(0:nlevels_tot+1),INTENT(INOUT) :: collim_down_cs
2326	REAL(r_std), DIMENSION(0:nlevels_tot+1),INTENT(INOUT) :: isotrop_down_cs
2327	REAL(r_std), DIMENSION(0:nlevels_tot+1),INTENT(INOUT) :: isotrop_up_cs
2328
2329	!! 0.4 Local variables
2330	INTEGER :: istep,ilevel
2331	REAL(r_std),DIMENSION(0:nlevels_tot+1) :: collim_down_ns,isotrop_down_ns,&
2332	isotrop_up_ns,Tran_Uncoll_Tot_temp
2333
2334	!_ ================================================================================================================================
2335
2336	Tran_Uncoll_Tot(:)=zero
2337	Tran_Coll_Tot(:)=zero
2338	Alb_Coll_Tot(:)=zero
2339
2340	Tran_Uncoll_Tot_temp(:)=un
2341
2342	istep=0
2343
2344	DO
2345
2346	istep=istep+1
2347
2348	lconverged=.TRUE.
2349
2350	! Zero out the counters for the next step
2351	collim_down_ns(:)=zero
2352	isotrop_down_ns(:)=zero
2353	isotrop_up_ns(:)=zero
2354
2355
2356	! Now we need to loop over all levels and see what light is entering the
2357	! level, and how it will propagate in the next step
2358	DO ilevel=1,nlevels_tot
2359
2360	! For collimated downwelling light into the level, it can be scattered
2361	! up, down, or pass through uncollided
2362	IF(collim_down_cs(ilevel) .GT. zero)THEN
2363	collim_down_ns(ilevel-1)=collim_down_ns(ilevel-1)+&
2364	collim_down_cs(ilevel)*Collim_Tran_UnColl_Unscaled(ilevel)
2365
2366
2367	! This statement checks to see if this light has been previously
2368	! scattered or not. This term is only present in level nlevels_tot
2369	! for the first step, level nlevels_tot-1 for the second step,
2370	! level nlevels_tot-2 for the third step, etc., and it only happens
2371	! in the case of a collimated light source.
2372	IF(ilevel == nlevels_tot-istep+1 .AND. .NOT. lisotrop) THEN
2373	Tran_Uncoll_Tot_temp(ilevel)=Tran_Uncoll_Tot_temp(ilevel+1)*&
2374	Collim_Tran_UnColl_Unscaled(ilevel)
2375	ENDIF
2376
2377	isotrop_down_ns(ilevel-1)=isotrop_down_ns(ilevel-1)+&
2378	collim_down_cs(ilevel)*Collim_Tran_Coll_Unscaled(ilevel)
2379	isotrop_up_ns(ilevel+1)=isotrop_up_ns(ilevel+1)+&
2380	collim_down_cs(ilevel)*Collim_Alb_Coll_Unscaled(ilevel)
2381
2382	ENDIF ! collim_down_cs(ilevel) .GT. zero
2383
2384	! For isotropic downwelling light, it can also be scattered up, down,
2385	! or pass through uncollided
2386	IF(isotrop_down_cs(ilevel) .GT. zero)THEN
2387	isotrop_down_ns(ilevel-1)=isotrop_down_ns(ilevel-1)+&
2388	isotrop_down_cs(ilevel)*Isotrop_Tran_UnColl_Unscaled(ilevel)
2389
2390	! This is the same check as above, but this time the light has
2391	! an isotropic source and not a collimated source
2392	IF(ilevel == nlevels_tot-istep+1 .AND. lisotrop) THEN
2393	Tran_Uncoll_Tot_temp(ilevel)=Tran_Uncoll_Tot_temp(ilevel+1)&
2394	*Isotrop_Tran_UnColl_Unscaled(ilevel)
2395	ENDIF
2396
2397	isotrop_down_ns(ilevel-1)=isotrop_down_ns(ilevel-1)+&
2398	isotrop_down_cs(ilevel)*Isotrop_Tran_Coll_Unscaled(ilevel)
2399	isotrop_up_ns(ilevel+1)=isotrop_up_ns(ilevel+1)+&
2400	isotrop_down_cs(ilevel)*Isotrop_Alb_Coll_Unscaled(ilevel)
2401	ENDIF ! isotrop_down_cs(ilevel) .GT. zero
2402
2403	! Isotropic upwelling light can pass through upwards collided or
2404	! uncollided with vegetation in this level, or it can be reflected downwards
2405	IF(isotrop_up_cs(ilevel) .GT. zero)THEN
2406
2407	isotrop_up_ns(ilevel+1)=isotrop_up_ns(ilevel+1)+isotrop_up_cs(ilevel)&
2408	*Isotrop_Tran_UnColl_Unscaled(ilevel)
2409	isotrop_up_ns(ilevel+1)=isotrop_up_ns(ilevel+1)+isotrop_up_cs(ilevel)&
2410	*Isotrop_Tran_Coll_Unscaled(ilevel)
2411	isotrop_down_ns(ilevel-1)=isotrop_down_ns(ilevel-1)+&
2412	isotrop_up_cs(ilevel)*Isotrop_Alb_Coll_Unscaled(ilevel)
2413	ENDIF
2414
2415	! there can be no collimated upwards light, since upwards light must
2416	! always have reflected off something
2417
2418
2419	ENDDO ! ilevel=1,nlevels_tot
2420
2421	! The background is a bit special. there is no transmission, but there is
2422	! reflection, which leads to a source term to the bottom level from below.
2423	isotrop_up_ns(1)=isotrop_up_ns(1)+collim_down_cs(0)*br_base_temp_collim
2424	isotrop_up_ns(1)=isotrop_up_ns(1)+isotrop_down_cs(0)*br_base_temp_isotrop
2425
2426
2427	! Now we add the light generated here to our cumulative counters to track
2428	! the total amount that leaves the canopy, either through being absorbed
2429	! by the background or being reflected back into the atmosphere.
2430	! We keep track of the uncoll above, so here we track the total light that
2431	! is transmitted to the soil, and then at the end we taken the difference
2432	! to get the collided radation.
2433	Tran_Coll_Tot(1:nlevels_tot)=Tran_Coll_Tot(1:nlevels_tot)+&
2434	isotrop_down_ns(0:nlevels_tot-1)+collim_down_ns(0:nlevels_tot-1)
2435	Alb_Coll_Tot(1:nlevels_tot)=Alb_Coll_Tot(1:nlevels_tot)+isotrop_up_ns(2:nlevels_tot+1)
2436
2437	! now we update the light we are currently tracking
2438	collim_down_cs(:)=collim_down_ns(:)
2439	isotrop_down_cs(:)=isotrop_down_ns(:)
2440	isotrop_up_cs(:)=isotrop_up_ns(:)
2441
2442
2443	! check for convegence...if all the values of light are currently below a
2444	! threshold, we're probably good assuming the threshold is low enough.
2445	IF(MAXVAL(collim_down_cs,1) .GT. alb_threshold) lconverged=.FALSE.
2446	IF(MAXVAL(isotrop_down_cs,1) .GT. alb_threshold) lconverged=.FALSE.
2447	IF(MAXVAL(isotrop_up_cs,1) .GT. alb_threshold) lconverged=.FALSE.
2448
2449	IF(lconverged)THEN
2450	IF(lprint) WRITE(numout,*) 'Converged after this many steps: ',istep
2451	EXIT
2452	ENDIF
2453
2454	! This number here could also be externalized.
2455	IF(istep .GE. 1000)THEN
2456	WRITE(numout,) '********************************************************'
2457	WRITE(numout,'(A,I6,F14.10)') 'Albedo not converging!',istep,alb_threshold
2458	WRITE(numout,'(A)') ' collim_down_cs ' // &
2459	'isotrop_down_cs isotrop_up_cs'
2460	DO ilevel=0,nlevels_tot+1
2461	WRITE(numout,'(I4,3F14.10)') ilevel, &
2462	collim_down_cs(ilevel),isotrop_down_cs(ilevel),isotrop_up_cs(ilevel)
2463
2464	ENDDO
2465	WRITE(numout,*) 'You should increase either the number' // &
2466	'of steps or the alb_threshold.'
2467	WRITE(numout,) '********************************************************'
2468	EXIT
2469	ENDIF ! istep .GE. 1000
2470	ENDDO ! convergence loop
2471
2472
2473	! now separate the collided from the uncollided light. Notice that
2474	! this is not really needed for any purposes other than debugging, as the
2475	! important quantity is the total amount of light striking the ground.
2476	Tran_UnColl_Tot(1:nlevels_tot)=Tran_UnColl_Tot_temp(1:nlevels_tot)
2477	Tran_Coll_Tot(1:nlevels_tot)=Tran_Coll_Tot(1:nlevels_tot)-Tran_UnColl_Tot(1:nlevels_tot)
2478
2479	! Some debugging information
2480	IF(lprint)THEN
2481	WRITE(numout,'(7X,3(A15,3X))') 'Tran_Uncoll_Tot',' Tran_Coll_Tot',' Alb_Coll_Tot'
2482	DO ilevel=nlevels_tot,1,-1
2483	WRITE(numout,'(I4,3X,3(F15.6,3X))') ilevel, &
2484	Tran_Uncoll_Tot(ilevel),Tran_Coll_Tot(ilevel),Alb_Coll_Tot(ilevel)
2485	ENDDO
2486	ENDIF
2487
2488	END SUBROUTINE propagate_fluxes
2489
2490
2491	!! ==============================================================================================================================
2492	!! SUBROUTINE : read_background_albedo
2493	!!
2494	!>\BRIEF This subroutine reads the background albedo
2495	!!
2496	!! DESCRIPTION This subroutine reads the background albedo map in 0.5 x 0.5 deg resolution
2497	!! derived from JRCTIP product. These values are then interpolated to the resolution of the
2498	!! simulation. For deserts and fallow croplands, the background albedo will
2499	!! be similar to the bare soil albedo. For all other vegetated PFTs, the background albedo
2500	!! will be determined by the understory and/or the litter layer.
2501	!!
2502	!! RECENT CHANGE(S): None
2503	!!
2504	!! MAIN OUTPUT VARIABLE(S): bckgrnd_alb for visible and near-infrared range
2505	!!
2506	!! REFERENCES : None
2507	!!
2508	!! FLOWCHART : None
2509	!! \n
2510	!_ ==============================================================================================================================
2511
2512	SUBROUTINE read_background_albedo(nbpt, lalo, neighbours, resolution, contfrac)
2513
2514	USE interpweight
2515
2516	IMPLICIT NONE
2517
2518	!! 0. Variable and parameter declaration
2519
2520	!! 0.1 Input variables
2521
2522	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
2523	!! interpolated (unitless)
2524	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
2525	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
2526	!! (1=N, 2=E, 3=S, 4=W)
2527	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
2528	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
2529
2530	!! 0.4 Local variables
2531
2532	CHARACTER(LEN=80) :: filename !! Filename of background albedo
2533	REAL(r_std), DIMENSION(nbpt) :: aalb_bg !! Availability of the interpolation
2534	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
2535	REAL(r_std) :: vmin, vmax !! min/max values to use for the
2536	!! renormalization
2537	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2538	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2539	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2540	!! 'XYKindTime': Input values are kinds
2541	!! of something with a temporal
2542	!! evolution on the dx*dy matrix'
2543	LOGICAL :: nonegative !! whether negative values should be removed
2544	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2545	!! 'nomask': no-mask is applied
2546	!! 'mbelow': take values below maskvals(1)
2547	!! 'mabove': take values above maskvals(1)
2548	!! 'msumrange': take values within 2 ranges;
2549	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2550	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2551	!! (normalized by maskedvals(3))
2552	!! 'var': mask values are taken from a
2553	!! variable inside the file (>0)
2554	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2555	!! `maskingtype')
2556	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2557	REAL(r_std) :: albbg_norefinf !! No value
2558	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: albbg_default !! Default value
2559
2560	!_ ================================================================================================================================
2561
2562	!! 1. Open file and allocate memory
2563
2564	! Open file with background albedo
2565
2566	!Config Key = ALB_BG_FILE
2567	!Config Desc = Name of file from which the background albedo is read
2568	!Config Def = alb_bg.nc
2569	!Config If = ALB_BG_MODIS
2570	!Config Help = The name of the file to be opened to read background albedo
2571	!Config Units = [FILE]
2572	!
2573	filename = 'alb_bg.nc'
2574	CALL getin_p('ALB_BG_FILE',filename)
2575
2576
2577	IF (xios_interpolation) THEN
2578
2579	! Read and interpolation background albedo using XIOS
2580	CALL xios_orchidee_recv_field('bg_alb_vis_interp',bckgrnd_alb(:,ivis))
2581	CALL xios_orchidee_recv_field('bg_alb_nir_interp',bckgrnd_alb(:,inir))
2582
2583	aalb_bg(:)=1
2584
2585	ELSE
2586	ALLOCATE(albbg_default(2), STAT=ALLOC_ERR)
2587	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'read_background_albedo','Pb in allocation for albbg_default','','')
2588
2589	! For this case there are not types/categories. We have 'only' a continuous field
2590	! Assigning values to vmin, vmax
2591
2592	vmin = 0.
2593	vmax = 9999.
2594
2595	!! Variables for interpweight
2596	! Type of calculation of cell fractions (not used here)
2597	fractype = 'default'
2598	! Name of the longitude and latitude in the input file
2599	lonname = 'longitude'
2600	latname = 'latitude'
2601	! Default value when no value is get from input file
2602	albbg_default(ivis) = 0.129
2603	albbg_default(inir) = 0.247
2604	! Reference value when no value is get from input file (not used here)
2605	albbg_norefinf = undef_sechiba
2606	! Should negative values be set to zero from input file?
2607	nonegative = .FALSE.
2608	! Type of mask to apply to the input data (see header for more details)
2609	maskingtype = 'var'
2610	! Values to use for the masking (here not used)
2611	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
2612	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var')
2613	namemaskvar = 'mask'
2614
2615	! There is a variable for each chanel 'infrared' and 'visible'
2616	! Interpolate variable bg_alb_vis
2617	variablename = 'bg_alb_vis'
2618	IF (printlev_loc >= 2) WRITE(numout,*) "read_background_albedo: Start interpolate " &
2619	// TRIM(filename) // " for variable " // TRIM(variablename)
2620	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
2621	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2622	maskvals, namemaskvar, -1, fractype, albbg_default(ivis), albbg_norefinf, &
2623	bckgrnd_alb(:,ivis), aalb_bg)
2624	IF (printlev_loc >= 5) WRITE(numout,*)" read_background_albedo after InterpWeight2Dcont for '" // &
2625	TRIM(variablename) // "'"
2626
2627	! Interpolate variable bg_alb_nir in the same file
2628	variablename = 'bg_alb_nir'
2629	IF (printlev_loc >= 2) WRITE(numout,*) "read_background_albedo: Start interpolate " &
2630	// TRIM(filename) // " for variable " // TRIM(variablename)
2631	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
2632	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2633	maskvals, namemaskvar, -1, fractype, albbg_default(inir), albbg_norefinf, &
2634	bckgrnd_alb(:,inir), aalb_bg)
2635	IF (printlev_loc >= 5) WRITE(numout,*)" read_background_albedo after InterpWeight2Dcont for '" // &
2636	TRIM(variablename) // "'"
2637
2638	IF (ALLOCATED(albbg_default)) DEALLOCATE(albbg_default)
2639
2640	ENDIF
2641
2642	! Write diagnostics
2643	CALL xios_orchidee_send_field("interp_diag_alb_vis",bckgrnd_alb(:,ivis))
2644	CALL xios_orchidee_send_field("interp_diag_alb_nir",bckgrnd_alb(:,inir))
2645	CALL xios_orchidee_send_field("interp_avail_aalb_bg",aalb_bg)
2646
2647	IF (printlev_loc >= 3) WRITE(numout,*)' read_background_albedo ended'
2648
2649
2650	END SUBROUTINE read_background_albedo
2651
2652
2653	!! ==============================================================================================================================
2654	!! SUBROUTINE : albedo_surface_soilalb
2655	!!
2656	!>\BRIEF This subroutine calculates the albedo of soil (without snow).
2657	!!
2658	!! DESCRIPTION This subroutine reads the soil colour maps in 1 x 1 deg resolution
2659	!! from the Henderson-Sellers & Wilson database. These values are interpolated to
2660	!! the model's resolution and transformed into
2661	!! dry and wet albedos.\n
2662	!!
2663	!! If the soil albedo is calculated without the dependence of soil moisture, the
2664	!! soil colour values are transformed into mean soil albedo values.\n
2665	!!
2666	!! The calculations follow the assumption that the grid of the data is regular and
2667	!! it covers the globe. The calculation for the model grid are based on the borders
2668	!! of the grid of the resolution.
2669	!!
2670	!! RECENT CHANGE(S): None
2671	!!
2672	!! CALCULATED MODULE VARIABLE(S): soilalb_dry for visible and near-infrared range,
2673	!! soilalb_wet for visible and near-infrared range,
2674	!! soilalb_moy for visible and near-infrared range
2675	!!
2676	!! REFERENCE(S) :
2677	!! -Wilson, M.F., and A. Henderson-Sellers, 1985: A global archive of land cover and
2678	!! soils data for use in general circulation climate models. J. Clim., 5, 119-143.
2679	!!
2680	!! FLOWCHART : None
2681	!! \n
2682	!_ ================================================================================================================================
2683
2684	SUBROUTINE albedo_surface_soilalb(nbpt, lalo, neighbours, resolution, contfrac)
2685
2686	USE interpweight
2687
2688	IMPLICIT NONE
2689
2690
2691	!! 0. Variable and parameter declaration
2692
2693	!! 0.1 Input variables
2694
2695	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
2696	!! interpolated (unitless)
2697	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
2698	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
2699	!! (1=N, 2=E, 3=S, 4=W)
2700	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
2701	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
2702
2703	!! 0.4 Local variables
2704
2705	CHARACTER(LEN=80) :: filename !! Filename of soil colour map
2706	INTEGER(i_std) :: i, ib, ip, nbexp !! Indices
2707	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
2708	REAL(r_std), DIMENSION(nbpt) :: asoilcol !! Availability of the soilcol interpolation
2709	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
2710	!! (variabletypevals(1) = -un, not used)
2711	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: soilcolrefrac !! soilcol fractions re-dimensioned
2712	REAL(r_std) :: vmin, vmax !! min/max values to use for the
2713	!! renormalization
2714	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2715	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2716	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2717	!! 'XYKindTime': Input values are kinds
2718	!! of something with a temporal
2719	!! evolution on the dx*dy matrix'
2720	LOGICAL :: nonegative !! whether negative values should be removed
2721	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2722	!! 'nomask': no-mask is applied
2723	!! 'mbelow': take values below maskvals(1)
2724	!! 'mabove': take values above maskvals(1)
2725	!! 'msumrange': take values within 2 ranges;
2726	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2727	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2728	!! (normalized by maskvals(3))
2729	!! 'var': mask values are taken from a
2730	!! variable inside the file (>0)
2731	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2732	!! `maskingtype')
2733	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2734	CHARACTER(LEN=250) :: msg
2735	INTEGER :: fopt
2736	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: vecpos
2737	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: solt
2738
2739	!_ ================================================================================================================================
2740	!! 1. Open file and allocate memory
2741
2742	IF (grid_type==unstructured) THEN
2743	CALL ipslerr_p(3,'albedo_surface_soilalb','Reading of SOILALB_FILE must be implemented with XIOS to be used for unstructured grid.', &
2744	'Use option alb_bg_modis for unstructured grid for now.','')
2745	END IF
2746
2747	! Open file with soil colours
2748
2749	!Config Key = SOILALB_FILE
2750	!Config Desc = Name of file from which the bare soil albedo
2751	!Config Def = soils_param.nc
2752	!Config If = NOT(IMPOSE_AZE)
2753	!Config Help = The name of the file to be opened to read the soil types from
2754	!Config which we derive then the bare soil albedos. This file is 1x1
2755	!Config deg and based on the soil colors defined by Wilson and Henderson-Seller.
2756	!Config Units = [FILE]
2757	!
2758	filename = 'soils_param.nc'
2759	CALL getin_p('SOILALB_FILE',filename)
2760
2761
2762	ALLOCATE(soilcolrefrac(nbpt, classnb), STAT=ALLOC_ERR)
2763	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'albedo_surface_soilalb','Problem in allocation of variable soilcolrefrac','','')
2764	ALLOCATE(vecpos(classnb), STAT=ALLOC_ERR)
2765	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'albedo_surface_soilalb','Problem in allocation of variable vecpos','','')
2766	ALLOCATE(solt(classnb), STAT=ALLOC_ERR)
2767	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'albedo_surface_soilalb','Problem in allocation of variable solt','','')
2768
2769	! Assigning values to vmin, vmax
2770	vmin = 1.0
2771	vmax = classnb
2772
2773	ALLOCATE(variabletypevals(classnb),STAT=ALLOC_ERR)
2774	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variabletypevals','','')
2775	variabletypevals = -un
2776
2777	!! Variables for interpweight
2778	! Type of calculation of cell fractions
2779	fractype = 'default'
2780	! Name of the longitude and latitude in the input file
2781	lonname = 'nav_lon'
2782	latname = 'nav_lat'
2783	! Should negative values be set to zero from input file?
2784	nonegative = .FALSE.
2785	! Type of mask to apply to the input data (see header for more details)
2786	maskingtype = 'mabove'
2787	! Values to use for the masking
2788	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
2789	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2790	namemaskvar = ''
2791
2792	! Interpolate variable soilcolor
2793	variablename = 'soilcolor'
2794	IF (printlev_loc >= 2) WRITE(numout,*) "albedo_surface_soilalb: Start interpolate " &
2795	// TRIM(filename) // " for variable " // TRIM(variablename)
2796	CALL interpweight_2D(nbpt, classnb, variabletypevals, lalo, resolution, neighbours, &
2797	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2798	maskvals, namemaskvar, 0, 0, -1, fractype, &
2799	-1., -1., soilcolrefrac, asoilcol)
2800	IF (printlev_loc >= 5) WRITE(numout,*)' albedo_surface_soilalb after interpweight_2D'
2801
2802	! Check how many points with soil information are found
2803	nbexp = 0
2804
2805	soilalb_dry(:,:) = zero
2806	soilalb_wet(:,:) = zero
2807	soilalb_moy(:,:) = zero
2808	IF (printlev_loc >= 5) THEN
2809	WRITE(numout,*)' albedo_surface_soilalb before starting loop nbpt:', nbpt
2810	WRITE(numout,*)' albedo_surface_soilalb initial values classnb: ',classnb
2811	WRITE(numout,*)' albedo_surface_soilalb vis_dry. SUM:',SUM(vis_dry),' vis_dry= ',vis_dry
2812	WRITE(numout,*)' albedo_surface_soilalb nir_dry. SUM:',SUM(nir_dry),' nir_dry= ',nir_dry
2813	WRITE(numout,*)' albedo_surface_soilalb vis_wet. SUM:',SUM(vis_wet),' vis_wet= ',vis_wet
2814	WRITE(numout,*)' albedo_surface_soilalb nir_wet. SUM:',SUM(nir_wet),' nir_wet= ',nir_wet
2815	END IF
2816
2817	DO ib=1,nbpt ! Loop over domain size
2818
2819	! vecpos: List of positions where textures were not zero
2820	! vecpos(1): number of not null textures found
2821	vecpos = interpweight_ValVecR(soilcolrefrac(ib,:),classnb,zero,'neq')
2822	fopt = vecpos(1)
2823	IF (fopt == classnb) THEN
2824	! All textures are not zero
2825	solt(:) = (/(i,i=1,classnb)/)
2826	ELSE IF (fopt == 0) THEN
2827	WRITE(numout,*)' albedo_surface_soilalb: for point=', ib, ' no soil class!'
2828	ELSE
2829	DO ip = 1,fopt
2830	solt(ip) = vecpos(ip+1)
2831	END DO
2832	END IF
2833
2834	!! 3. Compute the average bare soil albedo parameters
2835
2836	IF ( (fopt .EQ. 0) .OR. (asoilcol(ib) .LT. min_sechiba)) THEN
2837	! Initialize with mean value if no points were interpolated or if no data was found
2838	nbexp = nbexp + 1
2839	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
2840	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
2841	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
2842	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
2843	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
2844	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
2845	ELSE
2846	! If points were interpolated
2847	DO ip=1, fopt
2848	IF ( solt(ip) .LE. classnb) THEN
2849	! Set to zero if the value is below min_sechiba
2850	IF (soilcolrefrac(ib,solt(ip)) < min_sechiba) soilcolrefrac(ib,solt(ip)) = zero
2851
2852	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(solt(ip))*soilcolrefrac(ib,solt(ip))
2853	soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(solt(ip))*soilcolrefrac(ib,solt(ip))
2854	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(solt(ip))*soilcolrefrac(ib,solt(ip))
2855	soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(solt(ip))*soilcolrefrac(ib,solt(ip))
2856	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(solt(ip))* &
2857	soilcolrefrac(ib,solt(ip))
2858	soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(solt(ip))* &
2859	soilcolrefrac(ib,solt(ip))
2860	ELSE
2861	msg = 'The file contains a soil color class which is incompatible with this program'
2862	CALL ipslerr_p(3,'albedo_surface_soilalb',TRIM(msg),'','')
2863	ENDIF
2864	ENDDO
2865	ENDIF
2866
2867	ENDDO
2868
2869	IF ( nbexp .GT. 0 ) THEN
2870	WRITE(numout,*) 'albedo_surface_soilalb _______'
2871	WRITE(numout,*) 'albedo_surface_soilalb: The interpolation of the bare soil albedo had ', nbexp
2872	WRITE(numout,*) 'albedo_surface_soilalb: points without data. This are either coastal points or'
2873	WRITE(numout,*) 'albedo_surface_soilalb: ice covered land.'
2874	WRITE(numout,*) 'albedo_surface_soilalb: The problem was solved by using the average of all soils'
2875	WRITE(numout,*) 'albedo_surface_soilalb: in dry and wet conditions'
2876	WRITE(numout,*) 'albedo_surface_soilalb: Use the diagnostic output field asoilcol to see location of these points'
2877	ENDIF
2878
2879	DEALLOCATE (soilcolrefrac)
2880	DEALLOCATE (variabletypevals)
2881
2882	! Write diagnostics
2883	CALL xios_orchidee_send_field("interp_avail_asoilcol",asoilcol)
2884
2885
2886	IF (printlev_loc >= 3) WRITE(numout,*)' albedo_surface_soilalb ended'
2887
2888	END SUBROUTINE albedo_surface_soilalb
2889
2890	! ----------------------------------------------------------------------------------------
2891
2892	END MODULE albedo_surface

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

Download in other formats: