Context Navigation

condveg.f90 @ 7108

Last change on this file since 7108 was 2680, checked in by sebastiaan.luyssaert, 9 years ago
resolved bug when running in parallel (versions 2620-2664 in the log will not run in parallel
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 94.3 KB

Line
1	! ===============================================================================================================================
2	! MODULE : condveg
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 Initialise, compute and update the surface parameters emissivity,
10	!! roughness and albedo.
11	!!
12	!! \n DESCRIPTION : The module uses 3 settings to control its flow:\n
13	!! 1. :: z0cdrag_ave to choose between two methods to calculate the grid average of
14	!! the roughness height. If set to true: the grid average is calculated by the drag coefficients
15	!! per PFT. If set to false: the grid average is calculated by the logarithm of the
16	!! roughness height per PFT.\n
17	!! 2. :: impaze for choosing surface parameters. If set to false, the values for the
18	!! soil albedo, emissivity and roughness height are set to default values which are read from
19	!! the run.def. If set to true, the user imposes its own values, fixed for the grid point.
20	!! This is useful if one performs site simulations, however,
21	!! it is not recommended to do so for spatialized simulations.
22	!! roughheight_scal imposes the roughness height in (m) ,
23	!! same for emis_scal (in %), albedo_scal (in %), zo_scal (in m)
24	!! Note that these values are only used if 'impaze' is true.\n
25	!! 3. :: alb_bare_model for choosing values of bare soil albedo. If set to TRUE bare
26	!! soil albedo depends on soil wetness. If set to FALSE bare soil albedo is the mean
27	!! values of wet and dry soil albedos.\n
28	!! The surface fluxes are calculated between two levels: the atmospheric level reference
29	!! and the effective roughness height defined as the difference between the mean height of
30	!! the vegetation and the displacement height (zero wind level). Over bare soils, the zero
31	!! wind level is equal to the soil roughness. Over vegetation, the zero wind level
32	!! is increased by the displacement height which depends on the height of the vegetation.
33	!! For a grid point composed of different types of vegetation, an effective surface roughness
34	!! has to be calculated
35	!!
36	!! +++CHECK+++
37	!! NOTE: The dofoco set-up calculates its own snow albedo for VIS and NIR. The explicit
38	!! scheme only calculates it for VIS. To avoid conflicts condveg_snow is not used yet when
39	!! running dofoco. Remove ok_dofoco by a flag that is read from the run.def
40	!! +++++++++++
41	!!
42	!! RECENT CHANGE(S) : None
43	!!
44	!! REFERENCES(S) : None
45	!!
46	!! SVN :
47	!! $HeadURL$
48	!! $Date$
49	!! $Revision$
50	!! \n
51	!_ ================================================================================================================================
52
53	MODULE condveg
54	!
55	USE ioipsl
56	USE xios_orchidee
57	!
58	! modules used :
59	USE constantes
60	USE constantes_soil
61	USE pft_parameters
62	USE qsat_moisture
63	USE interpol_help
64	USE mod_orchidee_para
65	USE ioipsl_para
66	USE albedo
67	USE sechiba_io
68	USE slowproc
69
70	IMPLICIT NONE
71
72	! public routines :
73	! condveg_main only
74	PRIVATE
75	PUBLIC :: condveg_main,condveg_clear
76
77	!
78	! variables used inside condveg module : declaration and initialisation
79	!
80	LOGICAL, SAVE :: l_first_condveg=.TRUE. !! To keep first call's trace
81	!$OMP THREADPRIVATE(l_first_condveg)
82	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_dry(:,:) !! Albedo values for the dry bare soil (unitless)
83	!$OMP THREADPRIVATE(soilalb_dry)
84	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_wet(:,:) !! Albedo values for the wet bare soil (unitless)
85	!$OMP THREADPRIVATE(soilalb_wet)
86	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_moy(:,:) !! Albedo values for the mean bare soil (unitless)
87	!$OMP THREADPRIVATE(soilalb_moy)
88	REAL(r_std), ALLOCATABLE, SAVE :: alb_bare(:,:) !! Mean bare soil albedo for visible and near-infrared
89	!! range (unitless)
90	!$OMP THREADPRIVATE(alb_bare)
91	REAL(r_std), ALLOCATABLE, SAVE :: alb_veget(:,:) !! Mean vegetation albedo for visible and
92	!! near-infrared range (unitless)
93	!$OMP THREADPRIVATE(alb_veget)
94	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: albedo_snow !! Mean snow albedo over visible and near-infrared
95	!! range (unitless)
96	!$OMP THREADPRIVATE(albedo_snow)
97	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: albedo_glob !! Mean surface albedo over visible and
98	!! near-infrared range (unitless)
99	!$OMP THREADPRIVATE(albedo_glob)
100	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: albedo_sveg !! Snow albedo on vegetation (unitless)
101	!$OMP THREADPRIVATE(albedo_sveg)
102
103
104	CONTAINS
105
106	!! ==============================================================================================================================
107	!! SUBROUTINE : condveg_main
108	!!
109	!>\BRIEF Calls the subroutines to initialise the variables, update the variables
110	!! and write out data and restart files.
111	!!
112	!!
113	!! MAIN OUTPUT VARIABLE(S): emis (emissivity), albedo (albedo of
114	!! vegetative PFTs in visible and near-infrared range), z0 (surface roughness height),
115	!! roughheight (grid effective roughness height), soil type (fraction of soil types)
116	!!
117	!!
118	!! REFERENCE(S) : None
119	!!
120	!! FLOWCHART : None
121	!!
122	!! REVISION(S) : None
123	!!
124	!_ ================================================================================================================================
125
126	SUBROUTINE condveg_main (kjit, kjpindex, dtradia, ldrestart_read,&
127	ldrestart_write, index, lalo, neighbours, resolution, contfrac, &
128	veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
129	snow, snow_age, snow_nobio, snow_nobio_age, drysoil_frac, &
130	height, emis, albedo, z0, roughheight, &
131	rest_id, hist_id, hist2_id, lai, sinang, &
132	circ_class_biomass, circ_class_n, lai_split, z0_veg, Isotrop_Abs_Tot_p, &
133	Isotrop_Tran_Tot_p, laieff_fit, albedo_pft, frac_snow_nobio, frac_snow_veg, &
134	! Explicit snow scheme
135	albedo_undersnow, snowdz, snowrho, snowliq, precip_snow, &
136	osfcmelt, temp_air, u, v, &
137	Collim_Abs_Tot, Collim_Alb_Tot, laieff_collim, laieff_isotrop)
138
139	!! 0. Variable and parameter declaration
140
141	!! 0.1 Input variables
142	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
143	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
144	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
145	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier
146	INTEGER(i_std), OPTIONAL, INTENT (in) :: hist2_id !! _History_ file 2 identifier
147	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
148	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for restart file to read
149	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for restart file to write
150	! input fields
151	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
152	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
153	INTEGER(i_std),DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land
154	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
155	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
156	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
157	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
158	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
159	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
160	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
161	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
162	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
163	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
164	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
165	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
166	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
167	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m^2 m^{-2})
168	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: sinang !! the sine of the solar angle from the horizon (unitless)
169	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nparts,nelements), &
170	INTENT(IN) :: circ_class_biomass !! Stem diameter
171	!! @tex $(m)$ @endtex
172	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc), &
173	INTENT(IN) :: circ_class_n !! Number of trees within each circumference
174	REAL(r_std),DIMENSION(nlevels),INTENT(IN) :: lai_split !! temp
175	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
176
177	! Explicit snow scheme
178	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer
179	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density at each snow layer
180	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content at each snow layer
181	LOGICAL, DIMENSION(kjpindex), INTENT(in) :: osfcmelt !! Flag to activate snow melting
182	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snowfall
183	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
184	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u !! Horizontal wind speed, u direction
185	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: v !! Horizontal wind speed, v direction
186
187
188	!! 0.2 Output variables
189	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
190	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
191	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Roughness
192	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
193	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0_veg !! Roughness height of vegetated part (m)
194	REAL(r_std),DIMENSION (:,:,:), &
195	INTENT (out) :: Isotrop_Abs_Tot_p !! Absorbed radiation per level for photosynthesis
196	REAL(r_std),DIMENSION (:,:,:), &
197	INTENT (out) :: Isotrop_Tran_Tot_p !! Transmitted radiation per level for photosynthesis
198	REAL(r_std), DIMENSION(kjpindex,nvm,n_spectralbands), &
199	INTENT (out) :: albedo_pft !! Albedo (two stream radiation transfer model)
200	!! for visible and near-infrared range
201	!! for each PFT (unitless)
202	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
203	!! (unitless ratio)
204	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: frac_snow_veg !! the fraction of ground covered with snow, between 0 and 1
205
206
207	!! 0.3 Modified variables
208	REAL(r_std),DIMENSION (kjpindex), INTENT(inout) :: albedo_undersnow !! Albedo under snow (unitless)
209	REAL(r_std),DIMENSION(:), INTENT(inout) :: Collim_Abs_Tot !! collimated total absorption for a given level
210	REAL(r_std),DIMENSION(:), INTENT(inout) :: Collim_Alb_Tot !! Collimated (direct) total albedo for a given level
211	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: laieff_collim !! Leaf Area Index Effective for direct light
212	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: laieff_isotrop !! Leaf Area Index Effective converts
213	!! 3D lai into 1D lai for two stream
214	!! radiation transfer model...this is for
215	!! isotropic light and only calculated once per day
216	!! @tex $(m^{2} m^{-2})$ @endtex
217
218	!! 0.4 Local variables
219	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
220
221	!_ ================================================================================================================================
222
223	!! 1. Call subroutines to start initialisation
224
225	! Is TRUE if condveg is called the first time
226	IF (l_first_condveg) THEN
227	! Document what the programm is doing
228	IF (long_print) WRITE (numout,*) ' l_first_condveg : call condveg_init '
229
230	! Call the subroutine 'condveg_init' to allocate local array
231	CALL condveg_init (kjit, kjpindex, index, veget, &
232	lalo, neighbours, resolution, contfrac, rest_id, &
233	! Explicit snow scheme
234	albedo_undersnow)
235
236	! Call the subroutine 'condveg_var_init' to initialise local array
237	! Reads in map for soil albedo
238	CALL condveg_var_init (kjpindex, veget, veget_max, frac_nobio, &
239	totfrac_nobio, drysoil_frac, zlev, height, &
240	circ_class_biomass, circ_class_n, emis, albedo, &
241	z0, roughheight, lai, sinang, &
242	snow, snow_age, snow_nobio, snow_nobio_age, albedo_snow, &
243	albedo_glob, lai_split, z0_veg, Isotrop_Abs_Tot_p,&
244	Isotrop_Tran_Tot_p, laieff_fit, albedo_pft, frac_snow_nobio, &
245	frac_snow_veg, &
246	Collim_Abs_Tot, Collim_Alb_Tot, laieff_collim, laieff_isotrop)
247
248	! Explicit snow scheme
249	! Call subroutine 'condveg snow'
250	IF (.NOT. control%ok_dofoco) THEN
251	CALL condveg_snow (kjpindex, veget, veget_max, frac_nobio, &
252	totfrac_nobio, snow, snow_age, snow_nobio, snow_nobio_age, &
253	z0, snowdz, snowrho, snowliq, dtradia, &
254	precip_snow, osfcmelt, temp_air, u, v, &
255	albedo, albedo_snow, albedo_glob, albedo_undersnow, &
256	frac_snow_veg, frac_snow_nobio)
257	ENDIF
258
259	RETURN
260
261	ENDIF
262
263	!! 2. Call subroutines to write restart files and data
264	! To write restart files for soil albedo variables
265	IF (ldrestart_write) THEN
266
267	var_name = 'soilalbedo_dry'
268	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_dry, 'scatter', nbp_glo, index_g)
269
270	var_name = 'soilalbedo_wet'
271	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_wet, 'scatter', nbp_glo, index_g)
272
273	var_name = 'soilalbedo_moy'
274	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_moy, 'scatter', nbp_glo, index_g)
275
276	IF (.NOT. control%ok_dofoco .AND. ok_explicitsnow) THEN
277	var_name = 'albedo_undersnow'
278	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, albedo_undersnow, 'scatter', &
279	nbp_glo, index_g)
280
281	var_name = 'albedo_sveg'
282	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, albedo_sveg, 'scatter', &
283	nbp_glo, index_g)
284
285	END IF
286
287	RETURN
288
289	ENDIF
290
291
292	!! 3. Call subroutines to update fields
293
294	! Call the routine 'condveg_var_update' to update the fields of albedo,
295	! emissivity and roughness
296	CALL condveg_var_update (kjpindex, veget, veget_max, frac_nobio, &
297	totfrac_nobio, drysoil_frac, zlev, height, &
298	circ_class_biomass, circ_class_n, emis, albedo, &
299	z0, roughheight, lai, sinang, &
300	snow, snow_age, snow_nobio, snow_nobio_age, albedo_snow, &
301	albedo_glob, lai_split,z0_veg, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
302	Collim_Abs_Tot, Collim_Alb_Tot, laieff_collim, laieff_isotrop, &
303	laieff_fit, albedo_pft, frac_snow_nobio, frac_snow_veg)
304
305	! Explicit snow scheme
306	! Update snow parameters by calling subroutine 'condveg_snow'
307	IF (.NOT. control%ok_dofoco) THEN
308	CALL condveg_snow (kjpindex, veget, veget_max, frac_nobio, &
309	totfrac_nobio, snow, snow_age, snow_nobio, snow_nobio_age, &
310	z0, snowdz, snowrho, snowliq, dtradia, &
311	precip_snow, osfcmelt, temp_air, u, v, &
312	albedo, albedo_snow, albedo_glob, albedo_undersnow, &
313	frac_snow_veg, frac_snow_nobio)
314	ENDIF
315
316	! If this logical flag is set to true, the model
317	! will output all its data according to the ALMA
318	! convention. To facilitate the exchange of forcing
319	! data for land-surface schemes and the results produced
320	! by these schemes, ALMA (Assistance for Land-surface Modelling
321	! activities) proposes to establish standards.
322	! http://www.lmd.jussieu.fr/~polcher/ALMA/
323	CALL xios_orchidee_send_field("soilalb_vis",alb_bare(:,1))
324	CALL xios_orchidee_send_field("soilalb_nir",alb_bare(:,2))
325	CALL xios_orchidee_send_field("vegalb_vis",alb_veget(:,1))
326	CALL xios_orchidee_send_field("vegalb_nir",alb_veget(:,2))
327	CALL xios_orchidee_send_field("albedo_alma",albedo_glob)
328	CALL xios_orchidee_send_field("SAlbedo",albedo_snow)
329
330	IF ( almaoutput ) THEN
331	CALL histwrite_p(hist_id, 'Albedo', kjit, albedo_glob, kjpindex, index)
332	CALL histwrite_p(hist_id, 'SAlbedo', kjit, albedo_snow, kjpindex, index)
333	IF ( hist2_id > 0 ) THEN
334	CALL histwrite_p(hist2_id, 'Albedo', kjit, albedo_glob, kjpindex, index)
335	CALL histwrite_p(hist2_id, 'SAlbedo', kjit, albedo_snow, kjpindex, index)
336	ENDIF
337	ELSE
338	CALL histwrite_p(hist_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
339	CALL histwrite_p(hist_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
340	CALL histwrite_p(hist_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
341	CALL histwrite_p(hist_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
342	IF ( hist2_id > 0 ) THEN
343	CALL histwrite_p(hist2_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
344	CALL histwrite_p(hist2_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
345	CALL histwrite_p(hist2_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
346	CALL histwrite_p(hist2_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
347	ENDIF
348	ENDIF
349
350	IF (long_print) WRITE (numout,*)' condveg_main done '
351
352	END SUBROUTINE condveg_main
353
354	!! ==============================================================================================================================
355	!! SUBROUTINE : condveg_init
356	!!
357	!>\BRIEF Dynamic allocation of the variables, the dimensions of the
358	!! variables are determined by user-specified settings.
359	!!
360	!! DESCRIPTION : The domain size (:: kjpindex) is used to allocate the correct
361	!! dimensions to all variables in condveg.
362	!!
363	!! RECENT CHANGE(S): None
364	!!
365	!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
366	!! variables. However, the routine allocates memory and builds new indexing
367	!! variables for later use.
368	!!
369	!! REFERENCE(S) : None
370	!!
371	!! FLOWCHART : None
372	!! \n
373	!_ ================================================================================================================================
374
375	SUBROUTINE condveg_init (kjit, kjpindex, index, veget, &
376	lalo, neighbours, resolution, contfrac, rest_id,&
377	!Explicit snow scheme
378	albedo_undersnow)
379
380	!! 0. Variable and parameter declaration
381
382	!! 0.1. Input variables
383
384	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
385	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
386	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Index for the points on the map (unitless)
387	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in):: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
388	!! (m^2 m^{-2})
389	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates (degree)
390	INTEGER(i_std),DIMENSION (kjpindex,8), INTENT(in):: neighbours !! Neighbouring land grid cell
391	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size of grid in x and y direction (m)
392	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: contfrac !! Fraction of land in each grid box
393	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier
394
395
396	!! 0.2. Output variables
397
398	REAL(r_std),DIMENSION(kjpindex),INTENT(out) :: albedo_undersnow !! Albedo under snow
399
400
401	!! 0.3 Modified variables
402
403	!! 0.4 Local variables
404
405	INTEGER(i_std) :: ji !! Index
406	INTEGER(i_std) :: ier !! Check errors in memory allocation
407	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
408
409	!_ ================================================================================================================================
410
411	!! 1. Choice of calculation of snow albedo and soil albedo
412
413	! Is TRUE if condveg is called the first time
414	IF (l_first_condveg) THEN
415	!
416	l_first_condveg=.FALSE.
417
418	!! 2. Allocate all albedo variables
419
420	! Dry soil albedo
421	ALLOCATE (soilalb_dry(kjpindex,2),stat=ier)
422	IF (ier.NE.0) THEN
423	WRITE (numout,*) ' error in soilalb_dry allocation'
424	CALL ipslerr_p (3,'condveg.f90','condveg_init','Pb in allocation for soilalb_dry','')
425	END IF
426	soilalb_dry(:,:) = val_exp
427
428	! Wet soil albedo
429	ALLOCATE (soilalb_wet(kjpindex,2),stat=ier)
430	IF (ier.NE.0) THEN
431	WRITE (numout,*) ' error in soilalb_wet allocation'
432	CALL ipslerr_p (3,'condveg.f90','condveg_init','Pb in allocation for soilalb_wet','')
433	END IF
434	soilalb_wet(:,:) = val_exp
435
436	! Mean soil albedo
437	ALLOCATE (soilalb_moy(kjpindex,2),stat=ier)
438	IF (ier.NE.0) THEN
439	WRITE (numout,*) ' error in soilalb_moy allocation'
440	CALL ipslerr_p (3,'condveg.f90','condveg_init','Pb in allocation for soilalb_moy','')
441	END IF
442	soilalb_moy(:,:) = val_exp
443
444	! Snow albedo of vegetative PFTs
445	ALLOCATE (albedo_snow(kjpindex),stat=ier)
446	IF (ier.NE.0) THEN
447	WRITE (numout,*) ' error in soilalb_snow allocation'
448	CALL ipslerr_p (3,'condveg.f90','condveg_init','Pb in allocation for albedo_snow','')
449	END IF
450
451	! Mean vegetation albedo over visible and near-infrared range
452	ALLOCATE (albedo_glob(kjpindex),stat=ier)
453	IF (ier.NE.0) THEN
454	WRITE (numout,*) ' error in soilalb_glob allocation'
455	CALL ipslerr_p (3,'condveg.f90','condveg_init','Pb in allocation for albedo_glob','')
456	END IF
457
458	! Mean bare soil albedo
459	ALLOCATE (alb_bare(kjpindex,2),stat=ier)
460	IF (ier.NE.0) THEN
461	WRITE (numout,*) ' error in soilalb_bar allocation'
462	CALL ipslerr_p (3,'condveg.f90','condveg_init','Pb in allocation for alb_bar','')
463	END IF
464	alb_bare(:,:) = val_exp
465
466	! Mean vegetation albedo
467	ALLOCATE (alb_veget(kjpindex,2),stat=ier)
468	IF (ier.NE.0) THEN
469	WRITE (numout,*) ' error in soilalb_veget allocation'
470	CALL ipslerr_p (3,'condveg.f90','condveg_init','Pb in allocation for alb_veget','')
471	END IF
472	alb_veget(:,:) = val_exp
473
474	IF (.NOT. control%ok_dofoco .AND. ok_explicitsnow) THEN
475	ALLOCATE (albedo_sveg(kjpindex,nvm),stat=ier)
476	IF (ier.NE.0) THEN
477	WRITE (numout,*) ' error in soilalb_sveg allocation'
478	CALL ipslerr_p(3,'condveg_init','Pb in allocation for albedo_sveg','','')
479	ENDIF
480	!+++CHECK+++
481	! Not initialized in the trunk
482	!!$ alb_sveget(:,:) = val_exp
483	!+++++++++++
484	END IF
485
486
487	!! 3. Initialise bare soil albedo
488
489	! dry soil albedo
490	var_name= 'soilalbedo_dry'
491	CALL ioconf_setatt_p('UNITS', '-')
492	CALL ioconf_setatt_p('LONG_NAME','Dry bare soil albedo')
493	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_dry, "gather", &
494	nbp_glo, index_g)
495
496	! wet soil albedo
497	var_name= 'soilalbedo_wet'
498	CALL ioconf_setatt_p('UNITS', '-')
499	CALL ioconf_setatt_p('LONG_NAME','Wet bare soil albedo')
500	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_wet, "gather", &
501	nbp_glo, index_g)
502
503	! mean soil albedo
504	var_name= 'soilalbedo_moy'
505	CALL ioconf_setatt_p('UNITS', '-')
506	CALL ioconf_setatt_p('LONG_NAME','Mean bare soil albedo')
507	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_moy, "gather", &
508	nbp_glo, index_g)
509
510	IF (.NOT.control%ok_dofoco .AND. ok_explicitsnow) THEN
511
512	! Albedo under snow
513	var_name= 'albedo_undersnow'
514	CALL ioconf_setatt('UNITS', '-')
515	CALL ioconf_setatt('LONG_NAME','albedo of the soil/snow interface')
516	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit,.TRUE.,albedo_undersnow, "gather", &
517	nbp_glo, index_g)
518	CALL setvar_p (albedo_undersnow, val_exp,'NOKEYWORD', 0.0)
519
520	! Snow albedo on vegetation
521	var_name = 'albedo_sveg'
522	CALL ioconf_setatt('UNITS', '-')
523	CALL ioconf_setatt('LONG_NAME','albedo_sveg')
524	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE.,albedo_sveg, "gather", &
525	nbp_glo, index_g)
526	CALL setvar_p (albedo_sveg, val_exp,'NOKEYWORD', 0.0)
527
528	END IF
529
530
531	! check if we have real values and not only missing values
532	IF ( MINVAL(soilalb_wet) .EQ. MAXVAL(soilalb_wet) .AND. MAXVAL(soilalb_wet) .EQ. val_exp .OR.&
533	& MINVAL(soilalb_dry) .EQ. MAXVAL(soilalb_dry) .AND. MAXVAL(soilalb_dry) .EQ. val_exp .OR.&
534	& MINVAL(soilalb_moy) .EQ. MAXVAL(soilalb_moy) .AND. MAXVAL(soilalb_moy) .EQ. val_exp) THEN
535	CALL condveg_soilalb(kjpindex, lalo, neighbours, resolution, contfrac, soilalb_dry,soilalb_wet)
536	WRITE(numout,*) '---> val_exp ', val_exp
537	WRITE(numout,*) '---> ALBEDO_wet VIS:', MINVAL(soilalb_wet(:,ivis)), MAXVAL(soilalb_wet(:,ivis))
538	WRITE(numout,*) '---> ALBEDO_wet NIR:', MINVAL(soilalb_wet(:,inir)), MAXVAL(soilalb_wet(:,inir))
539	WRITE(numout,*) '---> ALBEDO_dry VIS:', MINVAL(soilalb_dry(:,ivis)), MAXVAL(soilalb_dry(:,ivis))
540	WRITE(numout,*) '---> ALBEDO_dry NIR:', MINVAL(soilalb_dry(:,inir)), MAXVAL(soilalb_dry(:,inir))
541	WRITE(numout,*) '---> ALBEDO_moy VIS:', MINVAL(soilalb_moy(:,ivis)), MAXVAL(soilalb_moy(:,ivis))
542	WRITE(numout,*) '---> ALBEDO_moy NIR:', MINVAL(soilalb_moy(:,inir)), MAXVAL(soilalb_moy(:,inir))
543	ENDIF
544
545	ELSE
546	WRITE (numout,*) ' l_first_condveg false . we stop '
547	CALL ipslerr_p (3,'condveg.f90','condveg_init','','')
548	ENDIF
549
550	IF (long_print) WRITE (numout,*) ' condveg_init done '
551
552	END SUBROUTINE condveg_init
553
554	!! ==============================================================================================================================
555	!! SUBROUTINE : condveg_clear
556	!!
557	!>\BRIEF Deallocate albedo variables
558	!!
559	!! DESCRIPTION : None
560	!!
561	!! RECENT CHANGE(S): None
562	!!
563	!! MAIN OUTPUT VARIABLE(S): None
564	!!
565	!! REFERENCES : None
566	!!
567	!! FLOWCHART : None
568	!! \n
569	!_ ================================================================================================================================
570
571	SUBROUTINE condveg_clear ()
572
573	l_first_condveg=.TRUE.
574
575	! Dry soil albedo
576	IF (ALLOCATED (soilalb_dry)) DEALLOCATE (soilalb_dry)
577	! Wet soil albedo
578	IF (ALLOCATED(soilalb_wet)) DEALLOCATE (soilalb_wet)
579	! Mean soil albedo
580	IF (ALLOCATED(soilalb_moy)) DEALLOCATE (soilalb_moy)
581	! Mean snow albedo
582	IF (ALLOCATED(albedo_snow)) DEALLOCATE (albedo_snow)
583	! Mean albedo
584	IF (ALLOCATED(albedo_glob)) DEALLOCATE (albedo_glob)
585	! Mean albedo of bare soil
586	IF (ALLOCATED(alb_bare)) DEALLOCATE (alb_bare)
587	! Mean vegetation albedo
588	IF (ALLOCATED(alb_veget)) DEALLOCATE (alb_veget)
589	! Snow albedo on vegetation
590	IF (ALLOCATED(albedo_sveg)) DEALLOCATE (albedo_sveg)
591
592	END SUBROUTINE condveg_clear
593
594	!! ==============================================================================================================================
595	!! SUBROUTINE : condveg_var_init
596	!!
597	!>\BRIEF Initialisation of local array and calculation of emissivity,
598	!! albedo and roughness height.
599	!!
600	!! DESCRIPTION : None
601	!!
602	!! MAIN OUTPUT VARIABLE(S): None
603	!!
604	!! REFERENCE(S) : None
605	!!
606	!! FLOWCHART : None
607	!! \n
608	!_ ================================================================================================================================
609
610	SUBROUTINE condveg_var_init (kjpindex, veget, veget_max, frac_nobio, &
611	totfrac_nobio, drysoil_frac, zlev, height, &
612	circ_class_biomass, circ_class_n, emis, albedo, &
613	z0, roughheight, lai, sinang, &
614	snow, snow_age, snow_nobio, snow_nobio_age, albedo_snow, &
615	albedo_glob, lai_split,z0_veg, Isotrop_Abs_Tot_p, &
616	Isotrop_Tran_Tot_p, laieff_fit, albedo_pft, frac_snow_nobio,&
617	frac_snow_veg, &
618	Collim_Abs_Tot, Collim_Alb_Tot, laieff_collim, laieff_isotrop)
619
620	!! 0. Variable and parameter declaration
621
622	!! 0.1 Input variables
623
624	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
625	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
626	!! (m^2 m^{-2})
627	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
628	!! (= ind*cn_ind) (m^2/m^{-2})
629	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
630	!! continental ice, lakes, etc. (unitless)
631	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces, i.e.
632	!! continental ice, lakes, etc. (unitless)
633	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of dry soil in visible range (unitless)
634	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
635	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
636	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
637	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: sinang !! the sine of the solar angle from the horizon (unitless)
638	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
639	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
640	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age (days)
641	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
642	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nparts,nelements), &
643	INTENT(IN) :: circ_class_biomass !! Stem diameter
644	!! @tex $(m)$ @endtex
645	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc), INTENT(IN) :: circ_class_n !! Number of trees within each circumference
646	REAL(r_std),DIMENSION(nlevels),INTENT(IN) :: lai_split
647	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
648
649	!! 0.2 Output varialbes
650
651	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Surface emissivity (unitless)
652	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo of vegetative PFTs in visible and
653	!! near-infrared range
654	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Soil roughness height (m)
655	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Grid effective roughness height (m)
656	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
657	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_glob !! Mean albedo (unitless ratio)
658	REAL(r_std),DIMENSION (:,:,:), &
659	INTENT (out) :: Isotrop_Abs_Tot_p !! Absorbed radiation per level for photosynthesis
660	REAL(r_std),DIMENSION (:,:,:), &
661	INTENT (out) :: Isotrop_Tran_Tot_p !! Transmitted radiation per level for photosynthesis
662
663	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0_veg !! Roughness height of vegetated part (m)
664	REAL(r_std), DIMENSION(kjpindex,nvm,n_spectralbands), &
665	INTENT (out) :: albedo_pft !! Albedo (two stream radiation transfer model)
666	!! for visible and near-infrared range
667	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
668	!! (unitless ratio)
669	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: frac_snow_veg !! the fraction of ground covered with snow, between 0 and 1
670	!! for each PFT (unitless)
671
672	!! 0.3 Modified variables
673	REAL(r_std),DIMENSION(:), INTENT(inout) :: Collim_Abs_Tot !! collimated total absorption for a given level
674	REAL(r_std),DIMENSION(:), INTENT(inout) :: Collim_Alb_Tot !! Collimated (direct) total albedo for a given level
675	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: laieff_collim !! Leaf Area Index Effective for direct light
676	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: laieff_isotrop !! Leaf Area Index Effective converts
677	!! 3D lai into 1D lai for two stream
678	!! radiation transfer model...this is for
679	!! isotropic light and only calculated once per day
680	!! @tex $(m^{2} m^{-2})$ @endtex
681	!! 0.4 Local variables
682
683	INTEGER(i_std) :: ier !! Check errors
684	INTEGER(i_std) :: jv !! Index
685
686	!_ ================================================================================================================================
687
688
689	!! 1. Choice of parameter setting
690
691	! This switch is for choosing surface parameters.
692	! If it is set to true, the values for the soil roughness, emissivity
693	! and albedo are set to default values which are read in from the run.def.
694	! This is useful if one performs a single site simulations,
695	! it is not recommended to do so for global simulations.
696
697
698	!! 2. Calculate emissivity
699
700	! If TRUE read in default values for emissivity
701
702	IF ( impaze ) THEN
703	!
704	emis(:) = emis_scal
705	!
706	! If FALSE set emissivity to 1.
707	ELSE
708	emis_scal = un
709	emis(:) = emis_scal
710	ENDIF
711
712	!! 3. Calculate roughness height
713
714	! Chooses between two methods to calculate the grid average of the roughness.
715	! If set to true: The grid average is calculated by averaging the drag coefficients over PFT.
716	! If set to false: The grid average is calculated by averaging the logarithm of the roughness length per PFT.
717
718	! TRUE read in default values for roughness height
719	IF ( impaze ) THEN
720
721	z0(:) = z0_scal
722	roughheight(:) = roughheight_scal
723
724	! If FALSE calculate roughness height
725	ELSE
726
727	IF ( z0cdrag_ave ) THEN
728	CALL condveg_z0cdrag(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
729	& height, z0, roughheight, z0_veg)
730	ELSE
731	CALL condveg_z0logz(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
732	& z0, roughheight, z0_veg)
733	ENDIF
734
735	ENDIF
736
737	!! 4. Calculate albedo
738
739	! If TRUE read in default values for albedo
740	IF ( impaze ) THEN
741	!
742	albedo(:,ivis) = albedo_scal(ivis)
743	albedo(:,inir) = albedo_scal(inir)
744	!
745	ELSE
746
747	! check to see which type of albedo we are using
748	SELECTCASE (albedo_type)
749	CASE ('standard')
750	CALL albedo_calc (kjpindex, drysoil_frac, veget,&
751	soilalb_dry, soilalb_wet, soilalb_moy, tot_bare_soil, alb_veget, alb_bare, albedo)
752	! Initialize these variables here. They are only computed in the two_stream model.
753	Isotrop_Abs_Tot_p(:,:,:)=undef_sechiba
754	Isotrop_Tran_Tot_p(:,:,:)=undef_sechiba
755	frac_snow_nobio(:,:)=undef_sechiba
756	frac_snow_veg(:,:)=undef_sechiba
757
758	CASE ('pinty')
759	! We need to initialize albedo here to give a realistic guess to the coupled
760	! model.
761	CALL albedo_two_stream(kjpindex, nlevels_tot, drysoil_frac, lai, veget_max, &
762	sinang, soilalb_dry, soilalb_wet, frac_nobio, soilalb_moy, &
763	alb_bare, albedo, snow, snow_age, snow_nobio, &
764	snow_nobio_age, albedo_snow, albedo_glob, &
765	circ_class_biomass, circ_class_n,&
766	lai_split, z0_veg, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
767	laieff_fit, albedo_pft, frac_snow_nobio, frac_snow_veg, &
768	Collim_Abs_Tot, Collim_Alb_Tot, laieff_collim, laieff_isotrop)
769
770
771	CASE DEFAULT
772	WRITE(numout,*) "Unsupported albedo type. This should have been caught in intersurf!"
773	CALL ipslerr_p (3,'condveg.f90','condveg_var_init','','')
774	END SELECT
775
776
777	ENDIF
778
779	!! Calculate snow albedo. This is tricky because the new scheme requires information about snow to be
780	!! passed to it, while the old scheme uses existing albedo data
781	SELECTCASE (albedo_type)
782	CASE ('standard')
783	CALL calculate_surface_albedo_with_snow(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
784	snow, snow_age, snow_nobio, snow_nobio_age, tot_bare_soil, albedo, albedo_snow, albedo_glob)
785	CASE ('pinty')
786	CASE DEFAULT
787	WRITE(numout,*) "Unsupported albedo type. This should have been caught in intersurf!"
788	CALL ipslerr_p (3,'condveg.f90','condveg_var_init','','')
789	END SELECT
790
791	IF (long_print) WRITE (numout,*) ' condveg_var_init done '
792
793	END SUBROUTINE condveg_var_init
794
795
796
797	!! ==============================================================================================================================
798	!! SUBROUTINE : condveg_var_update
799	!!
800	!>\BRIEF Update emissivity, albedo and roughness height.
801	!!
802	!! DESCRIPTION : None
803	!!
804	!! MAIN OUTPUT VARIABLE(S): \n
805	!!
806	!! REFERENCE(S) : None
807	!!
808	!! FLOWCHART : None
809	!! \n
810	!_ ================================================================================================================================
811
812	SUBROUTINE condveg_var_update (kjpindex, veget, veget_max, frac_nobio, &
813	totfrac_nobio, drysoil_frac, zlev, height, &
814	circ_class_biomass, circ_class_n, emis, albedo, &
815	z0, roughheight, lai, sinang, &
816	snow, snow_age, snow_nobio, snow_nobio_age, albedo_snow, &
817	albedo_glob, lai_split, z0_veg, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
818	Collim_Abs_Tot, Collim_Alb_Tot, laieff_collim, laieff_isotrop, &
819	laieff_fit, albedo_pft, frac_snow_nobio, frac_snow_veg)
820
821	!! 0. Variable and parameter declaration
822
823	!! 0.1 Input variables
824
825	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
826	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
827	!! (m^2 m^{-2})
828	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
829	!! (= ind*cn_ind) (m^2 m^{-2})
830	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
831	!! continental ice, lakes, etc. (unitless)
832	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces, i.e.
833	!! continental ice, lakes, etc. (unitless)
834	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of dry soil in visible range (unitless)
835	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
836	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
837	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
838	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: sinang !! the sine of the solar angle from the horizon (unitless)
839	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
840	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
841	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age (days)
842	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
843	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nparts,nelements), &
844	INTENT(IN) :: circ_class_biomass !! Stem diameter
845	!! @tex $(m)$ @endtex
846	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc), INTENT(IN) :: circ_class_n !! Number of trees within each circumference
847	REAL(r_std),DIMENSION(nlevels),INTENT(IN) :: lai_split
848	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
849
850
851	!! 0.2 Output variables
852
853	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity (unitless)
854	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo of vegetative PFTs in visible and
855	!! near-infrared range
856	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Roughness height (m)
857	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Grid effective roughness height (m)
858	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
859	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_glob !! Mean albedo (unitless ratio)
860	REAL(r_std),DIMENSION (:,:,:), &
861	INTENT (out) :: Isotrop_Abs_Tot_p !!Absorbed radiation per level for photosynthesis
862	REAL(r_std),DIMENSION (:,:,:), &
863	INTENT (out) :: Isotrop_Tran_Tot_p !!Transmitted radiation per level for photosynthesis
864
865	REAL(r_std), DIMENSION(kjpindex), INTENT (out) :: z0_veg !! Roughness height of vegetated part (m)
866	REAL(r_std), DIMENSION(kjpindex,nvm,n_spectralbands), &
867	INTENT (out) :: albedo_pft !! Albedo (two stream radiation transfer model)
868	!! for visible and near-infrared range
869	!! for each PFT (unitless)
870	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out) :: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
871	!! (unitless ratio)
872	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: frac_snow_veg !! the fraction of ground covered with snow, between 0 and 1
873
874	!! 0.3 Modified variables
875	REAL(r_std),DIMENSION(:), INTENT (inout) :: Collim_Abs_Tot !! collimated total absorption for a given level * for modified enerbil
876	REAL(r_std),DIMENSION(:), INTENT (inout) :: Collim_Alb_Tot !! Collimated (direct) total albedo for a given level * for modified enerbil
877	REAL(r_std), DIMENSION(:,:,:), INTENT (inout) :: laieff_collim !! Leaf Area Index Effective for direct light
878	REAL(r_std), DIMENSION(:,:,:), INTENT (inout) :: laieff_isotrop !! Leaf Area Index Effective converts
879	!! 3D lai into 1D lai for two stream
880	!! radiation transfer model...this is for
881	!! isotropic light and only calculated once per day
882	!! @tex $(m^{2} m^{-2})$ @endtex
883
884	!! 0.4 Local variables
885	INTEGER(i_std) :: ji, jv !! Indeces
886	!_ ================================================================================================================================
887
888	!! 1. Calculate emissivity
889
890	emis(:) = emis_scal
891
892	!! 2. Calculate roughness height
893
894	! If TRUE read in prescribed values for roughness height
895	IF ( impaze ) THEN
896
897	DO ji = 1, kjpindex
898	z0(ji) = z0_scal
899	roughheight(ji) = roughheight_scal
900	ENDDO
901
902	! Calculate roughness height
903	ELSE
904
905	IF ( z0cdrag_ave ) THEN
906	CALL condveg_z0cdrag (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, height, &
907	& z0, roughheight, z0_veg)
908	ELSE
909	CALL condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
910	& z0, roughheight, z0_veg)
911	ENDIF
912
913	ENDIF
914
915	!! 3. Calculate albedo
916	! This has to be done after the calculation of the roughness height, since it uses
917	! the roughness height for calculating the snow coverage fraction.
918
919	! If TRUE read in prescribed values for albedo
920	IF ( impaze ) THEN
921
922	albedo(:,ivis) = albedo_scal(ivis)
923	albedo(:,inir) = albedo_scal(inir)
924
925	! If FALSE calculate albedo from ORCHIDEE default values
926	ELSE
927
928	! check to see which type of albedo we are using
929	SELECTCASE (albedo_type)
930	CASE ('standard')
931	CALL albedo_calc (kjpindex, drysoil_frac, veget,&
932	soilalb_dry, soilalb_wet, soilalb_moy, tot_bare_soil,alb_veget, alb_bare, albedo)
933	! Initialize these variables here. They are only computed in the two_stream model.
934	Isotrop_Abs_Tot_p(:,:,:)=undef_sechiba
935	Isotrop_Tran_Tot_p(:,:,:)=undef_sechiba
936	frac_snow_nobio(:,:)=undef_sechiba
937	frac_snow_veg(:,:)=undef_sechiba
938	CASE ('pinty')
939
940
941	CALL albedo_two_stream(kjpindex, nlevels_tot, drysoil_frac, lai, veget_max, &
942	sinang, soilalb_dry, soilalb_wet, frac_nobio, soilalb_moy, &
943	alb_bare, albedo, snow, snow_age, snow_nobio, &
944	snow_nobio_age, albedo_snow, albedo_glob, &
945	circ_class_biomass, circ_class_n, lai_split, z0_veg, Isotrop_Abs_Tot_p,&
946	Isotrop_Tran_Tot_p, laieff_fit, albedo_pft, frac_snow_nobio,&
947	frac_snow_veg, Collim_Abs_Tot, Collim_Alb_Tot, &
948	laieff_collim, laieff_isotrop)
949
950
951	CASE DEFAULT
952	WRITE(numout,*) "Unsupported albedo type. This should have been caught in intersurf!"
953	CALL ipslerr_p (3,'condveg.f90','condveg_var_update','','')
954	END SELECT
955
956	ENDIF
957
958	!! Calculate snow albedo.
959	! The scheme of Pinty et al calculates the snow albedo in the albedo_two_stream routine,
960	! since it changes the background albedo. The old scheme calculates it separately.
961	SELECTCASE (albedo_type)
962	CASE ('standard')
963	CALL calculate_surface_albedo_with_snow(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
964	snow, snow_age, snow_nobio, snow_nobio_age, tot_bare_soil, albedo, albedo_snow, albedo_glob)
965	CASE ('pinty')
966	CASE DEFAULT
967	WRITE(numout,*) "Unsupported albedo type. This should have been caught in intersurf!"
968	CALL ipslerr_p (3,'condveg.f90','condveg_var_update','','')
969	END SELECT
970
971	IF (long_print) WRITE (numout,*) ' condveg_var_update done '
972
973	END SUBROUTINE condveg_var_update
974
975	!! ==============================================================================================================================\n
976	!! SUBROUTINE : condveg_snow
977	!!
978	!>\BRIEF Calcuating snow albedo and snow cover fraction, which are then used to
979	!! update the gridbox surface albedo following Chalita and Treut (1994).
980	!!
981	!! DESCRIPTION : The snow albedo scheme presented below belongs to prognostic albedo
982	!! category, i.e. the snow albedo value at a time step depends on the snow albedo value
983	!! at the previous time step.
984	!!
985	!! First, the following formula (described in Chalita and Treut 1994) is used to describe
986	!! the change in snow albedo with snow age on each PFT and each non-vegetative surfaces,
987	!! i.e. continental ice, lakes, etc.: \n
988	!! \latexonly
989	!! \input{SnowAlbedo.tex}
990	!! \endlatexonly
991	!! \n
992	!! Where snowAge is snow age, tcstSnowa is a critical aging time (tcstSnowa=5 days)
993	!! snowaIni and snowaIni+snowaDec corresponds to albedos measured for aged and
994	!! fresh snow respectively, and their values for each PFT and each non-vegetative surfaces
995	!! is precribed in in constantes_veg.f90.\n
996	!! In order to estimate gridbox snow albedo, snow albedo values for each PFT and
997	!! each non-vegetative surfaces with a grid box are weightedly summed up by their
998	!! respective fractions.\n
999	!! Secondly, the snow cover fraction is computed as:
1000	!! \latexonly
1001	!! \input{SnowFraction.tex}
1002	!! \endlatexonly
1003	!! \n
1004	!! Where fracSnow is the fraction of snow on total vegetative or total non-vegetative
1005	!! surfaces, snow is snow mass (kg/m^2) on total vegetated or total nobio surfaces.\n
1006	!! Finally, the surface albedo is then updated as the weighted sum of fracSnow, total
1007	!! vegetated fraction, total nobio fraction, gridbox snow albedo, and previous
1008	!! time step surface albedo.
1009	!!
1010	!! RECENT CHANGE(S): None
1011	!!
1012	!! MAIN OUTPUT VARIABLE(S): :: albedo; surface albedo. :: albedo_snow; snow
1013	!! albedo
1014	!!
1015	!! REFERENCE(S) :
1016	!! Chalita, S. and H Le Treut (1994), The albedo of temperate and boreal forest and
1017	!! the Northern Hemisphere climate: a sensitivity experiment using the LMD GCM,
1018	!! Climate Dynamics, 10 231-240.
1019	!!
1020	!! FLOWCHART : None
1021	!! \n
1022	!_ ================================================================================================================================
1023
1024	SUBROUTINE condveg_snow (kjpindex, veget, veget_max, frac_nobio, &
1025	totfrac_nobio, snow, snow_age, snow_nobio, snow_nobio_age, &
1026	z0, snowdz, snowrho, snowliq, dtradia, &
1027	precip_snow, osfcmelt, temp_air, u, v, &
1028	albedo, albedo_snow, albedo_glob, albedo_undersnow, &
1029	frac_snow_veg, frac_snow_nobio)
1030
1031
1032	!! 0. Variable and parameter declaration
1033
1034	!! 0.1 Input variables
1035
1036	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1037	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1038	!! (m^2 m^{-2})
1039	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max
1040	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
1041	!! continental ice, lakes, etc. (unitless)
1042	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces, i.e.
1043	!! continental ice, lakes, etc. (unitless)
1044	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
1045	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
1046	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age (days)
1047	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
1048	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: z0 !! Roughness length
1049	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer
1050	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density at each snow layer
1051	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowliq !! Liquid water content at each snow layer
1052	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
1053	REAL(r_std), DIMENSION(kjpindex),INTENT(in) :: precip_snow !! Snowfall
1054	LOGICAL, DIMENSION(kjpindex), INTENT(IN) :: osfcmelt !! Flag to activate snow melting
1055	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
1056	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u !! Horizontal wind speed
1057	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: v !! Horizontal wind speed
1058
1059	!! 0.2 Output variables
1060
1061	REAL(r_std),DIMENSION (kjpindex,2), INTENT (inout) :: albedo !! Albedo (unitless ratio)
1062	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
1063	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_glob !! Mean albedo (unitless ratio)
1064	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: albedo_undersnow!! Under snow albedo (unitless)
1065	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: frac_snow_veg !! Fraction of snow on vegetation (unitless ratio)
1066	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
1067	!! (unitless ratio)
1068
1069	!! 0.3 Modified variables
1070
1071	!! 0.4 Local variables
1072
1073	INTEGER(i_std) :: ji, jv, jb,iv !! indices (unitless)
1074	REAL(r_std), DIMENSION(kjpindex) :: snowa_veg !! Albedo of snow covered area on vegetation
1075	!! (unitless ratio)
1076	REAL(r_std), DIMENSION(kjpindex,nnobio) :: snowa_nobio !! Albedo of snow covered area on continental ice,
1077	!! lakes, etc. (unitless ratio)
1078	REAL(r_std), DIMENSION(kjpindex) :: fraction_veg !! Total vegetation fraction (unitless ratio)
1079	REAL(r_std), DIMENSION(kjpindex) :: agefunc_veg !! Age dependency of snow albedo on vegetation
1080	!! (unitless)
1081	REAL(r_std), DIMENSION(kjpindex,nnobio) :: agefunc_nobio !! Age dependency of snow albedo on ice,
1082	!! lakes, .. (unitless)
1083	REAL(r_std) :: alb_nobio !! Albedo of continental ice, lakes, etc.
1084	!!(unitless ratio)
1085	REAL(r_std),DIMENSION(kjpindex,nvm) :: zpcpalb !! Increase of snow albedo due to snowfall
1086	REAL(r_std),DIMENSION(kjpindex,nvm) :: zalbdry !! Dry snow albedo rate of change
1087	REAL(r_std),DIMENSION(kjpindex,nvm) :: zalbwet !! Wet snow albedo rate of change
1088	REAL(r_std),DIMENSION(kjpindex,nvm) :: zwholdmax !! Water holding capacity in snowpack
1089	REAL(r_std),DIMENSION(kjpindex,nvm) :: zfracliq !! Fraction of liquid water content in snowpack
1090
1091	!_ ================================================================================================================================
1092
1093	!! 1.Reset output variable (::albedo_snow ::albedo_glob)
1094	albedo_snow(:) = 0.0
1095	albedo_glob(:) = 0.0
1096	zfracliq(:,:) = 0.0
1097	zalbwet(:,:) = 0.0
1098	zalbdry(:,:) = 0.0
1099
1100	!! 2. Calculate snow albedos on both total vegetated and total nobio surfaces
1101
1102	! The snow albedo could be either prescribed (in condveg_init.f90) or
1103	! calculated following Chalita and Treut (1994).
1104	! Check if the precribed value fixed_snow_albedo exists
1105	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
1106	snowa_veg(:) = fixed_snow_albedo
1107	snowa_nobio(:,:) = fixed_snow_albedo
1108	ELSE ! calculated following Chalita and Treut (1994)
1109
1110	!! 2.1 Calculate age dependence
1111
1112	! On vegetated surfaces
1113	DO ji = 1, kjpindex
1114	agefunc_veg(ji) = EXP(-snow_age(ji)/tcst_snowa)
1115	ENDDO
1116
1117	! On non-vegtative surfaces
1118	DO jv = 1, nnobio ! Loop over # nobio types
1119	DO ji = 1, kjpindex
1120	agefunc_nobio(ji,jv) = EXP(-snow_nobio_age(ji,jv)/tcst_snowa)
1121	ENDDO
1122	ENDDO
1123
1124	!! 2.1 Calculate snow albedo
1125	IF (.NOT. control%ok_dofoco .AND. ok_explicitsnow) THEN
1126
1127	DO ji=1,kjpindex
1128	DO iv=1,nvm
1129	! The albedo scheme for the bare soil and short vegetation types
1130	IF (veget_max(ji,iv) .GT. 0.0 .AND. &
1131	(iv .EQ. 1 .OR. iv .GE. 10)) THEN
1132
1133	! 1. Albedo change (increase) due to snowfall:
1134	zpcpalb(ji,iv) = (precip_snow(ji)/snowcri_alb)*(xansmax-xansmin)
1135
1136	! 2. Dry snow albedo rate of change:
1137	zalbdry(ji,iv) = albedo_sveg(ji,iv) - xans_todry*(dtradia/86400.)
1138	IF (osfcmelt(ji)) THEN
1139	! 3. wet snow albedo rate of change:
1140	zalbwet(ji,iv) = (albedo_sveg(ji,iv) -xansmin)EXP(-xans_tdtradia/86400.) &
1141	+ xansmin
1142	! Calculate albedo rate based on liquid water
1143	! content of the first snow layer: During melt, if liquid
1144	! content is at the maximum, then snow albedo
1145	! decreases the fastest:
1146	zwholdmax(ji,iv) = snow3lhold_0d(snowrho(ji,1),snowdz(ji,1))
1147	zfracliq(ji,iv) = MIN(1.0,snowliq(ji,1)/MAX(xsnowdmin,zwholdmax(ji,iv)))
1148	END IF
1149
1150	! 4. Updated Albedo:
1151	albedo_sveg(ji,iv) = zfracliq(ji,iv)zalbwet(ji,iv) + (1.-zfracliq(ji,iv))&
1152	zalbdry(ji,iv) + zpcpalb(ji,iv)
1153	albedo_sveg(ji,iv) = MIN(xansmax,MAX(xansmin,albedo_sveg(ji,iv)))
1154	END IF
1155
1156	! The albedo scheme for tall trees
1157	IF (veget_max(ji,iv) .GT. 0.0 .AND. &
1158	(iv .GT. 1 .AND. iv .LT. 10)) THEN
1159	albedo_sveg(ji,iv) = snowa_aged(iv)+snowa_dec(iv)*agefunc_veg(ji)
1160	END IF
1161	END DO
1162	END DO
1163
1164	fraction_veg(:) = 1. - totfrac_nobio(:)
1165	snowa_veg(:) = 0.
1166
1167	DO jv = 1, nvm
1168	DO ji = 1, kjpindex
1169	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
1170	snowa_veg(ji) = snowa_veg(ji) + &
1171	veget_max(ji,jv)/fraction_veg(ji) * albedo_sveg(ji,jv)
1172	END IF
1173	END DO
1174	END DO
1175
1176	ELSE
1177
1178	! For vegetated surfaces
1179	fraction_veg(:) = un - totfrac_nobio(:)
1180	snowa_veg(:) = zero
1181	IF (control%ok_dgvm) THEN
1182	DO ji = 1, kjpindex
1183	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
1184	snowa_veg(ji) = snowa_veg(ji) + &
1185	tot_bare_soil(ji)/fraction_veg(ji) * ( snowa_aged(1)+snowa_dec(1)* &
1186	agefunc_veg(ji) )
1187	END IF
1188	END DO
1189
1190	DO jv = 2, nvm
1191	DO ji = 1, kjpindex
1192	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
1193	snowa_veg(ji) = snowa_veg(ji) + &
1194	veget(ji,jv)/fraction_veg(ji) * ( snowa_aged(jv)+snowa_dec(jv)* &
1195	agefunc_veg(ji) )
1196	ENDIF
1197	ENDDO
1198	ENDDO
1199	ELSE
1200	DO jv = 1, nvm
1201	DO ji = 1, kjpindex
1202	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
1203	snowa_veg(ji) = snowa_veg(ji) + &
1204	veget_max(ji,jv)/fraction_veg(ji) * ( snowa_aged(jv)+snowa_dec(jv)* &
1205	agefunc_veg(ji) )
1206	ENDIF
1207	ENDDO
1208	ENDDO
1209	ENDIF
1210	ENDIF
1211	!
1212	! snow albedo on other surfaces
1213	!
1214	DO jv = 1, nnobio
1215	DO ji = 1, kjpindex
1216	snowa_nobio(ji,jv) = ( snowa_aged(1) + snowa_dec(1) * agefunc_nobio(ji,jv) )
1217	ENDDO
1218	ENDDO
1219	ENDIF
1220
1221	!! 3. Calculate snow cover fraction for both total vegetated and total non-vegtative surfaces.\n
1222	frac_snow_veg(:) = MIN(MAX(snow(:),zero)/(MAX(snow(:),zero)+snowcri_alb*sn_dens/100.0),un)
1223	DO jv = 1, nnobio
1224	frac_snow_nobio(:,jv) = MIN(MAX(snow_nobio(:,jv),zero)/(MAX(snow_nobio(:,jv),zero)+ &
1225	snowcri_alb*sn_dens/100.0),un)
1226	ENDDO
1227
1228	!! 4. Update surface albedo
1229
1230	! Update surface albedo using the weighted sum of previous time step surface albedo,
1231	! total vegetated fraction, total nobio fraction, snow cover fraction (both vegetated and
1232	! non-vegetative surfaces), and snow albedo (both vegetated and non-vegetative surfaces).
1233	! Although both visible and near-infrared surface albedo are presented, their calculations
1234	! are the same.
1235	DO jb = 1, 2
1236	IF (.NOT. control%ok_dofoco .AND. ok_explicitsnow) THEN
1237	albedo_undersnow(:) = albedo_undersnow(:) + fraction_veg(:)*albedo(:,jb)
1238	ENDIF
1239
1240
1241	albedo(:,jb) = ( fraction_veg(:) ) * &
1242	( (un-frac_snow_veg(:)) * albedo(:,jb) + &
1243	( frac_snow_veg(:) ) * snowa_veg(:) )
1244	albedo_snow(:) = albedo_snow(:) + (fraction_veg(:)) * (frac_snow_veg(:)) * snowa_veg(:)
1245	!
1246	DO jv = 1, nnobio ! Loop over # nobio surfaces
1247	!
1248	IF ( jv .EQ. iice ) THEN
1249	alb_nobio = alb_ice(jb)
1250	ELSE
1251	WRITE(numout,*) 'jv=',jv
1252	WRITE(numout,*) 'DO NOT KNOW ALBEDO OF THIS SURFACE TYPE'
1253	CALL ipslerr_p(3,'condveg_snow','DO NOT KNOW ALBEDO OF THIS SURFACE TYPE','','')
1254	ENDIF
1255	!
1256	albedo(:,jb) = albedo(:,jb) + &
1257	( frac_nobio(:,jv) ) * &
1258	( (un-frac_snow_nobio(:,jv)) * alb_nobio + &
1259	( frac_snow_nobio(:,jv) ) * snowa_nobio(:,jv) )
1260	albedo_snow(:) = albedo_snow(:) + &
1261	( frac_nobio(:,jv) ) * ( frac_snow_nobio(:,jv) ) * &
1262	snowa_nobio(:,jv)
1263	albedo_glob(:) = albedo_glob(:) + albedo(:,jb)
1264
1265	IF (.NOT. control%ok_dofoco .AND. ok_explicitsnow) THEN
1266	albedo_undersnow(:)=albedo_undersnow(:)+frac_nobio(:,jv) * alb_nobio
1267	ENDIF
1268
1269	ENDDO
1270
1271	END DO
1272
1273	DO ji = 1, kjpindex
1274	albedo_snow(ji) = (albedo_snow(ji))/2.
1275	albedo_glob(ji) = (albedo_glob(ji))/2.
1276	IF (.NOT. control%ok_dofoco .AND. ok_explicitsnow) THEN
1277	albedo_undersnow(ji) = (albedo_undersnow(ji))/2.
1278	ENDIF
1279	ENDDO
1280
1281	IF (long_print) WRITE (numout,*) ' condveg_snow done '
1282
1283	END SUBROUTINE condveg_snow
1284
1285
1286	!! ==============================================================================================================================
1287	!! SUBROUTINE : condveg_soilalb
1288	!!
1289	!>\BRIEF This subroutine calculates the albedo of soil (without snow).
1290	!!
1291	!! DESCRIPTION This subroutine reads the soil colour maps in 1 x 1 deg resolution
1292	!! from the Henderson-Sellers & Wilson database. These values are interpolated to
1293	!! the model's resolution and transformed into
1294	!! dry and wet albedos.\n
1295	!!
1296	!! If the soil albedo is calculated without the dependence of soil moisture, the
1297	!! soil colour values are transformed into mean soil albedo values.\n
1298	!!
1299	!! The calculations follow the assumption that the grid of the data is regular and
1300	!! it covers the globe. The calculation for the model grid are based on the borders
1301	!! of the grid of the resolution.
1302	!!
1303	!! RECENT CHANGE(S): None
1304	!!
1305	!! MAIN OUTPUT VARIABLE(S): soilalb_dry for visible and near-infrared range,
1306	!! soilalb_wet for visible and near-infrared range,
1307	!! soilalb_moy for visible and near-infrared range
1308	!!
1309	!! REFERENCE(S) :
1310	!! -Wilson, M.F., and A. Henderson-Sellers, 1985: A global archive of land cover and
1311	!! soils data for use in general circulation climate models. J. Clim., 5, 119-143.
1312	!!
1313	!! FLOWCHART : None
1314	!! \n
1315	!_ ================================================================================================================================
1316
1317	SUBROUTINE condveg_soilalb(nbpt, lalo, neighbours, resolution, contfrac, soilalb_dry,soilalb_wet)
1318
1319	!! 0. Variable and parameter declaration
1320
1321	!! 0.1 Input variables
1322
1323	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
1324	!! interpolated (unitless)
1325	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
1326	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point
1327	!! (1=N, 2=E, 3=S, 4=W)
1328	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (m)
1329	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
1330	REAL(r_std), INTENT(inout) :: soilalb_dry(nbpt,2) !! Albedo of the dry bare soil (unitless)
1331	REAL(r_std), INTENT(inout) :: soilalb_wet(nbpt,2) !! Albedo of the wet dry soil (unitless)
1332
1333	!! 0.2 Output variables
1334
1335	!! 0.3 Modified variables
1336
1337	!! 0.4 Local variables
1338
1339	INTEGER(i_std) :: nbvmax !! nbvmax for interpolation (unitless). It is the
1340	!! dimension of the variables in which we store the list
1341	!! of points of the source grid which fit into one grid
1342	!! box of the target.
1343	CHARACTER(LEN=80) :: filename !! Filename of soil colour map
1344	INTEGER(i_std) :: iml, jml, lml, &
1345	&tml, fid, ib, ip, jp, fopt, ilf, lastjp, nbexp !! Indices
1346	REAL(r_std) :: lev(1), date, dt !! Help variables to read in file data
1347	INTEGER(i_std) :: itau(1) !! Help variables to read in file data
1348	REAL(r_std) :: sgn !! Help variable to compute average bare soil albedo
1349	REAL(r_std) :: coslat !! [DISPENSABLE]
1350	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel !! Help variable to read file data and allocate memory
1351	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon_rel !! Help variable to read file data and allocate memory
1352	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: soilcol !! Help variable to read file data and allocate memory
1353	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! Help variable to read file data and allocate memory
1354	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! Help variable to read file data and allocate memory
1355	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask !! Help variable to read file data and allocate memory
1356	CHARACTER(LEN=30) :: callsign !! Help variable to read file data and allocate memory
1357	LOGICAL :: ok_interpol = .FALSE. !! Optional return of aggregate_2d
1358	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
1359	!_ ================================================================================================================================
1360
1361	!! 1. Open file and allocate memory
1362
1363	! Open file with soil colours
1364
1365	!Config Key = SOILALB_FILE
1366	!Config Desc = Name of file from which the bare soil albedo
1367	!Config Def = soils_param.nc
1368	!Config If = NOT(IMPOSE_AZE)
1369	!Config Help = The name of the file to be opened to read the soil types from
1370	!Config which we derive then the bare soil albedos. This file is 1x1
1371	!Config deg and based on the soil colors defined by Wilson and Henderson-Seller.
1372	!Config Units = [FILE]
1373	!
1374	filename = 'soils_param.nc'
1375	CALL getin_p('SOILALB_FILE',filename)
1376
1377	! Read data from file
1378	IF (is_root_prc) CALL flininfo(filename,iml, jml, lml, tml, fid)
1379	CALL bcast(iml)
1380	CALL bcast(jml)
1381	CALL bcast(lml)
1382	CALL bcast(tml)
1383
1384	! Allocate memory for latitudes
1385	ALLOC_ERR=-1
1386	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
1387	IF (ALLOC_ERR/=0) THEN
1388	WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
1389	CALL ipslerr_p (3,'condveg.f90','condveg_soilalb','','')
1390	ENDIF
1391
1392	! Allcoate memory for longitude
1393	ALLOC_ERR=-1
1394	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
1395	IF (ALLOC_ERR/=0) THEN
1396	WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
1397	CALL ipslerr_p (3,'condveg.f90','condveg_soilalb','','')
1398	ENDIF
1399
1400	! Allocate memory for mask
1401	ALLOC_ERR=-1
1402	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
1403	IF (ALLOC_ERR/=0) THEN
1404	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
1405	CALL ipslerr_p (3,'condveg.f90','condveg_soilalb','','')
1406	ENDIF
1407
1408	! Allocate memory for soil data
1409	ALLOC_ERR=-1
1410	ALLOCATE(soilcol(iml,jml), STAT=ALLOC_ERR)
1411	IF (ALLOC_ERR/=0) THEN
1412	WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
1413	CALL ipslerr_p (3,'condveg.f90','condveg_soilalb','','')
1414	ENDIF
1415
1416	! Set values
1417	IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
1418	CALL bcast(lon_rel)
1419	CALL bcast(lat_rel)
1420	CALL bcast(lev)
1421	CALL bcast(itau)
1422	CALL bcast(date)
1423	CALL bcast(dt)
1424
1425	IF (is_root_prc) CALL flinget(fid, 'soilcolor', iml, jml, lml, tml, 1, 1, soilcol)
1426	CALL bcast(soilcol)
1427
1428	IF (is_root_prc) CALL flinclo(fid)
1429
1430	! Create mask with values of soil colour
1431	mask(:,:) = zero
1432	DO ip=1,iml
1433	DO jp=1,jml
1434	IF (soilcol(ip,jp) > min_sechiba) THEN
1435	mask(ip,jp) = un
1436	ENDIF
1437	ENDDO
1438	ENDDO
1439
1440	! Set nbvmax to 200 for interpolation
1441	! This number is the dimension of the variables in which we store
1442	! the list of points of the source grid which fit into one grid box of the target.
1443	nbvmax = 200
1444
1445	callsign = 'Soil color map'
1446
1447	! Start with interpolation
1448	DO WHILE ( .NOT. ok_interpol )
1449	WRITE(numout,*) "Projection arrays for ",callsign," : "
1450	WRITE(numout,*) "nbvmax = ",nbvmax
1451
1452	ALLOC_ERR=-1
1453	ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
1454	IF (ALLOC_ERR/=0) THEN
1455	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
1456	CALL ipslerr_p (3,'condveg.f90','condveg_soilalb','','')
1457	ENDIF
1458	sub_area(:,:)=zero
1459	ALLOC_ERR=-1
1460	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
1461	IF (ALLOC_ERR/=0) THEN
1462	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
1463	CALL ipslerr_p (3,'condveg.f90','condveg_soilalb','','')
1464	ENDIF
1465	sub_index(:,:,:)=0
1466
1467	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
1468	& iml, jml, lon_rel, lat_rel, mask, callsign, &
1469	& nbvmax, sub_index, sub_area, ok_interpol)
1470
1471	IF ( .NOT. ok_interpol ) THEN
1472	DEALLOCATE(sub_area)
1473	DEALLOCATE(sub_index)
1474	ENDIF
1475
1476	nbvmax = nbvmax * 2
1477	ENDDO
1478
1479	! Check how many points with soil information are found
1480	nbexp = 0
1481
1482	soilalb_dry(:,:) = zero
1483	soilalb_wet(:,:) = zero
1484	soilalb_moy(:,:) = zero
1485
1486	DO ib=1,nbpt ! Loop over domain size
1487
1488	fopt = COUNT(sub_area(ib,:) > zero)
1489
1490	!! 3. Compute the average bare soil albedo parameters
1491
1492	IF ( fopt .EQ. 0) THEN ! If no points were interpolated
1493	nbexp = nbexp + 1
1494	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1495	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1496	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1497	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1498	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
1499	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
1500	ELSE
1501	sgn = zero
1502
1503	DO ilf = 1,fopt ! If points were interpolated
1504
1505	ip = sub_index(ib,ilf,1)
1506	jp = sub_index(ib,ilf,2)
1507
1508	! Weighted albedo values by interpolation area
1509	IF ( NINT(soilcol(ip,jp)) .LE. classnb) THEN
1510	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1511	soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1512	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1513	soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1514	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1515	soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1516	sgn = sgn + sub_area(ib,ilf)
1517	ELSE
1518	WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
1519	CALL ipslerr_p (3,'condveg.f90','condveg_soilalb','','')
1520	ENDIF
1521
1522	ENDDO
1523
1524	! Normalize the surface
1525	IF ( sgn .LT. min_sechiba) THEN
1526	nbexp = nbexp + 1
1527	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1528	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1529	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1530	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1531	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
1532	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
1533	ELSE
1534	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis)/sgn
1535	soilalb_dry(ib,inir) = soilalb_dry(ib,inir)/sgn
1536	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis)/sgn
1537	soilalb_wet(ib,inir) = soilalb_wet(ib,inir)/sgn
1538	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis)/sgn
1539	soilalb_moy(ib,inir) = soilalb_moy(ib,inir)/sgn
1540	ENDIF
1541
1542	ENDIF
1543
1544	ENDDO
1545
1546	IF ( nbexp .GT. 0 ) THEN
1547	WRITE(numout,*) 'CONDVEG_soilalb : The interpolation of the bare soil albedo had ', nbexp
1548	WRITE(numout,*) 'CONDVEG_soilalb : points without data. This are either coastal points or'
1549	WRITE(numout,*) 'CONDVEG_soilalb : ice covered land.'
1550	WRITE(numout,*) 'CONDVEG_soilalb : The problem was solved by using the average of all soils'
1551	WRITE(numout,*) 'CONDVEG_soilalb : in dry and wet conditions'
1552	ENDIF
1553
1554	DEALLOCATE (lat_rel)
1555	DEALLOCATE (lon_rel)
1556	DEALLOCATE (mask)
1557	DEALLOCATE (sub_index)
1558	DEALLOCATE (sub_area)
1559	DEALLOCATE (soilcol)
1560
1561	RETURN
1562
1563	END SUBROUTINE condveg_soilalb
1564
1565
1566	!! ==============================================================================================================================
1567	!! SUBROUTINE : condveg_z0logz
1568	!!
1569	!>\BRIEF Computation of grid average of roughness height by averaging the
1570	!! logarithm of the roughness height of each grid box components fracbio and fracnobio.
1571	!!
1572	!! DESCRIPTION : Calculates mean roughness height
1573	!! over the grid cell. The mean roughness height is derived from the vegetation
1574	!! height which is scaled by the roughness parameter. The sum of the logarithm of the
1575	!! roughness times the fraction per grid cell gives the average roughness height per
1576	!! grid cell for the vegetative PFTs. The roughness height for the non-vegetative PFTs
1577	!! is calculated in a second step. \n
1578	!!
1579	!! To compute the fluxes,
1580	!! the difference between the height of the vegetation and the zero plane displacement height
1581	!! is needed and called roughheight .\n
1582	!!
1583	!! RECENT CHANGE(S): None
1584	!!
1585	!! MAIN OUTPUT VARIABLE(S): roughness height (z0), grid effective roughness height (roughheight)
1586	!!
1587	!! REFERENCE(S) :
1588	!!
1589	!! FLOWCHART : None
1590	!! \n
1591	!_ ================================================================================================================================
1592
1593	SUBROUTINE condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
1594	& z0, roughheight, z0_veg)
1595
1596	!! 0. Variable and parameter declaration
1597
1598	!! 0.1 Input variables
1599
1600	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1601	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1602	!! (m^2 m^{-2})
1603	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1604	!! (= ind*cn_ind) (m^2 m^{-2})
1605	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1606	!! i.e. continental ice, lakes, etc. (unitless)
1607	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1608	!! i.e. continental ice, lakes, etc. (unitless)
1609	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1610
1611	!! 0.2 Output variables
1612
1613	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0 !! Soil roughness height (m)
1614	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1615	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0_veg !! Roughness height of vegetated part (m)
1616
1617	!! 0.3 Modified variables
1618
1619	!! 0.4 Local variables
1620
1621	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1622	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1623	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1624	REAL(r_std), DIMENSION(kjpindex) :: d_veg !! PFT coverage of vegetative PFTs
1625	!! (= ind*cn_ind) (m^2 m^{-2})
1626	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1627	REAL(r_std) :: z0_nobio !! Roughness of non-vegetative fraction (m),
1628	!! i.e. continental ice, lakes, etc.
1629	!_ ================================================================================================================================
1630
1631	!! 1. Preliminary calculation
1632
1633	! Calculate the roughness (m) of bare soil, z0_bare
1634	! taken from constantes_veg.f90
1635	z0(:) = tot_bare_soil(:) * LOG(z0_bare)
1636
1637	! Define fraction of bare soil
1638	sumveg(:) = tot_bare_soil(:)
1639
1640	! Set average vegetation height to zero
1641	ave_height(:) = zero
1642
1643	!! 2. Calculate the mean roughness length
1644
1645	! Calculate the mean roughness height of
1646	! vegetative PFTs over the grid cell
1647	DO jv = 2, nvm !Loop over # vegetative PFTs
1648
1649	! In the case of forest, use parameter veget_max because
1650	! tree trunks influence the roughness even when there are no leaves
1651	IF ( is_tree(jv) ) THEN
1652	d_veg(:) = veget_max(:,jv)
1653	ELSE
1654
1655	! In the case of grass, use parameter veget because grasses
1656	! only influence the roughness during the growing season
1657	d_veg(:) = veget(:,jv)
1658	ENDIF
1659
1660	! Calculate the average roughness over the grid cell:
1661	! The roughness for vegetative PFTs is calculated by
1662	! the vegetation height per PFT multiplied by the roughness
1663	! parameter 'z0_over_height= 1/16'. If this scaled value is
1664	! lower than 0.01 than the value for the roughness length
1665	! of bare soil (0.01) is used. The sum of the logarithm of
1666	! the roughness times the fraction per grid cell gives the
1667	! logarithm of roughness length per grid cell for the vegetative
1668	! PFTs.
1669	z0(:) = z0(:) + d_veg(:) * &
1670	LOG( MAX(height(:,jv)*z0_over_height,z0_bare) )
1671	! Sum of bare soil and fraction vegetated fraction
1672	sumveg(:) = sumveg(:) + d_veg(:)
1673
1674	! Weighted height of vegetation with maximal cover fraction
1675	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1676
1677	ENDDO !Loop over # vegetative PFTs
1678
1679	!! 3. Calculate the mean roughness length of non-vegetative surfaces \n
1680
1681	! Search for pixels with vegetated part to normalise
1682	! roughness length
1683	WHERE ( sumveg(:) > zero ) z0(:) = z0(:) / sumveg(:)
1684
1685	! Calculate fraction of roughness for vegetated part
1686	z0(:) = (un - totfrac_nobio(:)) * z0(:)
1687	! Save roughness of vegetated part for calculation of snow fraction
1688	z0_veg(:) = z0(:)
1689
1690	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1691
1692	! Set rougness for ice
1693	IF ( jv .EQ. iice ) THEN
1694	z0_nobio = z0_ice
1695	ELSE
1696	WRITE(numout,*) 'jv=',jv
1697	CALL ipslerr_p (3,'condveg.f90',&
1698	'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1699	ENDIF
1700
1701	! Sum of vegetative roughness length and non-vegetative
1702	! roughness length
1703	z0(:) = z0(:) + frac_nobio(:,jv) * LOG(z0_nobio)
1704
1705
1706
1707	ENDDO ! loop over # of non-vegative surfaces
1708
1709	!! 4. Calculate the zero plane displacement height and effective roughness length
1710
1711	! Take the exponential of the roughness length
1712	z0(:) = EXP( z0(:) )
1713
1714	! Compute the zero plane displacement height which
1715	! is an equivalent height of the vegetation for the absorption of momentum
1716	zhdispl(:) = ave_height(:) * height_displacement
1717
1718	! Then we compute what we call the grid effective roughness height.
1719	! This is the height over which the roughness acts. It combines the
1720	! zero plane displacement height and the vegetation height. This
1721	! effective value is the difference between the height of the
1722	! vegetation and the zero plane displacement height.
1723	roughheight(:) = ave_height(:) - zhdispl(:)
1724
1725
1726	END SUBROUTINE condveg_z0logz
1727
1728
1729	!! ==============================================================================================================================
1730	!! SUBROUTINE : condveg_z0cdrag
1731	!!
1732	!>\BRIEF Computation of grid average of roughness length by calculating
1733	!! the drag coefficient.
1734	!!
1735	!! DESCRIPTION : This routine calculates the mean roughness height and mean
1736	!! effective roughness height over the grid cell. The mean roughness height (z0)
1737	!! is computed by averaging the drag coefficients \n
1738	!!
1739	!! \latexonly
1740	!! \input{z0cdrag1.tex}
1741	!! \endlatexonly
1742	!! \n
1743	!!
1744	!! where C is the drag coefficient at the height of the vegetation, kappa is the
1745	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
1746	!! The reference level for z needs to be high enough above the canopy to avoid
1747	!! singularities of the LOG. This height is set to minimum 10m above ground.
1748	!! The drag coefficient increases with roughness height to represent the greater
1749	!! turbulence generated by rougher surfaces.
1750	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
1751	!!
1752	!! The roughness height for the non-vegetative surfaces is calculated in a second step.
1753	!! In order to calculate the transfer coefficients the
1754	!! effective roughness height is calculated. This effective value is the difference
1755	!! between the height of the vegetation and the zero plane displacement height.\nn
1756	!!
1757	!! RECENT CHANGE(S): None
1758	!!
1759	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
1760	!!
1761	!! REFERENCE(S) : None
1762	!!
1763	!! FLOWCHART : None
1764	!! \n
1765	!_ ================================================================================================================================
1766
1767	SUBROUTINE condveg_z0cdrag (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, &
1768	& z0, roughheight,z0_veg)
1769
1770	!! 0. Variable and parameter declaration
1771
1772	!! 0.1 Input variables
1773
1774	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1775	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1776	!! (m^2 m^{-2})
1777	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1778	!! (= ind*cn_ind) (m^2 m^{-2})
1779	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1780	!! i.e. continental ice, lakes, etc. (unitless)
1781	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1782	!! i.e. continental ice, lakes, etc. (unitless)
1783	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
1784	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1785
1786	!! 0.2 Output variables
1787
1788	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0 !! Roughness height (m)
1789	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1790	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0_veg !! Roughness height of vegetated part (m)
1791
1792
1793	!! 0.3 Modified variables
1794
1795	!! 0.4 Local variables
1796
1797	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1798	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1799	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
1800	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1801	REAL(r_std), DIMENSION(kjpindex) :: d_veg !! PFT coverage of vegetative PFTs
1802	!! (= ind*cn_ind) (m^2 m^{-2})
1803	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1804	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
1805	!! i.e. continental ice, lakes, etc.
1806	!_ ================================================================================================================================
1807
1808	!! 1. Preliminary calculation
1809
1810	! Set maximal height of first layer
1811	ztmp(:) = MAX(10., zlev(:))
1812
1813	! Calculate roughness for non-vegetative surfaces
1814	! with the von Karman constant
1815	z0(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/z0_bare))**2
1816
1817	! Fraction of bare soil
1818	sumveg(:) = tot_bare_soil(:)
1819
1820	! Set average vegetation height to zero
1821	ave_height(:) = zero
1822
1823	!! 2. Calculate the mean roughness height
1824
1825	! Calculate the mean roughness height of
1826	! vegetative PFTs over the grid cell
1827	DO jv = 2, nvm
1828
1829	! In the case of forest, use parameter veget_max because
1830	! tree trunks influence the roughness even when there are no leaves
1831	IF ( is_tree(jv) ) THEN
1832	! In the case of grass, use parameter veget because grasses
1833	! only influence the roughness during the growing season
1834	d_veg(:) = veget_max(:,jv)
1835	ELSE
1836	! grasses only have an influence if they are really there!
1837	d_veg(:) = veget(:,jv)
1838	ENDIF
1839
1840	! Calculate the average roughness over the grid cell:
1841	! The unitless drag coefficient is per vegetative PFT
1842	! calculated by use of the von Karman constant, the height
1843	! of the first layer and the roughness. The roughness
1844	! is calculated as the vegetation height per PFT
1845	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
1846	! If this scaled value is lower than 0.01 then the value for
1847	! the roughness of bare soil (0.01) is used.
1848	! The sum over all PFTs gives the average roughness
1849	! per grid cell for the vegetative PFTs.
1850	z0(:) = z0(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)z0_over_height,z0_bare)))*2
1851
1852	! Sum of bare soil and fraction vegetated fraction
1853	sumveg(:) = sumveg(:) + d_veg(:)
1854
1855	! Weigh height of vegetation with maximal cover fraction
1856	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1857
1858	ENDDO
1859
1860	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1861
1862	! Search for pixels with vegetated part to normalise
1863	! roughness height
1864	WHERE ( sumveg(:) .GT. zero ) z0(:) = z0(:) / sumveg(:)
1865
1866	! Calculate fraction of roughness for vegetated part
1867	z0(:) = (un - totfrac_nobio(:)) * z0(:)
1868	! Save roughness of vegetated part for calculation of snow fraction
1869	z0_veg(:) = z0(:)
1870
1871
1872	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1873
1874	! Set rougness for ice
1875	IF ( jv .EQ. iice ) THEN
1876	z0_nobio = z0_ice
1877
1878
1879	ELSE
1880	WRITE(numout,*) 'jv=',jv
1881	CALL ipslerr_p (3,'condveg.f90',&
1882	'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1883	ENDIF
1884
1885	! Sum of vegetative roughness length and non-vegetative
1886	! roughness length
1887	z0(:) = z0(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1888
1889	ENDDO ! Loop over # of non-vegative surfaces
1890
1891
1892	!! 4. Calculate the zero plane displacement height and effective roughness length
1893
1894	! Take the exponential of the roughness
1895	z0(:) = ztmp(:) / EXP(ct_karman/SQRT(z0(:)))
1896
1897	! Compute the zero plane displacement height which
1898	! is an equivalent height for the absorption of momentum
1899	zhdispl(:) = ave_height(:) * height_displacement
1900
1901	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1902	! This is the height over which the roughness acts. It combines the
1903	! zero plane displacement height and the vegetation height.
1904	roughheight(:) = ave_height(:) - zhdispl(:)
1905
1906
1907	END SUBROUTINE condveg_z0cdrag
1908
1909
1910	END MODULE condveg

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

Context Navigation

source: branches/publications/ORCHIDEE_CAN_r3069/src_sechiba/condveg.f90 @ 7108

Download in other formats: