Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

condveg.f90

Last change on this file was 5476, checked in by albert.jornet, 6 years ago
Merge: revisions from [5320:5363/trunk/ORCHIDEE]
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 92.8 KB

Line
1	! ===============================================================================================================================
2	! MODULE : condveg
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Initialise, compute and update the surface parameters emissivity,
10	!! roughness and albedo.
11	!!
12	!! \n DESCRIPTION : The module uses 3 settings to control its flow:\n
13	!! 1. :: rough_dyn to choose between two methods to calculate
14	!! the roughness height. If set to false: the roughness height is calculated by the old formulation
15	!! which does not distinguish between z0m and z0h and which does not vary with LAI
16	!! If set to true: the grid average is calculated by the formulation proposed by Su et al. (2001)
17	!! 2. :: impaze for choosing surface parameters. If set to false, the values for the
18	!! soil albedo, emissivity and roughness height are set to default values which are read from
19	!! the run.def. If set to true, the user imposes its own values, fixed for the grid point. This is useful if
20	!! one performs site simulations, however,
21	!! it is not recommended to do so for spatialized simulations.
22	!! roughheight_scal imposes the roughness height in (m) ,
23	!! same for emis_scal (in %), albedo_scal (in %), zo_scal (in m)
24	!! Note that these values are only used if 'impaze' is true.\n
25	!! 3. :: alb_bare_model for choosing values of bare soil albedo. If set to TRUE bare
26	!! soil albedo depends on soil wetness. If set to FALSE bare soil albedo is the mean
27	!! values of wet and dry soil albedos.\n
28	!! The surface fluxes are calculated between two levels: the atmospheric level reference and the effective roughness height
29	!! defined as the difference between the mean height of the vegetation and the displacement height (zero wind
30	!! level). Over bare soils, the zero wind level is equal to the soil roughness. Over vegetation, the zero wind level
31	!! is increased by the displacement height
32	!! which depends on the height of the vegetation. For a grid point composed of different types of vegetation,
33	!! an effective surface roughness has to be calculated
34	!!
35	!! RECENT CHANGE(S): Added option rough_dyn and subroutine condveg_z0cdrag_dyn. Removed subroutine condveg_z0logz. June 2016.
36	!!
37	!! REFERENCES(S) : None
38	!!
39	!! SVN :
40	!! $HeadURL$
41	!! $Date$
42	!! $Revision$
43	!! \n
44	!_ ================================================================================================================================
45
46	MODULE condveg
47
48	USE ioipsl
49	USE xios_orchidee
50	!
51	! modules used :
52	USE constantes
53	USE constantes_soil
54	USE pft_parameters
55	USE qsat_moisture
56	USE interpol_help
57	USE mod_orchidee_para
58	USE ioipsl_para
59	USE sechiba_io_p
60	USE grid
61
62	IMPLICIT NONE
63
64	PRIVATE
65	PUBLIC :: condveg_main, condveg_initialize, condveg_finalize, condveg_clear
66
67	!
68	! variables used inside condveg module
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 :: soilalb_bg(:,:) !! Albedo values for the background bare soil (unitless)
79	!$OMP THREADPRIVATE(soilalb_bg)
80	INTEGER, SAVE :: printlev_loc !! Output debug level
81	!$OMP THREADPRIVATE(printlev_loc)
82
83	CONTAINS
84
85	!! =============================================================================================================================
86	!! SUBROUTINE : condveg_initialize
87	!!
88	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
89	!!
90	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values.
91	!! condveg_snow is called to initialize corresponding variables.
92	!!
93	!! RECENT CHANGE(S) : None
94	!!
95	!! MAIN OUTPUT VARIABLE(S)
96	!!
97	!! REFERENCE(S) : None
98	!!
99	!! FLOWCHART : None
100	!! \n
101	!_ ==============================================================================================================================
102	SUBROUTINE condveg_initialize (kjit, kjpindex, index, veget, &
103	veget_max, frac_nobio, totfrac_nobio, &
104	lalo, neighbours, resolution, contfrac, rest_id, &
105	zlev, drysoil_frac, height, snowdz, snowrho, tot_bare_soil, &
106	snow, snow_age, snow_nobio, snow_nobio_age, &
107	temp_air, pb, u, v, lai, &
108	emis, albedo, z0m, z0h, roughheight, roughheight_pft, &
109	frac_snow_veg, frac_snow_nobio)
110
111	!! 0. Variable and parameter declaration
112
113	!! 0.1. Input variables
114
115	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
116	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
117	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Index for the points on the map (unitless)
118	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in):: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
119	!! (m^2 m^{-2})
120	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates (degree)
121	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours !! Neighbouring land grid cell
122	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size of grid in x and y direction (m)
123	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: contfrac !! Fraction of land in each grid box
124	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier
125
126	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
127	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
128	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
129	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
130	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
131	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowdz !! Snow depth at each snow layer
132	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowrho !! Snow density at each snow layer
133	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
134
135	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
136	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
137	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
138	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
139	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
140	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
141	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
142	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u !! Horizontal wind speed, u direction
143	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: v !! Horizontal wind speed, v direction
144	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
145
146	!! 0.2. Output variables
147	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
148	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
149	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m !! Roughness for momentum (m)
150	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h !! Roughness for heat (m)
151	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
152	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (out) :: roughheight_pft
153	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: frac_snow_veg !! Snow cover fraction on vegeted area
154	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Snow cover fraction on non-vegeted area
155
156	!! 0.3 Modified variables
157
158	!! 0.4 Local variables
159	INTEGER(i_std) :: ji, jv !! Index
160	INTEGER(i_std) :: ier !! Check errors in memory allocation
161	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
162
163	REAL(r_std), DIMENSION(kjpindex,2) :: albedo_snow !! Snow albedo for visible and near-infrared range(unitless)
164	REAL(r_std), DIMENSION(kjpindex,2) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
165	!! range (unitless)
166	REAL(r_std), DIMENSION(kjpindex,2) :: alb_veget !! Mean vegetation albedo for visible and near-infrared
167	! !! range (unitless)
168	!_ ================================================================================================================================
169
170	!! 1. Choice of calculation of snow albedo and soil albedo
171	IF (.NOT. l_first_condveg) CALL ipslerr_p(3,'condveg_initialize','Error: initialization already done','','')
172	l_first_condveg=.FALSE.
173
174	!! Initialize local printlev
175	printlev_loc=get_printlev('condveg')
176
177	IF (printlev>=3) WRITE (numout,*) 'Start condveg_initialize'
178
179	!! 1. Allocate module variables and read from restart or initialize
180	IF (alb_bg_modis) THEN
181	! Allocate background soil albedo
182	ALLOCATE (soilalb_bg(kjpindex,2),stat=ier)
183	IF (ier /= 0) CALL ipslerr_p(3,'condveg_initialize','Pb in allocation for soilalb_bg','','')
184
185	! Read background albedo from restart file
186	CALL ioconf_setatt_p('UNITS', '-')
187	CALL ioconf_setatt_p('LONG_NAME','Background soil albedo for visible and near-infrared range')
188	CALL restget_p (rest_id, 'soilalbedo_bg', nbp_glo, 2, 1, kjit, .TRUE., soilalb_bg, "gather", nbp_glo, index_g)
189
190	! Initialize by interpolating from file if the variable was not in restart file
191	IF ( ALL(soilalb_bg(:,:) == val_exp) ) THEN
192	CALL condveg_background_soilalb(kjpindex, lalo, neighbours, resolution, contfrac)
193	END IF
194	CALL xios_orchidee_send_field("soilalb_bg",soilalb_bg)
195
196	ELSE
197	! Allocate
198	! Dry soil albedo
199	ALLOCATE (soilalb_dry(kjpindex,2),stat=ier)
200	IF (ier /= 0) CALL ipslerr_p(3,'condveg_initialize','Pb in allocation for soilalb_dry','','')
201
202	! Wet soil albedo
203	ALLOCATE (soilalb_wet(kjpindex,2),stat=ier)
204	IF (ier /= 0) CALL ipslerr_p(3,'condveg_initialize','Pb in allocation for soilalb_wet','','')
205
206	! Mean soil albedo
207	ALLOCATE (soilalb_moy(kjpindex,2),stat=ier)
208	IF (ier /= 0) CALL ipslerr_p(3,'condveg_initialize','Pb in allocation for soilalb_moy','','')
209
210	! Read variables from restart file
211	! dry soil albedo
212	CALL ioconf_setatt_p('UNITS', '-')
213	CALL ioconf_setatt_p('LONG_NAME','Dry bare soil albedo')
214	CALL restget_p (rest_id,'soilalbedo_dry' , nbp_glo, 2, 1, kjit, .TRUE., soilalb_dry, "gather", nbp_glo, index_g)
215
216	! wet soil albedo
217	CALL ioconf_setatt_p('UNITS', '-')
218	CALL ioconf_setatt_p('LONG_NAME','Wet bare soil albedo')
219	CALL restget_p (rest_id, 'soilalbedo_wet', nbp_glo, 2, 1, kjit, .TRUE., soilalb_wet, "gather", nbp_glo, index_g)
220
221	! mean soil aledo
222	CALL ioconf_setatt_p('UNITS', '-')
223	CALL ioconf_setatt_p('LONG_NAME','Mean bare soil albedo')
224	CALL restget_p (rest_id, 'soilalbedo_moy', nbp_glo, 2, 1, kjit, .TRUE., soilalb_moy, "gather", nbp_glo, index_g)
225
226
227	! Initialize the variables if not found in restart file
228	IF ( ALL(soilalb_wet(:,:) == val_exp) .OR. &
229	ALL(soilalb_dry(:,:) == val_exp) .OR. &
230	ALL(soilalb_moy(:,:) == val_exp)) THEN
231	! One or more of the variables were not in the restart file.
232	! Call routine condveg_soilalb to calculate them.
233	CALL condveg_soilalb(kjpindex, lalo, neighbours, resolution, contfrac)
234	WRITE(numout,*) '---> val_exp ', val_exp
235	WRITE(numout,*) '---> ALBEDO_wet VIS:', MINVAL(soilalb_wet(:,ivis)), MAXVAL(soilalb_wet(:,ivis))
236	WRITE(numout,*) '---> ALBEDO_wet NIR:', MINVAL(soilalb_wet(:,inir)), MAXVAL(soilalb_wet(:,inir))
237	WRITE(numout,*) '---> ALBEDO_dry VIS:', MINVAL(soilalb_dry(:,ivis)), MAXVAL(soilalb_dry(:,ivis))
238	WRITE(numout,*) '---> ALBEDO_dry NIR:', MINVAL(soilalb_dry(:,inir)), MAXVAL(soilalb_dry(:,inir))
239	WRITE(numout,*) '---> ALBEDO_moy VIS:', MINVAL(soilalb_moy(:,ivis)), MAXVAL(soilalb_moy(:,ivis))
240	WRITE(numout,*) '---> ALBEDO_moy NIR:', MINVAL(soilalb_moy(:,inir)), MAXVAL(soilalb_moy(:,inir))
241	ENDIF
242	END IF
243
244	! z0m
245	CALL ioconf_setatt_p('UNITS', '-')
246	CALL ioconf_setatt_p('LONG_NAME','Roughness for momentum')
247	CALL restget_p (rest_id, 'z0m', nbp_glo, 1, 1, kjit, .TRUE., z0m, "gather", nbp_glo, index_g)
248
249	! z0h
250	CALL ioconf_setatt_p('UNITS', '-')
251	CALL ioconf_setatt_p('LONG_NAME','Roughness for heat')
252	CALL restget_p (rest_id, 'z0h', nbp_glo, 1, 1, kjit, .TRUE., z0h, "gather", nbp_glo, index_g)
253
254	! roughness height
255	CALL ioconf_setatt_p('UNITS', '-')
256	CALL ioconf_setatt_p('LONG_NAME','Roughness height')
257	CALL restget_p (rest_id, 'roughheight', nbp_glo, 1, 1, kjit, .TRUE., roughheight, "gather", nbp_glo, index_g)
258	CALL restget_p (rest_id, 'roughheight_pft', nbp_glo, nvm, 1, kjit, .TRUE., roughheight_pft, "gather", nbp_glo, index_g)
259
260
261	!! 2. Calculate emissivity
262	! If TRUE read in default values for emissivity
263	IF ( impaze ) THEN
264	!
265	emis(:) = emis_scal
266	!
267	! If FALSE set emissivity to 1.
268	ELSE
269	emis_scal = un
270	emis(:) = emis_scal
271	ENDIF
272
273
274	!! 3. Calculate the fraction of snow on vegetation and nobio
275	CALL condveg_frac_snow(kjpindex, snow, snow_nobio, snowrho, snowdz, &
276	frac_snow_veg, frac_snow_nobio)
277
278	!! 4. Calculate roughness height if it was not found in the restart file
279	IF ( ALL(z0m(:) == val_exp) .OR. ALL(z0h(:) == val_exp) .OR. ALL(roughheight(:) == val_exp)) THEN
280	!! Calculate roughness height
281	! Chooses between two methods to calculate the grid average of the roughness.
282	! If impaze set to true: The grid average is calculated by averaging the drag coefficients over PFT.
283	! If impaze set to false: The grid average is calculated by averaging the logarithm of the roughness length per PFT.
284	IF ( impaze ) THEN
285	! Use parameter CONDVEG_Z0 and ROUGHHEIGHT from run.def
286	z0m(:) = z0_scal
287	z0h(:) = z0_scal
288	roughheight(:) = roughheight_scal
289	roughheight_pft(:,:) = roughheight_scal
290	ELSE
291	! Caluculate roughness height
292	IF( rough_dyn ) THEN
293	CALL condveg_z0cdrag_dyn(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
294	& height, temp_air, pb, u, v, lai, frac_snow_veg, z0m, z0h, roughheight, roughheight_pft)
295	ELSE
296	CALL condveg_z0cdrag(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
297	height, tot_bare_soil, frac_snow_veg, z0m, z0h, roughheight, roughheight_pft)
298	ENDIF
299	END IF
300	END IF
301
302	!! 5. Calculate albedo
303	CALL condveg_albedo (kjpindex, veget, veget_max, drysoil_frac, frac_nobio, &
304	totfrac_nobio, snow, snow_age, snow_nobio, &
305	snow_nobio_age, snowdz, snowrho, &
306	tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
307	albedo, albedo_snow, alb_bare, alb_veget)
308
309	IF (printlev>=3) WRITE (numout,*) 'condveg_initialize done '
310
311	END SUBROUTINE condveg_initialize
312
313	!! ==============================================================================================================================
314	!! SUBROUTINE : condveg_main
315	!!
316	!>\BRIEF Calls the subroutines to initialise the variables, update the variables
317	!! and write out data and restart files.
318	!!
319	!!
320	!! MAIN OUTPUT VARIABLE(S): emis (emissivity), albedo (albedo of
321	!! vegetative PFTs in visible and near-infrared range), z0 (surface roughness height),
322	!! roughheight (grid effective roughness height), soil type (fraction of soil types)
323	!!
324	!!
325	!! REFERENCE(S) : None
326	!!
327	!! FLOWCHART : None
328	!!
329	!! REVISION(S) : None
330	!!
331	!_ ================================================================================================================================
332
333	SUBROUTINE condveg_main (kjit, kjpindex, index,&
334	& lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
335	& zlev, snow, snow_age, snow_nobio, snow_nobio_age, tot_bare_soil, &
336	temp_air, pb, u, v, lai, &
337	& drysoil_frac, height, snowdz, snowrho, emis, albedo, &
338	& frac_snow_veg, frac_snow_nobio, &
339	& z0m, z0h, roughheight, roughheight_pft, rest_id, hist_id, hist2_id)
340
341	!! 0. Variable and parameter declaration
342
343	!! 0.1 Input variables
344
345	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
346	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
347	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
348	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier
349	INTEGER(i_std), OPTIONAL, INTENT (in) :: hist2_id !! _History_ file 2 identifier
350	! input fields
351	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
352	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
353	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours!! neighoring grid points if land
354	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
355	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
356	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
357	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
358	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
359	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
360	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
361	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
362	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
363	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
364	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
365	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
366	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
367	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowdz !! Snow depth at each snow layer
368	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in):: snowrho !! Snow density at each snow layer
369	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
370	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: frac_snow_veg !! Snow cover fraction on vegeted area
371	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Snow cover fraction on non-vegeted area
372	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp_air !! Air temperature
373	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
374	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: u !! Horizontal wind speed, u direction
375	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: v !! Horizontal wind speed, v direction
376	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
377
378	!! 0.2 Output variables
379
380	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
381	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
382	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0m !! Roughness for momentum (m)
383	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0h !! Roughness for heat (m)
384	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
385	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: roughheight_pft
386
387	!! 0.3 Modified variables
388
389	!! 0.4 Local variables
390	REAL(r_std), DIMENSION(kjpindex,2) :: albedo_snow !! Snow albedo (unitless ratio)
391	REAL(r_std), DIMENSION(kjpindex) :: albedo_snow_mean !! Mean snow albedo over all wave length, for diag (unitless ratio)
392	REAL(r_std), DIMENSION(kjpindex,2) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
393	!! range (unitless)
394	REAL(r_std), DIMENSION(kjpindex,2) :: alb_veget !! Mean vegetation albedo for visible and near-infrared
395	!! range (unitless)
396	INTEGER(i_std) :: ji
397	!_ ================================================================================================================================
398
399	!! 1. Calculate the fraction of snow on vegetation and nobio
400	CALL condveg_frac_snow(kjpindex, snow, snow_nobio, snowrho, snowdz, &
401	frac_snow_veg, frac_snow_nobio)
402
403	!! 2. Calculate emissivity
404	emis(:) = emis_scal
405
406	!! 3. Calculate roughness height
407
408	! If TRUE read in prescribed values for roughness height
409	IF ( impaze ) THEN
410
411	DO ji = 1, kjpindex
412	z0m(ji) = z0_scal
413	z0h(ji) = z0_scal
414	roughheight(ji) = roughheight_scal
415	roughheight_pft(ji,:) = roughheight_scal
416	ENDDO
417
418	! Calculate roughness height
419	ELSE
420
421	IF ( rough_dyn ) THEN
422	CALL condveg_z0cdrag_dyn (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, height, temp_air, pb, u, v, lai, &
423	frac_snow_veg, z0m, z0h, roughheight, roughheight_pft)
424	ELSE
425	CALL condveg_z0cdrag (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, height, tot_bare_soil, &
426	frac_snow_veg, z0m, z0h, roughheight, roughheight_pft)
427	ENDIF
428
429	ENDIF
430
431	!! 4. Calculate albedo
432	CALL condveg_albedo (kjpindex, veget, veget_max, drysoil_frac, frac_nobio, &
433	totfrac_nobio, snow, snow_age, snow_nobio, &
434	snow_nobio_age, snowdz, snowrho, &
435	tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
436	albedo, albedo_snow, alb_bare, alb_veget)
437
438
439
440	!! 5. Output diagnostics
441	CALL xios_orchidee_send_field("soilalb_vis",alb_bare(:,1))
442	CALL xios_orchidee_send_field("soilalb_nir",alb_bare(:,2))
443	CALL xios_orchidee_send_field("vegalb_vis",alb_veget(:,1))
444	CALL xios_orchidee_send_field("vegalb_nir",alb_veget(:,2))
445	CALL xios_orchidee_send_field("albedo_vis",albedo(:,1))
446	CALL xios_orchidee_send_field("albedo_nir",albedo(:,2))
447
448	! Calculcate albedo_snow mean over wave length, setting xios_default_val when there is no snow
449	DO ji=1,kjpindex
450	IF (snow(ji) > 0) THEN
451	albedo_snow_mean(ji) = (albedo_snow(ji,1) + albedo_snow(ji,2))/2
452	ELSE
453	albedo_snow_mean(ji) = xios_default_val
454	END IF
455	END DO
456	CALL xios_orchidee_send_field("albedo_snow", albedo_snow_mean)
457
458	IF ( almaoutput ) THEN
459	CALL histwrite_p(hist_id, 'Albedo', kjit, (albedo(:,1) + albedo(:,2))/2, kjpindex, index)
460	CALL histwrite_p(hist_id, 'SAlbedo', kjit, (albedo_snow(:,1) + albedo_snow(:,2))/2, kjpindex, index)
461	IF ( hist2_id > 0 ) THEN
462	CALL histwrite_p(hist2_id, 'Albedo', kjit, (albedo(:,1) + albedo(:,2))/2, kjpindex, index)
463	CALL histwrite_p(hist2_id, 'SAlbedo', kjit, (albedo_snow(:,1) + albedo_snow(:,2))/2, kjpindex, index)
464	ENDIF
465	ELSE
466	CALL histwrite_p(hist_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
467	CALL histwrite_p(hist_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
468	CALL histwrite_p(hist_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
469	CALL histwrite_p(hist_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
470	IF ( hist2_id > 0 ) THEN
471	CALL histwrite_p(hist2_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
472	CALL histwrite_p(hist2_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
473	CALL histwrite_p(hist2_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
474	CALL histwrite_p(hist2_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
475	ENDIF
476	ENDIF
477
478	IF (printlev>=3) WRITE (numout,*)' condveg_main done '
479
480	END SUBROUTINE condveg_main
481
482	!!
483	!=============================================================================================================================
484	!! SUBROUTINE : condveg_finalize
485	!!
486	!>\BRIEF Write to restart file
487	!!
488	!! DESCRIPTION : This subroutine writes the module
489	!variables and variables calculated in condveg
490	!! to restart file
491	!!
492	!! RECENT CHANGE(S) : None
493	!!
494	!! REFERENCE(S) : None
495	!!
496	!! FLOWCHART : None
497	!! \n
498	!_
499	!==============================================================================================================================
500	SUBROUTINE condveg_finalize (kjit, kjpindex, rest_id, z0m, z0h, roughheight, roughheight_pft)
501
502	!! 0. Variable and parameter declaration
503	!! 0.1 Input variables
504	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
505	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
506	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
507	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: z0m !! Roughness for momentum
508	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: z0h !! Roughness for heat
509	REAL(r_std),DIMENSION(kjpindex), INTENT(in) :: roughheight !! Grid effective roughness height (m)
510	REAL(r_std),DIMENSION(kjpindex,nvm), INTENT(in) :: roughheight_pft
511
512	!_ ================================================================================================================================
513
514	CALL restput_p (rest_id, 'z0m', nbp_glo, 1, 1, kjit, z0m, 'scatter', nbp_glo, index_g)
515	CALL restput_p (rest_id, 'z0h', nbp_glo, 1, 1, kjit, z0h, 'scatter', nbp_glo, index_g)
516	CALL restput_p (rest_id, 'roughheight', nbp_glo, 1, 1, kjit, roughheight, 'scatter', nbp_glo, index_g)
517	CALL restput_p (rest_id, 'roughheight_pft', nbp_glo, nvm, 1, kjit, roughheight_pft, 'scatter', nbp_glo, index_g)
518
519	IF ( alb_bg_modis ) THEN
520	CALL restput_p (rest_id, 'soilalbedo_bg', nbp_glo, 2, 1, kjit, soilalb_bg, 'scatter', nbp_glo, index_g)
521	ELSE
522	CALL restput_p (rest_id, 'soilalbedo_dry', nbp_glo, 2, 1, kjit, soilalb_dry, 'scatter', nbp_glo, index_g)
523	CALL restput_p (rest_id, 'soilalbedo_wet', nbp_glo, 2, 1, kjit, soilalb_wet, 'scatter', nbp_glo, index_g)
524	CALL restput_p (rest_id, 'soilalbedo_moy', nbp_glo, 2, 1, kjit, soilalb_moy, 'scatter', nbp_glo, index_g)
525	END IF
526	END SUBROUTINE condveg_finalize
527
528	!! ==============================================================================================================================
529	!! SUBROUTINE : condveg_clear
530	!!
531	!>\BRIEF Deallocate albedo variables
532	!!
533	!! DESCRIPTION : None
534	!!
535	!! RECENT CHANGE(S): None
536	!!
537	!! MAIN OUTPUT VARIABLE(S): None
538	!!
539	!! REFERENCES : None
540	!!
541	!! FLOWCHART : None
542	!! \n
543	!_ ================================================================================================================================
544
545	SUBROUTINE condveg_clear ()
546
547	l_first_condveg=.TRUE.
548
549	! Dry soil albedo
550	IF (ALLOCATED (soilalb_dry)) DEALLOCATE (soilalb_dry)
551	! Wet soil albedo
552	IF (ALLOCATED(soilalb_wet)) DEALLOCATE (soilalb_wet)
553	! Mean soil albedo
554	IF (ALLOCATED(soilalb_moy)) DEALLOCATE (soilalb_moy)
555	! BG soil albedo
556	IF (ALLOCATED(soilalb_bg)) DEALLOCATE (soilalb_bg)
557
558	END SUBROUTINE condveg_clear
559
560	!! ==============================================================================================================================\n
561	!! SUBROUTINE : condveg_albedo
562	!!
563	!>\BRIEF Calculate albedo
564	!!
565	!! DESCRIPTION : The albedo is calculated for both the visible and near-infrared
566	!! domain. First the mean albedo of the bare soil is calculated. Two options exist:
567	!! either the soil albedo depends on soil wetness (drysoil_frac variable), or the soil albedo
568	!! is set to a mean soil albedo value.
569	!! The snow albedo scheme presented below belongs to prognostic albedo
570	!! category, i.e. the snow albedo value at a time step depends on the snow albedo value
571	!! at the previous time step.
572	!!
573	!! First, the following formula (described in Chalita and Treut 1994) is used to describe
574	!! the change in snow albedo with snow age on each PFT and each non-vegetative surfaces,
575	!! i.e. continental ice, lakes, etc.: \n
576	!! \latexonly
577	!! \input{SnowAlbedo.tex}
578	!! \endlatexonly
579	!! \n
580	!! Where snowAge is snow age, tcstSnowa is a critical aging time (tcstSnowa=5 days)
581	!! snowaIni and snowaIni+snowaDec corresponds to albedos measured for aged and
582	!! fresh snow respectively, and their values for each PFT and each non-vegetative surfaces
583	!! is precribed in in constantes_veg.f90.\n
584	!! In order to estimate gridbox snow albedo, snow albedo values for each PFT and
585	!! each non-vegetative surfaces with a grid box are weightedly summed up by their
586	!! respective fractions.\n
587	!! Secondly, the snow cover fraction is computed as:
588	!! \latexonly
589	!! \input{SnowFraction.tex}
590	!! \endlatexonly
591	!! \n
592	!! Where fracSnow is the fraction of snow on total vegetative or total non-vegetative
593	!! surfaces, snow is snow mass (kg/m^2) on total vegetated or total nobio surfaces.\n
594	!! Finally, the surface albedo is then updated as the weighted sum of fracSnow, total
595	!! vegetated fraction, total nobio fraction, gridbox snow albedo, and previous
596	!! time step surface albedo.
597	!!
598	!! RECENT CHANGE(S): These calculations were previously done in condveg_albcalc and condveg_snow
599	!!
600	!! MAIN OUTPUT VARIABLE(S): :: albedo; surface albedo. :: albedo_snow; snow
601	!! albedo
602	!!
603	!! REFERENCE(S) :
604	!! Chalita, S. and H Le Treut (1994), The albedo of temperate and boreal forest and
605	!! the Northern Hemisphere climate: a sensitivity experiment using the LMD GCM,
606	!! Climate Dynamics, 10 231-240.
607	!!
608	!! FLOWCHART : None
609	!! \n
610	!_ ================================================================================================================================
611
612	SUBROUTINE condveg_albedo (kjpindex, veget, veget_max, drysoil_frac, frac_nobio, &
613	totfrac_nobio, snow, snow_age, snow_nobio, &
614	snow_nobio_age, snowdz, snowrho, &
615	tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
616	albedo, albedo_snow, alb_bare, alb_veget)
617
618	!! 0. Variable and parameter declarations
619
620	!! 0.1 Input variables
621
622	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
623	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
624	!! (m^2 m^{-2})
625	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max
626	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
627	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces, i.e.
628	!! continental ice, lakes, etc. (unitless)
629	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces, i.e.
630	!! continental ice, lakes, etc. (unitless)
631	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
632	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
633	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age (days)
634	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on continental ice, lakes, etc. (days)
635	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer
636	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density at each snow layer
637	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
638	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: frac_snow_veg !! Fraction of snow on vegetation (unitless ratio)
639	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc. (unitless ratio)
640
641	!! 0.2 Output variables
642	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo (unitless ratio)
643	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo_snow !! Snow albedo (unitless ratio)
644	REAL(r_std), DIMENSION(kjpindex,2), INTENT(out) :: alb_bare !! Mean bare soil albedo for visible and near-infrared
645	!! range (unitless). Only calculated for .NOT. impaze
646	REAL(r_std), DIMENSION(kjpindex,2), INTENT(out) :: alb_veget !! Mean vegetation albedo for visible and near-infrared
647	!! range (unitless). Only calculated for .NOT. impaze
648
649	!! 0.3 Local variables
650	INTEGER(i_std) :: ji, jv, jb,ks !! indices (unitless)
651	REAL(r_std), DIMENSION(kjpindex,2) :: snowa_veg !! Albedo of snow covered area on vegetation
652	!! (unitless ratio)
653	REAL(r_std), DIMENSION(kjpindex,nnobio,2) :: snowa_nobio !! Albedo of snow covered area on continental ice,
654	!! lakes, etc. (unitless ratio)
655	REAL(r_std), DIMENSION(kjpindex) :: fraction_veg !! Total vegetation fraction (unitless ratio)
656	REAL(r_std), DIMENSION(kjpindex) :: agefunc_veg !! Age dependency of snow albedo on vegetation
657	!! (unitless)
658	REAL(r_std), DIMENSION(kjpindex,nnobio) :: agefunc_nobio !! Age dependency of snow albedo on ice,
659	!! lakes, .. (unitless)
660	REAL(r_std) :: alb_nobio !! Albedo of continental ice, lakes, etc.
661	!!(unitless ratio)
662	REAL(r_std),DIMENSION (nvm,2) :: alb_leaf_tmp !! Variables for albedo values for all PFTs and
663	REAL(r_std),DIMENSION (nvm,2) :: snowa_aged_tmp !! spectral domains (unitless)
664	REAL(r_std),DIMENSION (nvm,2) :: snowa_dec_tmp
665	!_ ================================================================================================================================
666
667
668	snowa_aged_tmp(:,ivis) = snowa_aged_vis(:)
669	snowa_aged_tmp(:,inir) = snowa_aged_nir(:)
670	snowa_dec_tmp(:,ivis) = snowa_dec_vis(:)
671	snowa_dec_tmp(:,inir) = snowa_dec_nir(:)
672
673	!! 1. Preliminary calculation without considering snow
674	snowa_aged_tmp(:,ivis) = snowa_aged_vis(:)
675	snowa_aged_tmp(:,inir) = snowa_aged_nir(:)
676	snowa_dec_tmp(:,ivis) = snowa_dec_vis(:)
677	snowa_dec_tmp(:,inir) = snowa_dec_nir(:)
678
679	IF ( impaze ) THEN
680	!! No caluculation, set default value
681	albedo(:,ivis) = albedo_scal(ivis)
682	albedo(:,inir) = albedo_scal(inir)
683	!! need these variables for snow albedo and for diagnostic output
684	alb_veget(:,ivis) = albedo_scal(ivis)
685	alb_veget(:,inir) = albedo_scal(inir)
686	alb_bare(:,ivis) = albedo_scal(ivis)
687	alb_bare(:,inir) = albedo_scal(inir)
688
689	! These variables are needed for snow albedo and for diagnostic output
690	alb_veget(:,ivis) = albedo_scal(ivis)
691	alb_veget(:,inir) = albedo_scal(inir)
692	alb_bare(:,ivis) = albedo_scal(ivis)
693	alb_bare(:,inir) = albedo_scal(inir)
694	ELSE
695	!! Preliminary calculation without considering snow (previously done in condveg_albcalc)
696	! Assign values of leaf and snow albedo for visible and near-infrared range
697	! to local variable (constantes_veg.f90)
698	alb_leaf_tmp(:,ivis) = alb_leaf_vis(:)
699	alb_leaf_tmp(:,inir) = alb_leaf_nir(:)
700
701	!! 1.1 Calculation and assignment of soil albedo
702
703	DO ks = 1, 2! Loop over # of spectra
704
705	! If alb_bg_modis=TRUE, the background soil albedo map for the current simulated month is used
706	! If alb_bg_modis=FALSE and alb_bare_model=TRUE, the soil albedo calculation depends on soil moisture
707	! If alb_bg_modis=FALSE and alb_bare_model=FALSE, the mean soil albedo is used without the dependance on soil moisture
708	! see subroutines 'condveg_soilalb' and 'condveg_background_soilalb'
709	IF ( alb_bg_modis ) THEN
710	alb_bare(:,ks) = soilalb_bg(:,ks)
711	ELSE
712	IF ( alb_bare_model ) THEN
713	alb_bare(:,ks) = soilalb_wet(:,ks) + drysoil_frac(:) * (soilalb_dry(:,ks) - soilalb_wet(:,ks))
714	ELSE
715	alb_bare(:,ks) = soilalb_moy(:,ks)
716	ENDIF
717	ENDIF
718
719	! Soil albedo is weighed by fraction of bare soil
720	albedo(:,ks) = tot_bare_soil(:) * alb_bare(:,ks)
721
722	!! 1.2 Calculation of mean albedo of over the grid cell
723
724	! Calculation of mean albedo of over the grid cell and
725	! mean albedo of only vegetative PFTs over the grid cell
726	alb_veget(:,ks) = zero
727
728	DO jv = 2, nvm ! Loop over # of PFTs
729
730	! Mean albedo of grid cell for visible and near-infrared range
731	albedo(:,ks) = albedo(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
732
733	! Mean albedo of vegetation for visible and near-infrared range
734	alb_veget(:,ks) = alb_veget(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
735	ENDDO ! Loop over # of PFTs
736
737	ENDDO
738	END IF
739
740
741	!! 2. Calculate snow albedos on both total vegetated and total nobio surfaces
742
743	! The snow albedo could be either prescribed (in condveg_init.f90) or
744	! calculated following Chalita and Treut (1994).
745	! Check if the precribed value fixed_snow_albedo exists
746	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
747	snowa_veg(:,:) = fixed_snow_albedo
748	snowa_nobio(:,:,:) = fixed_snow_albedo
749	fraction_veg(:) = un - totfrac_nobio(:)
750	ELSE ! calculated following Chalita and Treut (1994)
751
752	!! 2.1 Calculate age dependence
753
754	! On vegetated surfaces
755	DO ji = 1, kjpindex
756	agefunc_veg(ji) = EXP(-snow_age(ji)/tcst_snowa)
757	ENDDO
758
759	! On non-vegtative surfaces
760	DO jv = 1, nnobio ! Loop over # nobio types
761	DO ji = 1, kjpindex
762	agefunc_nobio(ji,jv) = EXP(-snow_nobio_age(ji,jv)/tcst_snowa)
763	ENDDO
764	ENDDO
765
766	!! 2.1 Calculate snow albedo
767	! For vegetated surfaces
768	fraction_veg(:) = un - totfrac_nobio(:)
769	snowa_veg(:,:) = zero
770	!
771	! Alternative formulation based on veget and not veget_max that needs to be tested
772	! See ticket 223
773	!
774	! IF (ok_dgvm) THEN
775	! DO jb = 1, 2
776	! DO ji = 1, kjpindex
777	! IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
778	! snowa_veg(ji,jb) = snowa_veg(ji,jb) + &
779	! tot_bare_soil(ji)/fraction_veg(ji) * ( snowa_aged_tmp(1,jb)+snowa_dec_tmp(1,jb)*agefunc_veg(ji) )
780	! END IF
781	! END DO
782	! END DO
783	!
784	! DO jb = 1, 2
785	! DO jv = 2, nvm
786	! DO ji = 1, kjpindex
787	! IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
788	! snowa_veg(ji,jb) = snowa_veg(ji,jb) + &
789	! veget(ji,jv)/fraction_veg(ji) * ( snowa_aged_tmp(jv,jb)+snowa_dec_tmp(jv,jb)*agefunc_veg(ji) )
790	! ENDIF
791	! ENDDO
792	! ENDDO
793	! ENDDO
794	! ELSE
795	DO jb = 1, 2
796	DO jv = 1, nvm
797	DO ji = 1, kjpindex
798	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
799	snowa_veg(ji,jb) = snowa_veg(ji,jb) + &
800	veget_max(ji,jv)/fraction_veg(ji) * ( snowa_aged_tmp(jv,jb)+snowa_dec_tmp(jv,jb)*agefunc_veg(ji) )
801	ENDIF
802	ENDDO
803	ENDDO
804	ENDDO
805	! ENDIF
806	!
807	! snow albedo on other surfaces
808	!
809	DO jb = 1, 2
810	DO jv = 1, nnobio
811	DO ji = 1, kjpindex
812	snowa_nobio(ji,jv,jb) = ( snowa_aged_tmp(1,jb) + snowa_dec_tmp(1,jb) * agefunc_nobio(ji,jv) )
813	ENDDO
814	ENDDO
815	ENDDO
816	ENDIF
817
818	!! 3. Update surface albedo
819
820	! Update surface albedo using the weighted sum of previous time step surface albedo,
821	! total vegetated fraction, total nobio fraction, snow cover fraction (both vegetated and
822	! non-vegetative surfaces), and snow albedo (both vegetated and non-vegetative surfaces).
823	! Although both visible and near-infrared surface albedo are presented, their calculations
824	! are the same.
825	DO jb = 1, 2
826
827	albedo(:,jb) = ( fraction_veg(:) ) * &
828	( (un-frac_snow_veg(:)) * albedo(:,jb) + &
829	( frac_snow_veg(:) ) * snowa_veg(:,jb) )
830	DO jv = 1, nnobio ! Loop over # nobio surfaces
831
832	IF ( jv .EQ. iice ) THEN
833	alb_nobio = alb_ice(jb)
834	ELSE
835	WRITE(numout,*) 'jv=',jv
836	WRITE(numout,*) 'DO NOT KNOW ALBEDO OF THIS SURFACE TYPE'
837	CALL ipslerr_p(3,'condveg_snow','DO NOT KNOW ALBEDO OF THIS SURFACE TYPE','','')
838	ENDIF
839
840	albedo(:,jb) = albedo(:,jb) + &
841	( frac_nobio(:,jv) ) * &
842	( (un-frac_snow_nobio(:,jv)) * alb_nobio + &
843	( frac_snow_nobio(:,jv) ) * snowa_nobio(:,jv,jb) )
844	ENDDO
845
846	END DO
847
848	! Calculate snow albedo
849	DO jb = 1, 2
850	albedo_snow(:,jb) = fraction_veg(:) * frac_snow_veg(:) * snowa_veg(:,jb)
851	DO jv = 1, nnobio
852	albedo_snow(:,jb) = albedo_snow(:,jb) + &
853	frac_nobio(:,jv) * frac_snow_nobio(:,jv) * snowa_nobio(:,jv,jb)
854	ENDDO
855	ENDDO
856
857	IF (printlev>=3) WRITE (numout,*) ' condveg_albedo done '
858
859	END SUBROUTINE condveg_albedo
860
861
862
863	!! ==============================================================================================================================
864	!! SUBROUTINE : condveg_frac_snow
865	!!
866	!>\BRIEF This subroutine calculates the fraction of snow on vegetation and nobio
867	!!
868	!! DESCRIPTION
869	!!
870	!! RECENT CHANGE(S): These calculations were previously done in condveg_snow.
871	!!
872	!! REFERENCE(S) :
873	!!
874	!! FLOWCHART : None
875	!! \n
876	!_ ================================================================================================================================
877
878	SUBROUTINE condveg_frac_snow(kjpindex, snow, snow_nobio, snowrho, snowdz, &
879	frac_snow_veg, frac_snow_nobio)
880	!! 0. Variable and parameter declaration
881	!! 0.1 Input variables
882	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
883	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass in vegetation (kg m^{-2})
884	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass on continental ice, lakes, etc. (kg m^{-2})
885	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowrho !! Snow density at each snow layer
886	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(in) :: snowdz !! Snow depth at each snow layer
887
888	!! 0.2 Output variables
889	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: frac_snow_veg !! Fraction of snow on vegetation (unitless ratio)
890	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(out):: frac_snow_nobio !! Fraction of snow on continental ice, lakes, etc.
891
892	!! 0.3 Local variables
893	REAL(r_std), DIMENSION(kjpindex) :: snowrho_ave !! Average snow density
894	REAL(r_std), DIMENSION(kjpindex) :: snowdepth !! Snow depth
895	REAL(r_std), DIMENSION(kjpindex) :: snowrho_snowdz !! Snow rho time snowdz
896	INTEGER(i_std) :: jv
897
898	!! Calculate snow cover fraction for both total vegetated and total non-vegetative surfaces.
899	IF (ok_explicitsnow) THEN
900	snowdepth=sum(snowdz,2)
901	snowrho_snowdz=sum(snowrho*snowdz,2)
902	WHERE(snowdepth(:) .LT. min_sechiba)
903	frac_snow_veg(:) = 0.
904	ELSEWHERE
905	snowrho_ave(:)=snowrho_snowdz(:)/snowdepth(:)
906	frac_snow_veg(:) = tanh(snowdepth(:)/(0.025*(snowrho_ave(:)/50.)))
907	END WHERE
908	ELSE
909	frac_snow_veg(:) = MIN(MAX(snow(:),zero)/(MAX(snow(:),zero)+snowcri_alb*sn_dens/100.0),un)
910	END IF
911
912	DO jv = 1, nnobio
913	frac_snow_nobio(:,jv) = MIN(MAX(snow_nobio(:,jv),zero)/(MAX(snow_nobio(:,jv),zero)+snowcri_alb),un)
914	ENDDO
915
916	IF (printlev>=3) WRITE (numout,*) ' condveg_frac_snow done '
917
918	END SUBROUTINE condveg_frac_snow
919
920
921	!! ==============================================================================================================================
922	!! SUBROUTINE : condveg_soilalb
923	!!
924	!>\BRIEF This subroutine calculates the albedo of soil (without snow).
925	!!
926	!! DESCRIPTION This subroutine reads the soil colour maps in 1 x 1 deg resolution
927	!! from the Henderson-Sellers & Wilson database. These values are interpolated to
928	!! the model's resolution and transformed into
929	!! dry and wet albedos.\n
930	!!
931	!! If the soil albedo is calculated without the dependence of soil moisture, the
932	!! soil colour values are transformed into mean soil albedo values.\n
933	!!
934	!! The calculations follow the assumption that the grid of the data is regular and
935	!! it covers the globe. The calculation for the model grid are based on the borders
936	!! of the grid of the resolution.
937	!!
938	!! RECENT CHANGE(S): None
939	!!
940	!! MAIN OUTPUT VARIABLE(S): soilalb_dry for visible and near-infrared range,
941	!! soilalb_wet for visible and near-infrared range,
942	!! soilalb_moy for visible and near-infrared range
943	!!
944	!! REFERENCE(S) :
945	!! -Wilson, M.F., and A. Henderson-Sellers, 1985: A global archive of land cover and
946	!! soils data for use in general circulation climate models. J. Clim., 5, 119-143.
947	!!
948	!! FLOWCHART : None
949	!! \n
950	!_ ================================================================================================================================
951
952	SUBROUTINE condveg_soilalb(nbpt, lalo, neighbours, resolution, contfrac)
953
954	USE interpweight
955
956	IMPLICIT NONE
957
958
959	!! 0. Variable and parameter declaration
960
961	!! 0.1 Input variables
962
963	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
964	!! interpolated (unitless)
965	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
966	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
967	!! (1=N, 2=E, 3=S, 4=W)
968	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
969	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
970
971	!! 0.4 Local variables
972
973	CHARACTER(LEN=80) :: filename !! Filename of soil colour map
974	INTEGER(i_std) :: i, ib, ip, nbexp !! Indices
975	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
976	REAL(r_std), DIMENSION(nbpt) :: asoilcol !! Availability of the soilcol interpolation
977	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
978	!! (variabletypevals(1) = -un, not used)
979	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: soilcolrefrac !! soilcol fractions re-dimensioned
980	REAL(r_std) :: vmin, vmax !! min/max values to use for the
981	!! renormalization
982	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
983	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
984	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
985	!! 'XYKindTime': Input values are kinds
986	!! of something with a temporal
987	!! evolution on the dx*dy matrix'
988	LOGICAL :: nonegative !! whether negative values should be removed
989	CHARACTER(LEN=50) :: maskingtype !! Type of masking
990	!! 'nomask': no-mask is applied
991	!! 'mbelow': take values below maskvals(1)
992	!! 'mabove': take values above maskvals(1)
993	!! 'msumrange': take values within 2 ranges;
994	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
995	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
996	!! (normalized by maskvals(3))
997	!! 'var': mask values are taken from a
998	!! variable inside the file (>0)
999	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
1000	!! `maskingtype')
1001	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
1002	CHARACTER(LEN=250) :: msg
1003	INTEGER :: fopt
1004	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: vecpos
1005	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: solt
1006
1007	!_ ================================================================================================================================
1008	!! 1. Open file and allocate memory
1009
1010	! Open file with soil colours
1011
1012	!Config Key = SOILALB_FILE
1013	!Config Desc = Name of file from which the bare soil albedo
1014	!Config Def = soils_param.nc
1015	!Config If = NOT(IMPOSE_AZE)
1016	!Config Help = The name of the file to be opened to read the soil types from
1017	!Config which we derive then the bare soil albedos. This file is 1x1
1018	!Config deg and based on the soil colors defined by Wilson and Henderson-Seller.
1019	!Config Units = [FILE]
1020	!
1021	filename = 'soils_param.nc'
1022	CALL getin_p('SOILALB_FILE',filename)
1023
1024
1025	ALLOCATE(soilcolrefrac(nbpt, classnb), STAT=ALLOC_ERR)
1026	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable soilcolrefrac','','')
1027	ALLOCATE(vecpos(classnb), STAT=ALLOC_ERR)
1028	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable vecpos','','')
1029	ALLOCATE(solt(classnb), STAT=ALLOC_ERR)
1030	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable solt','','')
1031
1032	! Assigning values to vmin, vmax
1033	vmin = 1.0
1034	vmax = classnb
1035
1036	ALLOCATE(variabletypevals(classnb),STAT=ALLOC_ERR)
1037	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variabletypevals','','')
1038	variabletypevals = -un
1039
1040	!! Variables for interpweight
1041	! Type of calculation of cell fractions
1042	fractype = 'default'
1043	! Name of the longitude and latitude in the input file
1044	lonname = 'nav_lon'
1045	latname = 'nav_lat'
1046	! Should negative values be set to zero from input file?
1047	nonegative = .FALSE.
1048	! Type of mask to apply to the input data (see header for more details)
1049	maskingtype = 'mabove'
1050	! Values to use for the masking
1051	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
1052	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
1053	namemaskvar = ''
1054
1055	! Interpolate variable soilcolor
1056	variablename = 'soilcolor'
1057	IF (printlev_loc >= 1) WRITE(numout,*) "condveg_soilalb: Read and interpolate " &
1058	// TRIM(filename) // " for variable " // TRIM(variablename)
1059	CALL interpweight_2D(nbpt, classnb, variabletypevals, lalo, resolution, neighbours, &
1060	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1061	maskvals, namemaskvar, 0, 0, -1, fractype, &
1062	-1., -1., soilcolrefrac, asoilcol)
1063	IF (printlev_loc >= 5) WRITE(numout,*)' condveg_soilalb after interpweight_2D'
1064
1065	! Check how many points with soil information are found
1066	nbexp = 0
1067
1068	soilalb_dry(:,:) = zero
1069	soilalb_wet(:,:) = zero
1070	soilalb_moy(:,:) = zero
1071	IF (printlev_loc >= 5) THEN
1072	WRITE(numout,*)' condveg_soilalb before starting loop nbpt:', nbpt
1073	WRITE(numout,*)' condveg_soilalb initial values classnb: ',classnb
1074	WRITE(numout,*)' condveg_soilalb vis_dry. SUM:',SUM(vis_dry),' vis_dry= ',vis_dry
1075	WRITE(numout,*)' condveg_soilalb nir_dry. SUM:',SUM(nir_dry),' nir_dry= ',nir_dry
1076	WRITE(numout,*)' condveg_soilalb vis_wet. SUM:',SUM(vis_wet),' vis_wet= ',vis_wet
1077	WRITE(numout,*)' condveg_soilalb nir_wet. SUM:',SUM(nir_wet),' nir_wet= ',nir_wet
1078	END IF
1079
1080	DO ib=1,nbpt ! Loop over domain size
1081
1082	! vecpos: List of positions where textures were not zero
1083	! vecpos(1): number of not null textures found
1084	vecpos = interpweight_ValVecR(soilcolrefrac(ib,:),classnb,zero,'neq')
1085	fopt = vecpos(1)
1086	IF (fopt == classnb) THEN
1087	! All textures are not zero
1088	solt(:) = (/(i,i=1,classnb)/)
1089	ELSE IF (fopt == 0) THEN
1090	WRITE(numout,*)' condveg_soilalb: for point=', ib, ' no soil class!'
1091	ELSE
1092	DO ip = 1,fopt
1093	solt(ip) = vecpos(ip+1)
1094	END DO
1095	END IF
1096
1097	!! 3. Compute the average bare soil albedo parameters
1098
1099	IF ( (fopt .EQ. 0) .OR. (asoilcol(ib) .LT. min_sechiba)) THEN
1100	! Initialize with mean value if no points were interpolated or if no data was found
1101	nbexp = nbexp + 1
1102	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1103	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1104	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
1105	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
1106	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
1107	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
1108	ELSE
1109	! If points were interpolated
1110	DO ip=1, fopt
1111	IF ( solt(ip) .LE. classnb) THEN
1112	! Set to zero if the value is below min_sechiba
1113	IF (soilcolrefrac(ib,solt(ip)) < min_sechiba) soilcolrefrac(ib,solt(ip)) = zero
1114
1115	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(solt(ip))*soilcolrefrac(ib,solt(ip))
1116	soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(solt(ip))*soilcolrefrac(ib,solt(ip))
1117	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(solt(ip))*soilcolrefrac(ib,solt(ip))
1118	soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(solt(ip))*soilcolrefrac(ib,solt(ip))
1119	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(solt(ip))* &
1120	soilcolrefrac(ib,solt(ip))
1121	soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(solt(ip))* &
1122	soilcolrefrac(ib,solt(ip))
1123	ELSE
1124	msg = 'The file contains a soil color class which is incompatible with this program'
1125	CALL ipslerr_p(3,'condveg_soilalb',TRIM(msg),'','')
1126	ENDIF
1127	ENDDO
1128	ENDIF
1129
1130	ENDDO
1131
1132	IF ( nbexp .GT. 0 ) THEN
1133	WRITE(numout,*) 'condveg_soilalb _______'
1134	WRITE(numout,*) 'condveg_soilalb: The interpolation of the bare soil albedo had ', nbexp
1135	WRITE(numout,*) 'condveg_soilalb: points without data. This are either coastal points or'
1136	WRITE(numout,*) 'condveg_soilalb: ice covered land.'
1137	WRITE(numout,*) 'condveg_soilalb: The problem was solved by using the average of all soils'
1138	WRITE(numout,*) 'condveg_soilalb: in dry and wet conditions'
1139	WRITE(numout,*) 'condveg_soilalb: Use the diagnostic output field asoilcol to see location of these points'
1140	ENDIF
1141
1142	DEALLOCATE (soilcolrefrac)
1143	DEALLOCATE (variabletypevals)
1144
1145	! Write diagnostics
1146	CALL xios_orchidee_send_field("asoilcol",asoilcol)
1147
1148
1149	IF (printlev_loc >= 3) WRITE(numout,*)' condveg_soilalb ended'
1150
1151	END SUBROUTINE condveg_soilalb
1152
1153
1154	!! ==============================================================================================================================
1155	!! SUBROUTINE : condveg_background_soilalb
1156	!!
1157	!>\BRIEF This subroutine reads the albedo of bare soil
1158	!!
1159	!! DESCRIPTION This subroutine reads the background albedo map in 0.5 x 0.5 deg resolution
1160	!! derived from JRCTIP product to be used as bare soil albedo. These values are then interpolated
1161	!! to the model's resolution.\n
1162	!!
1163	!! RECENT CHANGE(S): None
1164	!!
1165	!! MAIN OUTPUT VARIABLE(S): soilalb_bg for visible and near-infrared range
1166	!!
1167	!! REFERENCES : None
1168	!!
1169	!! FLOWCHART : None
1170	!! \n
1171	!_ ================================================================================================================================
1172
1173	SUBROUTINE condveg_background_soilalb(nbpt, lalo, neighbours, resolution, contfrac)
1174
1175	USE interpweight
1176
1177	IMPLICIT NONE
1178
1179	!! 0. Variable and parameter declaration
1180
1181	!! 0.1 Input variables
1182
1183	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be
1184	!! interpolated (unitless)
1185	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (degree)
1186	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
1187	!! (1=N, 2=E, 3=S, 4=W)
1188	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid cell in X and Y (km)
1189	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid cell (unitless)
1190
1191	!! 0.4 Local variables
1192
1193	CHARACTER(LEN=80) :: filename !! Filename of background albedo
1194	REAL(r_std), DIMENSION(nbpt) :: aalb_bg !! Availability of the interpolation
1195	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! Latitudes and longitudes read from input file
1196	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel !! Help variable to read file data and allocate memory
1197	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu !! Help variable to read file data and allocate memory
1198	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
1199	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: soilalbedo_bg !! Help variable to read file data and allocate memory
1200	INTEGER :: ALLOC_ERR !! Help varialbe to count allocation error
1201	REAL(r_std) :: vmin, vmax !! min/max values to use for the
1202	!! renormalization
1203	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
1204	CHARACTER(LEN=250) :: maskvname !! Variable to read the mask from
1205	!! the file
1206	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
1207	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
1208	!! 'XYKindTime': Input values are kinds
1209	!! of something with a temporal
1210	!! evolution on the dx*dy matrix'
1211	LOGICAL :: nonegative !! whether negative values should be removed
1212	CHARACTER(LEN=50) :: maskingtype !! Type of masking
1213	!! 'nomask': no-mask is applied
1214	!! 'mbelow': take values below maskvals(1)
1215	!! 'mabove': take values above maskvals(1)
1216	!! 'msumrange': take values within 2 ranges;
1217	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
1218	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
1219	!! (normalized by maskedvals(3))
1220	!! 'var': mask values are taken from a
1221	!! variable inside the file (>0)
1222	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
1223	!! `maskingtype')
1224	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
1225	REAL(r_std) :: albbg_norefinf !! No value
1226	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: albbg_default !! Default value
1227
1228	!_ ================================================================================================================================
1229
1230	!! 1. Open file and allocate memory
1231
1232	! Open file with background albedo
1233
1234	!Config Key = ALB_BG_FILE
1235	!Config Desc = Name of file from which the background albedo is read
1236	!Config Def = alb_bg.nc
1237	!Config If =
1238	!Config Help = The name of the file to be opened to read background albedo
1239	!Config Units = [FILE]
1240	!
1241	filename = 'alb_bg.nc'
1242	CALL getin_p('ALB_BG_FILE',filename)
1243
1244
1245	ALLOCATE(albbg_default(2), STAT=ALLOC_ERR)
1246	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'condveg_background_soilalb','Pb in allocation for albbg_default','','')
1247
1248	! For this case there are not types/categories. We have 'only' a continuos field
1249	! Assigning values to vmin, vmax
1250
1251	vmin = 0.
1252	vmax = 9999.
1253
1254	!! Variables for interpweight
1255	! Type of calculation of cell fractions (not used here)
1256	fractype = 'default'
1257	! Name of the longitude and latitude in the input file
1258	lonname = 'longitude'
1259	latname = 'latitude'
1260	! Default value when no value is get from input file
1261	albbg_default(ivis) = 0.129
1262	albbg_default(inir) = 0.247
1263	! Reference value when no value is get from input file (not used here)
1264	albbg_norefinf = undef_sechiba
1265	! Should negative values be set to zero from input file?
1266	nonegative = .FALSE.
1267	! Type of mask to apply to the input data (see header for more details)
1268	maskingtype = 'var'
1269	! Values to use for the masking (here not used)
1270	maskvals = (/ undef_sechiba, undef_sechiba, undef_sechiba /)
1271	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var')
1272	namemaskvar = 'mask'
1273
1274	! There is a variable for each chanel 'infrared' and 'visible'
1275	! Interpolate variable bg_alb_vis
1276	variablename = 'bg_alb_vis'
1277	IF (printlev_loc >= 1) WRITE(numout,*) "condveg_background_soilalb: Read and interpolate " &
1278	// TRIM(filename) // " for variable " // TRIM(variablename)
1279	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
1280	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1281	maskvals, namemaskvar, -1, fractype, albbg_default(ivis), albbg_norefinf, &
1282	soilalb_bg(:,ivis), aalb_bg)
1283	IF (printlev_loc >= 5) WRITE(numout,*)" condveg_background_soilalb after InterpWeight2Dcont for '" // &
1284	TRIM(variablename) // "'"
1285
1286	! Interpolate variable bg_alb_nir in the same file
1287	variablename = 'bg_alb_nir'
1288	IF (printlev_loc >= 1) WRITE(numout,*) "condveg_background_soilalb: Read and interpolate " &
1289	// TRIM(filename) // " for variable " // TRIM(variablename)
1290	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
1291	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
1292	maskvals, namemaskvar, -1, fractype, albbg_default(inir), albbg_norefinf, &
1293	soilalb_bg(:,inir), aalb_bg)
1294	IF (printlev_loc >= 5) WRITE(numout,*)" condveg_background_soilalb after InterpWeight2Dcont for '" // &
1295	TRIM(variablename) // "'"
1296
1297	IF (ALLOCATED(albbg_default)) DEALLOCATE(albbg_default)
1298
1299	IF (printlev_loc >= 3) WRITE(numout,*)' condveg_background_soilalb ended'
1300
1301	! Write diagnostics
1302	CALL xios_orchidee_send_field("aalb_bg",aalb_bg)
1303
1304	END SUBROUTINE condveg_background_soilalb
1305
1306
1307	!! ==============================================================================================================================
1308	!! SUBROUTINE : condveg_z0cdrag
1309	!!
1310	!>\BRIEF Computation of grid average of roughness length by calculating
1311	!! the drag coefficient.
1312	!!
1313	!! DESCRIPTION : This routine calculates the mean roughness height and mean
1314	!! effective roughness height over the grid cell. The mean roughness height (z0)
1315	!! is computed by averaging the drag coefficients \n
1316	!!
1317	!! \latexonly
1318	!! \input{z0cdrag1.tex}
1319	!! \endlatexonly
1320	!! \n
1321	!!
1322	!! where C is the drag coefficient at the height of the vegetation, kappa is the
1323	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
1324	!! The reference level for z needs to be high enough above the canopy to avoid
1325	!! singularities of the LOG. This height is set to minimum 10m above ground.
1326	!! The drag coefficient increases with roughness height to represent the greater
1327	!! turbulence generated by rougher surfaces.
1328	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
1329	!!
1330	!! The roughness height for the non-vegetative surfaces is calculated in a second step.
1331	!! In order to calculate the transfer coefficients the
1332	!! effective roughness height is calculated. This effective value is the difference
1333	!! between the height of the vegetation and the zero plane displacement height.\nn
1334	!!
1335	!! RECENT CHANGE(S): None
1336	!!
1337	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
1338	!!
1339	!! REFERENCE(S) : None
1340	!!
1341	!! FLOWCHART : None
1342	!! \n
1343	!_ ================================================================================================================================
1344
1345	SUBROUTINE condveg_z0cdrag (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, tot_bare_soil, frac_snow_veg, &
1346	& z0m, z0h, roughheight, roughheight_pft)
1347
1348	!! 0. Variable and parameter declaration
1349
1350	!! 0.1 Input variables
1351
1352	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1353	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1354	!! (m^2 m^{-2})
1355	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1356	!! (= ind*cn_ind) (m^2 m^{-2})
1357	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1358	!! i.e. continental ice, lakes, etc. (unitless)
1359	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1360	!! i.e. continental ice, lakes, etc. (unitless)
1361	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
1362	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1363	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
1364	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1365
1366	!! 0.2 Output variables
1367
1368	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0m !! Roughness height for momentum (m)
1369	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0h !! Roughness height for heat (m)
1370	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1371	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: roughheight_pft
1372
1373	!! 0.3 Modified variables
1374
1375	!! 0.4 Local variables
1376
1377	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1378	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1379	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
1380	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1381	REAL(r_std), DIMENSION(kjpindex) :: d_veg !! PFT coverage of vegetative PFTs
1382	!! (= ind*cn_ind) (m^2 m^{-2})
1383	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1384	REAL(r_std), DIMENSION(kjpindex,nvm) :: zhdispl_pft
1385	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
1386	!! i.e. continental ice, lakes, etc.
1387	REAL(r_std), DIMENSION(kjpindex) :: dragm !! Drag coefficient for momentum
1388	REAL(r_std), DIMENSION(kjpindex) :: dragh !! Drag coefficient for heat
1389	REAL(r_std), DIMENSION(kjpindex) :: z0_ground !! z0m value used for ground surface
1390	!_ ================================================================================================================================
1391
1392	!! 1. Preliminary calculation
1393
1394	! Set maximal height of first layer
1395	ztmp(:) = MAX(10., zlev(:))
1396
1397	z0_ground(:) = (1.-frac_snow_veg(:))z0_bare + frac_snow_veg(:)z0_bare/10.
1398
1399	! Calculate roughness for non-vegetative surfaces
1400	! with the von Karman constant
1401	dragm(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/z0_ground))**2
1402	dragh(:) = tot_bare_soil(:) * (ct_karman/LOG(ztmp(:)/(z0_ground/ratio_z0m_z0h(1))))*(ct_karman/LOG(ztmp(:)/z0_ground))
1403	! Fraction of bare soil
1404	sumveg(:) = tot_bare_soil(:)
1405
1406	! Set average vegetation height to zero
1407	ave_height(:) = zero
1408
1409	!! 2. Calculate the mean roughness height
1410
1411	! Calculate the mean roughness height of
1412	! vegetative PFTs over the grid cell
1413	DO jv = 2, nvm
1414
1415	! In the case of forest, use parameter veget_max because
1416	! tree trunks influence the roughness even when there are no leaves
1417	IF ( is_tree(jv) ) THEN
1418	! In the case of grass, use parameter veget because grasses
1419	! only influence the roughness during the growing season
1420	d_veg(:) = veget_max(:,jv)
1421	ELSE
1422	! grasses only have an influence if they are really there!
1423	d_veg(:) = veget(:,jv)
1424	ENDIF
1425
1426	! Calculate the average roughness over the grid cell:
1427	! The unitless drag coefficient is per vegetative PFT
1428	! calculated by use of the von Karman constant, the height
1429	! of the first layer and the roughness. The roughness
1430	! is calculated as the vegetation height per PFT
1431	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
1432	! If this scaled value is lower than 0.01 then the value for
1433	! the roughness of bare soil (0.01) is used.
1434	! The sum over all PFTs gives the average roughness
1435	! per grid cell for the vegetative PFTs.
1436	dragm(:) = dragm(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)z0_over_height(jv),z0_ground)))*2
1437	dragh(:) = dragh(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/(MAX(height(:,jv)*z0_over_height(jv),z0_ground) / &
1438	ratio_z0m_z0h(jv)))) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)*z0_over_height(jv),z0_ground)))
1439
1440	! Sum of bare soil and fraction vegetated fraction
1441	sumveg(:) = sumveg(:) + d_veg(:)
1442
1443	! Weigh height of vegetation with maximal cover fraction
1444	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1445
1446	ENDDO
1447
1448	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1449
1450	! Search for pixels with vegetated part to normalise
1451	! roughness height
1452	WHERE ( sumveg(:) .GT. min_sechiba )
1453	dragm(:) = dragm(:) / sumveg(:)
1454	dragh(:) = dragh(:) / sumveg(:)
1455	ENDWHERE
1456	! Calculate fraction of roughness for vegetated part
1457	dragm(:) = (un - totfrac_nobio(:)) * dragm(:)
1458	dragh(:) = (un - totfrac_nobio(:)) * dragh(:)
1459
1460	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1461
1462	! Set rougness for ice
1463	IF ( jv .EQ. iice ) THEN
1464	z0_nobio = z0_ice
1465	ELSE
1466	WRITE(numout,*) 'jv=',jv
1467	WRITE(numout,*) 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1468	CALL ipslerr_p(3,'condveg_z0cdrag','DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1469	ENDIF
1470
1471	! Sum of vegetative roughness length and non-vegetative
1472	! roughness length
1473	dragm(:) = dragm(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1474	dragh(:) = dragh(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio/ratio_z0m_z0h(1)))*(ct_karman/LOG(ztmp(:)/z0_nobio))
1475
1476	ENDDO ! Loop over # of non-vegative surfaces
1477
1478	!! 4. Calculate the zero plane displacement height and effective roughness length
1479
1480	! Take the exponential of the roughness
1481	z0m(:) = ztmp(:) / EXP(ct_karman/SQRT(dragm(:)))
1482	z0h(:) = ztmp(:) / EXP((ct_karman*2.)/(dragh(:)LOG(ztmp(:)/z0m(:))))
1483
1484	! Compute the zero plane displacement height which
1485	! is an equivalent height for the absorption of momentum
1486	zhdispl(:) = ave_height(:) * height_displacement
1487	DO jv = 2,nvm
1488	zhdispl_pft(:,jv) = height(:,jv) * height_displacement
1489	ENDDO
1490
1491	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1492	! This is the height over which the roughness acts. It combines the
1493	! zero plane displacement height and the vegetation height.
1494	roughheight(:) = ave_height(:) - zhdispl(:)
1495	DO jv = 2,nvm
1496	roughheight_pft(:,jv) = height(:,jv) - zhdispl_pft(:,jv)
1497	!!! note that we only use roughheight_pft for croplands
1498	ENDDO
1499
1500	END SUBROUTINE condveg_z0cdrag
1501
1502
1503	!! ==============================================================================================================================
1504	!! SUBROUTINE : condveg_z0cdrag_dyn
1505	!!
1506	!>\BRIEF Computation of grid average of roughness length by calculating
1507	!! the drag coefficient based on formulation proposed by Su et al. (2001).
1508	!!
1509	!! DESCRIPTION : This routine calculates the mean roughness height and mean
1510	!! effective roughness height over the grid cell. The mean roughness height (z0)
1511	!! is computed by averaging the drag coefficients \n
1512	!!
1513	!! \latexonly
1514	!! \input{z0cdrag1.tex}
1515	!! \endlatexonly
1516	!! \n
1517	!!
1518	!! where C is the drag coefficient at the height of the vegetation, kappa is the
1519	!! von Karman constant, z (Ztmp) is the height at which the fluxes are estimated and z0 the roughness height.
1520	!! The reference level for z needs to be high enough above the canopy to avoid
1521	!! singularities of the LOG. This height is set to minimum 10m above ground.
1522	!! The drag coefficient increases with roughness height to represent the greater
1523	!! turbulence generated by rougher surfaces.
1524	!! The roughenss height is obtained by the inversion of the drag coefficient equation.\n
1525	!! In the formulation of Su et al. (2001), one distinguishes the roughness height for
1526	!! momentum (z0m) and the one for heat (z0h).
1527	!! z0m is computed as a function of LAI (z0m increases with LAI) and z0h is computed
1528	!! with a so-called kB-1 term (z0m/z0h=exp(kB-1))
1529	!!
1530	!! RECENT CHANGE(S): Written by N. Vuichard (2016)
1531	!!
1532	!! MAIN OUTPUT VARIABLE(S): :: roughness height(z0) and grid effective roughness height(roughheight)
1533	!!
1534	!! REFERENCE(S) :
1535	!! - Su, Z., Schmugge, T., Kustas, W.P., Massman, W.J., 2001. An Evaluation of Two Models for
1536	!! Estimation of the Roughness Height for Heat Transfer between the Land Surface and the Atmosphere. J. Appl.
1537	!! Meteorol. 40, 1933â1951. doi:10.1175/1520-0450(2001)
1538	!! - Ershadi, A., McCabe, M.F., Evans, J.P., Wood, E.F., 2015. Impact of model structure and parameterization
1539	!! on Penman-Monteith type evaporation models. J. Hydrol. 525, 521â535. doi:10.1016/j.jhydrol.2015.04.008
1540	!!
1541	!! FLOWCHART : None
1542	!! \n
1543	!_ ================================================================================================================================
1544
1545	SUBROUTINE condveg_z0cdrag_dyn (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, &
1546	& temp_air, pb, u, v, lai, frac_snow_veg, z0m, z0h, roughheight, roughheight_pft)
1547
1548	!! 0. Variable and parameter declaration
1549
1550	!! 0.1 Input variables
1551
1552	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - Number of land pixels (unitless)
1553	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! PFT coverage fraction of a PFT (= ind*cn_ind)
1554	!! (m^2 m^{-2})
1555	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
1556	!! (= ind*cn_ind) (m^2 m^{-2})
1557	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of non-vegetative surfaces,
1558	!! i.e. continental ice, lakes, etc. (unitless)
1559	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of non-vegetative surfaces,
1560	!! i.e. continental ice, lakes, etc. (unitless)
1561	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
1562	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height !! Vegetation height (m)
1563	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: temp_air !! 2m air temperature (K)
1564	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: pb !! Surface pressure (hPa)
1565	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed in direction u
1566	!! @tex $(m.s^{-1})$ @endtex
1567	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed in direction v
1568	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m2[leaf]/m2[ground])
1569	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
1570	!! 0.2 Output variables
1571
1572	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0m !! Roughness height for momentum (m)
1573	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0h !! Roughness height for heat (m)
1574	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight !! Grid effective roughness height (m)
1575	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: roughheight_pft
1576
1577	!! 0.3 Modified variables
1578
1579	!! 0.4 Local variables
1580
1581	INTEGER(i_std) :: jv !! Loop index over PFTs (unitless)
1582	REAL(r_std), DIMENSION(kjpindex) :: sumveg !! Fraction of bare soil (unitless)
1583	REAL(r_std), DIMENSION(kjpindex) :: ztmp !! Max height of the atmospheric level (m)
1584	REAL(r_std), DIMENSION(kjpindex) :: ave_height !! Average vegetation height (m)
1585	REAL(r_std), DIMENSION(kjpindex) :: zhdispl !! Zero plane displacement height (m)
1586	REAL(r_std), DIMENSION(kjpindex,nvm) :: zhdispl_pft
1587	REAL(r_std) :: z0_nobio !! Roughness height of non-vegetative fraction (m),
1588	!! i.e. continental ice, lakes, etc.
1589	REAL(r_std), DIMENSION(kjpindex) :: z0m_pft !! Roughness height for momentum for a specific PFT
1590	REAL(r_std), DIMENSION(kjpindex) :: z0h_pft !! Roughness height for heat for a specific PFT
1591	REAL(r_std), DIMENSION(kjpindex) :: dragm !! Drag coefficient for momentum
1592	REAL(r_std), DIMENSION(kjpindex) :: dragh !! Drag coefficient for heat
1593	REAL(r_std), DIMENSION(kjpindex) :: eta !! Ratio of friction velocity to the wind speed at the canopy top - See Ershadi et al. (2015)
1594	REAL(r_std), DIMENSION(kjpindex) :: eta_ec !! Within-canopy wind speed profile estimation coefficient - See Ershadi et al. (2015)
1595	REAL(r_std), DIMENSION(kjpindex) :: Ct_star !! Heat transfer coefficient of the soil - see Su et al. (2001)
1596	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))
1597	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))
1598	REAL(r_std), DIMENSION(kjpindex) :: fc !! fractional canopy coverage
1599	REAL(r_std), DIMENSION(kjpindex) :: fs !! fractional soil coverage
1600	REAL(r_std), DIMENSION(kjpindex) :: Reynolds !! Reynolds number
1601	REAL(r_std), DIMENSION(kjpindex) :: wind !! wind Speed (m)
1602	REAL(r_std), DIMENSION(kjpindex) :: u_star !! friction velocity
1603	REAL(r_std), DIMENSION(kjpindex) :: z0_ground !! z0m value used for ground surface
1604	!_ ================================================================================================================================
1605
1606	!! 1. Preliminary calculation
1607
1608	! Set maximal height of first layer
1609	ztmp(:) = MAX(10., zlev(:))
1610
1611	z0_ground(:) = (1.-frac_snow_veg(:))z0_bare + frac_snow_veg(:)z0_bare/10.
1612
1613	! Calculate roughness for non-vegetative surfaces
1614	! with the von Karman constant
1615	dragm(:) = veget_max(:,1) * (ct_karman/LOG(ztmp(:)/z0_ground(:)))**2
1616
1617	wind(:) = SQRT(u(:)u(:)+v(:)v(:))
1618	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG(zlev(:)/z0_ground(:))
1619	Reynolds(:) = z0_ground(:) * u_star(:) &
1620	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1621
1622	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1623
1624	dragh(:) = veget_max(:,1) * (ct_karman/LOG(ztmp(:)/z0_ground(:)))*(ct_karman/LOG(ztmp(:)/(z0_ground(:)/ exp(kBs_m1(:))) ))
1625
1626	! Fraction of bare soil
1627	sumveg(:) = veget_max(:,1)
1628
1629	! Set average vegetation height to zero
1630	ave_height(:) = zero
1631
1632	!! 2. Calculate the mean roughness height
1633
1634	! Calculate the mean roughness height of
1635	! vegetative PFTs over the grid cell
1636	DO jv = 2, nvm
1637
1638	WHERE(veget_max(:,jv) .GT. zero)
1639	! Calculate the average roughness over the grid cell:
1640	! The unitless drag coefficient is per vegetative PFT
1641	! calculated by use of the von Karman constant, the height
1642	! of the first layer and the roughness. The roughness
1643	! is calculated as the vegetation height per PFT
1644	! multiplied by the roughness parameter 'z0_over_height= 1/16'.
1645	! If this scaled value is lower than 0.01 then the value for
1646	! the roughness of bare soil (0.01) is used.
1647	! The sum over all PFTs gives the average roughness
1648	! per grid cell for the vegetative PFTs.
1649	eta(:) = c1 - c2 * exp(-c3 * Cdrag_foliage * lai(:,jv))
1650
1651	z0m_pft(:) = (height(:,jv)(1-height_displacement)(exp(-ct_karman/eta(:))-exp(-ct_karman/(c1-c2)))) &
1652	+ z0_ground(:)
1653
1654	dragm(:) = dragm(:) + veget_max(:,jv) * (ct_karman/LOG(ztmp(:)/z0m_pft(:)))**2
1655
1656	fc(:) = veget(:,jv)/veget_max(:,jv)
1657	fs(:) = 1. - fc(:)
1658
1659	eta_ec(:) = ( Cdrag_foliage * lai(:,jv)) / (2 * eta(:)*eta(:))
1660	wind(:) = SQRT(u(:)u(:)+v(:)v(:))
1661	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG((zlev(:)+(height(:,jv)*(1-height_displacement)))/z0m_pft(:))
1662	Reynolds(:) = z0_ground(:) * u_star(:) &
1663	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1664
1665	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1666	Ct_star(:) = Prandtl*(-2./3.) SQRT(1./Reynolds(:))
1667
1668	WHERE(lai(:,jv) .GT. min_sechiba)
1669	kB_m1(:) = (ct_karman * Cdrag_foliage) / (4 * Ct * eta(:) * (1 - exp(-eta_ec(:)/2.))) * fc(:)**2. &
1670	+ 2fc(:)fs(:) * (ct_karman * eta(:) * z0m_pft(:) / height(:,jv)) / Ct_star(:) &
1671	+ kBs_m1(:) * fs(:)**2.
1672	ELSEWHERE
1673	kB_m1(:) = kBs_m1(:) * fs(:)**2.
1674	ENDWHERE
1675
1676	z0h_pft(:) = z0m_pft(:) / exp(kB_m1(:))
1677
1678	dragh(:) = dragh(:) + veget_max(:,jv) * (ct_karman/LOG(ztmp(:)/z0m_pft(:)))*(ct_karman/LOG(ztmp(:)/z0h_pft(:)))
1679
1680	! Sum of bare soil and fraction vegetated fraction
1681	sumveg(:) = sumveg(:) + veget_max(:,jv)
1682
1683	! Weigh height of vegetation with maximal cover fraction
1684	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1685
1686	ENDWHERE
1687	ENDDO
1688
1689	!! 3. Calculate the mean roughness height of vegetative PFTs over the grid cell
1690
1691	! Search for pixels with vegetated part to normalise
1692	! roughness height
1693	WHERE ( sumveg(:) .GT. min_sechiba )
1694	dragh(:) = dragh(:) / sumveg(:)
1695	dragm(:) = dragm(:) / sumveg(:)
1696	ENDWHERE
1697
1698	! Calculate fraction of roughness for vegetated part
1699	dragh(:) = (un - totfrac_nobio(:)) * dragh(:)
1700	dragm(:) = (un - totfrac_nobio(:)) * dragm(:)
1701
1702	DO jv = 1, nnobio ! Loop over # of non-vegative surfaces
1703
1704	! Set rougness for ice
1705	IF ( jv .EQ. iice ) THEN
1706	z0_nobio = z0_ice
1707	ELSE
1708	WRITE(numout,*) 'jv=',jv
1709	WRITE(numout,*) 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1710	CALL ipslerr_p(3,'condveg_z0cdrag_dyn','DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE','','')
1711	ENDIF
1712
1713	! Sum of vegetative roughness length and non-vegetative roughness length
1714	! Note that z0m could be made dependent of frac_snow_nobio
1715	dragm(:) = dragm(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1716
1717	u_star(:)= ct_karman * MAX(min_wind,wind(:)) / LOG(zlev(:)/z0_nobio)
1718	Reynolds(:) = z0_nobio * u_star(:) &
1719	/ (1.3271e-5 (pb_std/pb(:)) * (temp_air(:)/ZeroCelsius)**(1.81))
1720
1721	kBs_m1(:) = 2.46 * reynolds**(1./4.) - LOG(7.4)
1722
1723	dragh(:) = dragh(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio)) * &
1724	(ct_karman/LOG(ztmp(:)/(z0_nobio/ exp(kBs_m1(:))) ))
1725	ENDDO ! Loop over # of non-vegative surfaces
1726
1727	!! 4. Calculate the zero plane displacement height and effective roughness length
1728	! Take the exponential of the roughness
1729	z0m(:) = ztmp(:) / EXP(ct_karman/SQRT(dragm(:)))
1730	z0h(:) = ztmp(:) / EXP((ct_karman*2.)/(dragh(:)LOG(ztmp(:)/z0m(:))))
1731
1732	! Compute the zero plane displacement height which
1733	! is an equivalent height for the absorption of momentum
1734	zhdispl(:) = ave_height(:) * height_displacement
1735	DO jv = 2,nvm
1736	zhdispl_pft(:,jv) = height(:,jv) * height_displacement
1737	ENDDO
1738
1739	! In order to calculate the fluxes we compute what we call the grid effective roughness height.
1740	! This is the height over which the roughness acts. It combines the
1741	! zero plane displacement height and the vegetation height.
1742	roughheight(:) = ave_height(:) - zhdispl(:)
1743	DO jv = 2,nvm
1744	roughheight_pft(:,jv) = height(:,jv) - zhdispl_pft(:,jv)
1745	!!!! roughheight_pft only used for croplands, xuhui
1746	ENDDO
1747
1748	END SUBROUTINE condveg_z0cdrag_dyn
1749
1750
1751	END MODULE condveg

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-OP-r6850/src_sechiba/condveg.f90

Download in other formats: