Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_sechiba/condveg.f90 @ 7346

Last change on this file since 7346 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

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

File size: 92.0 KB

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

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

Download in other formats: