Context Navigation

source: branches/publications/ORCHIDEE-MUSLE-r6129/src_sechiba/condveg.f90 @ 6591

Last change on this file since 6591 was 3171, checked in by josefine.ghattas, 9 years ago

Added option ALB_BG_MODIS : read background albedo from file. See ticket #227

Vladislav Bastrikov, Philippe Peylin

Property svn:keywords set to Revision Date HeadURL Date Author Revision

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

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

Download in other formats: