Context Navigation

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

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

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

Download in other formats: