Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/condveg.f90 @ 7266

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

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

Download in other formats: