Context Navigation

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_sechiba/condveg.f90 @ 8398

Last change on this file since 8398 was 6957, checked in by fabienne.maignan, 4 years ago
Update before parameters optimization
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 89.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. :: rough_dyn to choose between two methods to calculate
14	!! the roughness height. If set to false: the roughness height is calculated by the old formulation
15	!! which does not distinguish between z0m and z0h and which does not vary with LAI
16	!! If set to true: the grid average is calculated by the formulation proposed by Su et al. (2001)
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. This is useful if
20	!! 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 and the effective roughness height
29	!! defined as the difference between the mean height of the vegetation and the displacement height (zero wind
30	!! level). Over bare soils, the zero wind level is equal to the soil roughness. Over vegetation, the zero wind level
31	!! is increased by the displacement height
32	!! which depends on the height of the vegetation. For a grid point composed of different types of vegetation,
33	!! an effective surface roughness has to be calculated
34	!!
35	!! RECENT CHANGE(S): Added option rough_dyn and subroutine condveg_z0cdrag_dyn. Removed subroutine condveg_z0logz. June 2016.
36	!!
37	!! REFERENCES(S) : None
38	!!
39
40	!! SVN :
41	!! $HeadURL$
42	!! $Date$
43	!! $Revision$
44	!! \n
45	!_ ================================================================================================================================
46
47	MODULE condveg
48
49	USE ioipsl
50	USE xios_orchidee
51	USE constantes
52	USE constantes_soil
53	USE pft_parameters
54	USE qsat_moisture
55	USE interpol_help
56	USE mod_orchidee_para
57	USE ioipsl_para
58	USE sechiba_io_p
59	USE grid
60
61	IMPLICIT NONE
62
63	PRIVATE
64	PUBLIC :: condveg_main, condveg_initialize, condveg_finalize, condveg_clear
65
66	!
67	! Variables used inside condveg module
68	!
69	LOGICAL, SAVE :: l_first_condveg=.TRUE. !! To keep first call's trace
70	!$OMP THREADPRIVATE(l_first_condveg)
71	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_dry(:,:) !! Albedo values for the dry bare soil (unitless)
72	!$OMP THREADPRIVATE(soilalb_dry)
73	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_wet(:,:) !! Albedo values for the wet bare soil (unitless)
74	!$OMP THREADPRIVATE(soilalb_wet)
75	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_moy(:,:) !! Albedo values for the mean bare soil (unitless)
76	!$OMP THREADPRIVATE(soilalb_moy)
77	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_bg(:,:) !! Albedo values for the background bare soil (unitless)
78	!$OMP THREADPRIVATE(soilalb_bg)
79	INTEGER, SAVE :: printlev_loc !! Output debug level
80	!$OMP THREADPRIVATE(printlev_loc)
81
82	CONTAINS
83
84	!! =============================================================================================================================
85	!! SUBROUTINE : condveg_initialize
86	!!
87	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
88	!!
89	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
90	!! condveg_snow is called to initialize corresponding variables.
91	!!
92	!! RECENT CHANGE(S) : None
93	!!
94	!! MAIN OUTPUT VARIABLE(S)
95	!!
96	!! REFERENCE(S) : None
97	!!
98	!! FLOWCHART : None
99	!! \n
100	!_ ==============================================================================================================================
101	SUBROUTINE condveg_initialize (kjit, kjpindex, index, rest_id, &
102	lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
103	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
104	drysoil_frac, height, snowdz, snowrho, tot_bare_soil, &
105	temp_air, pb, u, v, lai, &
106	emis, albedo, z0m, z0h, roughheight, &
107	frac_snow_veg,frac_snow_nobio)
108
109	!! 0. Variable and parameter declaration
110	!! 0.1 Input variables
111	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
112	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
113	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
114	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
115	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
116	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours!! neighoring grid points if land
117	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
118	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
119	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
120	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
121	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
122	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
123	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
124	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
125	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
126	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
127	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
128	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
129	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
130	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowdz !! Snow depth at each snow layer
131	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowrho !! Snow density at each snow layer
132	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
133	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
134	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
135	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u !! Horizontal wind speed, u direction
136	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: v !! Horizontal wind speed, v direction
137	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
138
139	!! 0.2 Output variables
140	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
141	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
142	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m !! Roughness for momentum (m)
143	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h !! Roughness for heat (m)
144	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
145	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: frac_snow_veg !! Snow cover fraction on vegeted area
146	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Snow cover fraction on non-vegeted area
147
148	!! 0.4 Local variables
149	INTEGER :: ier
150	REAL(r_std), DIMENSION(kjpindex,2) :: albedo_snow !! Snow albedo for visible and near-infrared range(unitless)
151	REAL(r_std), DIMENSION(kjpindex,2) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
152	!! range (unitless)
153	REAL(r_std), DIMENSION(kjpindex,2) :: alb_veget !! Mean vegetation albedo for visible and near-infrared
154	! !! range (unitless)
155	!_ ================================================================================================================================
156
157	IF (.NOT. l_first_condveg) CALL ipslerr_p(3,'condveg_initialize','Error: initialization already done','','')
158	l_first_condveg=.FALSE.
159
160	!! Initialize local printlev
161	printlev_loc=get_printlev('condveg')
162
163	IF (printlev>=3) WRITE (numout,*) 'Start condveg_initialize'
164
165	!! 1. Allocate module variables and read from restart or initialize
166
167	IF (alb_bg_modis) THEN
168	! Allocate background soil albedo
169	ALLOCATE (soilalb_bg(kjpindex,2),stat=ier)
170	IF (ier /= 0) CALL ipslerr_p(3,'condveg_initialize','Pb in allocation for soilalb_bg','','')
171
172	! Read background albedo from restart file
173	CALL ioconf_setatt_p('UNITS', '-')
174	CALL ioconf_setatt_p('LONG_NAME','Background soil albedo for visible and near-infrared range')
175	CALL restget_p (rest_id, 'soilalbedo_bg', nbp_glo, 2, 1, kjit, .TRUE., soilalb_bg, "gather", nbp_glo, index_g)
176
177	! Initialize by interpolating from file if the variable was not in restart file
178	IF ( ALL(soilalb_bg(:,:) == val_exp) ) THEN
179	CALL condveg_background_soilalb(kjpindex, lalo, neighbours, resolution, contfrac)
180	END IF
181	CALL xios_orchidee_send_field("soilalb_bg",soilalb_bg)
182
183	ELSE
184	! Allocate
185	! Dry soil albedo
186	ALLOCATE (soilalb_dry(kjpindex,2),stat=ier)
187	IF (ier /= 0) CALL ipslerr_p(3,'condveg_initialize','Pb in allocation for soilalb_dry','','')
188
189	! Wet soil albedo
190	ALLOCATE (soilalb_wet(kjpindex,2),stat=ier)
191	IF (ier /= 0) CALL ipslerr_p(3,'condveg_initialize','Pb in allocation for soilalb_wet','','')
192
193	! Mean soil albedo
194	ALLOCATE (soilalb_moy(kjpindex,2),stat=ier)
195	IF (ier /= 0) CALL ipslerr_p(3,'condveg_initialize','Pb in allocation for soilalb_moy','','')
196
197	! Read variables from restart file
198	! dry soil albedo
199	CALL ioconf_setatt_p('UNITS', '-')
200	CALL ioconf_setatt_p('LONG_NAME','Dry bare soil albedo')
201	CALL restget_p (rest_id,'soilalbedo_dry' , nbp_glo, 2, 1, kjit, .TRUE., soilalb_dry, "gather", nbp_glo, index_g)
202
203	! wet soil albedo
204	CALL ioconf_setatt_p('UNITS', '-')
205	CALL ioconf_setatt_p('LONG_NAME','Wet bare soil albedo')
206	CALL restget_p (rest_id, 'soilalbedo_wet', nbp_glo, 2, 1, kjit, .TRUE., soilalb_wet, "gather", nbp_glo, index_g)
207
208	! mean soil aledo
209	CALL ioconf_setatt_p('UNITS', '-')
210	CALL ioconf_setatt_p('LONG_NAME','Mean bare soil albedo')
211	CALL restget_p (rest_id, 'soilalbedo_moy', nbp_glo, 2, 1, kjit, .TRUE., soilalb_moy, "gather", nbp_glo, index_g)
212
213
214	! Initialize the variables if not found in restart file
215	IF ( ALL(soilalb_wet(:,:) == val_exp) .OR. &
216	ALL(soilalb_dry(:,:) == val_exp) .OR. &
217	ALL(soilalb_moy(:,:) == val_exp)) THEN
218	! One or more of the variables were not in the restart file.
219	! Call routine condveg_soilalb to calculate them.
220	CALL condveg_soilalb(kjpindex, lalo, neighbours, resolution, contfrac)
221	WRITE(numout,*) '---> val_exp ', val_exp
222	WRITE(numout,*) '---> ALBEDO_wet VIS:', MINVAL(soilalb_wet(:,ivis)), MAXVAL(soilalb_wet(:,ivis))
223	WRITE(numout,*) '---> ALBEDO_wet NIR:', MINVAL(soilalb_wet(:,inir)), MAXVAL(soilalb_wet(:,inir))
224	WRITE(numout,*) '---> ALBEDO_dry VIS:', MINVAL(soilalb_dry(:,ivis)), MAXVAL(soilalb_dry(:,ivis))
225	WRITE(numout,*) '---> ALBEDO_dry NIR:', MINVAL(soilalb_dry(:,inir)), MAXVAL(soilalb_dry(:,inir))
226	WRITE(numout,*) '---> ALBEDO_moy VIS:', MINVAL(soilalb_moy(:,ivis)), MAXVAL(soilalb_moy(:,ivis))
227	WRITE(numout,*) '---> ALBEDO_moy NIR:', MINVAL(soilalb_moy(:,inir)), MAXVAL(soilalb_moy(:,inir))
228	ENDIF
229	END IF
230
231	! z0m
232	CALL ioconf_setatt_p('UNITS', '-')
233	CALL ioconf_setatt_p('LONG_NAME','Roughness for momentum')
234	CALL restget_p (rest_id, 'z0m', nbp_glo, 1, 1, kjit, .TRUE., z0m, "gather", nbp_glo, index_g)
235
236	! z0h
237	CALL ioconf_setatt_p('UNITS', '-')
238	CALL ioconf_setatt_p('LONG_NAME','Roughness for heat')
239	CALL restget_p (rest_id, 'z0h', nbp_glo, 1, 1, kjit, .TRUE., z0h, "gather", nbp_glo, index_g)
240
241	! roughness height
242	CALL ioconf_setatt_p('UNITS', '-')
243	CALL ioconf_setatt_p('LONG_NAME','Roughness height')
244	CALL restget_p (rest_id, 'roughheight', nbp_glo, 1, 1, kjit, .TRUE., roughheight, "gather", nbp_glo, index_g)
245
246
247	!! Initialize emissivity
248	IF ( impaze ) THEN
249	! Use parameter CONDVEG_EMIS from run.def
250	emis(:) = emis_scal
251	ELSE
252	! Set emissivity to 1.
253	emis_scal = un
254	emis(:) = emis_scal
255	ENDIF
256
257
258	!! 3. Calculate the fraction of snow on vegetation and nobio
259	CALL condveg_frac_snow(kjpindex, snow, snow_nobio, snowrho, snowdz, &
260	frac_snow_veg, frac_snow_nobio, totfrac_nobio)
261
262	!! 4. Calculate roughness height if it was not found in the restart file
263	IF ( ALL(z0m(:) == val_exp) .OR. ALL(z0h(:) == val_exp) .OR. ALL(roughheight(:) == val_exp)) THEN
264	!! Calculate roughness height
265	! Chooses between two methods to calculate the grid average of the roughness.
266	! If impaze set to true: The grid average is calculated by averaging the drag coefficients over PFT.
267	! If impaze set to false: The grid average is calculated by averaging the logarithm of the roughness length per PFT.
268	IF ( impaze ) THEN
269	! Use parameter CONDVEG_Z0 and ROUGHHEIGHT from run.def
270	z0m(:) = z0_scal
271	z0h(:) = z0_scal
272	roughheight(:) = roughheight_scal
273	ELSE
274	! Caluculate roughness height
275	IF( rough_dyn ) THEN
276	CALL condveg_z0cdrag_dyn(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
277	& height, temp_air, pb, u, v, lai, frac_snow_veg, z0m, z0h, roughheight)
278	ELSE
279	CALL condveg_z0cdrag(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
280	height, tot_bare_soil, z0m, z0h, roughheight)
281	ENDIF
282	END IF
283	END IF
284
285	!! 5. Calculate albedo
286	CALL condveg_albedo (kjpindex, veget, veget_max, drysoil_frac, frac_nobio, &
287	totfrac_nobio, snow, snow_age, snow_nobio, &
288	snow_nobio_age, snowdz, snowrho, &
289	tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
290	albedo, albedo_snow, alb_bare, alb_veget)
291
292	IF (printlev>=3) WRITE (numout,*) 'condveg_initialize done '
293
294	END SUBROUTINE condveg_initialize
295
296
297
298	!! ==============================================================================================================================
299	!! SUBROUTINE : condveg_main
300	!!
301	!>\BRIEF Calls the subroutines update the variables for current time step
302	!!
303	!!
304	!! MAIN OUTPUT VARIABLE(S): emis (emissivity), albedo (albedo of
305	!! vegetative PFTs in visible and near-infrared range), z0 (surface roughness height),
306	!! roughheight (grid effective roughness height), soil type (fraction of soil types)
307	!!
308	!!
309	!! REFERENCE(S) : None
310	!!
311	!! FLOWCHART : None
312	!!
313	!! REVISION(S) : None
314	!!
315	!_ ================================================================================================================================
316
317	SUBROUTINE condveg_main (kjit, kjpindex, index, rest_id, hist_id, hist2_id, &
318	lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
319	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
320	drysoil_frac, height, snowdz, snowrho, tot_bare_soil, &
321	temp_air, pb, u, v, lai, &
322	emis, albedo, z0m, z0h, roughheight, &
323	frac_snow_veg, frac_snow_nobio )
324
325	!! 0. Variable and parameter declaration
326
327	!! 0.1 Input variables
328
329	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
330	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
331	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
332	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier
333	INTEGER(i_std), OPTIONAL, INTENT (in) :: hist2_id !! _History_ file 2 identifier
334	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
335	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
336	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours!! neighoring grid points if land
337	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
338	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
339	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
340	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
341	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
342	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
343	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
344	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
345	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
346	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
347	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
348	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
349	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
350	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowdz !! Snow depth at each snow layer
351	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowrho !! Snow density at each snow layer
352	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
353	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
354	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
355	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u !! Horizontal wind speed, u direction
356	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: v !! Horizontal wind speed, v direction
357	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
358
359	!! 0.2 Output variables
360
361	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
362	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
363	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m !! Roughness for momentum (m)
364	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h !! Roughness for heat (m)
365	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
366	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: frac_snow_veg !! Snow cover fraction on vegeted area
367	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Snow cover fraction on non-vegeted area
368
369	!! 0.3 Modified variables
370
371	!! 0.4 Local variables
372	REAL(r_std), DIMENSION(kjpindex,2) :: albedo_snow !! Snow albedo (unitless ratio)
373	REAL(r_std), DIMENSION(kjpindex,2) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
374	!! range (unitless)
375	REAL(r_std), DIMENSION(kjpindex,2) :: alb_veget !! Mean vegetation albedo for visible and near-infrared
376	!! range (unitless)
377	INTEGER(i_std) :: ji
378	!_ ================================================================================================================================
379
380	!! 1. Calculate the fraction of snow on vegetation and nobio
381	CALL condveg_frac_snow(kjpindex, snow, snow_nobio, snowrho, snowdz, &
382	frac_snow_veg, frac_snow_nobio, totfrac_nobio)
383
384	!! 2. Calculate emissivity
385	emis(:) = emis_scal
386
387	!! 3. Calculate roughness height
388	! If TRUE read in prescribed values for roughness height
389	IF ( impaze ) THEN
390
391	DO ji = 1, kjpindex
392	z0m(ji) = z0_scal
393	z0h(ji) = z0_scal
394	roughheight(ji) = roughheight_scal
395	ENDDO
396
397	ELSE
398	! Calculate roughness height
399	IF ( rough_dyn ) THEN
400	CALL condveg_z0cdrag_dyn (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, height, &
401	& temp_air, pb, u, v, lai, frac_snow_veg, z0m, z0h, roughheight)
402	ELSE
403	CALL condveg_z0cdrag (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
404	height, tot_bare_soil, z0m, z0h, roughheight)
405	ENDIF
406
407	ENDIF
408
409
410	!! 4. Calculate albedo
411	CALL condveg_albedo (kjpindex, veget, veget_max, drysoil_frac, frac_nobio, &
412	totfrac_nobio, snow, snow_age, snow_nobio, &
413	snow_nobio_age, snowdz, snowrho, &
414	tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
415	albedo, albedo_snow, alb_bare, alb_veget)
416
417
418
419	!! 5. Output diagnostics
420	IF (.NOT. impaze) THEN
421	CALL xios_orchidee_send_field("soilalb_vis",alb_bare(:,1))
422	CALL xios_orchidee_send_field("soilalb_nir",alb_bare(:,2))
423	CALL xios_orchidee_send_field("vegalb_vis",alb_veget(:,1))
424	CALL xios_orchidee_send_field("vegalb_nir",alb_veget(:,2))
425	END IF
426	CALL xios_orchidee_send_field("albedo_vis",albedo(:,1))
427	CALL xios_orchidee_send_field("albedo_nir",albedo(:,2))
428	CALL xios_orchidee_send_field("albedo_snow",(albedo_snow(:,1) + albedo_snow(:,2))/2)
429
430	IF ( almaoutput ) THEN
431	CALL histwrite_p(hist_id, 'Albedo', kjit, (albedo(:,1) + albedo(:,2))/2, kjpindex, index)
432	CALL histwrite_p(hist_id, 'SAlbedo', kjit, (albedo_snow(:,1) + albedo_snow(:,2))/2, kjpindex, index)
433	IF ( hist2_id > 0 ) THEN
434	CALL histwrite_p(hist2_id, 'Albedo', kjit, (albedo(:,1) + albedo(:,2))/2, kjpindex, index)
435	CALL histwrite_p(hist2_id, 'SAlbedo', kjit, (albedo_snow(:,1) + albedo_snow(:,2))/2, kjpindex, index)
436	ENDIF
437	ELSE
438	IF (.NOT. impaze) THEN
439	CALL histwrite_p(hist_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
440	CALL histwrite_p(hist_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
441	CALL histwrite_p(hist_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
442	CALL histwrite_p(hist_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
443	IF ( hist2_id > 0 ) THEN
444	CALL histwrite_p(hist2_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
445	CALL histwrite_p(hist2_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
446	CALL histwrite_p(hist2_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
447	CALL histwrite_p(hist2_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
448	ENDIF
449	END IF
450	ENDIF
451
452	IF (printlev>=3) WRITE (numout,*)' condveg_main done '
453
454	END SUBROUTINE condveg_main
455
456	!! =============================================================================================================================
457	!! SUBROUTINE : condveg_finalize
458	!!
459	!>\BRIEF Write to restart file
460	!!
461	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in condveg
462	!! to restart file
463	!!
464	!! RECENT CHANGE(S) : None
465	!!
466	!! REFERENCE(S) : None
467	!!
468	!! FLOWCHART : None
469	!! \n
470	!_ ==============================================================================================================================
471	SUBROUTINE condveg_finalize (kjit, kjpindex, rest_id, z0m, z0h, roughheight)
472
473	!! 0. Variable and parameter declaration
474	!! 0.1 Input variables
475	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
476	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
477	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
478	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: z0m !! Roughness for momentum
479	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: z0h !! Roughness for heat
480	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: roughheight !! Grid effective roughness height (m)
481
482	!_ ================================================================================================================================
483
484	CALL restput_p (rest_id, 'z0m', nbp_glo, 1, 1, kjit, z0m, 'scatter', nbp_glo, index_g)
485	CALL restput_p (rest_id, 'z0h', nbp_glo, 1, 1, kjit, z0h, 'scatter', nbp_glo, index_g)
486	CALL restput_p (rest_id, 'roughheight', nbp_glo, 1, 1, kjit, roughheight, 'scatter', nbp_glo, index_g)
487
488	IF ( alb_bg_modis ) THEN
489	CALL restput_p (rest_id, 'soilalbedo_bg', nbp_glo, 2, 1, kjit, soilalb_bg, 'scatter', nbp_glo, index_g)
490	ELSE
491	CALL restput_p (rest_id, 'soilalbedo_dry', nbp_glo, 2, 1, kjit, soilalb_dry, 'scatter', nbp_glo, index_g)
492	CALL restput_p (rest_id, 'soilalbedo_wet', nbp_glo, 2, 1, kjit, soilalb_wet, 'scatter', nbp_glo, index_g)
493	CALL restput_p (rest_id, 'soilalbedo_moy', nbp_glo, 2, 1, kjit, soilalb_moy, 'scatter', nbp_glo, index_g)
494	END IF
495	END SUBROUTINE condveg_finalize
496
497	!! ==============================================================================================================================
498	!! SUBROUTINE : condveg_clear
499	!!
500	!>\BRIEF Deallocate albedo variables
501	!!
502	!! DESCRIPTION : None
503	!!
504	!! RECENT CHANGE(S): None
505	!!
506	!! MAIN OUTPUT VARIABLE(S): None
507	!!
508	!! REFERENCES : None
509	!!
510	!! FLOWCHART : None
511	!! \n
512	!_ ================================================================================================================================
513
514	SUBROUTINE condveg_clear ()
515
516	l_first_condveg=.TRUE.
517
518	! Dry soil albedo
519	IF (ALLOCATED (soilalb_dry)) DEALLOCATE (soilalb_dry)
520	! Wet soil albedo
521	IF (ALLOCATED(soilalb_wet)) DEALLOCATE (soilalb_wet)
522	! Mean soil albedo
523	IF (ALLOCATED(soilalb_moy)) DEALLOCATE (soilalb_moy)
524	! BG soil albedo
525	IF (ALLOCATED(soilalb_bg)) DEALLOCATE (soilalb_bg)
526
527	END SUBROUTINE condveg_clear
528
529	!! ==============================================================================================================================\n
530	!! SUBROUTINE : condveg_albedo
531	!!
532	!>\BRIEF Calculate albedo
533	!!
534	!! DESCRIPTION : The albedo is calculated for both the visible and near-infrared
535	!! domain. First the mean albedo of the bare soil is calculated. Two options exist:
536	!! either the soil albedo depends on soil wetness (drysoil_frac variable), or the soil albedo
537	!! is set to a mean soil albedo value.
538	!! The snow albedo scheme presented below belongs to prognostic albedo
539	!! category, i.e. the snow albedo value at a time step depends on the snow albedo value
540	!! at the previous time step.
541	!!
542	!! First, the following formula (described in Chalita and Treut 1994) is used to describe
543	!! the change in snow albedo with snow age on each PFT and each non-vegetative surfaces,
544	!! i.e. continental ice, lakes, etc.: \n
545	!! \latexonly
546	!! \input{SnowAlbedo.tex}
547	!! \endlatexonly
548	!! \n
549	!! Where snowAge is snow age, tcstSnowa is a critical aging time (tcstSnowa=5 days)
550	!! snowaIni and snowaIni+snowaDec corresponds to albedos measured for aged and
551	!! fresh snow respectively, and their values for each PFT and each non-vegetative surfaces
552	!! is precribed in in constantes_veg.f90.\n
553	!! In order to estimate gridbox snow albedo, snow albedo values for each PFT and
554	!! each non-vegetative surfaces with a grid box are weightedly summed up by their
555	!! respective fractions.\n
556	!! Secondly, the snow cover fraction is computed as:
557	!! \latexonly
558	!! \input{SnowFraction.tex}
559	!! \endlatexonly
560	!! \n
561	!! Where fracSnow is the fraction of snow on total vegetative or total non-vegetative
562	!! surfaces, snow is snow mass (kg/m^2) on total vegetated or total nobio surfaces.\n
563	!! Finally, the surface albedo is then updated as the weighted sum of fracSnow, total
564	!! vegetated fraction, total nobio fraction, gridbox snow albedo, and previous
565	!! time step surface albedo.
566	!!
567	!! RECENT CHANGE(S): These calculations were previously done in condveg_albcalc and condveg_snow
568	!!
569	!! MAIN OUTPUT VARIABLE(S): :: albedo; surface albedo. :: albedo_snow; snow
570	!! albedo
571	!!
572	!! REFERENCE(S) :
573	!! Chalita, S. and H Le Treut (1994), The albedo of temperate and boreal forest and
574	!! the Northern Hemisphere climate: a sensitivity experiment using the LMD GCM,
575	!! Climate Dynamics, 10 231-240.
576	!!
577	!! FLOWCHART : None
578	!! \n
579	!_ ================================================================================================================================
580
581	SUBROUTINE condveg_albedo (kjpindex, veget, veget_max, drysoil_frac, frac_nobio, &
582	totfrac_nobio, snow, snow_age, snow_nobio, &
583	snow_nobio_age, snowdz, snowrho, &
584	tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
585	albedo, albedo_snow, alb_bare, alb_veget)
586
587
588	!! 0. Variable and parameter declarations
589
590	!! 0.1 Input variables
591
592	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
593	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
594	!! (m^2 m^{-2})
595	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max
596	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
597	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
598	!! continental ice, lakes, etc. (unitless)
599	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces, i.e.
600	!! continental ice, lakes, etc. (unitless)
601	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
602	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
603	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age (days)
604	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
605	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer
606	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density at each snow layer
607	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
608	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: frac_snow_veg !! Fraction of snow on vegetation (unitless ratio)
609	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc. (unitless ratio)
610
611	!! 0.2 Output variables
612	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo (unitless ratio)
613	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
614	REAL(r_std), DIMENSION(kjpindex,2), INTENT(out) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
615	!! range (unitless). Only calculated for .NOT. impaze
616	REAL(r_std), DIMENSION(kjpindex,2), INTENT(out) :: alb_veget !! Mean vegetation albedo for visible and near-infrared
617	!! range (unitless). Only calculated for .NOT. impaze
618
619	!! 0.3 Local variables
620	INTEGER(i_std) :: ji, jv, jb,ks !! indices (unitless)
621	REAL(r_std), DIMENSION(kjpindex,2) :: snowa_veg !! Albedo of snow covered area on vegetation
622	!! (unitless ratio)
623	REAL(r_std), DIMENSION(kjpindex,nnobio,2) :: snowa_nobio !! Albedo of snow covered area on continental ice,
624	!! lakes, etc. (unitless ratio)
625	REAL(r_std), DIMENSION(kjpindex) :: fraction_veg !! Total vegetation fraction (unitless ratio)
626	REAL(r_std), DIMENSION(kjpindex) :: agefunc_veg !! Age dependency of snow albedo on vegetation
627	!! (unitless)
628	REAL(r_std), DIMENSION(kjpindex,nnobio) :: agefunc_nobio !! Age dependency of snow albedo on ice,
629	!! lakes, .. (unitless)
630	REAL(r_std) :: alb_nobio !! Albedo of continental ice, lakes, etc.
631	!!(unitless ratio)
632	REAL(r_std),DIMENSION (nvm,2) :: alb_leaf_tmp !! Variables for albedo values for all PFTs and
633	REAL(r_std),DIMENSION (nvm,2) :: snowa_aged_tmp !! spectral domains (unitless)
634	REAL(r_std),DIMENSION (nvm,2) :: snowa_dec_tmp
635	!_ ================================================================================================================================
636
637
638
639	!! 1. Preliminary calculation without considering snow
640	IF ( impaze ) THEN
641	!! No caluculation, set default value
642	albedo(:,ivis) = albedo_scal(ivis)
643	albedo(:,inir) = albedo_scal(inir)
644
645	ELSE
646	!! Preliminary calculation without considering snow (previously done in condveg_albcalc)
647	! Assign values of leaf and snow albedo for visible and near-infrared range
648	! to local variable (constantes_veg.f90)
649	alb_leaf_tmp(:,ivis) = alb_leaf_vis(:)
650	alb_leaf_tmp(:,inir) = alb_leaf_nir(:)
651
652	!!$!!-SC- AlbÃ©do moyen
653
654	snowa_aged_tmp(:,ivis) = snowa_aged_vis(:)
655	snowa_aged_tmp(:,inir) = snowa_aged_nir(:)
656	snowa_dec_tmp(:,ivis) = snowa_dec_vis(:)
657	snowa_dec_tmp(:,inir) = snowa_dec_nir(:)
658
659	!!$!!-SC ! Attention snowa_aged_vis/ snowa_dec_vis est la moyenne des valeurs vis + nir (Ã modifier plus tard)
660	!!$ snowa_aged_tmp(:,ivis) = snowa_aged_vis(:)
661	!!$ snowa_dec_tmp(:,ivis) = snowa_dec_vis(:)
662	!!$
663	!!$ snowa_aged_tmp(:,inir) = snowa_aged_tmp(:,ivis)
664	!!$ snowa_dec_tmp(:,inir) = snowa_dec_tmp(:,ivis)
665	!!$!!-SC_ Fin albÃ©do moyen
666	!!$
667
668	!! 1.1 Calculation and assignment of soil albedo
669
670	DO ks = 1, 2! Loop over # of spectra
671
672	! If alb_bg_modis=TRUE, the background soil albedo map for the current simulated month is used
673	! If alb_bg_modis=FALSE and alb_bare_model=TRUE, the soil albedo calculation depends on soil moisture
674	! If alb_bg_modis=FALSE and alb_bare_model=FALSE, the mean soil albedo is used without the dependance on soil moisture
675	! see subroutines 'condveg_soilalb' and 'condveg_background_soilalb'
676	IF ( alb_bg_modis ) THEN
677	alb_bare(:,ks) = soilalb_bg(:,ks)
678	ELSE
679	IF ( alb_bare_model ) THEN
680	alb_bare(:,ks) = soilalb_wet(:,ks) + drysoil_frac(:) * (soilalb_dry(:,ks) - soilalb_wet(:,ks))
681	ELSE
682	alb_bare(:,ks) = soilalb_moy(:,ks)
683	ENDIF
684	ENDIF
685
686	! Soil albedo is weighed by fraction of bare soil
687	albedo(:,ks) = tot_bare_soil(:) * alb_bare(:,ks)
688
689	!! 1.2 Calculation of mean albedo of over the grid cell
690
691	! Calculation of mean albedo of over the grid cell and
692	! mean albedo of only vegetative PFTs over the grid cell
693	alb_veget(:,ks) = zero
694
695	DO jv = 2, nvm ! Loop over # of PFTs
696
697	! Mean albedo of grid cell for visible and near-infrared range
698	albedo(:,ks) = albedo(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
699
700	! Mean albedo of vegetation for visible and near-infrared range
701	alb_veget(:,ks) = alb_veget(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
702	ENDDO ! Loop over # of PFTs
703
704	ENDDO
705	END IF
706
707
708	!! 2. Calculate snow albedos on both total vegetated and total nobio surfaces
709
710	! The snow albedo could be either prescribed (in condveg_init.f90) or
711	! calculated following Chalita and Treut (1994).
712	! Check if the precribed value fixed_snow_albedo exists
713	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
714	snowa_veg(:,:) = fixed_snow_albedo
715	snowa_nobio(:,:,:) = fixed_snow_albedo
716	fraction_veg(:) = un - totfrac_nobio(:)
717	ELSE ! calculated following Chalita and Treut (1994)
718
719	!! 2.1 Calculate age dependence
720
721	! On vegetated surfaces
722	DO ji = 1, kjpindex
723	agefunc_veg(ji) = EXP(-snow_age(ji)/tcst_snowa)
724	ENDDO
725
726	! On non-vegtative surfaces
727	DO jv = 1, nnobio ! Loop over # nobio types
728	DO ji = 1, kjpindex
729	agefunc_nobio(ji,jv) = EXP(-snow_nobio_age(ji,jv)/tcst_snowa_nobio)
730	ENDDO
731	ENDDO
732
733	!! 2.1 Calculate snow albedo
734	! For vegetated surfaces
735	fraction_veg(:) = un - totfrac_nobio(:)
736	snowa_veg(:,:) = zero
737	IF (ok_dgvm) THEN
738	DO jb = 1, 2
739	DO ji = 1, kjpindex
740	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
741	snowa_veg(ji,jb) = snowa_veg(ji,jb) + &
742	tot_bare_soil(ji)/fraction_veg(ji) * ( snowa_aged_tmp(1,jb)+snowa_dec_tmp(1,jb)*agefunc_veg(ji) )
743	END IF
744	END DO
745	END DO
746
747	DO jb = 1, 2
748	DO jv = 2, nvm
749	DO ji = 1, kjpindex
750	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
751	snowa_veg(ji,jb) = snowa_veg(ji,jb) + &
752	veget(ji,jv)/fraction_veg(ji) * ( snowa_aged_tmp(jv,jb)+snowa_dec_tmp(jv,jb)*agefunc_veg(ji) )
753	ENDIF
754	ENDDO
755	ENDDO
756	ENDDO
757	ELSE
758	DO jb = 1, 2
759	DO jv = 1, nvm
760	DO ji = 1, kjpindex
761	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
762	snowa_veg(ji,jb) = snowa_veg(ji,jb) + &
763	veget_max(ji,jv)/fraction_veg(ji) * ( snowa_aged_tmp(jv,jb)+snowa_dec_tmp(jv,jb)*agefunc_veg(ji) )
764	ENDIF
765	ENDDO
766	ENDDO
767	ENDDO
768	ENDIF
769	!
770	! snow albedo on other surfaces
771	!
772	DO jb = 1, 2
773	DO jv = 1, nnobio
774	DO ji = 1, kjpindex
775	snowa_nobio(ji,jv,jb) = ( snowa_aged_tmp(1,jb) + snowa_dec_tmp(1,jb) * agefunc_nobio(ji,jv) )
776	ENDDO
777	ENDDO
778	ENDDO
779	ENDIF
780
781	!! 3. Update surface albedo
782
783	! Update surface albedo using the weighted sum of previous time step surface albedo,
784	! total vegetated fraction, total nobio fraction, snow cover fraction (both vegetated and
785	! non-vegetative surfaces), and snow albedo (both vegetated and non-vegetative surfaces).
786	! Although both visible and near-infrared surface albedo are presented, their calculations
787	! are the same.
788	DO jb = 1, 2
789
790	albedo(:,jb) = ( fraction_veg(:) ) * &
791	( (un-frac_snow_veg(:)) * albedo(:,jb) + &
792	( frac_snow_veg(:) ) * snowa_veg(:,jb) )
793	DO jv = 1, nnobio ! Loop over # nobio surfaces
794
795	IF ( jv .EQ. iice ) THEN
796	alb_nobio = alb_ice(jb)
797	ELSE
798	WRITE(numout,*) 'jv=',jv
799	WRITE(numout,*) 'DO NOT KNOW ALBEDO OF THIS SURFACE TYPE'
800	CALL ipslerr_p(3,'condveg_snow','DO NOT KNOW ALBEDO OF THIS SURFACE TYPE','','')
801	ENDIF
802
803	albedo(:,jb) = albedo(:,jb) + &
804	( frac_nobio(:,jv) ) * &
805	( (un-frac_snow_nobio(:,jv)) * alb_nobio + &
806	( frac_snow_nobio(:,jv) ) * snowa_nobio(:,jv,jb) )
807	ENDDO
808
809	END DO
810
811	! Calculate snow albedo
812	DO jb = 1, 2
813	albedo_snow(:,jb) = fraction_veg(:) * frac_snow_veg(:) * snowa_veg(:,jb)
814	DO jv = 1, nnobio
815	albedo_snow(:,jb) = albedo_snow(:,jb) + &
816	frac_nobio(:,jv) * frac_snow_nobio(:,jv) * snowa_nobio(:,jv,jb)
817	ENDDO
818	ENDDO
819
820	IF (printlev>=3) WRITE (numout,*) ' condveg_albedo done '
821
822	END SUBROUTINE condveg_albedo
823
824
825
826	!! ==============================================================================================================================
827	!! SUBROUTINE : condveg_frac_snow
828	!!
829	!>\BRIEF This subroutine calculates the fraction of snow on vegetation and nobio
830	!!
831	!! DESCRIPTION
832	!!
833	!! RECENT CHANGE(S): These calculations were previously done in condveg_snow.
834	!!
835	!! REFERENCE(S) :
836	!!
837	!! FLOWCHART : None
838	!! \n
839	!_ ================================================================================================================================
840
841	SUBROUTINE condveg_frac_snow(kjpindex, snow, snow_nobio, snowrho, snowdz, &
842	frac_snow_veg, frac_snow_nobio, totfrac_nobio)
843	!! 0. Variable and parameter declaration
844	!! 0.1 Input variables
845	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
846	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
847	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
848	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density at each snow layer
849	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer
850	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
851
852	!! 0.2 Output variables
853	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: frac_snow_veg !! Fraction of snow on vegetation (unitless ratio)
854	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
855
856	!! 0.3 Local variables
857	REAL(r_std), DIMENSION(kjpindex) :: snowrho_ave !! Average snow density
858	REAL(r_std), DIMENSION(kjpindex) :: snowdepth !! Snow depth
859	REAL(r_std), DIMENSION(kjpindex) :: snowrho_snowdz !! Snow rho time snowdz
860	INTEGER(i_std) :: jv
861
862	!! Calculate snow cover fraction for both total vegetated and total non-vegetative surfaces.
863	IF (ok_explicitsnow) THEN
864	snowdepth=sum(snowdz,2)
865	snowrho_snowdz=sum(snowrho*snowdz,2)
866	WHERE(snowdepth(:) .LT. min_sechiba)
867	frac_snow_veg(:) = 0.
868	ELSEWHERE
869	snowrho_ave(:)=snowrho_snowdz(:)/snowdepth(:)
870	frac_snow_veg(:) = tanh(snowdepth(:)/(0.025*(snowrho_ave(:)/50.)))
871	END WHERE
872	ELSE
873	frac_snow_veg(:) = MIN(MAX(snow(:),zero)/(MAX(snow(:),zero)+snowcri_alb*sn_dens/100.0),un)
874	END IF
875
876	WHERE (totfrac_nobio(:) .EQ. 1)
877	frac_snow_veg(:)=0.
878	END WHERE
879
880	frac_snow_nobio(:,:)=0.
881	!cdc frac_snow_nobio(:,iice) = MIN(MAX(snow(:),zero)/(MAX(snow(:),zero)+snowcri_alb*sn_dens/100.0),un)
882	!cdc MAR-ice80 => version frac_snow_nobio en tanh ( + bas si peu de neige dense)
883	snowdepth=SUM(snowdz,2)
884	snowrho_snowdz=sum(snowrho*snowdz,2)
885	WHERE(snowdepth(:) .LT. min_sechiba)
886	frac_snow_nobio(:,iice) = 0.
887	ELSEWHERE
888	snowrho_ave(:)=snowrho_snowdz(:)/snowdepth(:)
889	frac_snow_nobio(:,iice) = tanh(snowdepth(:)/(0.025*(snowrho_ave(:)/50.)))
890	ENDWHERE
891	!cdc MAR-ice80 fin modif
892
893	IF (printlev>=3) WRITE (numout,*) ' condveg_frac_snow done '
894
895	END SUBROUTINE condveg_frac_snow
896
897
898	!! ==============================================================================================================================
899	!! SUBROUTINE : condveg_soilalb
900	!!
901	!>\BRIEF This subroutine calculates the albedo of soil (without snow).
902	!!
903	!! DESCRIPTION This subroutine reads the soil colour maps in 1 x 1 deg resolution
904	!! from the Henderson-Sellers & Wilson database. These values are interpolated to
905	!! the model's resolution and transformed into
906	!! dry and wet albedos.\n
907	!!
908	!! If the soil albedo is calculated without the dependence of soil moisture, the
909	!! soil colour values are transformed into mean soil albedo values.\n
910	!!
911	!! The calculations follow the assumption that the grid of the data is regular and
912	!! it covers the globe. The calculation for the model grid are based on the borders
913	!! of the grid of the resolution.
914	!!
915	!! RECENT CHANGE(S): None
916	!!
917	!! CALCULATED MODULE VARIABLE(S): soilalb_dry for visible and near-infrared range,
918	!! soilalb_wet for visible and near-infrared range,
919	!! soilalb_moy for visible and near-infrared range
920	!!
921	!! REFERENCE(S) :
922	!! -Wilson, M.F., and A. Henderson-Sellers, 1985: A global archive of land cover and
923	!! soils data for use in general circulation climate models. J. Clim., 5, 119-143.
924	!!
925	!! FLOWCHART : None
926	!! \n
927	!_ ================================================================================================================================
928
929	SUBROUTINE condveg_soilalb(nbpt, lalo, neighbours, resolution, contfrac)
930
931	USE interpweight
932
933	IMPLICIT NONE
934
935
936	!! 0. Variable and parameter declaration
937
938	!! 0.1 Input variables
939
940	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
941	!! interpolated (unitless)
942	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
943	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
944	!! (1=N, 2=E, 3=S, 4=W)
945	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
946	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
947
948	!! 0.4 Local variables
949
950	CHARACTER(LEN=80) :: filename !! Filename of soil colour map
951	INTEGER(i_std) :: i, ib, ip, nbexp !! Indices
952	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
953	REAL(r_std), DIMENSION(nbpt) :: asoilcol !! Availability of the soilcol interpolation
954	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
955	!! (variabletypevals(1) = -un, not used)
956	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: soilcolrefrac !! soilcol fractions re-dimensioned
957	REAL(r_std) :: vmin, vmax !! min/max values to use for the
958	!! renormalization
959	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
960	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
961	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
962	!! 'XYKindTime': Input values are kinds
963	!! of something with a temporal
964	!! evolution on the dx*dy matrix'
965	LOGICAL :: nonegative !! whether negative values should be removed
966	CHARACTER(LEN=50) :: maskingtype !! Type of masking
967	!! 'nomask': no-mask is applied
968	!! 'mbelow': take values below maskvals(1)
969	!! 'mabove': take values above maskvals(1)
970	!! 'msumrange': take values within 2 ranges;
971	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
972	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
973	!! (normalized by maskvals(3))
974	!! 'var': mask values are taken from a
975	!! variable inside the file (>0)
976	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
977	!! `maskingtype')
978	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
979	CHARACTER(LEN=250) :: msg
980	INTEGER :: fopt
981	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: vecpos
982	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: solt
983
984	!_ ================================================================================================================================
985	!! 1. Open file and allocate memory
986
987	! Open file with soil colours
988
989	!Config Key = SOILALB_FILE
990	!Config Desc = Name of file from which the bare soil albedo
991	!Config Def = soils_param.nc
992	!Config If = NOT(IMPOSE_AZE)
993	!Config Help = The name of the file to be opened to read the soil types from
994	!Config which we derive then the bare soil albedos. This file is 1x1
995	!Config deg and based on the soil colors defined by Wilson and Henderson-Seller.
996	!Config Units = [FILE]
997	!
998	filename = 'soils_param.nc'
999	CALL getin_p('SOILALB_FILE',filename)
1000
1001
1002	ALLOCATE(soilcolrefrac(nbpt, classnb), STAT=ALLOC_ERR)
1003	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable soilcolrefrac','','')
1004	ALLOCATE(vecpos(classnb), STAT=ALLOC_ERR)
1005	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable vecpos','','')
1006	ALLOCATE(solt(classnb), STAT=ALLOC_ERR)
1007	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable solt','','')
1008
1009	! Assigning values to vmin, vmax
1010	vmin = 1.0
1011	vmax = classnb
1012
1013	ALLOCATE(variabletypevals(classnb),STAT=ALLOC_ERR)
1014	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variabletypevals','','')
1015	variabletypevals = -un
1016
1017	!! Variables for interpweight
1018	! Type of calculation of cell fractions
1019	fractype = 'default'
1020	! Name of the longitude and latitude in the input file
1021	lonname = 'nav_lon'
1022	latname = 'nav_lat'
1023	! Should negative values be set to zero from input file?
1024	nonegative = .FALSE.
1025	! Type of mask to apply to the input data (see header for more details)
1026	maskingtype = 'mabove'
1027	! Values to use for the masking
1028	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
1029	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
1030	namemaskvar = ''
1031
1032	! Interpolate variable soilcolor
1033	variablename = 'soilcolor'
1034	IF (printlev_loc >= 2) WRITE(numout,*) "condveg_soilalb: Start interpolate " &
1035	// TRIM(filename) // " for variable " // TRIM(variablename)
1036	CALL interpweight_2D(nbpt, classnb, variabletypevals, lalo, resolution, neighbours, &
1037	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1038	maskvals, namemaskvar, 0, 0, -1, fractype, &
1039	-1., -1., soilcolrefrac, asoilcol)
1040	IF (printlev_loc >= 5) WRITE(numout,*)' condveg_soilalb after interpweight_2D'
1041
1042	! Check how many points with soil information are found
1043	nbexp = 0
1044
1045	soilalb_dry(:,:) = zero
1046	soilalb_wet(:,:) = zero
1047	soilalb_moy(:,:) = zero
1048	IF (printlev_loc >= 5) THEN
1049	WRITE(numout,*)' condveg_soilalb before starting loop nbpt:', nbpt
1050	WRITE(numout,*)' condveg_soilalb initial values classnb: ',classnb
1051	WRITE(numout,*)' condveg_soilalb vis_dry. SUM:',SUM(vis_dry),' vis_dry= ',vis_dry
1052	WRITE(numout,*)' condveg_soilalb nir_dry. SUM:',SUM(nir_dry),' nir_dry= ',nir_dry
1053	WRITE(numout,*)' condveg_soilalb vis_wet. SUM:',SUM(vis_wet),' vis_wet= ',vis_wet
1054	WRITE(numout,*)' condveg_soilalb nir_wet. SUM:',SUM(nir_wet),' nir_wet= ',nir_wet
1055	END IF
1056
1057	DO ib=1,nbpt ! Loop over domain size
1058
1059	! vecpos: List of positions where textures were not zero
1060	! vecpos(1): number of not null textures found
1061	vecpos = interpweight_ValVecR(soilcolrefrac(ib,:),classnb,zero,'neq')
1062	fopt = vecpos(1)
1063	IF (fopt == classnb) THEN
1064	! All textures are not zero
1065	solt(:) = (/(i,i=1,classnb)/)
1066	ELSE IF (fopt == 0) THEN
1067	WRITE(numout,*)' condveg_soilalb: for point=', ib, ' no soil class!'
1068	ELSE
1069	DO ip = 1,fopt
1070	solt(ip) = vecpos(ip+1)
1071	END DO
1072	END IF
1073
1074	!! 3. Compute the average bare soil albedo parameters
1075
1076	IF ( (fopt .EQ. 0) .OR. (asoilcol(ib) .LT. min_sechiba)) THEN
1077	! Initialize with mean value if no points were interpolated or if no data was found
1078	nbexp = nbexp + 1
1079	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1080	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1081	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1082	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1083	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
1084	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
1085	ELSE
1086	! If points were interpolated
1087	DO ip=1, fopt
1088	IF ( solt(ip) .LE. classnb) THEN
1089	! Set to zero if the value is below min_sechiba
1090	IF (soilcolrefrac(ib,solt(ip)) < min_sechiba) soilcolrefrac(ib,solt(ip)) = zero
1091
1092	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(solt(ip))*soilcolrefrac(ib,solt(ip))
1093	soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(solt(ip))*soilcolrefrac(ib,solt(ip))
1094	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(solt(ip))*soilcolrefrac(ib,solt(ip))
1095	soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(solt(ip))*soilcolrefrac(ib,solt(ip))
1096	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(solt(ip))* &
1097	soilcolrefrac(ib,solt(ip))
1098	soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(solt(ip))* &
1099	soilcolrefrac(ib,solt(ip))
1100	ELSE
1101	msg = 'The file contains a soil color class which is incompatible with this program'
1102	CALL ipslerr_p(3,'condveg_soilalb',TRIM(msg),'','')
1103	ENDIF
1104	ENDDO
1105	ENDIF
1106
1107	ENDDO
1108
1109	IF ( nbexp .GT. 0 ) THEN
1110	WRITE(numout,*) 'condveg_soilalb _______'
1111	WRITE(numout,*) 'condveg_soilalb: The interpolation of the bare soil albedo had ', nbexp
1112	WRITE(numout,*) 'condveg_soilalb: points without data. This are either coastal points or'
1113	WRITE(numout,*) 'condveg_soilalb: ice covered land.'
1114	WRITE(numout,*) 'condveg_soilalb: The problem was solved by using the average of all soils'
1115	WRITE(numout,*) 'condveg_soilalb: in dry and wet conditions'
1116	WRITE(numout,*) 'condveg_soilalb: Use the diagnostic output field asoilcol to see location of these points'
1117	ENDIF
1118
1119	DEALLOCATE (soilcolrefrac)
1120	DEALLOCATE (variabletypevals)
1121
1122	! Write diagnostics
1123	CALL xios_orchidee_send_field("asoilcol",asoilcol)
1124
1125
1126	IF (printlev_loc >= 3) WRITE(numout,*)' condveg_soilalb ended'
1127
1128	END SUBROUTINE condveg_soilalb
1129
1130
1131	!! ==============================================================================================================================
1132	!! SUBROUTINE : condveg_background_soilalb
1133	!!
1134	!>\BRIEF This subroutine reads the albedo of bare soil
1135	!!
1136	!! DESCRIPTION This subroutine reads the background albedo map in 0.5 x 0.5 deg resolution
1137	!! derived from JRCTIP product to be used as bare soil albedo. These values are then interpolated
1138	!! to the model's resolution.\n
1139	!!
1140	!! RECENT CHANGE(S): None
1141	!!
1142	!! MAIN OUTPUT VARIABLE(S): soilalb_bg for visible and near-infrared range
1143	!!
1144	!! REFERENCES : None
1145	!!
1146	!! FLOWCHART : None
1147	!! \n
1148	!_ ================================================================================================================================
1149
1150	SUBROUTINE condveg_background_soilalb(nbpt, lalo, neighbours, resolution, contfrac)
1151
1152	USE interpweight
1153
1154	IMPLICIT NONE
1155
1156	!! 0. Variable and parameter declaration
1157
1158	!! 0.1 Input variables
1159
1160	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
1161	!! interpolated (unitless)
1162	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
1163	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
1164	!! (1=N, 2=E, 3=S, 4=W)
1165	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
1166	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
1167
1168	!! 0.4 Local variables
1169
1170	CHARACTER(LEN=80) :: filename !! Filename of background albedo
1171	REAL(r_std), DIMENSION(nbpt) :: aalb_bg !! Availability of the interpolation
1172	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! Latitudes and longitudes read from input file
1173	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel !! Help variable to read file data and allocate memory
1174	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu !! Help variable to read file data and allocate memory
1175	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
1176	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: soilalbedo_bg !! Help variable to read file data and allocate memory
1177	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
1178	REAL(r_std) :: vmin, vmax !! min/max values to use for the
1179	!! renormalization
1180	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
1181	CHARACTER(LEN=250) :: maskvname !! Variable to read the mask from
1182	!! the file
1183	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
1184	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
1185	!! 'XYKindTime': Input values are kinds
1186	!! of something with a temporal
1187	!! evolution on the dx*dy matrix'
1188	LOGICAL :: nonegative !! whether negative values should be removed
1189	CHARACTER(LEN=50) :: maskingtype !! Type of masking
1190	!! 'nomask': no-mask is applied
1191	!! 'mbelow': take values below maskvals(1)
1192	!! 'mabove': take values above maskvals(1)
1193	!! 'msumrange': take values within 2 ranges;
1194	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
1195	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
1196	!! (normalized by maskedvals(3))
1197	!! 'var': mask values are taken from a
1198	!! variable inside the file (>0)
1199	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
1200	!! `maskingtype')
1201	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
1202	REAL(r_std) :: albbg_norefinf !! No value
1203	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: albbg_default !! Default value
1204
1205	!_ ================================================================================================================================
1206
1207	!! 1. Open file and allocate memory
1208
1209	! Open file with background albedo
1210
1211	!Config Key = ALB_BG_FILE
1212	!Config Desc = Name of file from which the background albedo is read
1213	!Config Def = alb_bg.nc
1214	!Config If =
1215	!Config Help = The name of the file to be opened to read background albedo
1216	!Config Units = [FILE]
1217	!
1218	filename = 'alb_bg.nc'
1219	CALL getin_p('ALB_BG_FILE',filename)
1220
1221
1222	ALLOCATE(albbg_default(2), STAT=ALLOC_ERR)
1223	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'condveg_background_soilalb','Pb in allocation for albbg_default','','')
1224
1225	! For this case there are not types/categories. We have 'only' a continuos field
1226	! Assigning values to vmin, vmax
1227
1228	vmin = 0.
1229	vmax = 9999.
1230
1231	!! Variables for interpweight
1232	! Type of calculation of cell fractions (not used here)
1233	fractype = 'default'
1234	! Name of the longitude and latitude in the input file
1235	lonname = 'longitude'
1236	latname = 'latitude'
1237	! Default value when no value is get from input file
1238	albbg_default(ivis) = 0.129
1239	albbg_default(inir) = 0.247
1240	! Reference value when no value is get from input file (not used here)
1241	albbg_norefinf = undef_sechiba
1242	! Should negative values be set to zero from input file?
1243	nonegative = .FALSE.
1244	! Type of mask to apply to the input data (see header for more details)
1245	maskingtype = 'var'
1246	! Values to use for the masking (here not used)
1247	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
1248	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var')
1249	namemaskvar = 'mask'
1250
1251	! There is a variable for each chanel 'infrared' and 'visible'
1252	! Interpolate variable bg_alb_vis
1253	variablename = 'bg_alb_vis'
1254	IF (printlev_loc >= 2) WRITE(numout,*) "condveg_background_soilalb: Start interpolate " &
1255	// TRIM(filename) // " for variable " // TRIM(variablename)
1256	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
1257	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1258	maskvals, namemaskvar, -1, fractype, albbg_default(inir), albbg_norefinf, &
1259	soilalb_bg(:,ivis), aalb_bg)
1260	IF (printlev_loc >= 5) WRITE(numout,*)" condveg_background_soilalb after InterpWeight2Dcont for '" // &
1261	TRIM(variablename) // "'"
1262
1263	! Interpolate variable bg_alb_nir in the same file
1264	variablename = 'bg_alb_nir'
1265	IF (printlev_loc >= 2) WRITE(numout,*) "condveg_background_soilalb: Start interpolate " &
1266	// TRIM(filename) // " for variable " // TRIM(variablename)
1267	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
1268	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1269	maskvals, namemaskvar, -1, fractype, albbg_default(inir), albbg_norefinf, &
1270	soilalb_bg(:,inir), aalb_bg)
1271	IF (printlev_loc >= 5) WRITE(numout,*)" condveg_background_soilalb after InterpWeight2Dcont for '" // &
1272	TRIM(variablename) // "'"
1273
1274	IF (ALLOCATED(albbg_default)) DEALLOCATE(albbg_default)
1275
1276	IF (printlev_loc >= 3) WRITE(numout,*)' condveg_background_soilalb ended'
1277
1278	! Write diagnostics
1279	CALL xios_orchidee_send_field("aalb_bg",aalb_bg)
1280
1281	END SUBROUTINE condveg_background_soilalb
1282
1283
1284	!! ==============================================================================================================================
1285	!! SUBROUTINE : condveg_z0cdrag
1286	!!
1287	!>\BRIEF Computation of grid average of roughness length by calculating
1288	!! the drag coefficient.
1289	!!
1290	!! DESCRIPTION : This routine calculates the mean roughness height and mean
1291	!! effective roughness height over the grid cell. The mean roughness height (z0)
1292	!! is computed by averaging the drag coefficients \n
1293	!!
1294	!! \latexonly
1295	!! \input{z0cdrag1.tex}
1296	!! \endlatexonly
1297	!! \n
1298	!!
1299	!! where C is the drag coefficient at the height of the vegetation, kappa is the
1300	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
1301	!! The reference level for z needs to be high enough above the canopy to avoid
1302	!! singularities of the LOG. This height is set to minimum 10m above ground.
1303	!! The drag coefficient increases with roughness height to represent the greater
1304	!! turbulence generated by rougher surfaces.
1305	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
1306	!!
1307	!! The roughness height for the non-vegetative surfaces is calculated in a second step.
1308	!! In order to calculate the transfer coefficients the
1309	!! effective roughness height is calculated. This effective value is the difference
1310	!! between the height of the vegetation and the zero plane displacement height.\nn
1311	!!
1312	!! RECENT CHANGE(S): None
1313	!!
1314	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
1315	!!
1316	!! REFERENCE(S) : None
1317	!!
1318	!! FLOWCHART : None
1319	!! \n
1320	!_ ================================================================================================================================
1321
1322	SUBROUTINE condveg_z0cdrag (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, tot_bare_soil, &
1323	& z0m, z0h, roughheight)
1324
1325	!! 0. Variable and parameter declaration
1326
1327	!! 0.1 Input variables
1328
1329	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1330	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1331	!! (m^2 m^{-2})
1332	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1333	!! (= ind*cn_ind) (m^2 m^{-2})
1334	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1335	!! i.e. continental ice, lakes, etc. (unitless)
1336	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1337	!! i.e. continental ice, lakes, etc. (unitless)
1338	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
1339	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1340	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
1341
1342	!! 0.2 Output variables
1343
1344	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0m !! Roughness height for momentum (m)
1345	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0h !! Roughness height for heat (m)
1346	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1347
1348	!! 0.3 Modified variables
1349
1350	!! 0.4 Local variables
1351
1352	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1353	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1354	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
1355	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1356	REAL(r_std), DIMENSION(kjpindex) :: d_veg !! PFT coverage of vegetative PFTs
1357	!! (= ind*cn_ind) (m^2 m^{-2})
1358	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1359	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
1360	!! i.e. continental ice, lakes, etc.
1361
1362	!_ ================================================================================================================================
1363
1364	!! 1. Preliminary calculation
1365
1366	! Set maximal height of first layer
1367	ztmp(:) = MAX(10., zlev(:))
1368
1369	! Calculate roughness for non-vegetative surfaces
1370	! with the von Karman constant
1371	z0m(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/z0_bare))**2
1372	z0h(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/(z0_bare/ratio_z0m_z0h(1))))**2
1373	! Fraction of bare soil
1374	sumveg(:) = tot_bare_soil(:)
1375
1376	! Set average vegetation height to zero
1377	ave_height(:) = zero
1378
1379	!! 2. Calculate the mean roughness height
1380
1381	! Calculate the mean roughness height of
1382	! vegetative PFTs over the grid cell
1383	DO jv = 2, nvm
1384
1385	! In the case of forest, use parameter veget_max because
1386	! tree trunks influence the roughness even when there are no leaves
1387	IF ( is_tree(jv) ) THEN
1388	! In the case of grass, use parameter veget because grasses
1389	! only influence the roughness during the growing season
1390	d_veg(:) = veget_max(:,jv)
1391	ELSE
1392	! grasses only have an influence if they are really there!
1393	d_veg(:) = veget(:,jv)
1394	ENDIF
1395
1396	! Calculate the average roughness over the grid cell:
1397	! The unitless drag coefficient is per vegetative PFT
1398	! calculated by use of the von Karman constant, the height
1399	! of the first layer and the roughness. The roughness
1400	! is calculated as the vegetation height per PFT
1401	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
1402	! If this scaled value is lower than 0.01 then the value for
1403	! the roughness of bare soil (0.01) is used.
1404	! The sum over all PFTs gives the average roughness
1405	! per grid cell for the vegetative PFTs.
1406	z0m(:) = z0m(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)z0_over_height(jv),z0_bare)))*2
1407	z0h(:) = z0h(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/(MAX(height(:,jv)z0_over_height(jv),z0_bare)/ratio_z0m_z0h(jv))))*2
1408
1409	! Sum of bare soil and fraction vegetated fraction
1410	sumveg(:) = sumveg(:) + d_veg(:)
1411
1412	! Weigh height of vegetation with maximal cover fraction
1413	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1414
1415	ENDDO
1416
1417	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1418
1419	! Search for pixels with vegetated part to normalise
1420	! roughness height
1421	WHERE ( sumveg(:) .GT. zero )
1422	z0m(:) = z0m(:) / sumveg(:)
1423	z0h(:) = z0h(:) / sumveg(:)
1424	ENDWHERE
1425	! Calculate fraction of roughness for vegetated part
1426	z0m(:) = (un - totfrac_nobio(:)) * z0m(:)
1427	z0h(:) = (un - totfrac_nobio(:)) * z0h(:)
1428
1429	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1430
1431	! Set rougness for ice
1432	IF ( jv .EQ. iice ) THEN
1433	z0_nobio = z0_ice
1434	ELSE
1435	WRITE(numout,*) 'jv=',jv
1436	WRITE(numout,*) 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1437	CALL ipslerr_p(3,'condveg_z0cdrag','DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1438	ENDIF
1439
1440	! Sum of vegetative roughness length and non-vegetative
1441	! roughness length
1442	z0m(:) = z0m(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1443	z0h(:) = z0h(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio/ratio_z0m_z0h(1)))**2
1444
1445	ENDDO ! Loop over # of non-vegative surfaces
1446
1447	!! 4. Calculate the zero plane displacement height and effective roughness length
1448
1449	! Take the exponential of the roughness
1450	z0m(:) = ztmp(:) / EXP(ct_karman/SQRT(z0m(:)))
1451	z0h(:) = ztmp(:) / EXP(ct_karman/SQRT(z0h(:)))
1452
1453	! Compute the zero plane displacement height which
1454	! is an equivalent height for the absorption of momentum
1455	zhdispl(:) = ave_height(:) * height_displacement
1456
1457	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1458	! This is the height over which the roughness acts. It combines the
1459	! zero plane displacement height and the vegetation height.
1460	roughheight(:) = ave_height(:) - zhdispl(:)
1461
1462	END SUBROUTINE condveg_z0cdrag
1463
1464
1465	!! ==============================================================================================================================
1466	!! SUBROUTINE : condveg_z0cdrag_dyn
1467	!!
1468	!>\BRIEF Computation of grid average of roughness length by calculating
1469	!! the drag coefficient based on formulation proposed by Su et al. (2001).
1470	!!
1471	!! DESCRIPTION : This routine calculates the mean roughness height and mean
1472	!! effective roughness height over the grid cell. The mean roughness height (z0)
1473	!! is computed by averaging the drag coefficients \n
1474	!!
1475	!! \latexonly
1476	!! \input{z0cdrag1.tex}
1477	!! \endlatexonly
1478	!! \n
1479	!!
1480	!! where C is the drag coefficient at the height of the vegetation, kappa is the
1481	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
1482	!! The reference level for z needs to be high enough above the canopy to avoid
1483	!! singularities of the LOG. This height is set to minimum 10m above ground.
1484	!! The drag coefficient increases with roughness height to represent the greater
1485	!! turbulence generated by rougher surfaces.
1486	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
1487	!! In the formulation of Su et al. (2001), one distinguishes the roughness height for
1488	!! momentum (z0m) and the one for heat (z0h).
1489	!! z0m is computed as a function of LAI (z0m increases with LAI) and z0h is computed
1490	!! with a so-called kB-1 term (z0m/z0h=exp(kB-1))
1491	!!
1492	!! RECENT CHANGE(S): Written by N. Vuichard (2016)
1493	!!
1494	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
1495	!!
1496	!! REFERENCE(S) :
1497	!! - Su, Z., Schmugge, T., Kustas, W.P., Massman, W.J., 2001. An Evaluation of Two Models for
1498	!! Estimation of the Roughness Height for Heat Transfer between the Land Surface and the Atmosphere. J. Appl.
1499	!! Meteorol. 40, 1933â1951. doi:10.1175/1520-0450(2001)
1500	!! - Ershadi, A., McCabe, M.F., Evans, J.P., Wood, E.F., 2015. Impact of model structure and parameterization
1501	!! on Penman-Monteith type evaporation models. J. Hydrol. 525, 521â535. doi:10.1016/j.jhydrol.2015.04.008
1502	!!
1503	!! FLOWCHART : None
1504	!! \n
1505	!_ ================================================================================================================================
1506
1507	SUBROUTINE condveg_z0cdrag_dyn (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, &
1508	& temp_air, pb, u, v, lai, frac_snow_veg, z0m, z0h, roughheight)
1509
1510	!! 0. Variable and parameter declaration
1511
1512	!! 0.1 Input variables
1513
1514	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1515	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1516	!! (m^2 m^{-2})
1517	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1518	!! (= ind*cn_ind) (m^2 m^{-2})
1519	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1520	!! i.e. continental ice, lakes, etc. (unitless)
1521	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1522	!! i.e. continental ice, lakes, etc. (unitless)
1523	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
1524	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1525	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: temp_air !! 2m air temperature (K)
1526	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
1527	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
1528	!! @tex $(m.s^{-1})$ @endtex
1529	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
1530	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
1531	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1532	!! 0.2 Output variables
1533
1534	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0m !! Roughness height for momentum (m)
1535	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0h !! Roughness height for heat (m)
1536	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1537
1538	!! 0.3 Modified variables
1539
1540	!! 0.4 Local variables
1541
1542	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1543	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1544	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
1545	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1546	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1547	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
1548	REAL(r_std), DIMENSION(kjpindex) :: z0m_pft !! Roughness height for momentum for a specific PFT
1549	REAL(r_std), DIMENSION(kjpindex) :: z0h_pft !! Roughness height for heat for a specific PFT
1550	REAL(r_std), DIMENSION(kjpindex) :: eta !! Ratio of friction velocity to the wind speed at the canopy top - See Ershadi et al. (2015)
1551	REAL(r_std), DIMENSION(kjpindex) :: eta_ec !! Within-canopy wind speed profile estimation coefficient - See Ershadi et al. (2015)
1552	REAL(r_std), DIMENSION(kjpindex) :: Ct_star !! Heat transfer coefficient of the soil - see Su et al. (2001)
1553	REAL(r_std), DIMENSION(kjpindex) :: kBs_m1 !! Canopy model of Brutsaert (1982) for a bare soil surface - used in the calculation of kB_m1 (see Ershadi et al. (2015))
1554	REAL(r_std), DIMENSION(kjpindex) :: kB_m1 !! kB**-1: Term used in the calculation of z0h where B-1 is the inverse Stanton number (see Ershadi et al. (2015))
1555	REAL(r_std), DIMENSION(kjpindex) :: fc !! fractional canopy coverage
1556	REAL(r_std), DIMENSION(kjpindex) :: fs !! fractional soil coverage
1557	REAL(r_std), DIMENSION(kjpindex) :: Reynolds !! Reynolds number
1558	REAL(r_std), DIMENSION(kjpindex) :: wind !! wind Speed (m)
1559	REAL(r_std), DIMENSION(kjpindex) :: u_star !! friction velocity
1560	REAL(r_std), DIMENSION(kjpindex) :: z0_ground !! z0m value used for ground surface
1561	!_ ================================================================================================================================
1562
1563	!! 1. Preliminary calculation
1564
1565	! Set maximal height of first layer
1566	ztmp(:) = MAX(10., zlev(:))
1567
1568	z0_ground(:) = (1.-frac_snow_veg(:))z0_bare + frac_snow_veg(:)z0_bare/10.
1569
1570	! Calculate roughness for non-vegetative surfaces
1571	! with the von Karman constant
1572	z0m(:) = veget_max(:,1) * (ct_karman/LOG(ztmp(:)/z0_ground(:)))**2
1573
1574	wind(:) = SQRT(u(:)u(:)+v(:)v(:))
1575	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG(zlev(:)/z0_ground(:))
1576	Reynolds(:) = z0_ground(:) * u_star(:) &
1577	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1578
1579	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1580
1581	z0h(:) = veget_max(:,1) * (ct_karman/LOG(ztmp(:)/(z0_ground(:)/ exp(kBs_m1(:))) ))**2
1582
1583
1584	! Fraction of bare soil
1585	sumveg(:) = veget_max(:,1)
1586
1587	! Set average vegetation height to zero
1588	ave_height(:) = zero
1589
1590	!! 2. Calculate the mean roughness height
1591
1592	! Calculate the mean roughness height of
1593	! vegetative PFTs over the grid cell
1594	DO jv = 2, nvm
1595
1596	WHERE(veget_max(:,jv) .GT. zero)
1597	! Calculate the average roughness over the grid cell:
1598	! The unitless drag coefficient is per vegetative PFT
1599	! calculated by use of the von Karman constant, the height
1600	! of the first layer and the roughness. The roughness
1601	! is calculated as the vegetation height per PFT
1602	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
1603	! If this scaled value is lower than 0.01 then the value for
1604	! the roughness of bare soil (0.01) is used.
1605	! The sum over all PFTs gives the average roughness
1606	! per grid cell for the vegetative PFTs.
1607	eta(:) = c1 - c2 * exp(-c3 * Cdrag_foliage * lai(:,jv))
1608
1609	z0m_pft(:) = (height(:,jv)(1-height_displacement)(exp(-ct_karman/eta(:))-exp(-ct_karman/(c1-c2)))) &
1610	+ z0_ground(:)
1611
1612	z0m(:) = z0m(:) + veget_max(:,jv) * (ct_karman/LOG(ztmp(:)/z0m_pft(:)))**2
1613
1614	fc(:) = veget(:,jv)/veget_max(:,jv)
1615	fs(:) = 1. - fc(:)
1616
1617	eta_ec(:) = ( Cdrag_foliage * lai(:,jv)) / (2 * eta(:)*eta(:))
1618	wind(:) = SQRT(u(:)u(:)+v(:)v(:))
1619	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG((zlev(:)+(height(:,jv)*(1-height_displacement)))/z0m_pft(:))
1620	Reynolds(:) = z0_ground(:) * u_star(:) &
1621	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1622
1623	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1624	Ct_star(:) = Prandtl*(-2./3.) SQRT(1./Reynolds(:))
1625
1626	WHERE(lai(:,jv) .GT. min_sechiba)
1627	kB_m1(:) = (ct_karman * Cdrag_foliage) / (4 * Ct * eta(:) * (1 - exp(-eta_ec(:)/2.))) * fc(:)**2. &
1628	+ 2fc(:)fs(:) * (ct_karman * eta(:) * z0m_pft(:) / height(:,jv)) / Ct_star(:) &
1629	+ kBs_m1(:) * fs(:)**2.
1630	ELSEWHERE
1631	kB_m1(:) = kBs_m1(:) * fs(:)**2.
1632	ENDWHERE
1633
1634	z0h_pft(:) = z0m_pft(:) / exp(kB_m1(:))
1635
1636	z0h(:) = z0h(:) + veget_max(:,jv) * (ct_karman/LOG(ztmp(:)/z0h_pft(:)))**2
1637
1638	! Sum of bare soil and fraction vegetated fraction
1639	sumveg(:) = sumveg(:) + veget_max(:,jv)
1640
1641	! Weigh height of vegetation with maximal cover fraction
1642	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1643
1644	ENDWHERE
1645	ENDDO
1646
1647	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1648
1649	! Search for pixels with vegetated part to normalise
1650	! roughness height
1651	WHERE ( sumveg(:) .GT. zero )
1652	z0h(:) = z0h(:) / sumveg(:)
1653	z0m(:) = z0m(:) / sumveg(:)
1654	ENDWHERE
1655
1656	! Calculate fraction of roughness for vegetated part
1657	z0h(:) = (un - totfrac_nobio(:)) * z0h(:)
1658	z0m(:) = (un - totfrac_nobio(:)) * z0m(:)
1659
1660	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1661
1662	! Set rougness for ice
1663	IF ( jv .EQ. iice ) THEN
1664	z0_nobio = z0_ice
1665	ELSE
1666	WRITE(numout,*) 'jv=',jv
1667	WRITE(numout,*) 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1668	CALL ipslerr_p(3,'condveg_z0cdrag_dyn','DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1669	ENDIF
1670
1671	! Sum of vegetative roughness length and non-vegetative roughness length
1672	! Note that z0m could be made dependent of frac_snow_nobio
1673	z0m(:) = z0m(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1674
1675	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG(zlev(:)/z0_nobio)
1676	Reynolds(:) = z0_nobio * u_star(:) &
1677	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1678
1679	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1680
1681	z0h(:) = z0h(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/(z0_nobio/ exp(kBs_m1(:))) ))**2
1682	ENDDO ! Loop over # of non-vegative surfaces
1683
1684	!! 4. Calculate the zero plane displacement height and effective roughness length
1685	! Take the exponential of the roughness
1686	z0h(:) = ztmp(:) / EXP(ct_karman/SQRT(z0h(:)))
1687	z0m(:) = ztmp(:) / EXP(ct_karman/SQRT(z0m(:)))
1688
1689	! Compute the zero plane displacement height which
1690	! is an equivalent height for the absorption of momentum
1691	zhdispl(:) = ave_height(:) * height_displacement
1692
1693	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1694	! This is the height over which the roughness acts. It combines the
1695	! zero plane displacement height and the vegetation height.
1696	roughheight(:) = ave_height(:) - zhdispl(:)
1697
1698	END SUBROUTINE condveg_z0cdrag_dyn
1699
1700
1701	END MODULE condveg

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

Download in other formats: