Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_sechiba/condveg.f90 @ 7475

Last change on this file since 7475 was 2189, checked in by matthew.mcgrath, 10 years ago
DEV: First stage cleaning up histwrites...still need to add a few variables
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 69.6 KB

Line
1	! ===============================================================================================================================
2	! MODULE : condveg
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.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. :: z0cdrag_ave to choose between two methods to calculate the grid average of
14	!! the roughness height. If set to true: the grid average is calculated by the drag coefficients
15	!! per PFT. If set to false: the grid average is calculated by the logarithm of the
16	!! roughness height per PFT.\n
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.
20	!! This is useful if 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) : None
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	!
51	! modules used :
52	USE constantes
53	USE constantes_soil
54	USE pft_parameters
55	USE interpol_help
56	USE ioipsl_para
57	USE albedo
58	USE slowproc
59
60	IMPLICIT NONE
61
62	! public routines :
63	! condveg_main only
64	PRIVATE
65	PUBLIC :: condveg_main,condveg_clear
66
67	!
68	! variables used inside condveg module : declaration and initialisation
69	!
70	LOGICAL, SAVE :: l_first_condveg=.TRUE. !! To keep first call's trace
71	!$OMP THREADPRIVATE(l_first_condveg)
72	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_dry(:,:) !! Albedo values for the dry bare soil (unitless)
73	!$OMP THREADPRIVATE(soilalb_dry)
74	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_wet(:,:) !! Albedo values for the wet bare soil (unitless)
75	!$OMP THREADPRIVATE(soilalb_wet)
76	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_moy(:,:) !! Albedo values for the mean bare soil (unitless)
77	!$OMP THREADPRIVATE(soilalb_moy)
78	REAL(r_std), ALLOCATABLE, SAVE :: alb_bare(:,:) !! Mean bare soil albedo for visible and near-infrared
79	!! range (unitless)
80	!$OMP THREADPRIVATE(alb_bare)
81	REAL(r_std), ALLOCATABLE, SAVE :: alb_veget(:,:) !! Mean vegetation albedo for visible and
82	!! near-infrared range (unitless)
83	!$OMP THREADPRIVATE(alb_veget)
84	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: albedo_snow !! Mean snow albedo over visible and near-infrared
85	!! range (unitless)
86	!$OMP THREADPRIVATE(albedo_snow)
87	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: albedo_glob !! Mean surface albedo over visible and
88	!! near-infrared range (unitless)
89	!$OMP THREADPRIVATE(albedo_glob)
90	CONTAINS
91
92	!! ==============================================================================================================================
93	!! SUBROUTINE : condveg_main
94	!!
95	!>\BRIEF Calls the subroutines to initialise the variables, update the variables
96	!! and write out data and restart files.
97	!!
98	!!
99	!! MAIN OUTPUT VARIABLE(S): emis (emissivity), albedo (albedo of
100	!! vegetative PFTs in visible and near-infrared range), z0 (surface roughness height),
101	!! roughheight (grid effective roughness height), soil type (fraction of soil types)
102	!!
103	!!
104	!! REFERENCE(S) : None
105	!!
106	!! FLOWCHART : None
107	!!
108	!! REVISION(S) : None
109	!!
110	!_ ================================================================================================================================
111
112	SUBROUTINE condveg_main (kjit, kjpindex, dtradia, ldrestart_write, &
113	index, lalo, neighbours, resolution, contfrac, &
114	veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
115	snow, snow_age, snow_nobio, snow_nobio_age, drysoil_frac, &
116	height, emis, albedo, z0, roughheight, &
117	rest_id, hist_id, hist2_id, lai, &
118	sinang, circ_class_biomass, circ_class_n, &
119	lai_split,z0_veg, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
120	laieff_fit, albedo_pft, frac_snow_nobio, frac_snow_veg)
121
122	!! 0. Variable and parameter declaration
123
124	!! 0.1 Input variables
125
126	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
127	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
128	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
129	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier
130	INTEGER(i_std), OPTIONAL, INTENT (in) :: hist2_id !! _History_ file 2 identifier
131	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
132	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for restart file to write
133	! input fields
134	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
135	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
136	INTEGER(i_std),DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land
137	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
138	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
139	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
140	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
141	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
142	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
143	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
144	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
145	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
146	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
147	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
148	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
149	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
150	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m^2 m^{-2})
151	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: sinang !! the sine of the solar angle from the horizon (unitless)
152	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nparts,nelements), &
153	INTENT(IN) :: circ_class_biomass !! Stem diameter
154	!! @tex $(m)$ @endtex
155	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc), &
156	INTENT(IN) :: circ_class_n !! Number of trees within each circumference
157	REAL(r_std),DIMENSION(nlevels),INTENT(IN) :: lai_split !!!!temp
158	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
159
160	!! 0.2 Output variables
161
162	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
163	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
164	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Roughness
165	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
166	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0_veg !! Roughness height of vegetated part (m)
167	REAL(r_std),DIMENSION (:,:,:), &
168	INTENT (out) :: Isotrop_Abs_Tot_p !! Absorbed radiation per level for photosynthesis
169	REAL(r_std),DIMENSION (:,:,:), &
170	INTENT (out) :: Isotrop_Tran_Tot_p !! Transmitted radiation per level for photosynthesis
171	REAL(r_std), DIMENSION(kjpindex,nvm,n_spectralbands), &
172	INTENT (out) :: albedo_pft !! Albedo (two stream radiation transfer model)
173	!! for visible and near-infrared range
174	!! for each PFT (unitless)
175	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
176	!! (unitless ratio)
177	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: frac_snow_veg !! the fraction of ground covered with snow, between 0 and 1
178
179	!! 0.3 Modified variables
180
181	!! 0.4 Local variables
182
183	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
184	!_ ================================================================================================================================
185
186	!! 1. Call subroutines to start initialisation
187	! Is TRUE if condveg is called the first time
188	IF (l_first_condveg) THEN
189	! Document what the programm is doing
190	IF (long_print) WRITE (numout,*) ' l_first_condveg : call condveg_init '
191
192	! Call the subroutine 'condveg_init' to allocate local array
193	CALL condveg_init (kjit, kjpindex, index, veget, &
194	lalo, neighbours, resolution, contfrac, rest_id)
195
196	! Call the subroutine 'condveg_var_init' to initialise local array
197	! Reads in map for soil albedo
198	CALL condveg_var_init (kjpindex, veget, veget_max, frac_nobio, &
199	totfrac_nobio, drysoil_frac, zlev, height, &
200	circ_class_biomass, circ_class_n, emis, albedo, &
201	z0, roughheight, lai, sinang, &
202	snow, snow_age, snow_nobio, snow_nobio_age, albedo_snow, &
203	albedo_glob, lai_split, z0_veg, Isotrop_Abs_Tot_p,&
204	Isotrop_Tran_Tot_p, laieff_fit, albedo_pft, frac_snow_nobio,&
205	frac_snow_veg)
206
207
208	RETURN
209
210	ENDIF
211
212	!! 2. Call subroutines to update fields
213
214	! Call the routine 'condveg_var_update' to update the fields of albedo, emissivity and roughness
215	CALL condveg_var_update (kjpindex, veget, veget_max, frac_nobio, &
216	totfrac_nobio, drysoil_frac, zlev, height, &
217	circ_class_biomass, circ_class_n, emis, albedo, &
218	z0, roughheight, lai, sinang, &
219	snow, snow_age, snow_nobio, snow_nobio_age, albedo_snow, &
220	albedo_glob, lai_split,z0_veg, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
221	laieff_fit, albedo_pft, frac_snow_nobio, frac_snow_veg)
222
223	!! 3. Call subroutines to write restart files and data
224
225	! To write restart files for soil albedo variables
226	IF (ldrestart_write) THEN
227	!
228	var_name = 'soilalbedo_dry'
229	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_dry, 'scatter', nbp_glo, index_g)
230	!
231	var_name = 'soilalbedo_wet'
232	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_wet, 'scatter', nbp_glo, index_g)
233	!
234	var_name = 'soilalbedo_moy'
235	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_moy, 'scatter', nbp_glo, index_g)
236	!
237	RETURN
238	!
239	ENDIF
240
241	! If this logical flag is set to true, the model
242	! will output all its data according to the ALMA
243	! convention.
244	! To facilitate the exchange of forcing data for land-surface schemes and the results produced by these schemes,
245	! ALMA (Assistance for Land-surface Modelling activities) proposes to establish standards.
246	! http://www.lmd.jussieu.fr/~polcher/ALMA/
247	CALL xios_orchidee_send_field("soilalb_vis",alb_bare(:,1))
248	CALL xios_orchidee_send_field("soilalb_nir",alb_bare(:,2))
249	CALL xios_orchidee_send_field("vegalb_vis",alb_veget(:,1))
250	CALL xios_orchidee_send_field("vegalb_nir",alb_veget(:,2))
251	CALL xios_orchidee_send_field("albedo_alma",albedo_glob)
252	CALL xios_orchidee_send_field("SAlbedo",albedo_snow)
253
254	! These variables are NOT written out in the ALMA convention.
255	! At least, they are not declared this way in intersurf.
256	IF( .NOT. almaoutput)THEN
257	CALL histwrite_p(hist_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
258	CALL histwrite_p(hist_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
259	CALL histwrite_p(hist_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
260	CALL histwrite_p(hist_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
261	CALL histwrite_p(hist2_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
262	CALL histwrite_p(hist2_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
263	CALL histwrite_p(hist2_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
264	CALL histwrite_p(hist2_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
265	ENDIF
266
267	IF (long_print) WRITE (numout,*)' condveg_main done '
268
269
270	END SUBROUTINE condveg_main
271
272	!! ==============================================================================================================================
273	!! SUBROUTINE : condveg_init
274	!!
275	!>\BRIEF Dynamic allocation of the variables, the dimensions of the
276	!! variables are determined by user-specified settings.
277	!!
278	!! DESCRIPTION : The domain size (:: kjpindex) is used to allocate the correct
279	!! dimensions to all variables in condveg.
280	!!
281	!! RECENT CHANGE(S): None
282	!!
283	!! MAIN OUTPUT VARIABLE(S): Strictly speaking the subroutine has no output
284	!! variables. However, the routine allocates memory and builds new indexing
285	!! variables for later use.
286	!!
287	!! REFERENCE(S) : None
288	!!
289	!! FLOWCHART : None
290	!! \n
291	!_ ================================================================================================================================
292
293	SUBROUTINE condveg_init (kjit, kjpindex, index, veget, &
294	lalo, neighbours, resolution, contfrac, rest_id)
295
296	!! 0. Variable and parameter declaration
297
298	!! 0.1. Input variables
299
300	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
301	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
302	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Index for the points on the map (unitless)
303	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in):: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
304	!! (m^2 m^{-2})
305	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates (degree)
306	INTEGER(i_std),DIMENSION (kjpindex,8), INTENT(in):: neighbours !! Neighbouring land grid cell
307	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size of grid in x and y direction (m)
308	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: contfrac !! Fraction of land in each grid box
309	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier
310
311	!! 0.2. Output variables
312
313	!! 0.3 Modified variables
314
315	!! 0.4 Local variables
316
317	INTEGER(i_std) :: ji !! Index
318	INTEGER(i_std) :: ier !! Check errors in memory allocation
319	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
320
321	!_ ================================================================================================================================
322
323	!! 1. Choice of calculation of snow albedo and soil albedo
324
325	! Is TRUE if condveg is called the first time
326	IF (l_first_condveg) THEN
327	!
328	l_first_condveg=.FALSE.
329
330	!! 2. Allocate all albedo variables
331
332	! Dry soil albedo
333	ALLOCATE (soilalb_dry(kjpindex,2),stat=ier)
334	IF (ier.NE.0) THEN
335	WRITE (numout,) ' error in soilalb_dry allocation. We stop.We need kjpindex2 words = ',kjpindex*2
336	STOP 'condveg_init'
337	END IF
338	soilalb_dry(:,:) = val_exp
339
340	! Wet soil albedo
341	ALLOCATE (soilalb_wet(kjpindex,2),stat=ier)
342	IF (ier.NE.0) THEN
343	WRITE (numout,) ' error in soilalb_wet allocation. We stop.We need kjpindex2 words = ',kjpindex*2
344	STOP 'condveg_init'
345	END IF
346	soilalb_wet(:,:) = val_exp
347
348	! Mean soil albedo
349	ALLOCATE (soilalb_moy(kjpindex,2),stat=ier)
350	IF (ier.NE.0) THEN
351	WRITE (numout,) ' error in soilalb_moy allocation. We stop.We need kjpindex2 words = ',kjpindex*2
352	STOP 'condveg_init'
353	END IF
354	soilalb_moy(:,:) = val_exp
355
356	! Snow albedo of vegetative PFTs
357	ALLOCATE (albedo_snow(kjpindex),stat=ier)
358	IF (ier.NE.0) THEN
359	WRITE (numout,*) ' error in albedo_snow allocation. We stop.We need kjpindex words = ',kjpindex
360	STOP 'condveg_init'
361	END IF
362
363	! Mean vegetation albedo over visible and near-infrared range
364	ALLOCATE (albedo_glob(kjpindex),stat=ier)
365	IF (ier.NE.0) THEN
366	WRITE (numout,*) ' error in albedo_glob allocation. We stop.We need kjpindex words = ',kjpindex
367	STOP 'condveg_init'
368	END IF
369
370	! Mean bare soil albedo
371	ALLOCATE (alb_bare(kjpindex,2),stat=ier)
372	IF (ier.NE.0) THEN
373	WRITE (numout,) ' error in alb_bare allocation. We stop.We need kjpindex2 words = ',kjpindex*2
374	STOP 'condveg_init'
375	END IF
376	alb_bare(:,:) = val_exp
377
378	! Mean vegetation albedo
379	ALLOCATE (alb_veget(kjpindex,2),stat=ier)
380	IF (ier.NE.0) THEN
381	WRITE (numout,) ' error in alb_veget allocation. We stop.We need kjpindex2 words = ',kjpindex*2
382	STOP 'condveg_init'
383	END IF
384	alb_veget(:,:) = val_exp
385
386
387	!! 3. Initialise bare soil albedo
388
389	! dry soil albedo
390	var_name= 'soilalbedo_dry'
391	CALL ioconf_setatt_p('UNITS', '-')
392	CALL ioconf_setatt_p('LONG_NAME','Dry bare soil albedo')
393	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_dry, "gather", nbp_glo, index_g)
394
395	! wet soil albedo
396	var_name= 'soilalbedo_wet'
397	CALL ioconf_setatt_p('UNITS', '-')
398	CALL ioconf_setatt_p('LONG_NAME','Wet bare soil albedo')
399	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_wet, "gather", nbp_glo, index_g)
400
401	! mean soil aledo
402	var_name= 'soilalbedo_moy'
403	CALL ioconf_setatt_p('UNITS', '-')
404	CALL ioconf_setatt_p('LONG_NAME','Mean bare soil albedo')
405	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_moy, "gather", nbp_glo, index_g)
406
407	! check if we have real values and not only missing values
408	IF ( MINVAL(soilalb_wet) .EQ. MAXVAL(soilalb_wet) .AND. MAXVAL(soilalb_wet) .EQ. val_exp .OR.&
409	& MINVAL(soilalb_dry) .EQ. MAXVAL(soilalb_dry) .AND. MAXVAL(soilalb_dry) .EQ. val_exp .OR.&
410	& MINVAL(soilalb_moy) .EQ. MAXVAL(soilalb_moy) .AND. MAXVAL(soilalb_moy) .EQ. val_exp) THEN
411	CALL condveg_soilalb(kjpindex, lalo, neighbours, resolution, contfrac, soilalb_dry,soilalb_wet)
412	WRITE(numout,*) '---> val_exp ', val_exp
413	WRITE(numout,*) '---> ALBEDO_wet VIS:', MINVAL(soilalb_wet(:,ivis)), MAXVAL(soilalb_wet(:,ivis))
414	WRITE(numout,*) '---> ALBEDO_wet NIR:', MINVAL(soilalb_wet(:,inir)), MAXVAL(soilalb_wet(:,inir))
415	WRITE(numout,*) '---> ALBEDO_dry VIS:', MINVAL(soilalb_dry(:,ivis)), MAXVAL(soilalb_dry(:,ivis))
416	WRITE(numout,*) '---> ALBEDO_dry NIR:', MINVAL(soilalb_dry(:,inir)), MAXVAL(soilalb_dry(:,inir))
417	WRITE(numout,*) '---> ALBEDO_moy VIS:', MINVAL(soilalb_moy(:,ivis)), MAXVAL(soilalb_moy(:,ivis))
418	WRITE(numout,*) '---> ALBEDO_moy NIR:', MINVAL(soilalb_moy(:,inir)), MAXVAL(soilalb_moy(:,inir))
419	ENDIF
420
421	ELSE
422	WRITE (numout,*) ' l_first_condveg false . we stop '
423	STOP 'condveg_init'
424	ENDIF
425
426	IF (long_print) WRITE (numout,*) ' condveg_init done '
427
428	END SUBROUTINE condveg_init
429
430	!! ==============================================================================================================================
431	!! SUBROUTINE : condveg_clear
432	!!
433	!>\BRIEF Deallocate albedo variables
434	!!
435	!! DESCRIPTION : None
436	!!
437	!! RECENT CHANGE(S): None
438	!!
439	!! MAIN OUTPUT VARIABLE(S): None
440	!!
441	!! REFERENCES : None
442	!!
443	!! FLOWCHART : None
444	!! \n
445	!_ ================================================================================================================================
446
447	SUBROUTINE condveg_clear ()
448
449	l_first_condveg=.TRUE.
450
451	! Dry soil albedo
452	IF (ALLOCATED (soilalb_dry)) DEALLOCATE (soilalb_dry)
453	! Wet soil albedo
454	IF (ALLOCATED(soilalb_wet)) DEALLOCATE (soilalb_wet)
455	! Mean soil albedo
456	IF (ALLOCATED(soilalb_moy)) DEALLOCATE (soilalb_moy)
457	! Mean snow albedo
458	IF (ALLOCATED(albedo_snow)) DEALLOCATE (albedo_snow)
459	! Mean albedo
460	IF (ALLOCATED(albedo_glob)) DEALLOCATE (albedo_glob)
461	! Mean albedo of bare soil
462	IF (ALLOCATED(alb_bare)) DEALLOCATE (alb_bare)
463	! Mean vegetation albedo
464	IF (ALLOCATED(alb_veget)) DEALLOCATE (alb_veget)
465
466	END SUBROUTINE condveg_clear
467
468	!! ==============================================================================================================================
469	!! SUBROUTINE : condveg_var_init
470	!!
471	!>\BRIEF Initialisation of local array and calculation of emissivity,
472	!! albedo and roughness height.
473	!!
474	!! DESCRIPTION : None
475	!!
476	!! MAIN OUTPUT VARIABLE(S): None
477	!!
478	!! REFERENCE(S) : None
479	!!
480	!! FLOWCHART : None
481	!! \n
482	!_ ================================================================================================================================
483
484	SUBROUTINE condveg_var_init (kjpindex, veget, veget_max, frac_nobio, &
485	totfrac_nobio, drysoil_frac, zlev, height, &
486	circ_class_biomass, circ_class_n, emis, albedo, &
487	z0, roughheight, lai, sinang, &
488	snow, snow_age, snow_nobio, snow_nobio_age, albedo_snow, &
489	albedo_glob, lai_split,z0_veg, Isotrop_Abs_Tot_p, &
490	Isotrop_Tran_Tot_p, laieff_fit, albedo_pft, frac_snow_nobio,&
491	frac_snow_veg)
492
493	!! 0. Variable and parameter declaration
494
495	!! 0.1 Input variables
496
497	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
498	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
499	!! (m^2 m^{-2})
500	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
501	!! (= ind*cn_ind) (m^2/m^{-2})
502	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
503	!! continental ice, lakes, etc. (unitless)
504	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces, i.e.
505	!! continental ice, lakes, etc. (unitless)
506	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of dry soil in visible range (unitless)
507	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
508	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
509	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
510	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: sinang !! the sine of the solar angle from the horizon (unitless)
511	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
512	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
513	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age (days)
514	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
515	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nparts,nelements), &
516	INTENT(IN) :: circ_class_biomass !! Stem diameter
517	!! @tex $(m)$ @endtex
518	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc), INTENT(IN) :: circ_class_n !! Number of trees within each circumference
519	REAL(r_std),DIMENSION(nlevels),INTENT(IN) :: lai_split
520	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
521
522	!! 0.2 Output varialbes
523
524	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Surface emissivity (unitless)
525	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo of vegetative PFTs in visible and
526	!! near-infrared range
527	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Soil roughness height (m)
528	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Grid effective roughness height (m)
529	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
530	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_glob !! Mean albedo (unitless ratio)
531	REAL(r_std),DIMENSION (:,:,:), &
532	INTENT (out) :: Isotrop_Abs_Tot_p !! Absorbed radiation per level for photosynthesis
533	REAL(r_std),DIMENSION (:,:,:), &
534	INTENT (out) :: Isotrop_Tran_Tot_p !! Transmitted radiation per level for photosynthesis
535
536	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0_veg !! Roughness height of vegetated part (m)
537	REAL(r_std), DIMENSION(kjpindex,nvm,n_spectralbands), &
538	INTENT (out) :: albedo_pft !! Albedo (two stream radiation transfer model)
539	!! for visible and near-infrared range
540	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
541	!! (unitless ratio)
542	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: frac_snow_veg !! the fraction of ground covered with snow, between 0 and 1
543	!! for each PFT (unitless)
544	!! 0.3 Modified variables
545
546	!! 0.4 Local variables
547
548	INTEGER(i_std) :: ier !! Check errors
549	INTEGER(i_std) :: jv !! Index
550
551	!_ ================================================================================================================================
552
553
554	!! 1. Choice of parameter setting
555
556	! This switch is for choosing surface parameters.
557	! If it is set to true, the values for the soil roughness, emissivity
558	! and albedo are set to default values which are read in from the run.def.
559	! This is useful if one performs a single site simulations,
560	! it is not recommended to do so for global simulations.
561
562
563	!! 2. Calculate emissivity
564
565	! If TRUE read in default values for emissivity
566
567	IF ( impaze ) THEN
568	!
569	emis(:) = emis_scal
570	!
571	! If FALSE set emissivity to 1.
572	ELSE
573	emis_scal = un
574	emis(:) = emis_scal
575	ENDIF
576
577	!! 3. Calculate roughness height
578
579	! Chooses between two methods to calculate the grid average of the roughness.
580	! If set to true: The grid average is calculated by averaging the drag coefficients over PFT.
581	! If set to false: The grid average is calculated by averaging the logarithm of the roughness length per PFT.
582
583	! TRUE read in default values for roughness height
584	IF ( impaze ) THEN
585
586	z0(:) = z0_scal
587	roughheight(:) = roughheight_scal
588
589	! If FALSE calculate roughness height
590	ELSE
591
592	IF ( z0cdrag_ave ) THEN
593	CALL condveg_z0cdrag(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
594	& height, z0, roughheight, z0_veg)
595	ELSE
596	CALL condveg_z0logz(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
597	& z0, roughheight, z0_veg)
598	ENDIF
599
600	ENDIF
601
602	!! 4. Calculate albedo
603
604	! If TRUE read in default values for albedo
605	IF ( impaze ) THEN
606	!
607	albedo(:,ivis) = albedo_scal(ivis)
608	albedo(:,inir) = albedo_scal(inir)
609	!
610	ELSE
611
612	! check to see which type of albedo we are using
613	SELECTCASE (albedo_type)
614	CASE ('standard')
615	CALL albedo_calc (kjpindex, drysoil_frac, veget,&
616	soilalb_dry, soilalb_wet, soilalb_moy, tot_bare_soil, alb_veget, alb_bare, albedo)
617	! Initialize these variables here. They are only computed in the two_stream model.
618	Isotrop_Abs_Tot_p(:,:,:)=undef_sechiba
619	Isotrop_Tran_Tot_p(:,:,:)=undef_sechiba
620	frac_snow_nobio(:,:)=undef_sechiba
621	frac_snow_veg(:,:)=undef_sechiba
622
623	CASE ('pinty')
624	! We need to initialize albedo here to give a realistic guess to the coupled
625	! model.
626	CALL albedo_two_stream(kjpindex, nlevels_tot, drysoil_frac, lai, veget_max, &
627	sinang, soilalb_dry, soilalb_wet, frac_nobio, soilalb_moy, &
628	alb_bare, albedo, snow, snow_age, snow_nobio, &
629	snow_nobio_age, albedo_snow, albedo_glob, &
630	circ_class_biomass, circ_class_n,&
631	lai_split, z0_veg, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
632	laieff_fit, albedo_pft, frac_snow_nobio, frac_snow_veg)
633
634
635	CASE DEFAULT
636	WRITE(numout,*) "Unsupported albedo type. This should have been caught in intersurf!"
637	STOP 'condveg'
638	END SELECT
639
640
641	ENDIF
642
643	!! Calculate snow albedo. This is tricky because the new scheme requires information about snow to be
644	!! passed to it, while the old scheme uses existing albedo data
645	SELECTCASE (albedo_type)
646	CASE ('standard')
647	CALL calculate_surface_albedo_with_snow(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
648	snow, snow_age, snow_nobio, snow_nobio_age, tot_bare_soil, albedo, albedo_snow, albedo_glob)
649	CASE ('pinty')
650	CASE DEFAULT
651	WRITE(numout,*) "Unsupported albedo type. This should have been caught in intersurf!"
652	STOP 'condveg'
653	END SELECT
654
655	IF (long_print) WRITE (numout,*) ' condveg_var_init done '
656
657	END SUBROUTINE condveg_var_init
658
659
660
661	!! ==============================================================================================================================
662	!! SUBROUTINE : condveg_var_update
663	!!
664	!>\BRIEF Update emissivity, albedo and roughness height.
665	!!
666	!! DESCRIPTION : None
667	!!
668	!! MAIN OUTPUT VARIABLE(S): \n
669	!!
670	!! REFERENCE(S) : None
671	!!
672	!! FLOWCHART : None
673	!! \n
674	!_ ================================================================================================================================
675
676	SUBROUTINE condveg_var_update (kjpindex, veget, veget_max, frac_nobio, &
677	totfrac_nobio, drysoil_frac, zlev, height, &
678	circ_class_biomass, circ_class_n, emis, albedo, &
679	z0, roughheight, lai, sinang, &
680	snow, snow_age, snow_nobio, snow_nobio_age, albedo_snow, &
681	albedo_glob, lai_split, z0_veg, Isotrop_Abs_Tot_p, Isotrop_Tran_Tot_p, &
682	laieff_fit, albedo_pft, frac_snow_nobio, frac_snow_veg)
683
684	!! 0. Variable and parameter declaration
685
686	!! 0.1 Input variables
687
688	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
689	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
690	!! (m^2 m^{-2})
691	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
692	!! (= ind*cn_ind) (m^2 m^{-2})
693	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
694	!! continental ice, lakes, etc. (unitless)
695	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces, i.e.
696	!! continental ice, lakes, etc. (unitless)
697	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of dry soil in visible range (unitless)
698	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
699	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
700	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
701	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: sinang !! the sine of the solar angle from the horizon (unitless)
702	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
703	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
704	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age (days)
705	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
706	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc,nparts,nelements), &
707	INTENT(IN) :: circ_class_biomass !! Stem diameter
708	!! @tex $(m)$ @endtex
709	REAL(r_std), DIMENSION(kjpindex,nvm,ncirc), INTENT(IN) :: circ_class_n !! Number of trees within each circumference
710	REAL(r_std),DIMENSION(nlevels),INTENT(IN) :: lai_split
711	TYPE(laieff_type),DIMENSION (:,:,:),INTENT(in) :: laieff_fit !! Fitted parameters for the effective LAI
712
713
714	!! 0.2 Output variables
715
716	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity (unitless)
717	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo of vegetative PFTs in visible and
718	!! near-infrared range
719	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Roughness height (m)
720	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Grid effective roughness height (m)
721	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
722	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_glob !! Mean albedo (unitless ratio)
723	REAL(r_std),DIMENSION (:,:,:), &
724	INTENT (out) :: Isotrop_Abs_Tot_p !!Absorbed radiation per level for photosynthesis
725	REAL(r_std),DIMENSION (:,:,:), &
726	INTENT (out) :: Isotrop_Tran_Tot_p !!Transmitted radiation per level for photosynthesis
727
728	REAL(r_std), DIMENSION(kjpindex), INTENT (out) :: z0_veg !! Roughness height of vegetated part (m)
729	REAL(r_std), DIMENSION(kjpindex,nvm,n_spectralbands), &
730	INTENT (out) :: albedo_pft !! Albedo (two stream radiation transfer model)
731	!! for visible and near-infrared range
732	!! for each PFT (unitless)
733	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out) :: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
734	!! (unitless ratio)
735	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: frac_snow_veg !! the fraction of ground covered with snow, between 0 and 1
736
737	!! 0.3 Modified variables
738
739	!! 0.4 Local variables
740	INTEGER(i_std) :: ji, jv !! Indeces
741	!_ ================================================================================================================================
742
743	!! 1. Calculate emissivity
744
745	emis(:) = emis_scal
746
747	!! 2. Calculate roughness height
748
749	! If TRUE read in prescribed values for roughness height
750	IF ( impaze ) THEN
751
752	DO ji = 1, kjpindex
753	z0(ji) = z0_scal
754	roughheight(ji) = roughheight_scal
755	ENDDO
756
757	! Calculate roughness height
758	ELSE
759
760	IF ( z0cdrag_ave ) THEN
761	CALL condveg_z0cdrag (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, height, &
762	& z0, roughheight, z0_veg)
763	ELSE
764	CALL condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
765	& z0, roughheight, z0_veg)
766	ENDIF
767
768	ENDIF
769
770	!! 3. Calculate albedo
771	! This has to be done after the calculation of the roughness height, since it uses
772	! the roughness height for calculating the snow coverage fraction.
773
774	! If TRUE read in prescribed values for albedo
775	IF ( impaze ) THEN
776
777	albedo(:,ivis) = albedo_scal(ivis)
778	albedo(:,inir) = albedo_scal(inir)
779
780	! If FALSE calculate albedo from ORCHIDEE default values
781	ELSE
782
783	! check to see which type of albedo we are using
784	SELECTCASE (albedo_type)
785	CASE ('standard')
786	CALL albedo_calc (kjpindex, drysoil_frac, veget,&
787	soilalb_dry, soilalb_wet, soilalb_moy, tot_bare_soil,alb_veget, alb_bare, albedo)
788	! Initialize these variables here. They are only computed in the two_stream model.
789	Isotrop_Abs_Tot_p(:,:,:)=undef_sechiba
790	Isotrop_Tran_Tot_p(:,:,:)=undef_sechiba
791	frac_snow_nobio(:,:)=undef_sechiba
792	frac_snow_veg(:,:)=undef_sechiba
793	CASE ('pinty')
794
795
796	CALL albedo_two_stream(kjpindex, nlevels_tot, drysoil_frac, lai, veget_max, &
797	sinang, soilalb_dry, soilalb_wet, frac_nobio, soilalb_moy, &
798	alb_bare, albedo, snow, snow_age, snow_nobio, &
799	snow_nobio_age, albedo_snow, albedo_glob, &
800	circ_class_biomass, circ_class_n, lai_split, z0_veg, Isotrop_Abs_Tot_p,&
801	Isotrop_Tran_Tot_p, laieff_fit, albedo_pft, frac_snow_nobio,&
802	frac_snow_veg)
803
804
805	CASE DEFAULT
806	WRITE(numout,*) "Unsupported albedo type. This should have been caught in intersurf!"
807	STOP 'condveg'
808	END SELECT
809
810	ENDIF
811
812	!! Calculate snow albedo.
813	! The scheme of Pinty et al calculates the snow albedo in the albedo_two_stream routine,
814	! since it changes the background albedo. The old scheme calculates it separately.
815	SELECTCASE (albedo_type)
816	CASE ('standard')
817	CALL calculate_surface_albedo_with_snow(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
818	snow, snow_age, snow_nobio, snow_nobio_age, tot_bare_soil, albedo, albedo_snow, albedo_glob)
819	CASE ('pinty')
820	CASE DEFAULT
821	WRITE(numout,*) "Unsupported albedo type. This should have been caught in intersurf!"
822	STOP 'condveg'
823	END SELECT
824
825	IF (long_print) WRITE (numout,*) ' condveg_var_update done '
826
827	END SUBROUTINE condveg_var_update
828
829
830	!! ==============================================================================================================================
831	!! SUBROUTINE : condveg_soilalb
832	!!
833	!>\BRIEF This subroutine calculates the albedo of soil (without snow).
834	!!
835	!! DESCRIPTION This subroutine reads the soil colour maps in 1 x 1 deg resolution
836	!! from the Henderson-Sellers & Wilson database. These values are interpolated to
837	!! the model's resolution and transformed into
838	!! dry and wet albedos.\n
839	!!
840	!! If the soil albedo is calculated without the dependence of soil moisture, the
841	!! soil colour values are transformed into mean soil albedo values.\n
842	!!
843	!! The calculations follow the assumption that the grid of the data is regular and
844	!! it covers the globe. The calculation for the model grid are based on the borders
845	!! of the grid of the resolution.
846	!!
847	!! RECENT CHANGE(S): None
848	!!
849	!! MAIN OUTPUT VARIABLE(S): soilalb_dry for visible and near-infrared range,
850	!! soilalb_wet for visible and near-infrared range,
851	!! soilalb_moy for visible and near-infrared range
852	!!
853	!! REFERENCE(S) :
854	!! -Wilson, M.F., and A. Henderson-Sellers, 1985: A global archive of land cover and
855	!! soils data for use in general circulation climate models. J. Clim., 5, 119-143.
856	!!
857	!! FLOWCHART : None
858	!! \n
859	!_ ================================================================================================================================
860
861	SUBROUTINE condveg_soilalb(nbpt, lalo, neighbours, resolution, contfrac, soilalb_dry,soilalb_wet)
862
863	!! 0. Variable and parameter declaration
864
865	!! 0.1 Input variables
866
867	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
868	!! interpolated (unitless)
869	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
870	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point
871	!! (1=N, 2=E, 3=S, 4=W)
872	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (m)
873	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
874	REAL(r_std), INTENT(inout) :: soilalb_dry(nbpt,2) !! Albedo of the dry bare soil (unitless)
875	REAL(r_std), INTENT(inout) :: soilalb_wet(nbpt,2) !! Albedo of the wet dry soil (unitless)
876
877	!! 0.2 Output variables
878
879	!! 0.3 Modified variables
880
881	!! 0.4 Local variables
882
883	INTEGER(i_std) :: nbvmax !! nbvmax for interpolation (unitless). It is the
884	!! dimension of the variables in which we store the list
885	!! of points of the source grid which fit into one grid
886	!! box of the target.
887	CHARACTER(LEN=80) :: filename !! Filename of soil colour map
888	INTEGER(i_std) :: iml, jml, lml, &
889	&tml, fid, ib, ip, jp, fopt, ilf, lastjp, nbexp !! Indices
890	REAL(r_std) :: lev(1), date, dt !! Help variables to read in file data
891	INTEGER(i_std) :: itau(1) !! Help variables to read in file data
892	REAL(r_std) :: sgn !! Help variable to compute average bare soil albedo
893	REAL(r_std) :: coslat !! [DISPENSABLE]
894	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel !! Help variable to read file data and allocate memory
895	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon_rel !! Help variable to read file data and allocate memory
896	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: soilcol !! Help variable to read file data and allocate memory
897	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! Help variable to read file data and allocate memory
898	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! Help variable to read file data and allocate memory
899	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask !! Help variable to read file data and allocate memory
900	CHARACTER(LEN=30) :: callsign !! Help variable to read file data and allocate memory
901	LOGICAL :: ok_interpol = .FALSE. !! Optional return of aggregate_2d
902	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
903	!_ ================================================================================================================================
904
905	!! 1. Open file and allocate memory
906
907	! Open file with soil colours
908
909	!Config Key = SOILALB_FILE
910	!Config Desc = Name of file from which the bare soil albedo
911	!Config Def = soils_param.nc
912	!Config If = NOT(IMPOSE_AZE)
913	!Config Help = The name of the file to be opened to read the soil types from
914	!Config which we derive then the bare soil albedos. This file is 1x1
915	!Config deg and based on the soil colors defined by Wilson and Henderson-Seller.
916	!Config Units = [FILE]
917	!
918	filename = 'soils_param.nc'
919	CALL getin_p('SOILALB_FILE',filename)
920
921	! Read data from file
922	IF (is_root_prc) CALL flininfo(filename,iml, jml, lml, tml, fid)
923	CALL bcast(iml)
924	CALL bcast(jml)
925	CALL bcast(lml)
926	CALL bcast(tml)
927
928	! Allocate memory for latitudes
929	ALLOC_ERR=-1
930	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
931	IF (ALLOC_ERR/=0) THEN
932	WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
933	STOP
934	ENDIF
935
936	! Allcoate memory for longitude
937	ALLOC_ERR=-1
938	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
939	IF (ALLOC_ERR/=0) THEN
940	WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
941	STOP
942	ENDIF
943
944	! Allocate memory for mask
945	ALLOC_ERR=-1
946	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
947	IF (ALLOC_ERR/=0) THEN
948	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
949	STOP
950	ENDIF
951
952	! Allocate memory for soil data
953	ALLOC_ERR=-1
954	ALLOCATE(soilcol(iml,jml), STAT=ALLOC_ERR)
955	IF (ALLOC_ERR/=0) THEN
956	WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
957	STOP
958	ENDIF
959
960	! Set values
961	IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
962	CALL bcast(lon_rel)
963	CALL bcast(lat_rel)
964	CALL bcast(lev)
965	CALL bcast(itau)
966	CALL bcast(date)
967	CALL bcast(dt)
968
969	IF (is_root_prc) CALL flinget(fid, 'soilcolor', iml, jml, lml, tml, 1, 1, soilcol)
970	CALL bcast(soilcol)
971
972	IF (is_root_prc) CALL flinclo(fid)
973
974	! Create mask with values of soil colour
975	mask(:,:) = zero
976	DO ip=1,iml
977	DO jp=1,jml
978	IF (soilcol(ip,jp) > min_sechiba) THEN
979	mask(ip,jp) = un
980	ENDIF
981	ENDDO
982	ENDDO
983
984	! Set nbvmax to 200 for interpolation
985	! This number is the dimension of the variables in which we store
986	! the list of points of the source grid which fit into one grid box of the target.
987	nbvmax = 200
988
989	callsign = 'Soil color map'
990
991	! Start with interpolation
992	DO WHILE ( .NOT. ok_interpol )
993	WRITE(numout,*) "Projection arrays for ",callsign," : "
994	WRITE(numout,*) "nbvmax = ",nbvmax
995
996	ALLOC_ERR=-1
997	ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
998	IF (ALLOC_ERR/=0) THEN
999	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
1000	STOP
1001	ENDIF
1002	sub_area(:,:)=zero
1003	ALLOC_ERR=-1
1004	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
1005	IF (ALLOC_ERR/=0) THEN
1006	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
1007	STOP
1008	ENDIF
1009	sub_index(:,:,:)=0
1010
1011	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
1012	& iml, jml, lon_rel, lat_rel, mask, callsign, &
1013	& nbvmax, sub_index, sub_area, ok_interpol)
1014
1015	IF ( .NOT. ok_interpol ) THEN
1016	DEALLOCATE(sub_area)
1017	DEALLOCATE(sub_index)
1018	ENDIF
1019
1020	nbvmax = nbvmax * 2
1021	ENDDO
1022
1023	! Check how many points with soil information are found
1024	nbexp = 0
1025
1026	soilalb_dry(:,:) = zero
1027	soilalb_wet(:,:) = zero
1028	soilalb_moy(:,:) = zero
1029
1030	DO ib=1,nbpt ! Loop over domain size
1031
1032	fopt = COUNT(sub_area(ib,:) > zero)
1033
1034	!! 3. Compute the average bare soil albedo parameters
1035
1036	IF ( fopt .EQ. 0) THEN ! If no points were interpolated
1037	nbexp = nbexp + 1
1038	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1039	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1040	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1041	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1042	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
1043	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
1044	ELSE
1045	sgn = zero
1046
1047	DO ilf = 1,fopt ! If points were interpolated
1048
1049	ip = sub_index(ib,ilf,1)
1050	jp = sub_index(ib,ilf,2)
1051
1052	! Weighted albedo values by interpolation area
1053	IF ( NINT(soilcol(ip,jp)) .LE. classnb) THEN
1054	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1055	soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1056	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1057	soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1058	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1059	soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
1060	sgn = sgn + sub_area(ib,ilf)
1061	ELSE
1062	WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
1063	STOP
1064	ENDIF
1065
1066	ENDDO
1067
1068	! Normalize the surface
1069	IF ( sgn .LT. min_sechiba) THEN
1070	nbexp = nbexp + 1
1071	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1072	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1073	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1074	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1075	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
1076	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
1077	ELSE
1078	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis)/sgn
1079	soilalb_dry(ib,inir) = soilalb_dry(ib,inir)/sgn
1080	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis)/sgn
1081	soilalb_wet(ib,inir) = soilalb_wet(ib,inir)/sgn
1082	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis)/sgn
1083	soilalb_moy(ib,inir) = soilalb_moy(ib,inir)/sgn
1084	ENDIF
1085
1086	ENDIF
1087
1088	ENDDO
1089
1090	IF ( nbexp .GT. 0 ) THEN
1091	WRITE(numout,*) 'CONDVEG_soilalb : The interpolation of the bare soil albedo had ', nbexp
1092	WRITE(numout,*) 'CONDVEG_soilalb : points without data. This are either coastal points or'
1093	WRITE(numout,*) 'CONDVEG_soilalb : ice covered land.'
1094	WRITE(numout,*) 'CONDVEG_soilalb : The problem was solved by using the average of all soils'
1095	WRITE(numout,*) 'CONDVEG_soilalb : in dry and wet conditions'
1096	ENDIF
1097
1098	DEALLOCATE (lat_rel)
1099	DEALLOCATE (lon_rel)
1100	DEALLOCATE (mask)
1101	DEALLOCATE (sub_index)
1102	DEALLOCATE (sub_area)
1103	DEALLOCATE (soilcol)
1104
1105	RETURN
1106
1107	END SUBROUTINE condveg_soilalb
1108
1109
1110	!! ==============================================================================================================================
1111	!! SUBROUTINE : condveg_z0logz
1112	!!
1113	!>\BRIEF Computation of grid average of roughness height by averaging the
1114	!! logarithm of the roughness height of each grid box components fracbio and fracnobio.
1115	!!
1116	!! DESCRIPTION : Calculates mean roughness height
1117	!! over the grid cell. The mean roughness height is derived from the vegetation
1118	!! height which is scaled by the roughness parameter. The sum of the logarithm of the
1119	!! roughness times the fraction per grid cell gives the average roughness height per
1120	!! grid cell for the vegetative PFTs. The roughness height for the non-vegetative PFTs
1121	!! is calculated in a second step. \n
1122	!!
1123	!! To compute the fluxes,
1124	!! the difference between the height of the vegetation and the zero plane displacement height
1125	!! is needed and called roughheight .\n
1126	!!
1127	!! RECENT CHANGE(S): None
1128	!!
1129	!! MAIN OUTPUT VARIABLE(S): roughness height (z0), grid effective roughness height (roughheight)
1130	!!
1131	!! REFERENCE(S) :
1132	!!
1133	!! FLOWCHART : None
1134	!! \n
1135	!_ ================================================================================================================================
1136
1137	SUBROUTINE condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
1138	& z0, roughheight, z0_veg)
1139
1140	!! 0. Variable and parameter declaration
1141
1142	!! 0.1 Input variables
1143
1144	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1145	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1146	!! (m^2 m^{-2})
1147	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1148	!! (= ind*cn_ind) (m^2 m^{-2})
1149	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1150	!! i.e. continental ice, lakes, etc. (unitless)
1151	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1152	!! i.e. continental ice, lakes, etc. (unitless)
1153	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1154
1155	!! 0.2 Output variables
1156
1157	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0 !! Soil roughness height (m)
1158	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1159	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0_veg !! Roughness height of vegetated part (m)
1160
1161	!! 0.3 Modified variables
1162
1163	!! 0.4 Local variables
1164
1165	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1166	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1167	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1168	REAL(r_std), DIMENSION(kjpindex) :: d_veg !! PFT coverage of vegetative PFTs
1169	!! (= ind*cn_ind) (m^2 m^{-2})
1170	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1171	REAL(r_std) :: z0_nobio !! Roughness of non-vegetative fraction (m),
1172	!! i.e. continental ice, lakes, etc.
1173	!_ ================================================================================================================================
1174
1175	!! 1. Preliminary calculation
1176
1177	! Calculate the roughness (m) of bare soil, z0_bare
1178	! taken from constantes_veg.f90
1179	z0(:) = tot_bare_soil(:) * LOG(z0_bare)
1180
1181	! Define fraction of bare soil
1182	sumveg(:) = tot_bare_soil(:)
1183
1184	! Set average vegetation height to zero
1185	ave_height(:) = zero
1186
1187	!! 2. Calculate the mean roughness length
1188
1189	! Calculate the mean roughness height of
1190	! vegetative PFTs over the grid cell
1191	DO jv = 2, nvm !Loop over # vegetative PFTs
1192
1193	! In the case of forest, use parameter veget_max because
1194	! tree trunks influence the roughness even when there are no leaves
1195	IF ( is_tree(jv) ) THEN
1196	d_veg(:) = veget_max(:,jv)
1197	ELSE
1198
1199	! In the case of grass, use parameter veget because grasses
1200	! only influence the roughness during the growing season
1201	d_veg(:) = veget(:,jv)
1202	ENDIF
1203
1204	! Calculate the average roughness over the grid cell:
1205	! The roughness for vegetative PFTs is calculated by
1206	! the vegetation height per PFT multiplied by the roughness
1207	! parameter 'z0_over_height= 1/16'. If this scaled value is
1208	! lower than 0.01 than the value for the roughness length
1209	! of bare soil (0.01) is used. The sum of the logarithm of
1210	! the roughness times the fraction per grid cell gives the
1211	! logarithm of roughness length per grid cell for the vegetative
1212	! PFTs.
1213	z0(:) = z0(:) + d_veg(:) * &
1214	LOG( MAX(height(:,jv)*z0_over_height,z0_bare) )
1215	! Sum of bare soil and fraction vegetated fraction
1216	sumveg(:) = sumveg(:) + d_veg(:)
1217
1218	! Weighted height of vegetation with maximal cover fraction
1219	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1220
1221	ENDDO !Loop over # vegetative PFTs
1222
1223	!! 3. Calculate the mean roughness length of non-vegetative surfaces \n
1224
1225	! Search for pixels with vegetated part to normalise
1226	! roughness length
1227	WHERE ( sumveg(:) > zero ) z0(:) = z0(:) / sumveg(:)
1228
1229	! Calculate fraction of roughness for vegetated part
1230	z0(:) = (un - totfrac_nobio(:)) * z0(:)
1231	! Save roughness of vegetated part for calculation of snow fraction
1232	z0_veg(:) = z0(:)
1233
1234	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1235
1236	! Set rougness for ice
1237	IF ( jv .EQ. iice ) THEN
1238	z0_nobio = z0_ice
1239	ELSE
1240	WRITE(numout,*) 'jv=',jv
1241	STOP 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1242	ENDIF
1243
1244	! Sum of vegetative roughness length and non-vegetative
1245	! roughness length
1246	z0(:) = z0(:) + frac_nobio(:,jv) * LOG(z0_nobio)
1247
1248
1249
1250	ENDDO ! loop over # of non-vegative surfaces
1251
1252	!! 4. Calculate the zero plane displacement height and effective roughness length
1253
1254	! Take the exponential of the roughness length
1255	z0(:) = EXP( z0(:) )
1256
1257	! Compute the zero plane displacement height which
1258	! is an equivalent height of the vegetation for the absorption of momentum
1259	zhdispl(:) = ave_height(:) * height_displacement
1260
1261	! Then we compute what we call the grid effective roughness height.
1262	! This is the height over which the roughness acts. It combines the
1263	! zero plane displacement height and the vegetation height. This
1264	! effective value is the difference between the height of the
1265	! vegetation and the zero plane displacement height.
1266	roughheight(:) = ave_height(:) - zhdispl(:)
1267
1268
1269	END SUBROUTINE condveg_z0logz
1270
1271
1272	!! ==============================================================================================================================
1273	!! SUBROUTINE : condveg_z0cdrag
1274	!!
1275	!>\BRIEF Computation of grid average of roughness length by calculating
1276	!! the drag coefficient.
1277	!!
1278	!! DESCRIPTION : This routine calculates the mean roughness height and mean
1279	!! effective roughness height over the grid cell. The mean roughness height (z0)
1280	!! is computed by averaging the drag coefficients \n
1281	!!
1282	!! \latexonly
1283	!! \input{z0cdrag1.tex}
1284	!! \endlatexonly
1285	!! \n
1286	!!
1287	!! where C is the drag coefficient at the height of the vegetation, kappa is the
1288	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
1289	!! The reference level for z needs to be high enough above the canopy to avoid
1290	!! singularities of the LOG. This height is set to minimum 10m above ground.
1291	!! The drag coefficient increases with roughness height to represent the greater
1292	!! turbulence generated by rougher surfaces.
1293	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
1294	!!
1295	!! The roughness height for the non-vegetative surfaces is calculated in a second step.
1296	!! In order to calculate the transfer coefficients the
1297	!! effective roughness height is calculated. This effective value is the difference
1298	!! between the height of the vegetation and the zero plane displacement height.\nn
1299	!!
1300	!! RECENT CHANGE(S): None
1301	!!
1302	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
1303	!!
1304	!! REFERENCE(S) : None
1305	!!
1306	!! FLOWCHART : None
1307	!! \n
1308	!_ ================================================================================================================================
1309
1310	SUBROUTINE condveg_z0cdrag (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, &
1311	& z0, roughheight,z0_veg)
1312
1313	!! 0. Variable and parameter declaration
1314
1315	!! 0.1 Input variables
1316
1317	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1318	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1319	!! (m^2 m^{-2})
1320	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1321	!! (= ind*cn_ind) (m^2 m^{-2})
1322	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1323	!! i.e. continental ice, lakes, etc. (unitless)
1324	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1325	!! i.e. continental ice, lakes, etc. (unitless)
1326	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
1327	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1328
1329	!! 0.2 Output variables
1330
1331	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0 !! Roughness height (m)
1332	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1333	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0_veg !! Roughness height of vegetated part (m)
1334
1335
1336	!! 0.3 Modified variables
1337
1338	!! 0.4 Local variables
1339
1340	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1341	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1342	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
1343	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1344	REAL(r_std), DIMENSION(kjpindex) :: d_veg !! PFT coverage of vegetative PFTs
1345	!! (= ind*cn_ind) (m^2 m^{-2})
1346	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1347	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
1348	!! i.e. continental ice, lakes, etc.
1349	!_ ================================================================================================================================
1350
1351	!! 1. Preliminary calculation
1352
1353	! Set maximal height of first layer
1354	ztmp(:) = MAX(10., zlev(:))
1355
1356	! Calculate roughness for non-vegetative surfaces
1357	! with the von Karman constant
1358	z0(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/z0_bare))**2
1359
1360	! Fraction of bare soil
1361	sumveg(:) = tot_bare_soil(:)
1362
1363	! Set average vegetation height to zero
1364	ave_height(:) = zero
1365
1366	!! 2. Calculate the mean roughness height
1367
1368	! Calculate the mean roughness height of
1369	! vegetative PFTs over the grid cell
1370	DO jv = 2, nvm
1371
1372	! In the case of forest, use parameter veget_max because
1373	! tree trunks influence the roughness even when there are no leaves
1374	IF ( is_tree(jv) ) THEN
1375	! In the case of grass, use parameter veget because grasses
1376	! only influence the roughness during the growing season
1377	d_veg(:) = veget_max(:,jv)
1378	ELSE
1379	! grasses only have an influence if they are really there!
1380	d_veg(:) = veget(:,jv)
1381	ENDIF
1382
1383	! Calculate the average roughness over the grid cell:
1384	! The unitless drag coefficient is per vegetative PFT
1385	! calculated by use of the von Karman constant, the height
1386	! of the first layer and the roughness. The roughness
1387	! is calculated as the vegetation height per PFT
1388	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
1389	! If this scaled value is lower than 0.01 then the value for
1390	! the roughness of bare soil (0.01) is used.
1391	! The sum over all PFTs gives the average roughness
1392	! per grid cell for the vegetative PFTs.
1393	z0(:) = z0(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)z0_over_height,z0_bare)))*2
1394
1395	! Sum of bare soil and fraction vegetated fraction
1396	sumveg(:) = sumveg(:) + d_veg(:)
1397
1398	! Weigh height of vegetation with maximal cover fraction
1399	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1400
1401	ENDDO
1402
1403	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1404
1405	! Search for pixels with vegetated part to normalise
1406	! roughness height
1407	WHERE ( sumveg(:) .GT. zero ) z0(:) = z0(:) / sumveg(:)
1408
1409	! Calculate fraction of roughness for vegetated part
1410	z0(:) = (un - totfrac_nobio(:)) * z0(:)
1411	! Save roughness of vegetated part for calculation of snow fraction
1412	z0_veg(:) = z0(:)
1413
1414
1415	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1416
1417	! Set rougness for ice
1418	IF ( jv .EQ. iice ) THEN
1419	z0_nobio = z0_ice
1420
1421
1422	ELSE
1423	WRITE(numout,*) 'jv=',jv
1424	STOP 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1425	ENDIF
1426
1427	! Sum of vegetative roughness length and non-vegetative
1428	! roughness length
1429	z0(:) = z0(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1430
1431	ENDDO ! Loop over # of non-vegative surfaces
1432
1433
1434	!! 4. Calculate the zero plane displacement height and effective roughness length
1435
1436	! Take the exponential of the roughness
1437	z0(:) = ztmp(:) / EXP(ct_karman/SQRT(z0(:)))
1438
1439	! Compute the zero plane displacement height which
1440	! is an equivalent height for the absorption of momentum
1441	zhdispl(:) = ave_height(:) * height_displacement
1442
1443	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1444	! This is the height over which the roughness acts. It combines the
1445	! zero plane displacement height and the vegetation height.
1446	roughheight(:) = ave_height(:) - zhdispl(:)
1447
1448
1449	END SUBROUTINE condveg_z0cdrag
1450
1451
1452	END MODULE condveg

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

Download in other formats: