Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_sechiba/condveg.f90

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

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

Download in other formats: