Context Navigation

condveg.f90

Last change on this file was 3952, checked in by albert.jornet, 8 years ago

Fix: for the thermics: (1) revise the calculation of thermal conductivity when SOC insulating effect is activated, which follows Lawrence&Slater,2008

(2) add the option to read prescribed SOC map only used in thermics for the SOC insulating effect.

for the permafrost carbon dynamics:
(1) corrected a bug related to time step
(2) change the condition for cryoturbation, and the default configuration for cryoturbation method

Done by dan.

Property svn:keywords set to Revision Date HeadURL Date Author Revision

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

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

Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_sechiba/condveg.f90

Download in other formats: