Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_sechiba/condveg.f90 @ 7704

Last change on this file since 7704 was 7615, checked in by sebastiaan.luyssaert, 2 years ago
Enhanced consistency of variable names: input has been changed in n_input throughout the code and the variable name vegstress introduced in sechiba is now also used in stomate. Enhnaced computational consistency: Pgap_cumul is used in stomate rather than recalculating it before calculating light_tran_to_floor_season. Edited getin_p while checking the code (but no real changes were made) and added several missing stomate and sechiba variables to age_class_distr to improve the 1+1=2 issue when comparing a model run with against a run without age classes. Finally: Write warning 10b in allocation to the history file rather than the out_orchidee file
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 58.5 KB

Line
1	! ===============================================================================================================================
2	! MODULE : condveg
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Initialise, compute and update the surface parameters emissivity and roughness.
10	!!
11	!! \n DESCRIPTION : This module calculates the emissivity and roughness. It calls albedo_surface_main for the albedo calculations.
12	!! The module uses 2 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
15	!! calculated by the old formulation which does not distinguish between
16	!! z0m and z0h and which does not vary with LAI. If set to true: the
17	!! grid average is calculated by the formulation proposed by Su et al. (2001)
18	!! 2. :: impaze for choosing surface parameters. If set to false, the
19	!! values for the soil albedo, emissivity and roughness height are
20	!! set to default values which are read from the run.def. If set to
21	!! true, the user imposes its own values, fixed for the grid point.
22	!! This is useful if one performs site simulations, however,
23	!! it is not recommended to do so for spatialized simulations.
24	!! roughheight_scal imposes the roughness height in (m) ,
25	!! same for emis_scal (in %), albedo_scal (in %), zo_scal (in m)
26	!! Note that these values are only used if 'impaze' is true.\n
27	!!
28	!! The surface fluxes are calculated between two levels: the
29	!! atmospheric level reference and the effective roughness height
30	!! defined as the difference between the mean height of the vegetation
31	!! and the displacement height (zero wind level). Over bare soils, the
32	!! zero wind level is equal to the soil roughness. Over vegetation,
33	!! the zero wind level is increased by the displacement height
34	!! which depends on the height of the vegetation. For a grid point
35	!! composed of different types of vegetation, an effective surface
36	!! roughness has to be calculated
37	!!
38	!! RECENT CHANGE(S): Added option rough_dyn and subroutine condveg_z0cdrag_dyn. Removed
39	!! subroutine condveg_z0logz. June 2016.
40	!! The condveg_albedo scheme has been replaced by albedo_surface_main,
41	!! developed for the DOFOCO branch. The albedo_surface_main is contained
42	!! in albedo.f90 module.
43	!! REFERENCES(S) : None
44	!!
45	!! SVN :
46	!! $HeadURL$
47	!! $Date$
48	!! $Revision$
49	!! \n
50	!_ ================================================================================================================================
51
52	MODULE condveg
53
54	USE ioipsl
55	USE xios_orchidee
56	USE constantes
57	USE constantes_soil
58	USE pft_parameters
59	USE qsat_moisture
60	USE interpol_help
61	USE mod_orchidee_para
62	USE ioipsl_para
63	USE sechiba_io_p
64	USE albedo_surface
65	USE structures
66	USE grid
67
68	IMPLICIT NONE
69
70	PRIVATE
71	PUBLIC :: condveg_xios_initialize, condveg_main, condveg_initialize, condveg_finalize, condveg_clear
72
73	!
74	! Variables used inside condveg module
75	!
76	LOGICAL, SAVE :: l_first_condveg=.TRUE. !! To keep first call's trace
77	!$OMP THREADPRIVATE(l_first_condveg)
78
79	INTEGER, SAVE :: printlev_loc !! Output debug level
80	!$OMP THREADPRIVATE(printlev_loc)
81
82	CONTAINS
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	CALL albedo_surface_xios_initialize()
101
102	END SUBROUTINE condveg_xios_initialize
103
104	!! =============================================================================================================================
105	!! SUBROUTINE : condveg_initialize
106	!!
107	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
108	!!
109	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
110	!! condveg_snow is called to initialize corresponding variables.
111	!!
112	!! RECENT CHANGE(S) : None
113	!!
114	!! MAIN OUTPUT VARIABLE(S)
115	!!
116	!! REFERENCE(S) : None
117	!!
118	!! FLOWCHART : None
119	!! \n
120	!_ ==============================================================================================================================
121	!+++CHECK+++
122	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
123	! when running a larger domain. Needs to be corrected when implementing
124	! a global use of the multi-layer energy budget.
125	SUBROUTINE condveg_initialize (kjit, kjpindex, index, rest_id, &
126	lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
127	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
128	drysoil_frac, height, height_dom, snowdz, snowrho, tot_bare_soil, &
129	temp_air, pb, u, v, &
130	lai, &
131	emis, albedo, z0m, z0h, roughheight, &
132	frac_snow_veg,frac_snow_nobio, &
133	coszang, &
134	Light_Abs_Tot, Light_Tran_Tot, laieff_fit, &
135	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
136	laieff_isotrop)
137	!+++++++++++
138
139	!! 0. Variable and parameter declaration
140	!! 0.1 Input variables
141	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
142	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
143	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
144	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
145	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
146	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours!! neighoring grid points if land
147	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
148	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
149	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
150	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
151	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
152	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
153	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
154	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
155	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
156	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
157	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
158	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
159	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
160	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height_dom !! Dominant vegetation height (m)
161	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowdz !! Snow depth at each snow layer (m)
162	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowrho !! Snow density at each snow layer (Kg/m^3)
163	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
164	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
165	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
166	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u !! Horizontal wind speed, u direction
167	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: v !! Horizontal wind speed, v direction
168	!! @tex $(gC m^{-2})$ @endtex
169	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
170	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: coszang !! the cosine of the solar zenith angle (unitless)
171	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
172
173	!! 0.2 Output variables
174	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
175	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
176	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m !! Roughness for momentum (m)
177	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h !! Roughness for heat (m)
178	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
179	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: frac_snow_veg !! Snow cover fraction on vegeted area
180	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Snow cover fraction on non-vegeted area
181	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), &
182	INTENT(out) :: Light_Abs_Tot !! Absorbed radiation per layer, averaged between light
183	!! from a collimated (direct) source and that from an
184	!! isotropic (diffuse) source. Expressed as the fraction
185	!! of overall light hitting the canopy.
186	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), &
187	INTENT(out) :: Light_Tran_Tot !! Transmitted radiation per layer, averaged between light
188	!! from a collimated (direct) source and that from an
189	!! isotropic (diffuse) source. Expressed as the fraction
190	!! of overall light hitting the canopy.
191	!+++CHECK+++
192	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
193	! when running a larger domain. Needs to be corrected when implementing
194	! a global use of the multi-layer energy budget.
195	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT (out) :: Light_Abs_Tot_mean !! total absorption for a given level
196	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT (out) :: Light_Alb_Tot_mean !! total albedo for a given level
197	!+++++++++++
198	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: laieff_isotrop !! Leaf Area Index Effective converts
199	!! 3D lai inot 1D lai for two steam
200	!! radiation transfer model...this is for
201	!! isotropic light and only calculated once per day
202	!! @tex $(m^{2} m^{-2})$ @endtex
203
204	!! 0.4 Local variables
205	INTEGER :: ier
206	INTEGER :: ji,jv
207
208	!_ ================================================================================================================================
209
210	IF (.NOT. l_first_condveg) CALL ipslerr_p(3,'condveg_initialize','Error: initialization already done','','')
211	l_first_condveg=.FALSE.
212
213	!! Initialize local printlev
214	printlev_loc=get_printlev('condveg')
215
216
217	IF (printlev>=3) WRITE (numout,*) 'Start condveg_initialize'
218
219	!! 1. Allocate module variables and read from restart or initialize
220
221	! z0m
222	CALL ioconf_setatt_p('UNITS', '-')
223	CALL ioconf_setatt_p('LONG_NAME','Roughness for momentum')
224	CALL restget_p (rest_id, 'z0m', nbp_glo, 1, 1, kjit, .TRUE., z0m, "gather", nbp_glo, index_g)
225
226	! z0h
227	CALL ioconf_setatt_p('UNITS', '-')
228	CALL ioconf_setatt_p('LONG_NAME','Roughness for heat')
229	CALL restget_p (rest_id, 'z0h', nbp_glo, 1, 1, kjit, .TRUE., z0h, "gather", nbp_glo, index_g)
230
231	! roughness height
232	CALL ioconf_setatt_p('UNITS', '-')
233	CALL ioconf_setatt_p('LONG_NAME','Roughness height')
234	CALL restget_p (rest_id, 'roughheight', nbp_glo, 1, 1, kjit, .TRUE., roughheight, "gather", nbp_glo, index_g)
235
236	!! Initialize emissivity
237	IF ( impaze ) THEN
238	! Use parameter CONDVEG_EMIS from run.def
239	emis(:) = emis_scal
240	ELSE
241	! Set emissivity to 1.
242	emis_scal = un
243	emis(:) = emis_scal
244	ENDIF
245
246
247	!! 3. Calculate the fraction of snow on vegetation and nobio
248	CALL condveg_frac_snow(kjpindex, snow_nobio, snowrho, snowdz, &
249	frac_snow_veg, frac_snow_nobio)
250
251	!! 4. Calculate roughness height if it was not found in the restart file
252	IF ( ALL(z0m(:) == val_exp) .OR. ALL(z0h(:) == val_exp) .OR. ALL(roughheight(:) == val_exp)) THEN
253	!! Calculate roughness height
254	! Chooses between two methods to calculate the grid average of the
255	! roughness. If impaze set to true: The grid average is calculated by
256	! averaging the drag coefficients over PFT. If impaze set to false:
257	! The grid average is calculated by averaging the logarithm of the
258	! roughness length per PFT.
259	IF ( impaze ) THEN
260	! Use parameter CONDVEG_Z0 and ROUGHHEIGHT from run.def
261	z0m(:) = z0_scal
262	z0h(:) = z0_scal
263	roughheight(:) = roughheight_scal
264	ELSE
265	! Caluculate roughness height
266	IF( rough_dyn ) THEN
267	CALL condveg_z0cdrag_dyn(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
268	& height, height_dom, temp_air, pb, u, v, lai, frac_snow_veg, z0m, z0h, roughheight)
269	ELSE
270	CALL condveg_z0cdrag(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
271	height, height_dom, tot_bare_soil, frac_snow_veg, z0m, z0h, roughheight)
272	ENDIF
273	END IF
274	END IF
275
276
277	!! 5. Initialze albedo module
278	!+++CHECK+++
279	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
280	! when running a larger domain. Needs to be corrected when implementing
281	! a global use of the multi-layer energy budget.
282	CALL albedo_surface_initialize(kjit, kjpindex, rest_id, &
283	lalo, neighbours, resolution, contfrac, &
284	drysoil_frac, veget_max, coszang, frac_nobio, &
285	snow, snow_age, snow_nobio, &
286	snow_nobio_age, frac_snow_veg, frac_snow_nobio,&
287	z0m, laieff_fit, &
288	albedo, &
289	Light_Abs_Tot, Light_Tran_Tot, &
290	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
291	laieff_isotrop, veget)
292	!+++++++++++
293
294	IF (printlev>=3) WRITE (numout,*) 'condveg_initialize done '
295
296	END SUBROUTINE condveg_initialize
297
298
299
300	!! ==============================================================================================================================
301	!! SUBROUTINE : condveg_main
302	!!
303	!>\BRIEF Calls the subroutines update the variables for current time step
304	!!
305	!!
306	!! MAIN OUTPUT VARIABLE(S): emis (emissivity), albedo (albedo of
307	!! vegetative PFTs in visible and near-infrared range), z0 (surface roughness height),
308	!! roughheight (grid effective roughness height), soil type (fraction of soil types)
309	!!
310	!!
311	!! REFERENCE(S) : None
312	!!
313	!! FLOWCHART : None
314	!!
315	!! REVISION(S) : None
316	!!
317	!_ ================================================================================================================================
318
319	!+++CHECK+++
320	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
321	! when running a larger domain. Needs to be corrected when implementing
322	! a global use of the multi-layer energy budget.
323	SUBROUTINE condveg_main (kjit, kjpindex, index, rest_id, hist_id, hist2_id, &
324	lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
325	zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
326	drysoil_frac, height, height_dom, snowdz, snowrho, tot_bare_soil, &
327	Temp_air, pb, u, v, lai,&
328	emis, albedo, z0m, z0h, roughheight, &
329	frac_snow_veg, frac_snow_nobio, coszang, &
330	Light_Abs_Tot, Light_Tran_Tot, laieff_fit, &
331	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
332	laieff_isotrop )
333	!+++++++++++
334
335	!! 0. Variable and parameter declaration
336
337	!! 0.1 Input variables
338
339	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
340	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
341	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
342	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier
343	INTEGER(i_std), OPTIONAL, INTENT (in) :: hist2_id !! _History_ file 2 identifier
344	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
345	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
346	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours!! neighoring grid points if land
347	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
348	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
349	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
350	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
351	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
352	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
353	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
354	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
355	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
356	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
357	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
358	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
359	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
360	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height_dom !! Dominant vegetation height (m)
361	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowdz !! Snow depth at each snow layer (m)
362	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowrho !! Snow density at each snow layer (Kg/m^3)
363	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
364	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
365	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
366	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u !! Horizontal wind speed, u direction
367	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: v !! Horizontal wind speed, v direction
368	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: coszang !! the cosine of the solar enith angle (unitless)
369	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
370	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
371
372
373	!! 0.2 Output variables
374
375	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
376	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
377	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m !! Roughness for momentum (m)
378	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h !! Roughness for heat (m)
379	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
380	REAL(r_std),DIMENSION (:,:,:), &
381	INTENT (out) :: Light_Abs_Tot !!Absorbed radiation per level for photosynthesis
382	REAL(r_std),DIMENSION (:,:,:), &
383	INTENT (out) :: Light_Tran_Tot !!Transmitted radiation per level for photosynthesis
384	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: frac_snow_veg !! Snow cover fraction on vegeted area
385	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Snow cover fraction on non-vegeted area
386
387	!! 0.3 Modified variables
388	!+++CHECK+++
389	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
390	! when running a larger domain. Needs to be corrected when implementing
391	! a global use of the multi-layer energy budget.
392	!!$ REAL(r_std),DIMENSION(:), INTENT(inout) :: Light_Abs_Tot_mean !! total light absorption for a given level
393	!!$ REAL(r_std),DIMENSION(:), INTENT(inout) :: Light_Alb_Tot_mean !! total albedo for a given level
394	!+++++++++++
395	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: laieff_isotrop !! Leaf Area Index Effective converts
396	!! 3D lai into 1D lai for two stream
397	!! radiation transfer model...this is for
398	!! isotropic light and only calculated once per day
399	!! @tex $(m^{2} m^{-2})$ @endtex
400
401	!! 0.4 Local variables
402	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
403	INTEGER(i_std) :: ji,ilevel,ivm
404
405	!_ ================================================================================================================================
406
407	!! 1. Calculate the fraction of snow on vegetation and nobio
408	CALL condveg_frac_snow(kjpindex, snow_nobio, snowrho, snowdz, &
409	frac_snow_veg, frac_snow_nobio)
410
411	!! 2. Calculate emissivity
412	emis(:) = emis_scal
413
414	!! 3. Calculate roughness height
415	! If TRUE read in prescribed values for roughness height
416	IF ( impaze ) THEN
417
418	DO ji = 1, kjpindex
419	z0m(ji) = z0_scal
420	z0h(ji) = z0_scal
421	roughheight(ji) = roughheight_scal
422	ENDDO
423
424	ELSE
425
426	! Calculate roughness height
427	IF ( rough_dyn ) THEN
428	CALL condveg_z0cdrag_dyn (kjpindex, veget, veget_max, &
429	frac_nobio, totfrac_nobio, zlev, height, height_dom, &
430	temp_air, pb, u, v, lai, frac_snow_veg, &
431	z0m, z0h, roughheight)
432	ELSE
433	CALL condveg_z0cdrag (kjpindex, veget, veget_max, &
434	frac_nobio, totfrac_nobio, zlev, &
435	height, height_dom, tot_bare_soil, frac_snow_veg, z0m, z0h, roughheight)
436	ENDIF
437
438	ENDIF ! impaze
439
440
441	!! 4. Calculate albedo
442	!+++CHECK+++
443	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
444	! when running a larger domain. Needs to be corrected when implementing
445	! a global use of the multi-layer energy budget.
446	CALL albedo_surface_main( &
447	kjit, kjpindex, hist_id, hist2_id, &
448	index, drysoil_frac, veget_max, coszang, &
449	frac_nobio, frac_snow_veg, frac_snow_nobio, &
450	snow, snow_age, snow_nobio, snow_nobio_age, &
451	z0m, laieff_fit, &
452	albedo, Light_Abs_Tot, Light_Tran_Tot, &
453	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
454	laieff_isotrop, veget)
455	!+++++++++++
456
457	! Debug
458	IF(printlev >= 4)THEN
459	DO ivm = 1,nvm
460	WRITE(numout,*) 'laieff_isotrop condveg_main for ivm', &
461	ivm,'is', laieff_isotrop(:,:,ivm)
462	ENDDO
463	ENDIF
464	!-
465
466
467	IF (printlev>=3) WRITE (numout,*)' condveg_main done '
468
469	END SUBROUTINE condveg_main
470
471	!! =============================================================================================================================
472	!! SUBROUTINE : condveg_finalize
473	!!
474	!>\BRIEF Write to restart file
475	!!
476	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in condveg
477	!! to restart file
478	!!
479	!! RECENT CHANGE(S) : None
480	!!
481	!! REFERENCE(S) : None
482	!!
483	!! FLOWCHART : None
484	!! \n
485	!_ ==============================================================================================================================
486	!+++CHECK+++
487	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
488	! when running a larger domain. Needs to be corrected when implementing
489	! a global use of the multi-layer energy budget.
490	SUBROUTINE condveg_finalize (kjit, kjpindex, rest_id, &
491	z0m, z0h, roughheight, &
492	albedo, Light_Abs_Tot, Light_Tran_Tot, &
493	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
494	laieff_isotrop)
495	!+++++++++++
496
497	!! 0. Variable and parameter declaration
498	!! 0.1 Input variables
499	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
500	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
501	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
502	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: z0m !! Roughness for momentum
503	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: z0h !! Roughness for heat
504	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: roughheight !! Grid effective roughness height (m)
505	REAL(r_std), DIMENSION(kjpindex,n_spectralbands), &
506	INTENT(in) :: albedo !! Albedo (two stream radiation transfer model)
507	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), &
508	INTENT(in) :: Light_Abs_Tot !! Absorbed radiation per layer, averaged between light
509	!! from a collimated (direct) source and that from an
510	!! isotropic (diffuse) source. Expressed as the fraction
511	!! of overall light hitting the canopy.
512	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), &
513	INTENT(in) :: Light_Tran_Tot !! Transmitted radiation per layer, averaged between light
514	!! from a collimated (direct) source and that from an
515	!! isotropic (diffuse) source. Expressed as the fraction
516	!! of overall light hitting the canopy.
517	!+++CHECK+++
518	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
519	! when running a larger domain. Needs to be corrected when implementing
520	! a global use of the multi-layer energy budget.
521	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT(in) :: Light_Abs_Tot_mean !! total light absorption for a given level
522	!!$ REAL(r_std),DIMENSION(nlevels_tot), INTENT(in) :: Light_Alb_Tot_mean !! total albedo for a given level
523	!+++++++++++
524	REAL(r_std), DIMENSION(kjpindex,nlevels_tot,nvm), INTENT(in) :: laieff_isotrop !! Leaf Area Index Effective converts
525	!! 3D lai into 1D lai for two stream
526
527	!_ ================================================================================================================================
528	CALL restput_p (rest_id, 'z0m', nbp_glo, 1, 1, kjit, z0m, 'scatter', nbp_glo, index_g)
529	!-
530	CALL restput_p (rest_id, 'z0h', nbp_glo, 1, 1, kjit, z0h, 'scatter', nbp_glo, index_g)
531	!-
532	CALL restput_p (rest_id, 'roughheight', nbp_glo, 1, 1, kjit, roughheight, 'scatter', nbp_glo, index_g)
533
534	! Finalize albedo module
535	!+++CHECK+++
536	! Variable dimensions xxx_Tot_mean are for a single pixel. Causing 1+1 problems
537	! when running a larger domain. Needs to be corrected when implementing
538	! a global use of the multi-layer energy budget.
539	CALL albedo_surface_finalize (kjit, kjpindex, rest_id, &
540	albedo, Light_Abs_Tot, Light_Tran_Tot, &
541	!!$ Light_Abs_Tot_mean, Light_Alb_Tot_mean, &
542	laieff_isotrop)
543	!+++++++++++
544
545	END SUBROUTINE condveg_finalize
546
547	!! ==============================================================================================================================
548	!! SUBROUTINE : condveg_clear
549	!!
550	!>\BRIEF Deallocate albedo variables
551	!!
552	!! DESCRIPTION : None
553	!!
554	!! RECENT CHANGE(S): None
555	!!
556	!! MAIN OUTPUT VARIABLE(S): None
557	!!
558	!! REFERENCES : None
559	!!
560	!! FLOWCHART : None
561	!! \n
562	!_ ================================================================================================================================
563
564	SUBROUTINE condveg_clear ()
565
566	l_first_condveg=.TRUE.
567
568	CALL albedo_surface_clear()
569
570	END SUBROUTINE condveg_clear
571
572	!! ==============================================================================================================================
573	!! SUBROUTINE : condveg_frac_snow
574	!!
575	!>\BRIEF This subroutine calculates the fraction of snow on vegetation and nobio
576	!!
577	!! DESCRIPTION
578	!!
579	!! RECENT CHANGE(S): These calculations were previously done in condveg_snow.
580	!!
581	!! REFERENCE(S) :
582	!!
583	!! FLOWCHART : None
584	!! \n
585	!_ ================================================================================================================================
586
587	SUBROUTINE condveg_frac_snow(kjpindex, snow_nobio, snowrho, snowdz, &
588	frac_snow_veg, frac_snow_nobio)
589
590	!! 0. Variable and parameter declaration
591	!! 0.1 Input variables
592	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
593	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
594	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density at each snow layer (Kg/m^3)
595	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer (m)
596
597	!! 0.2 Output variables
598	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: frac_snow_veg !! Fraction of snow on vegetation (unitless ratio)
599	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
600
601	!! 0.3 Local variables
602	REAL(r_std), DIMENSION(kjpindex) :: snowrho_ave !! Average snow density
603	REAL(r_std), DIMENSION(kjpindex) :: snowdepth !! Snow depth
604	REAL(r_std), DIMENSION(kjpindex) :: snowrho_snowdz !! Snow rho time snowdz
605	INTEGER(i_std) :: jv
606
607	!_ ================================================================================================================================
608
609	! Calculate snow cover fraction for both total vegetated and
610	! total non-vegetative surfaces.
611	snowdepth=sum(snowdz,2)
612	snowrho_snowdz=sum(snowrho*snowdz,2)
613	WHERE(snowdepth(:) .LT. min_sechiba)
614	frac_snow_veg(:) = 0.
615	ELSEWHERE
616	snowrho_ave(:)=snowrho_snowdz(:)/snowdepth(:)
617	frac_snow_veg(:) = tanh(snowdepth(:)/(0.025*(snowrho_ave(:)/50.)))
618	END WHERE
619
620	DO jv = 1, nnobio
621	frac_snow_nobio(:,jv) = MIN(MAX(snow_nobio(:,jv),zero)/&
622	(MAX(snow_nobio(:,jv),zero)+snowcri_alb*sn_dens/100.0),un)
623	ENDDO
624
625	IF (printlev>=3) WRITE (numout,*) ' condveg_frac_snow done '
626
627	END SUBROUTINE condveg_frac_snow
628
629
630
631	!! ==============================================================================================================================
632	!! SUBROUTINE : condveg_z0cdrag
633	!!
634	!>\BRIEF Computation of grid average of roughness length by calculating
635	!! the drag coefficient.
636	!!
637	!! DESCRIPTION : This routine calculates the mean roughness height and mean
638	!! effective roughness height over the grid cell. The mean roughness height (z0)
639	!! is computed by averaging the drag coefficients \n
640	!!
641	!! \latexonly
642	!! \input{z0cdrag1.tex}
643	!! \endlatexonly
644	!! \n
645	!!
646	!! where C is the drag coefficient at the height of the vegetation, kappa is the
647	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
648	!! The reference level for z needs to be high enough above the canopy to avoid
649	!! singularities of the LOG. This height is set to minimum 10m above ground.
650	!! The drag coefficient increases with roughness height to represent the greater
651	!! turbulence generated by rougher surfaces.
652	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
653	!!
654	!! The roughness height for the non-vegetative surfaces is calculated in a second step.
655	!! In order to calculate the transfer coefficients the
656	!! effective roughness height is calculated. This effective value is the difference
657	!! between the height of the vegetation and the zero plane displacement height.\nn
658	!!
659	!! RECENT CHANGE(S): None
660	!!
661	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
662	!!
663	!! REFERENCE(S) : None
664	!!
665	!! FLOWCHART : None
666	!! \n
667	!_ ================================================================================================================================
668
669	SUBROUTINE condveg_z0cdrag (kjpindex, veget, veget_max, frac_nobio,&
670	totfrac_nobio, zlev, height, height_dom, tot_bare_soil, frac_snow_veg,&
671	z0m, z0h, roughheight)
672
673	!! 0. Variable and parameter declaration
674
675	!! 0.1 Input variables
676
677	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
678	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
679	!! (m^2 m^{-2})
680	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
681	!! (= ind*cn_ind) (m^2 m^{-2})
682	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
683	!! i.e. continental ice, lakes, etc. (unitless)
684	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
685	!! i.e. continental ice, lakes, etc. (unitless)
686	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
687	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
688	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height_dom !! Dominant vegetation height (m)
689	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
690	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
691
692	!! 0.2 Output variables
693
694	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0m !! Roughness height for momentum (m)
695	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0h !! Roughness height for heat (m)
696	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
697
698	!! 0.3 Modified variables
699
700	!! 0.4 Local variables
701
702	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
703	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
704	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
705	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
706	REAL(r_std), DIMENSION(kjpindex) :: d_veg !! PFT coverage of vegetative PFTs
707	!! (= ind*cn_ind) (m^2 m^{-2})
708	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
709	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
710	!! i.e. continental ice, lakes, etc.
711	REAL(r_std), DIMENSION(kjpindex) :: dragm !! Dragcoefficient for momentum
712	REAL(r_std), DIMENSION(kjpindex) :: dragh !! Dragcoefficient for heat
713	REAL(r_std), DIMENSION(kjpindex) :: z0_ground !! z0m value used for ground surface
714
715	!_ ================================================================================================================================
716
717	!! 1. Preliminary calculation
718
719	! Set maximal height of first layer
720	ztmp(:) = MAX(10., zlev(:))
721
722	z0_ground(:) = (1.-frac_snow_veg(:))z0_bare + frac_snow_veg(:)z0_bare/10.
723
724	! Calculate roughness for non-vegetative surfaces
725	! with the von Karman constant
726	dragm(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/z0_ground))**2
727	dragh(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/(z0_ground/ratio_z0m_z0h(1))))*(ct_karman/LOG(ztmp(:)/z0_ground))
728
729	! Fraction of bare soil
730	sumveg(:) = tot_bare_soil(:)
731
732	! Set average vegetation height to zero
733	ave_height(:) = zero
734
735	!! 2. Calculate the mean roughness height
736
737	! Calculate the mean roughness height of
738	! vegetative PFTs over the grid cell
739	DO jv = 2, nvm
740
741	! In the case of forest, use parameter veget_max because
742	! tree trunks influence the roughness even when there are no leaves
743	IF ( is_tree(jv) ) THEN
744	! In the case of grass, use parameter veget because grasses
745	! only influence the roughness during the growing season
746	d_veg(:) = veget_max(:,jv)
747	ELSE
748	! grasses only have an influence if they are really there!
749	d_veg(:) = veget(:,jv)
750	ENDIF
751
752	! Calculate the average roughness over the grid cell:
753	! The unitless drag coefficient is per vegetative PFT
754	! calculated by use of the von Karman constant, the height
755	! of the first layer and the roughness. The roughness
756	! is calculated as the vegetation height per PFT
757	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
758	! If this scaled value is lower than 0.01 then the value for
759	! the roughness of bare soil (0.01) is used.
760	! The sum over all PFTs gives the average roughness
761	! per grid cell for the vegetative PFTs.
762	IF(use_height_dom)THEN
763	dragm(:) = dragm(:) + d_veg(:) (ct_karman/LOG(ztmp(:)/MAX(height_dom(:,jv)z0_over_height(jv),z0_ground)))**2
764	dragh(:) = dragh(:) + d_veg(:) (ct_karman/LOG(ztmp(:)/(MAX(height_dom(:,jv)z0_over_height(jv),z0_ground) / &
765	ratio_z0m_z0h(jv)))) * (ct_karman/LOG(ztmp(:)/MAX(height_dom(:,jv)*z0_over_height(jv),z0_ground)))
766	ELSE
767	dragm(:) = dragm(:) + d_veg(:) (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)z0_over_height(jv),z0_ground)))**2
768	dragh(:) = dragh(:) + d_veg(:) (ct_karman/LOG(ztmp(:)/(MAX(height(:,jv)z0_over_height(jv),z0_ground) / &
769	ratio_z0m_z0h(jv)))) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)*z0_over_height(jv),z0_ground)))
770	ENDIF
771
772	! Sum of bare soil and fraction vegetated fraction
773	sumveg(:) = sumveg(:) + d_veg(:)
774
775	! Weigh height of vegetation with maximal cover fraction
776	IF(use_height_dom)THEN
777	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
778	ELSE
779	ave_height(:) = ave_height(:) + veget_max(:,jv)*height_dom(:,jv)
780	ENDIF
781	ENDDO
782
783	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
784
785	! Search for pixels with vegetated part to normalise
786	! roughness height
787	WHERE ( sumveg(:) .GT. min_sechiba )
788	dragm(:) = dragm(:) / sumveg(:)
789	dragh(:) = dragh(:) / sumveg(:)
790	ENDWHERE
791	! Calculate fraction of roughness for vegetated part
792	dragm(:) = (un - totfrac_nobio(:)) * dragm(:)
793	dragh(:) = (un - totfrac_nobio(:)) * dragh(:)
794
795	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
796
797	! Set rougness for ice
798	IF ( jv .EQ. iice ) THEN
799	z0_nobio = z0_ice
800	ELSE
801	WRITE(numout,*) 'jv=',jv
802	WRITE(numout,*) 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
803	CALL ipslerr_p(3,'condveg_z0cdrag','DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
804	ENDIF
805
806	! Sum of vegetative roughness length and non-vegetative roughness length
807	dragm(:) = dragm(:) + frac_nobio(:,jv) (ct_karman/LOG(ztmp(:)/z0_nobio))*2
808	dragh(:) = dragh(:) + frac_nobio(:,jv) (ct_karman/LOG(ztmp(:)/z0_nobio/ratio_z0m_z0h(1)))(ct_karman/LOG(ztmp(:)/z0_nobio))
809
810	ENDDO ! Loop over # of non-vegative surfaces
811
812	!! 4. Calculate the zero plane displacement height and effective roughness length
813
814	! Take the exponential of the roughness
815	z0m(:) = ztmp(:) / EXP(ct_karman/SQRT(dragm(:)))
816	z0h(:) = ztmp(:) / EXP((ct_karman*2.)/(dragh(:)LOG(ztmp(:)/z0m(:))))
817
818	! Compute the zero plane displacement height which
819	! is an equivalent height for the absorption of momentum
820	zhdispl(:) = ave_height(:) * height_displacement
821
822	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
823	! This is the height over which the roughness acts. It combines the
824	! zero plane displacement height and the vegetation height.
825	roughheight(:) = ave_height(:) - zhdispl(:)
826
827	END SUBROUTINE condveg_z0cdrag
828
829
830	!! ==============================================================================================================================
831	!! SUBROUTINE : condveg_z0cdrag_dyn
832	!!
833	!>\BRIEF Computation of grid average of roughness length by calculating
834	!! the drag coefficient based on formulation proposed by Su et al. (2001).
835	!!
836	!! DESCRIPTION : This routine calculates the mean roughness height and mean
837	!! effective roughness height over the grid cell. The mean roughness height (z0)
838	!! is computed by averaging the drag coefficients \n
839	!!
840	!! \latexonly
841	!! \input{z0cdrag1.tex}
842	!! \endlatexonly
843	!! \n
844	!!
845	!! where C is the drag coefficient at the height of the vegetation, kappa is the
846	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
847	!! The reference level for z needs to be high enough above the canopy to avoid
848	!! singularities of the LOG. This height is set to minimum 10m above ground.
849	!! The drag coefficient increases with roughness height to represent the greater
850	!! turbulence generated by rougher surfaces.
851	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
852	!! In the formulation of Su et al. (2001), one distinguishes the roughness height for
853	!! momentum (z0m) and the one for heat (z0h).
854	!! z0m is computed as a function of LAI (z0m increases with LAI) and z0h is computed
855	!! with a so-called kB-1 term (z0m/z0h=exp(kB-1))
856	!!
857	!! RECENT CHANGE(S): Written by N. Vuichard (2016)
858	!!
859	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
860	!!
861	!! REFERENCE(S) :
862	!! - Su, Z., Schmugge, T., Kustas, W.P., Massman, W.J., 2001. An Evaluation of Two Models for
863	!! Estimation of the Roughness Height for Heat Transfer between the Land Surface and the Atmosphere. J. Appl.
864	!! Meteorol. 40, 1933â1951. doi:10.1175/1520-0450(2001)
865	!! - Ershadi, A., McCabe, M.F., Evans, J.P., Wood, E.F., 2015. Impact of model structure and parameterization
866	!! on Penman-Monteith type evaporation models. J. Hydrol. 525, 521â535. doi:10.1016/j.jhydrol.2015.04.008
867	!!
868	!! FLOWCHART : None
869	!! \n
870	!_ ================================================================================================================================
871
872	SUBROUTINE condveg_z0cdrag_dyn (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, height_dom, &
873	& temp_air, pb, u, v, lai, frac_snow_veg, z0m, z0h, roughheight)
874
875	!! 0. Variable and parameter declaration
876
877	!! 0.1 Input variables
878
879	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
880	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
881	!! (m^2 m^{-2})
882	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
883	!! (= ind*cn_ind) (m^2 m^{-2})
884	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
885	!! i.e. continental ice, lakes, etc. (unitless)
886	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
887	!! i.e. continental ice, lakes, etc. (unitless)
888	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
889	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
890	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height_dom !! Dominant vegetation height (m)
891	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: temp_air !! 2m air temperature (K)
892	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
893	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
894	!! @tex $(m.s^{-1})$ @endtex
895	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
896	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
897	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
898	!! 0.2 Output variables
899
900	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0m !! Roughness height for momentum (m)
901	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0h !! Roughness height for heat (m)
902	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
903
904	!! 0.3 Modified variables
905
906	!! 0.4 Local variables
907	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
908	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
909	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
910	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
911	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
912	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
913	!! i.e. continental ice, lakes, etc.
914	REAL(r_std), DIMENSION(kjpindex) :: z0m_pft !! Roughness height for momentum for a specific PFT
915	REAL(r_std), DIMENSION(kjpindex) :: z0h_pft !! Roughness height for heat for a specific PFT
916	REAL(r_std), DIMENSION(kjpindex) :: dragm !! Dragcoefficient for momentum
917	REAL(r_std), DIMENSION(kjpindex) :: dragh !! Dragcoefficient for heat
918	REAL(r_std), DIMENSION(kjpindex) :: eta !! Ratio of friction velocity to the wind speed at the canopy top - See Ershadi et al. (2015)
919	REAL(r_std), DIMENSION(kjpindex) :: eta_ec !! Within-canopy wind speed profile estimation coefficient - See Ershadi et al. (2015)
920	REAL(r_std), DIMENSION(kjpindex) :: Ct_star !! Heat transfer coefficient of the soil - see Su et al. (2001)
921	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))
922	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))
923	REAL(r_std), DIMENSION(kjpindex) :: fc !! fractional canopy coverage
924	REAL(r_std), DIMENSION(kjpindex) :: fs !! fractional soil coverage
925	REAL(r_std), DIMENSION(kjpindex) :: Reynolds !! Reynolds number
926	REAL(r_std), DIMENSION(kjpindex) :: wind !! wind Speed (m)
927	REAL(r_std), DIMENSION(kjpindex) :: u_star !! friction velocity
928	REAL(r_std), DIMENSION(kjpindex) :: z0_ground !! z0m value used for ground surface
929	REAL(r_std), DIMENSION(kjpindex,nvm) :: loc_height !! Vegetation height (m)
930
931	!_ ================================================================================================================================
932
933	!! 1. Preliminary calculation
934
935	! Set maximal height of first layer
936	ztmp(:) = MAX(10., zlev(:))
937
938	z0_ground(:) = (1.-frac_snow_veg(:))z0_bare + frac_snow_veg(:)z0_bare/10.
939
940	! Calculate roughness for non-vegetative surfaces
941	! with the von Karman constant
942	dragm(:) = veget_max(:,1) * (ct_karman/LOG(ztmp(:)/z0_ground(:)))**2
943
944	wind(:) = SQRT(u(:)u(:)+v(:)v(:))
945	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG(zlev(:)/z0_ground(:))
946	Reynolds(:) = z0_ground(:) * u_star(:) &
947	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
948
949	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
950
951	dragh(:) = veget_max(:,1) * (ct_karman/LOG(ztmp(:)/z0_ground(:)))*(ct_karman/LOG(ztmp(:)/(z0_ground(:)/exp(kBs_m1(:))) ))
952
953	! Fraction of bare soil
954	sumveg(:) = veget_max(:,1)
955
956	! Set average vegetation height to zero
957	ave_height(:) = zero
958
959	!! 2. Calculate the mean roughness height
960	! height is an input variable write it to a local variable
961	! such that the value can be changed where needed
962	IF(use_height_dom)THEN
963	loc_height(:,:) = height_dom(:,:)
964	ELSE
965	loc_height(:,:) = height(:,:)
966	ENDIF
967
968	! Calculate the mean roughness height of
969	! vegetative PFTs over the grid cell
970	DO jv = 2, nvm
971
972	WHERE(veget_max(:,jv) .GT. zero)
973
974	WHERE (lai(:,jv).GT.min_sechiba .AND. height(:,jv).LT.min_sechiba)
975
976	! This is possible when the trees are under coppicing. Coppicing removes
977	! the aboveground biomass but retains the belowground components. LAI
978	! is calculated making use of the sum of above and belowground biomass
979	! height only makes use of the aboverground biomass. Following
980	! coppincing we thus have an LAI but no height. I propose a quick fix
981	! by prescribing a height of 10 cm (= the stumps of the trees)
982	loc_height(:,jv) = 0.1
983
984	ENDWHERE
985
986	! Calculate the average roughness over the grid cell:
987	! The unitless drag coefficient is per vegetative PFT
988	! calculated by use of the von Karman constant, the height
989	! of the first layer and the roughness. The roughness
990	! is calculated as the vegetation height per PFT
991	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
992	! If this scaled value is lower than 0.01 then the value for
993	! the roughness of bare soil (0.01) is used.
994	! The sum over all PFTs gives the average roughness
995	! per grid cell for the vegetative PFTs.
996	eta(:) = c1 - c2 * exp(-c3 * Cdrag_foliage * lai(:,jv))
997
998	z0m_pft(:) = (loc_height(:,jv)(1-height_displacement)(exp(-ct_karman/eta(:))-exp(-ct_karman/(c1-c2)))) &
999	+ z0_ground(:)
1000
1001	dragm(:) = dragm(:) + veget_max(:,jv) (ct_karman/LOG(ztmp(:)/z0m_pft(:)))*2
1002
1003	fc(:) = veget(:,jv)/veget_max(:,jv)
1004	fs(:) = 1. - fc(:)
1005
1006	eta_ec(:) = ( Cdrag_foliage * lai(:,jv)) / (2 * eta(:)*eta(:))
1007	wind(:) = SQRT(u(:)u(:)+v(:)v(:))
1008	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG((zlev(:)+(loc_height(:,jv)*(1-height_displacement)))/z0m_pft(:))
1009	Reynolds(:) = z0_ground(:) * u_star(:) &
1010	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1011
1012	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1013	Ct_star(:) = Prandtl*(-2./3.) SQRT(1./Reynolds(:))
1014
1015	WHERE(lai(:,jv) .GT. min_sechiba .AND. loc_height(:,jv).GT.min_sechiba)
1016	kB_m1(:) = (ct_karman * Cdrag_foliage) / (4 * Ct * eta(:) * (1 - exp(-eta_ec(:)/2.))) * fc(:)**2. &
1017	+ 2fc(:)fs(:) * (ct_karman * eta(:) * z0m_pft(:) / loc_height(:,jv)) / Ct_star(:) &
1018	+ kBs_m1(:) * fs(:)**2.
1019	ELSEWHERE
1020	kB_m1(:) = kBs_m1(:) * fs(:)**2.
1021	ENDWHERE
1022
1023	z0h_pft(:) = z0m_pft(:) / exp(kB_m1(:))
1024
1025	dragh(:) = dragh(:) + veget_max(:,jv) * (ct_karman/LOG(ztmp(:)/z0m_pft(:)))*(ct_karman/LOG(ztmp(:)/z0h_pft(:)))
1026
1027	! Sum of bare soil and fraction vegetated fraction
1028	sumveg(:) = sumveg(:) + veget_max(:,jv)
1029
1030	! Weigh height of vegetation with maximal cover fraction
1031	ave_height(:) =ave_height(:) + veget_max(:,jv)*loc_height(:,jv)
1032
1033	ENDWHERE
1034	ENDDO
1035
1036	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1037
1038	! Search for pixels with vegetated part to normalise
1039	! roughness height
1040	WHERE ( sumveg(:) .GT. min_sechiba )
1041	dragh(:) = dragh(:) / sumveg(:)
1042	dragm(:) = dragm(:) / sumveg(:)
1043	ENDWHERE
1044
1045	! Calculate fraction of roughness for vegetated part
1046	dragh(:) = (un - totfrac_nobio(:)) * dragh(:)
1047	dragm(:) = (un - totfrac_nobio(:)) * dragm(:)
1048
1049	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1050
1051	! Set rougness for ice
1052	IF ( jv .EQ. iice ) THEN
1053	z0_nobio = z0_ice
1054	ELSE
1055	WRITE(numout,*) 'jv=',jv
1056	WRITE(numout,*) 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1057	CALL ipslerr_p(3,'condveg_z0cdrag_dyn','DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1058	ENDIF
1059
1060	! Sum of vegetative roughness length and non-vegetative roughness length
1061	dragm(:) = dragm(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1062
1063	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG(zlev(:)/z0_nobio)
1064	Reynolds(:) = z0_nobio * u_star(:) &
1065	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1066
1067	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1068
1069	dragh(:) = dragh(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio)) * &
1070	(ct_karman/LOG(ztmp(:)/(z0_nobio/ exp(kBs_m1(:))) ))
1071
1072	ENDDO ! Loop over # of non-vegative surfaces
1073
1074	!! 4. Calculate the zero plane displacement height and effective roughness length
1075	! Take the exponential of the roughness
1076	z0m(:) = ztmp(:) / EXP(ct_karman/SQRT(dragm(:)))
1077	z0h(:) = ztmp(:) / EXP((ct_karman*2.)/(dragh(:)LOG(ztmp(:)/z0m(:))))
1078
1079	! Compute the zero plane displacement height which
1080	! is an equivalent height for the absorption of momentum
1081	zhdispl(:) = ave_height(:) * height_displacement
1082
1083	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1084	! This is the height over which the roughness acts. It combines the
1085	! zero plane displacement height and the vegetation height.
1086	roughheight(:) = ave_height(:) - zhdispl(:)
1087
1088	END SUBROUTINE condveg_z0cdrag_dyn
1089
1090
1091	END MODULE condveg

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

Download in other formats: