Context Navigation

source: branches/publications/ORCHIDEE_DFv1.0_site/src_sechiba/condveg.f90 @ 7326

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

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

Download in other formats: