Context Navigation

source: branches/ORCHIDEE_Quest/ORCHIDEE/src_sechiba/condveg.f90 @ 7406

Last change on this file since 7406 was 5470, checked in by josefine.ghattas, 6 years ago

Ticket #431 : Remove old snow scheme and option OK_EXPLICITSNOW. Only explicitsnow scheme is now possible to use. Removed therefor the subroutines hydrol_snow, enerbil_fusion and thermosoil_getdiff_old_thermix_with_snow.

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

File size: 93.5 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_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_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_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,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
909	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density at each snow layer
910	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer
911
912	!! 0.2 Output variables
913	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: frac_snow_veg !! Fraction of snow on vegetation (unitless ratio)
914	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
915
916	!! 0.3 Local variables
917	REAL(r_std), DIMENSION(kjpindex) :: snowrho_ave !! Average snow density
918	REAL(r_std), DIMENSION(kjpindex) :: snowdepth !! Snow depth
919	REAL(r_std), DIMENSION(kjpindex) :: snowrho_snowdz !! Snow rho time snowdz
920	INTEGER(i_std) :: jv
921
922	!! Calculate snow cover fraction for both total vegetated and total non-vegetative surfaces.
923	snowdepth=sum(snowdz,2)
924	snowrho_snowdz=sum(snowrho*snowdz,2)
925	WHERE(snowdepth(:) .LT. min_sechiba)
926	frac_snow_veg(:) = 0.
927	ELSEWHERE
928	snowrho_ave(:)=snowrho_snowdz(:)/snowdepth(:)
929	frac_snow_veg(:) = tanh(snowdepth(:)/(0.025*(snowrho_ave(:)/50.)))
930	END WHERE
931
932	DO jv = 1, nnobio
933	frac_snow_nobio(:,jv) = MIN(MAX(snow_nobio(:,jv),zero)/(MAX(snow_nobio(:,jv),zero)+snowcri_alb*sn_dens/100.0),un)
934	ENDDO
935
936	IF (printlev>=3) WRITE (numout,*) ' condveg_frac_snow done '
937
938	END SUBROUTINE condveg_frac_snow
939
940
941	!! ==============================================================================================================================
942	!! SUBROUTINE : condveg_soilalb
943	!!
944	!>\BRIEF This subroutine calculates the albedo of soil (without snow).
945	!!
946	!! DESCRIPTION This subroutine reads the soil colour maps in 1 x 1 deg resolution
947	!! from the Henderson-Sellers & Wilson database. These values are interpolated to
948	!! the model's resolution and transformed into
949	!! dry and wet albedos.\n
950	!!
951	!! If the soil albedo is calculated without the dependence of soil moisture, the
952	!! soil colour values are transformed into mean soil albedo values.\n
953	!!
954	!! The calculations follow the assumption that the grid of the data is regular and
955	!! it covers the globe. The calculation for the model grid are based on the borders
956	!! of the grid of the resolution.
957	!!
958	!! RECENT CHANGE(S): None
959	!!
960	!! CALCULATED MODULE VARIABLE(S): soilalb_dry for visible and near-infrared range,
961	!! soilalb_wet for visible and near-infrared range,
962	!! soilalb_moy for visible and near-infrared range
963	!!
964	!! REFERENCE(S) :
965	!! -Wilson, M.F., and A. Henderson-Sellers, 1985: A global archive of land cover and
966	!! soils data for use in general circulation climate models. J. Clim., 5, 119-143.
967	!!
968	!! FLOWCHART : None
969	!! \n
970	!_ ================================================================================================================================
971
972	SUBROUTINE condveg_soilalb(nbpt, lalo, neighbours, resolution, contfrac)
973
974	USE interpweight
975
976	IMPLICIT NONE
977
978
979	!! 0. Variable and parameter declaration
980
981	!! 0.1 Input variables
982
983	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
984	!! interpolated (unitless)
985	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
986	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
987	!! (1=N, 2=E, 3=S, 4=W)
988	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
989	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
990
991	!! 0.4 Local variables
992
993	CHARACTER(LEN=80) :: filename !! Filename of soil colour map
994	INTEGER(i_std) :: i, ib, ip, nbexp !! Indices
995	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
996	REAL(r_std), DIMENSION(nbpt) :: asoilcol !! Availability of the soilcol interpolation
997	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
998	!! (variabletypevals(1) = -un, not used)
999	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: soilcolrefrac !! soilcol fractions re-dimensioned
1000	REAL(r_std) :: vmin, vmax !! min/max values to use for the
1001	!! renormalization
1002	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
1003	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
1004	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
1005	!! 'XYKindTime': Input values are kinds
1006	!! of something with a temporal
1007	!! evolution on the dx*dy matrix'
1008	LOGICAL :: nonegative !! whether negative values should be removed
1009	CHARACTER(LEN=50) :: maskingtype !! Type of masking
1010	!! 'nomask': no-mask is applied
1011	!! 'mbelow': take values below maskvals(1)
1012	!! 'mabove': take values above maskvals(1)
1013	!! 'msumrange': take values within 2 ranges;
1014	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
1015	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
1016	!! (normalized by maskvals(3))
1017	!! 'var': mask values are taken from a
1018	!! variable inside the file (>0)
1019	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
1020	!! `maskingtype')
1021	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
1022	CHARACTER(LEN=250) :: msg
1023	INTEGER :: fopt
1024	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: vecpos
1025	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: solt
1026
1027	!_ ================================================================================================================================
1028	!! 1. Open file and allocate memory
1029
1030	! Open file with soil colours
1031
1032	!Config Key = SOILALB_FILE
1033	!Config Desc = Name of file from which the bare soil albedo
1034	!Config Def = soils_param.nc
1035	!Config If = NOT(IMPOSE_AZE)
1036	!Config Help = The name of the file to be opened to read the soil types from
1037	!Config which we derive then the bare soil albedos. This file is 1x1
1038	!Config deg and based on the soil colors defined by Wilson and Henderson-Seller.
1039	!Config Units = [FILE]
1040	!
1041	filename = 'soils_param.nc'
1042	CALL getin_p('SOILALB_FILE',filename)
1043
1044
1045	ALLOCATE(soilcolrefrac(nbpt, classnb), STAT=ALLOC_ERR)
1046	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable soilcolrefrac','','')
1047	ALLOCATE(vecpos(classnb), STAT=ALLOC_ERR)
1048	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable vecpos','','')
1049	ALLOCATE(solt(classnb), STAT=ALLOC_ERR)
1050	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable solt','','')
1051
1052	! Assigning values to vmin, vmax
1053	vmin = 1.0
1054	vmax = classnb
1055
1056	ALLOCATE(variabletypevals(classnb),STAT=ALLOC_ERR)
1057	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variabletypevals','','')
1058	variabletypevals = -un
1059
1060	!! Variables for interpweight
1061	! Type of calculation of cell fractions
1062	fractype = 'default'
1063	! Name of the longitude and latitude in the input file
1064	lonname = 'nav_lon'
1065	latname = 'nav_lat'
1066	! Should negative values be set to zero from input file?
1067	nonegative = .FALSE.
1068	! Type of mask to apply to the input data (see header for more details)
1069	maskingtype = 'mabove'
1070	! Values to use for the masking
1071	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
1072	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
1073	namemaskvar = ''
1074
1075	! Interpolate variable soilcolor
1076	variablename = 'soilcolor'
1077	IF (printlev_loc >= 1) WRITE(numout,*) "condveg_soilalb: Read and interpolate " &
1078	// TRIM(filename) // " for variable " // TRIM(variablename)
1079	CALL interpweight_2D(nbpt, classnb, variabletypevals, lalo, resolution, neighbours, &
1080	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1081	maskvals, namemaskvar, 0, 0, -1, fractype, &
1082	-1., -1., soilcolrefrac, asoilcol)
1083	IF (printlev_loc >= 5) WRITE(numout,*)' condveg_soilalb after interpweight_2D'
1084
1085	! Check how many points with soil information are found
1086	nbexp = 0
1087
1088	soilalb_dry(:,:) = zero
1089	soilalb_wet(:,:) = zero
1090	soilalb_moy(:,:) = zero
1091	IF (printlev_loc >= 5) THEN
1092	WRITE(numout,*)' condveg_soilalb before starting loop nbpt:', nbpt
1093	WRITE(numout,*)' condveg_soilalb initial values classnb: ',classnb
1094	WRITE(numout,*)' condveg_soilalb vis_dry. SUM:',SUM(vis_dry),' vis_dry= ',vis_dry
1095	WRITE(numout,*)' condveg_soilalb nir_dry. SUM:',SUM(nir_dry),' nir_dry= ',nir_dry
1096	WRITE(numout,*)' condveg_soilalb vis_wet. SUM:',SUM(vis_wet),' vis_wet= ',vis_wet
1097	WRITE(numout,*)' condveg_soilalb nir_wet. SUM:',SUM(nir_wet),' nir_wet= ',nir_wet
1098	END IF
1099
1100	DO ib=1,nbpt ! Loop over domain size
1101
1102	! vecpos: List of positions where textures were not zero
1103	! vecpos(1): number of not null textures found
1104	vecpos = interpweight_ValVecR(soilcolrefrac(ib,:),classnb,zero,'neq')
1105	fopt = vecpos(1)
1106	IF (fopt == classnb) THEN
1107	! All textures are not zero
1108	solt(:) = (/(i,i=1,classnb)/)
1109	ELSE IF (fopt == 0) THEN
1110	WRITE(numout,*)' condveg_soilalb: for point=', ib, ' no soil class!'
1111	ELSE
1112	DO ip = 1,fopt
1113	solt(ip) = vecpos(ip+1)
1114	END DO
1115	END IF
1116
1117	!! 3. Compute the average bare soil albedo parameters
1118
1119	IF ( (fopt .EQ. 0) .OR. (asoilcol(ib) .LT. min_sechiba)) THEN
1120	! Initialize with mean value if no points were interpolated or if no data was found
1121	nbexp = nbexp + 1
1122	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1123	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1124	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1125	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1126	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
1127	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
1128	ELSE
1129	! If points were interpolated
1130	DO ip=1, fopt
1131	IF ( solt(ip) .LE. classnb) THEN
1132	! Set to zero if the value is below min_sechiba
1133	IF (soilcolrefrac(ib,solt(ip)) < min_sechiba) soilcolrefrac(ib,solt(ip)) = zero
1134
1135	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(solt(ip))*soilcolrefrac(ib,solt(ip))
1136	soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(solt(ip))*soilcolrefrac(ib,solt(ip))
1137	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(solt(ip))*soilcolrefrac(ib,solt(ip))
1138	soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(solt(ip))*soilcolrefrac(ib,solt(ip))
1139	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(solt(ip))* &
1140	soilcolrefrac(ib,solt(ip))
1141	soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(solt(ip))* &
1142	soilcolrefrac(ib,solt(ip))
1143	ELSE
1144	msg = 'The file contains a soil color class which is incompatible with this program'
1145	CALL ipslerr_p(3,'condveg_soilalb',TRIM(msg),'','')
1146	ENDIF
1147	ENDDO
1148	ENDIF
1149
1150	ENDDO
1151
1152	IF ( nbexp .GT. 0 ) THEN
1153	WRITE(numout,*) 'condveg_soilalb _______'
1154	WRITE(numout,*) 'condveg_soilalb: The interpolation of the bare soil albedo had ', nbexp
1155	WRITE(numout,*) 'condveg_soilalb: points without data. This are either coastal points or'
1156	WRITE(numout,*) 'condveg_soilalb: ice covered land.'
1157	WRITE(numout,*) 'condveg_soilalb: The problem was solved by using the average of all soils'
1158	WRITE(numout,*) 'condveg_soilalb: in dry and wet conditions'
1159	WRITE(numout,*) 'condveg_soilalb: Use the diagnostic output field asoilcol to see location of these points'
1160	ENDIF
1161
1162	DEALLOCATE (soilcolrefrac)
1163	DEALLOCATE (variabletypevals)
1164
1165	! Write diagnostics
1166	CALL xios_orchidee_send_field("asoilcol",asoilcol)
1167
1168
1169	IF (printlev_loc >= 3) WRITE(numout,*)' condveg_soilalb ended'
1170
1171	END SUBROUTINE condveg_soilalb
1172
1173
1174	!! ==============================================================================================================================
1175	!! SUBROUTINE : condveg_background_soilalb
1176	!!
1177	!>\BRIEF This subroutine reads the albedo of bare soil
1178	!!
1179	!! DESCRIPTION This subroutine reads the background albedo map in 0.5 x 0.5 deg resolution
1180	!! derived from JRCTIP product to be used as bare soil albedo. These values are then interpolated
1181	!! to the model's resolution.\n
1182	!!
1183	!! RECENT CHANGE(S): None
1184	!!
1185	!! MAIN OUTPUT VARIABLE(S): soilalb_bg for visible and near-infrared range
1186	!!
1187	!! REFERENCES : None
1188	!!
1189	!! FLOWCHART : None
1190	!! \n
1191	!_ ================================================================================================================================
1192
1193	SUBROUTINE condveg_background_soilalb(nbpt, lalo, neighbours, resolution, contfrac)
1194
1195	USE interpweight
1196
1197	IMPLICIT NONE
1198
1199	!! 0. Variable and parameter declaration
1200
1201	!! 0.1 Input variables
1202
1203	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
1204	!! interpolated (unitless)
1205	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
1206	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
1207	!! (1=N, 2=E, 3=S, 4=W)
1208	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
1209	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
1210
1211	!! 0.4 Local variables
1212
1213	CHARACTER(LEN=80) :: filename !! Filename of background albedo
1214	REAL(r_std), DIMENSION(nbpt) :: aalb_bg !! Availability of the interpolation
1215	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! Latitudes and longitudes read from input file
1216	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel !! Help variable to read file data and allocate memory
1217	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu !! Help variable to read file data and allocate memory
1218	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
1219	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: soilalbedo_bg !! Help variable to read file data and allocate memory
1220	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
1221	REAL(r_std) :: vmin, vmax !! min/max values to use for the
1222	!! renormalization
1223	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
1224	CHARACTER(LEN=250) :: maskvname !! Variable to read the mask from
1225	!! the file
1226	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
1227	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
1228	!! 'XYKindTime': Input values are kinds
1229	!! of something with a temporal
1230	!! evolution on the dx*dy matrix'
1231	LOGICAL :: nonegative !! whether negative values should be removed
1232	CHARACTER(LEN=50) :: maskingtype !! Type of masking
1233	!! 'nomask': no-mask is applied
1234	!! 'mbelow': take values below maskvals(1)
1235	!! 'mabove': take values above maskvals(1)
1236	!! 'msumrange': take values within 2 ranges;
1237	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
1238	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
1239	!! (normalized by maskedvals(3))
1240	!! 'var': mask values are taken from a
1241	!! variable inside the file (>0)
1242	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
1243	!! `maskingtype')
1244	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
1245	REAL(r_std) :: albbg_norefinf !! No value
1246	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: albbg_default !! Default value
1247
1248	!_ ================================================================================================================================
1249
1250	!! 1. Open file and allocate memory
1251
1252	! Open file with background albedo
1253
1254	!Config Key = ALB_BG_FILE
1255	!Config Desc = Name of file from which the background albedo is read
1256	!Config Def = alb_bg.nc
1257	!Config If = ALB_BG_MODIS
1258	!Config Help = The name of the file to be opened to read background albedo
1259	!Config Units = [FILE]
1260	!
1261	filename = 'alb_bg.nc'
1262	CALL getin_p('ALB_BG_FILE',filename)
1263
1264	IF (xios_interpolation) THEN
1265
1266	! Read and interpolation background albedo using XIOS
1267	CALL xios_orchidee_recv_field('bg_alb_vis_interp',soilalb_bg(:,ivis))
1268	CALL xios_orchidee_recv_field('bg_alb_nir_interp',soilalb_bg(:,inir))
1269
1270	aalb_bg(:)=1
1271
1272	ELSE
1273	! Read background albedo file using IOIPSL and interpolate using aggregate standard method
1274
1275	ALLOCATE(albbg_default(2), STAT=ALLOC_ERR)
1276	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'condveg_background_soilalb','Pb in allocation for albbg_default','','')
1277
1278	! For this case there are not types/categories. We have 'only' a continuos field
1279	! Assigning values to vmin, vmax
1280	vmin = 0.
1281	vmax = 9999.
1282
1283	!! Variables for interpweight
1284	! Type of calculation of cell fractions (not used here)
1285	fractype = 'default'
1286	! Name of the longitude and latitude in the input file
1287	lonname = 'longitude'
1288	latname = 'latitude'
1289	! Default value when no value is get from input file
1290	albbg_default(ivis) = 0.129
1291	albbg_default(inir) = 0.247
1292	! Reference value when no value is get from input file (not used here)
1293	albbg_norefinf = undef_sechiba
1294	! Should negative values be set to zero from input file?
1295	nonegative = .FALSE.
1296	! Type of mask to apply to the input data (see header for more details)
1297	maskingtype = 'var'
1298	! Values to use for the masking (here not used)
1299	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
1300	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var')
1301	namemaskvar = 'mask'
1302
1303	! There is a variable for each chanel 'infrared' and 'visible'
1304	! Interpolate variable bg_alb_vis
1305	variablename = 'bg_alb_vis'
1306	IF (printlev_loc >= 2) WRITE(numout,*) "condveg_background_soilalb: Start interpolate " &
1307	// TRIM(filename) // " for variable " // TRIM(variablename)
1308	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
1309	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1310	maskvals, namemaskvar, -1, fractype, albbg_default(ivis), albbg_norefinf, &
1311	soilalb_bg(:,ivis), aalb_bg)
1312	IF (printlev_loc >= 5) WRITE(numout,*)" condveg_background_soilalb after InterpWeight2Dcont for '" // &
1313	TRIM(variablename) // "'"
1314
1315	! Interpolate variable bg_alb_nir in the same file
1316	variablename = 'bg_alb_nir'
1317	IF (printlev_loc >= 2) WRITE(numout,*) "condveg_background_soilalb: Start interpolate " &
1318	// TRIM(filename) // " for variable " // TRIM(variablename)
1319	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
1320	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1321	maskvals, namemaskvar, -1, fractype, albbg_default(inir), albbg_norefinf, &
1322	soilalb_bg(:,inir), aalb_bg)
1323	IF (printlev_loc >= 5) WRITE(numout,*)" condveg_background_soilalb after InterpWeight2Dcont for '" // &
1324	TRIM(variablename) // "'"
1325
1326	IF (ALLOCATED(albbg_default)) DEALLOCATE(albbg_default)
1327
1328	IF (printlev_loc >= 3) WRITE(numout,*)' condveg_background_soilalb ended'
1329
1330	ENDIF
1331
1332	CALL xios_orchidee_send_field("interp_diag_alb_vis",soilalb_bg(:,ivis))
1333	CALL xios_orchidee_send_field("interp_diag_alb_nir",soilalb_bg(:,inir))
1334	CALL xios_orchidee_send_field("aalb_bg",aalb_bg)
1335
1336	END SUBROUTINE condveg_background_soilalb
1337
1338
1339	!! ==============================================================================================================================
1340	!! SUBROUTINE : condveg_z0cdrag
1341	!!
1342	!>\BRIEF Computation of grid average of roughness length by calculating
1343	!! the drag coefficient.
1344	!!
1345	!! DESCRIPTION : This routine calculates the mean roughness height and mean
1346	!! effective roughness height over the grid cell. The mean roughness height (z0)
1347	!! is computed by averaging the drag coefficients \n
1348	!!
1349	!! \latexonly
1350	!! \input{z0cdrag1.tex}
1351	!! \endlatexonly
1352	!! \n
1353	!!
1354	!! where C is the drag coefficient at the height of the vegetation, kappa is the
1355	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
1356	!! The reference level for z needs to be high enough above the canopy to avoid
1357	!! singularities of the LOG. This height is set to minimum 10m above ground.
1358	!! The drag coefficient increases with roughness height to represent the greater
1359	!! turbulence generated by rougher surfaces.
1360	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
1361	!!
1362	!! The roughness height for the non-vegetative surfaces is calculated in a second step.
1363	!! In order to calculate the transfer coefficients the
1364	!! effective roughness height is calculated. This effective value is the difference
1365	!! between the height of the vegetation and the zero plane displacement height.\nn
1366	!!
1367	!! RECENT CHANGE(S): None
1368	!!
1369	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
1370	!!
1371	!! REFERENCE(S) : None
1372	!!
1373	!! FLOWCHART : None
1374	!! \n
1375	!_ ================================================================================================================================
1376
1377	SUBROUTINE condveg_z0cdrag (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, tot_bare_soil, frac_snow_veg, &
1378	& z0m, z0h, roughheight)
1379
1380	!! 0. Variable and parameter declaration
1381
1382	!! 0.1 Input variables
1383
1384	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1385	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1386	!! (m^2 m^{-2})
1387	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1388	!! (= ind*cn_ind) (m^2 m^{-2})
1389	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1390	!! i.e. continental ice, lakes, etc. (unitless)
1391	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1392	!! i.e. continental ice, lakes, etc. (unitless)
1393	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
1394	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1395	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
1396	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1397
1398	!! 0.2 Output variables
1399
1400	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0m !! Roughness height for momentum (m)
1401	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0h !! Roughness height for heat (m)
1402	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1403
1404	!! 0.3 Modified variables
1405
1406	!! 0.4 Local variables
1407
1408	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1409	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1410	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
1411	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1412	REAL(r_std), DIMENSION(kjpindex) :: d_veg !! PFT coverage of vegetative PFTs
1413	!! (= ind*cn_ind) (m^2 m^{-2})
1414	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1415	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
1416	!! i.e. continental ice, lakes, etc.
1417	REAL(r_std), DIMENSION(kjpindex) :: dragm !! Drag coefficient for momentum
1418	REAL(r_std), DIMENSION(kjpindex) :: dragh !! Drag coefficient for heat
1419	REAL(r_std), DIMENSION(kjpindex) :: z0_ground !! z0m value used for ground surface
1420	!_ ================================================================================================================================
1421
1422	!! 1. Preliminary calculation
1423
1424	! Set maximal height of first layer
1425	ztmp(:) = MAX(10., zlev(:))
1426
1427	z0_ground(:) = (1.-frac_snow_veg(:))z0_bare + frac_snow_veg(:)z0_bare/10.
1428
1429	! Calculate roughness for non-vegetative surfaces
1430	! with the von Karman constant
1431	dragm(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/z0_ground))**2
1432	dragh(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/(z0_ground/ratio_z0m_z0h(1))))*(ct_karman/LOG(ztmp(:)/z0_ground))
1433	! Fraction of bare soil
1434	sumveg(:) = tot_bare_soil(:)
1435
1436	! Set average vegetation height to zero
1437	ave_height(:) = zero
1438
1439	!! 2. Calculate the mean roughness height
1440
1441	! Calculate the mean roughness height of
1442	! vegetative PFTs over the grid cell
1443	DO jv = 2, nvm
1444
1445	! In the case of forest, use parameter veget_max because
1446	! tree trunks influence the roughness even when there are no leaves
1447	IF ( is_tree(jv) ) THEN
1448	! In the case of grass, use parameter veget because grasses
1449	! only influence the roughness during the growing season
1450	d_veg(:) = veget_max(:,jv)
1451	ELSE
1452	! grasses only have an influence if they are really there!
1453	d_veg(:) = veget(:,jv)
1454	ENDIF
1455
1456	! Calculate the average roughness over the grid cell:
1457	! The unitless drag coefficient is per vegetative PFT
1458	! calculated by use of the von Karman constant, the height
1459	! of the first layer and the roughness. The roughness
1460	! is calculated as the vegetation height per PFT
1461	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
1462	! If this scaled value is lower than 0.01 then the value for
1463	! the roughness of bare soil (0.01) is used.
1464	! The sum over all PFTs gives the average roughness
1465	! per grid cell for the vegetative PFTs.
1466	dragm(:) = dragm(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)z0_over_height(jv),z0_ground)))*2
1467	dragh(:) = dragh(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/(MAX(height(:,jv)*z0_over_height(jv),z0_ground) / &
1468	ratio_z0m_z0h(jv)))) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)*z0_over_height(jv),z0_ground)))
1469
1470	! Sum of bare soil and fraction vegetated fraction
1471	sumveg(:) = sumveg(:) + d_veg(:)
1472
1473	! Weigh height of vegetation with maximal cover fraction
1474	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1475
1476	ENDDO
1477
1478	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1479
1480	! Search for pixels with vegetated part to normalise
1481	! roughness height
1482	WHERE ( sumveg(:) .GT. min_sechiba )
1483	dragm(:) = dragm(:) / sumveg(:)
1484	dragh(:) = dragh(:) / sumveg(:)
1485	ENDWHERE
1486	! Calculate fraction of roughness for vegetated part
1487	dragm(:) = (un - totfrac_nobio(:)) * dragm(:)
1488	dragh(:) = (un - totfrac_nobio(:)) * dragh(:)
1489
1490	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1491
1492	! Set rougness for ice
1493	IF ( jv .EQ. iice ) THEN
1494	z0_nobio = z0_ice
1495	ELSE
1496	WRITE(numout,*) 'jv=',jv
1497	WRITE(numout,*) 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1498	CALL ipslerr_p(3,'condveg_z0cdrag','DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1499	ENDIF
1500
1501	! Sum of vegetative roughness length and non-vegetative
1502	! roughness length
1503	dragm(:) = dragm(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1504	dragh(:) = dragh(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio/ratio_z0m_z0h(1)))*(ct_karman/LOG(ztmp(:)/z0_nobio))
1505
1506	ENDDO ! Loop over # of non-vegative surfaces
1507
1508	!! 4. Calculate the zero plane displacement height and effective roughness length
1509
1510	! Take the exponential of the roughness
1511	z0m(:) = ztmp(:) / EXP(ct_karman/SQRT(dragm(:)))
1512	z0h(:) = ztmp(:) / EXP((ct_karman*2.)/(dragh(:)LOG(ztmp(:)/z0m(:))))
1513
1514	! Compute the zero plane displacement height which
1515	! is an equivalent height for the absorption of momentum
1516	zhdispl(:) = ave_height(:) * height_displacement
1517
1518	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1519	! This is the height over which the roughness acts. It combines the
1520	! zero plane displacement height and the vegetation height.
1521	roughheight(:) = ave_height(:) - zhdispl(:)
1522
1523	END SUBROUTINE condveg_z0cdrag
1524
1525
1526	!! ==============================================================================================================================
1527	!! SUBROUTINE : condveg_z0cdrag_dyn
1528	!!
1529	!>\BRIEF Computation of grid average of roughness length by calculating
1530	!! the drag coefficient based on formulation proposed by Su et al. (2001).
1531	!!
1532	!! DESCRIPTION : This routine calculates the mean roughness height and mean
1533	!! effective roughness height over the grid cell. The mean roughness height (z0)
1534	!! is computed by averaging the drag coefficients \n
1535	!!
1536	!! \latexonly
1537	!! \input{z0cdrag1.tex}
1538	!! \endlatexonly
1539	!! \n
1540	!!
1541	!! where C is the drag coefficient at the height of the vegetation, kappa is the
1542	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
1543	!! The reference level for z needs to be high enough above the canopy to avoid
1544	!! singularities of the LOG. This height is set to minimum 10m above ground.
1545	!! The drag coefficient increases with roughness height to represent the greater
1546	!! turbulence generated by rougher surfaces.
1547	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
1548	!! In the formulation of Su et al. (2001), one distinguishes the roughness height for
1549	!! momentum (z0m) and the one for heat (z0h).
1550	!! z0m is computed as a function of LAI (z0m increases with LAI) and z0h is computed
1551	!! with a so-called kB-1 term (z0m/z0h=exp(kB-1))
1552	!!
1553	!! RECENT CHANGE(S): Written by N. Vuichard (2016)
1554	!!
1555	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
1556	!!
1557	!! REFERENCE(S) :
1558	!! - Su, Z., Schmugge, T., Kustas, W.P., Massman, W.J., 2001. An Evaluation of Two Models for
1559	!! Estimation of the Roughness Height for Heat Transfer between the Land Surface and the Atmosphere. J. Appl.
1560	!! Meteorol. 40, 1933â1951. doi:10.1175/1520-0450(2001)
1561	!! - Ershadi, A., McCabe, M.F., Evans, J.P., Wood, E.F., 2015. Impact of model structure and parameterization
1562	!! on Penman-Monteith type evaporation models. J. Hydrol. 525, 521â535. doi:10.1016/j.jhydrol.2015.04.008
1563	!!
1564	!! FLOWCHART : None
1565	!! \n
1566	!_ ================================================================================================================================
1567
1568	SUBROUTINE condveg_z0cdrag_dyn (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, &
1569	& temp_air, pb, u, v, lai, frac_snow_veg, z0m, z0h, roughheight)
1570
1571	!! 0. Variable and parameter declaration
1572
1573	!! 0.1 Input variables
1574
1575	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1576	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1577	!! (m^2 m^{-2})
1578	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1579	!! (= ind*cn_ind) (m^2 m^{-2})
1580	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1581	!! i.e. continental ice, lakes, etc. (unitless)
1582	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1583	!! i.e. continental ice, lakes, etc. (unitless)
1584	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
1585	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1586	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: temp_air !! 2m air temperature (K)
1587	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
1588	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
1589	!! @tex $(m.s^{-1})$ @endtex
1590	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
1591	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
1592	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1593	!! 0.2 Output variables
1594
1595	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0m !! Roughness height for momentum (m)
1596	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0h !! Roughness height for heat (m)
1597	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1598
1599	!! 0.3 Modified variables
1600
1601	!! 0.4 Local variables
1602
1603	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1604	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1605	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
1606	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1607	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1608	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
1609	!! i.e. continental ice, lakes, etc.
1610	REAL(r_std), DIMENSION(kjpindex) :: z0m_pft !! Roughness height for momentum for a specific PFT
1611	REAL(r_std), DIMENSION(kjpindex) :: z0h_pft !! Roughness height for heat for a specific PFT
1612	REAL(r_std), DIMENSION(kjpindex) :: dragm !! Drag coefficient for momentum
1613	REAL(r_std), DIMENSION(kjpindex) :: dragh !! Drag coefficient for heat
1614	REAL(r_std), DIMENSION(kjpindex) :: eta !! Ratio of friction velocity to the wind speed at the canopy top - See Ershadi et al. (2015)
1615	REAL(r_std), DIMENSION(kjpindex) :: eta_ec !! Within-canopy wind speed profile estimation coefficient - See Ershadi et al. (2015)
1616	REAL(r_std), DIMENSION(kjpindex) :: Ct_star !! Heat transfer coefficient of the soil - see Su et al. (2001)
1617	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))
1618	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))
1619	REAL(r_std), DIMENSION(kjpindex) :: fc !! fractional canopy coverage
1620	REAL(r_std), DIMENSION(kjpindex) :: fs !! fractional soil coverage
1621	REAL(r_std), DIMENSION(kjpindex) :: Reynolds !! Reynolds number
1622	REAL(r_std), DIMENSION(kjpindex) :: wind !! wind Speed (m)
1623	REAL(r_std), DIMENSION(kjpindex) :: u_star !! friction velocity
1624	REAL(r_std), DIMENSION(kjpindex) :: z0_ground !! z0m value used for ground surface
1625	!_ ================================================================================================================================
1626
1627	!! 1. Preliminary calculation
1628
1629	! Set maximal height of first layer
1630	ztmp(:) = MAX(10., zlev(:))
1631
1632	z0_ground(:) = (1.-frac_snow_veg(:))z0_bare + frac_snow_veg(:)z0_bare/10.
1633
1634	! Calculate roughness for non-vegetative surfaces
1635	! with the von Karman constant
1636	dragm(:) = veget_max(:,1) * (ct_karman/LOG(ztmp(:)/z0_ground(:)))**2
1637
1638	wind(:) = SQRT(u(:)u(:)+v(:)v(:))
1639	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG(zlev(:)/z0_ground(:))
1640	Reynolds(:) = z0_ground(:) * u_star(:) &
1641	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1642
1643	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1644
1645	dragh(:) = veget_max(:,1) * (ct_karman/LOG(ztmp(:)/z0_ground(:)))*(ct_karman/LOG(ztmp(:)/(z0_ground(:)/ exp(kBs_m1(:))) ))
1646
1647	! Fraction of bare soil
1648	sumveg(:) = veget_max(:,1)
1649
1650	! Set average vegetation height to zero
1651	ave_height(:) = zero
1652
1653	!! 2. Calculate the mean roughness height
1654
1655	! Calculate the mean roughness height of
1656	! vegetative PFTs over the grid cell
1657	DO jv = 2, nvm
1658
1659	WHERE(veget_max(:,jv) .GT. zero)
1660	! Calculate the average roughness over the grid cell:
1661	! The unitless drag coefficient is per vegetative PFT
1662	! calculated by use of the von Karman constant, the height
1663	! of the first layer and the roughness. The roughness
1664	! is calculated as the vegetation height per PFT
1665	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
1666	! If this scaled value is lower than 0.01 then the value for
1667	! the roughness of bare soil (0.01) is used.
1668	! The sum over all PFTs gives the average roughness
1669	! per grid cell for the vegetative PFTs.
1670	eta(:) = c1 - c2 * exp(-c3 * Cdrag_foliage * lai(:,jv))
1671
1672	z0m_pft(:) = (height(:,jv)(1-height_displacement)(exp(-ct_karman/eta(:))-exp(-ct_karman/(c1-c2)))) &
1673	+ z0_ground(:)
1674
1675	dragm(:) = dragm(:) + veget_max(:,jv) * (ct_karman/LOG(ztmp(:)/z0m_pft(:)))**2
1676
1677	fc(:) = veget(:,jv)/veget_max(:,jv)
1678	fs(:) = 1. - fc(:)
1679
1680	eta_ec(:) = ( Cdrag_foliage * lai(:,jv)) / (2 * eta(:)*eta(:))
1681	wind(:) = SQRT(u(:)u(:)+v(:)v(:))
1682	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG((zlev(:)+(height(:,jv)*(1-height_displacement)))/z0m_pft(:))
1683	Reynolds(:) = z0_ground(:) * u_star(:) &
1684	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1685
1686	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1687	Ct_star(:) = Prandtl*(-2./3.) SQRT(1./Reynolds(:))
1688
1689	WHERE(lai(:,jv) .GT. min_sechiba)
1690	kB_m1(:) = (ct_karman * Cdrag_foliage) / (4 * Ct * eta(:) * (1 - exp(-eta_ec(:)/2.))) * fc(:)**2. &
1691	+ 2fc(:)fs(:) * (ct_karman * eta(:) * z0m_pft(:) / height(:,jv)) / Ct_star(:) &
1692	+ kBs_m1(:) * fs(:)**2.
1693	ELSEWHERE
1694	kB_m1(:) = kBs_m1(:) * fs(:)**2.
1695	ENDWHERE
1696
1697	z0h_pft(:) = z0m_pft(:) / exp(kB_m1(:))
1698
1699	dragh(:) = dragh(:) + veget_max(:,jv) * (ct_karman/LOG(ztmp(:)/z0m_pft(:)))*(ct_karman/LOG(ztmp(:)/z0h_pft(:)))
1700
1701	! Sum of bare soil and fraction vegetated fraction
1702	sumveg(:) = sumveg(:) + veget_max(:,jv)
1703
1704	! Weigh height of vegetation with maximal cover fraction
1705	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1706
1707	ENDWHERE
1708	ENDDO
1709
1710	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1711
1712	! Search for pixels with vegetated part to normalise
1713	! roughness height
1714	WHERE ( sumveg(:) .GT. min_sechiba )
1715	dragh(:) = dragh(:) / sumveg(:)
1716	dragm(:) = dragm(:) / sumveg(:)
1717	ENDWHERE
1718
1719	! Calculate fraction of roughness for vegetated part
1720	dragh(:) = (un - totfrac_nobio(:)) * dragh(:)
1721	dragm(:) = (un - totfrac_nobio(:)) * dragm(:)
1722
1723	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1724
1725	! Set rougness for ice
1726	IF ( jv .EQ. iice ) THEN
1727	z0_nobio = z0_ice
1728	ELSE
1729	WRITE(numout,*) 'jv=',jv
1730	WRITE(numout,*) 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1731	CALL ipslerr_p(3,'condveg_z0cdrag_dyn','DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1732	ENDIF
1733
1734	! Sum of vegetative roughness length and non-vegetative roughness length
1735	! Note that z0m could be made dependent of frac_snow_nobio
1736	dragm(:) = dragm(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1737
1738	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG(zlev(:)/z0_nobio)
1739	Reynolds(:) = z0_nobio * u_star(:) &
1740	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1741
1742	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1743
1744	dragh(:) = dragh(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio)) * &
1745	(ct_karman/LOG(ztmp(:)/(z0_nobio/ exp(kBs_m1(:))) ))
1746	ENDDO ! Loop over # of non-vegative surfaces
1747
1748	!! 4. Calculate the zero plane displacement height and effective roughness length
1749	! Take the exponential of the roughness
1750	z0m(:) = ztmp(:) / EXP(ct_karman/SQRT(dragm(:)))
1751	z0h(:) = ztmp(:) / EXP((ct_karman*2.)/(dragh(:)LOG(ztmp(:)/z0m(:))))
1752
1753	! Compute the zero plane displacement height which
1754	! is an equivalent height for the absorption of momentum
1755	zhdispl(:) = ave_height(:) * height_displacement
1756
1757	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1758	! This is the height over which the roughness acts. It combines the
1759	! zero plane displacement height and the vegetation height.
1760	roughheight(:) = ave_height(:) - zhdispl(:)
1761
1762	END SUBROUTINE condveg_z0cdrag_dyn
1763
1764
1765	END MODULE condveg

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

Download in other formats: