Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/condveg.f90 @ 104

Last change on this file since 104 was 104, checked in by didier.solyga, 13 years ago
Put alb_leaf in order to have its components updated after the calling to getin
File size: 49.2 KB

Line
1	!!
2	!! This module computes surface conditions
3	!! - albedo
4	!! - roughness
5	!! - emissivity
6	!!
7	!! @author Marie-Alice Foujols and Jan Polcher
8	!! @Version : $Revision: 1.30 $, $Date: 2009/01/07 13:39:45 $
9	!!
10	!! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/condveg.f90,v 1.30 2009/01/07 13:39:45 ssipsl Exp $
11	!! IPSL (2006)
12	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
13	!!
14	MODULE condveg
15	!
16	USE ioipsl
17	!
18	! modules used :
19	USE constantes
20	USE pft_parameters
21	USE interpol_help
22	USE parallel
23
24	IMPLICIT NONE
25
26	! public routines :
27	! condveg_main only
28	PRIVATE
29	PUBLIC :: condveg_main,condveg_clear
30
31	!
32	! variables used inside condveg module : declaration and initialisation
33	!
34	LOGICAL, SAVE :: l_first_condveg=.TRUE. !! To keep first call's trace
35	LOGICAL, SAVE :: z0cdrag_ave=.FALSE. !! Chooses the method for the z0 average
36	!
37	REAL(r_std), SAVE :: fixed_snow_albedo !! In case we wish a fxed snow albedo
38	INTEGER(i_std), PARAMETER :: ivis = 1 !! index for visible albedo
39	INTEGER(i_std), PARAMETER :: inir = 2 !! index for near infrared albedo
40	LOGICAL, SAVE :: impaze !! Choice on the surface parameters
41	!
42	REAL(r_std), SAVE :: z0_scal !! Roughness used to initialize the scheme
43	REAL(r_std), SAVE :: roughheight_scal !! Height to displace the surface
44	!! from the zero wind height.
45	REAL(r_std), SAVE :: albedo_scal(2) !! Two albedos used to initialize the scheme
46	REAL(r_std), SAVE :: emis_scal !! Surface emissivity used to initialize the scheme
47	!
48	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_dry(:,:) !! albedo for the dry bare soil
49	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_wet(:,:) !! albedo for the wet bare soil
50	! Ajout Nathalie pour autre calcul soilalbedo
51	REAL(r_std), ALLOCATABLE, SAVE :: soilalb_moy(:,:) !! mean soil albedo.
52	! Ajouts Nathalie - Septembre 2008
53	REAL(r_std), ALLOCATABLE, SAVE :: alb_bare(:,:) !! Soil Albedo, vis(1) and nir(2)
54	REAL(r_std), ALLOCATABLE, SAVE :: alb_veget(:,:) !! Vegetation Albedo, vis(1) and nir(2)
55	!
56	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: albedo_snow !! Snow albedo
57	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: albedo_glob !! Mean albedo
58	!
59	LOGICAL, SAVE :: alb_bare_model !! Switch to old (albedo bare depend on soil wetness)
60	!! or new one (mean of soilalb)
61
62	CONTAINS
63
64	!!
65	!! Main routine for condveg module
66	!! - called only one time for initialisation
67	!! - called every time step
68	!! - called one more time at last time step for writing _restart_ file
69	!!
70	!! Algorithm:
71	!! - call condveg_init for initialisation
72	!! - call condveg_var_init for initialisation done every time step
73	!! - call condveg_snow for computing the modification of the albedo induced by snow cover
74	!!
75	!! @call condveg_init
76	!! @call condveg_var_init
77	!! @call condveg_snow
78	!!
79	!!
80	SUBROUTINE condveg_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index,&
81	& lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
82	& zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
83	& drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)
84
85	! interface description:
86	! input scalar
87	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
88	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
89	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
90	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier
91	INTEGER(i_std), OPTIONAL, INTENT (in) :: hist2_id !! _History_ file 2 identifier
92	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
93	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for restart file to read
94	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for restart file to write
95	! input fields
96	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
97	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
98	INTEGER(i_std),DIMENSION (kjpindex,8), INTENT(in):: neighbours !! neighoring grid points if land
99	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
100	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
101	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
102	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
103	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
104	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! total fraction of continental ice+lakes+...
105	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
106	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2]
107	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
108	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2] on ice, lakes, ...
109	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
110	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
111	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
112	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: deadleaf_cover !! Fraction of soil covered by dead leaves
113	! output scalar
114	! output fields
115	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
116	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
117	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Roughness
118	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
119	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! fraction of soil types
120	! local
121	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
122	!
123	IF (l_first_condveg) THEN
124
125	IF (long_print) WRITE (numout,*) ' l_first_condveg : call condveg_init '
126
127	CALL condveg_init (kjit, ldrestart_read, kjpindex, index, veget, &
128	lalo, neighbours, resolution, contfrac, rest_id)
129	CALL condveg_var_init (ldrestart_read, kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
130	& drysoil_frac, zlev, height, deadleaf_cover, emis, albedo, z0, roughheight)
131	! Nathalie - Fevrier 2007 - c'est veget_max qu'il faut passer
132	! Sauf en cas de DGVM .... pour le moment car la somme des maxvegetfrac depasse 1!
133	! Si DGVM alors on passe veget.
134	IF (control%ok_dgvm) THEN
135	CALL condveg_snow (ldrestart_read, kjpindex, veget, frac_nobio, totfrac_nobio, &
136	snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
137	ELSE
138	CALL condveg_snow (ldrestart_read, kjpindex, veget_max, frac_nobio, totfrac_nobio, &
139	snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
140	ENDIF
141
142	RETURN
143
144	ENDIF
145
146	CALL condveg_var_update (ldrestart_read, kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
147	& drysoil_frac, zlev, height, deadleaf_cover, emis, albedo, z0, roughheight)
148	! Nathalie - Fevrier 2007 - c'est veget_max qu'il faut passer
149	! Sauf en cas de DGVM .... pour le moment car la somme des maxvegetfrac depasse 1!
150	! Si DGVM alors on passe veget.
151	IF (control%ok_dgvm) THEN
152	CALL condveg_snow (ldrestart_read, kjpindex, veget, frac_nobio, totfrac_nobio, &
153	snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
154	ELSE
155	CALL condveg_snow (ldrestart_read, kjpindex, veget_max, frac_nobio, totfrac_nobio, &
156	snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
157	ENDIF
158
159	IF (ldrestart_write) THEN
160	!
161	var_name = 'soilalbedo_dry'
162	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_dry, 'scatter', nbp_glo, index_g)
163	!
164	var_name = 'soilalbedo_wet'
165	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_wet, 'scatter', nbp_glo, index_g)
166	!
167	var_name = 'soilalbedo_moy'
168	CALL restput_p (rest_id, var_name, nbp_glo, 2, 1, kjit, soilalb_moy, 'scatter', nbp_glo, index_g)
169	!
170	RETURN
171	!
172	ENDIF
173
174	IF ( almaoutput ) THEN
175	CALL histwrite(hist_id, 'Albedo', kjit, albedo_glob, kjpindex, index)
176	CALL histwrite(hist_id, 'SAlbedo', kjit, albedo_snow, kjpindex, index)
177	IF ( hist2_id > 0 ) THEN
178	CALL histwrite(hist2_id, 'Albedo', kjit, albedo_glob, kjpindex, index)
179	CALL histwrite(hist2_id, 'SAlbedo', kjit, albedo_snow, kjpindex, index)
180	ENDIF
181	ELSE
182	CALL histwrite(hist_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
183	CALL histwrite(hist_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
184	CALL histwrite(hist_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
185	CALL histwrite(hist_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
186	IF ( hist2_id > 0 ) THEN
187	CALL histwrite(hist2_id, 'soilalb_vis', kjit, alb_bare(:,1), kjpindex, index)
188	CALL histwrite(hist2_id, 'soilalb_nir', kjit, alb_bare(:,2), kjpindex, index)
189	CALL histwrite(hist2_id, 'vegalb_vis', kjit, alb_veget(:,1), kjpindex, index)
190	CALL histwrite(hist2_id, 'vegalb_nir', kjit, alb_veget(:,2), kjpindex, index)
191	ENDIF
192	ENDIF
193
194	IF (long_print) WRITE (numout,*)' condveg_main done '
195
196	END SUBROUTINE condveg_main
197
198	!! Algorithm:
199	!! - dynamic allocation for local array
200	!!
201	SUBROUTINE condveg_init (kjit, ldrestart_read, kjpindex, index, veget, &
202	lalo, neighbours, resolution, contfrac, rest_id)
203
204	! interface description
205	! input scalar
206	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
207	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for restart file to read
208	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
209	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
210	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in):: veget !! Vegetation distribution
211	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
212	INTEGER(i_std),DIMENSION (kjpindex,4), INTENT(in):: neighbours !! neighoring grid points if land
213	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
214	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: contfrac ! Fraction of land in each grid box.
215	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier
216	! input fields
217	! output scalar
218	! output fields
219
220	! local declaration
221	INTEGER(i_std) :: ji
222	INTEGER(i_std) :: ier
223	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
224	! initialisation
225	IF (l_first_condveg) THEN
226	!
227	! Get the fixed snow albedo if needed
228	!
229	!
230	!Config Key = CONDVEG_SNOWA
231	!Config Desc = The snow albedo used by SECHIBA
232	!Config Def = DEF
233	!Config Help = This option allows the user to impose a snow albedo.
234	!Config Default behaviour is to use the model of snow albedo
235	!Config developed by Chalita (1993).
236	!
237	fixed_snow_albedo = undef_sechiba
238	CALL getin_p('CONDVEG_SNOWA', fixed_snow_albedo)
239	!
240	!
241	!Config Key = ALB_BARE_MODEL
242	!Config Desc = Switch bare soil albedo dependent (if TRUE) on soil wetness
243	!Config Def = FALSE
244	!Config Help = If TRUE, the model for bare soil albedo is the old formulation.
245	!Config Then it depend on the soil dry or wetness. If FALSE, it is the
246	!Config new computation that is taken, it is the mean of soil albedo.
247	!
248	alb_bare_model=.FALSE.
249	CALL getin_p('ALB_BARE_MODEL', alb_bare_model)
250	!
251	l_first_condveg=.FALSE.
252	!
253	! Allocate variables which have to
254	!
255	ALLOCATE (soilalb_dry(kjpindex,2),stat=ier)
256	IF (ier.NE.0) THEN
257	WRITE (numout,) ' error in soilalb_dry allocation. We stop.We need kjpindex2 words = ',kjpindex*2
258	STOP 'condveg_init'
259	END IF
260	soilalb_dry(:,:) = val_exp
261
262	ALLOCATE (soilalb_wet(kjpindex,2),stat=ier)
263	IF (ier.NE.0) THEN
264	WRITE (numout,) ' error in soilalb_wet allocation. We stop.We need kjpindex2 words = ',kjpindex*2
265	STOP 'condveg_init'
266	END IF
267	soilalb_wet(:,:) = val_exp
268
269	! Ajout Nathalie
270	ALLOCATE (soilalb_moy(kjpindex,2),stat=ier)
271	IF (ier.NE.0) THEN
272	WRITE (numout,) ' error in soilalb_moy allocation. We stop.We need kjpindex2 words = ',kjpindex*2
273	STOP 'condveg_init'
274	END IF
275	soilalb_moy(:,:) = val_exp
276
277	ALLOCATE (albedo_snow(kjpindex),stat=ier)
278	IF (ier.NE.0) THEN
279	WRITE (numout,*) ' error in albedo_snow allocation. We stop.We need kjpindex words = ',kjpindex
280	STOP 'condveg_init'
281	END IF
282
283	ALLOCATE (albedo_glob(kjpindex),stat=ier)
284	IF (ier.NE.0) THEN
285	WRITE (numout,*) ' error in albedo_glob allocation. We stop.We need kjpindex words = ',kjpindex
286	STOP 'condveg_init'
287	END IF
288
289	ALLOCATE (alb_bare(kjpindex,2),stat=ier)
290	IF (ier.NE.0) THEN
291	WRITE (numout,) ' error in alb_bare allocation. We stop.We need kjpindex2 words = ',kjpindex*2
292	STOP 'condveg_init'
293	END IF
294	alb_bare(:,:) = val_exp
295
296	ALLOCATE (alb_veget(kjpindex,2),stat=ier)
297	IF (ier.NE.0) THEN
298	WRITE (numout,) ' error in alb_veget allocation. We stop.We need kjpindex2 words = ',kjpindex*2
299	STOP 'condveg_init'
300	END IF
301	alb_veget(:,:) = val_exp
302	!
303	! Get the bare soil albedo
304	!
305	!
306	var_name= 'soilalbedo_dry'
307	CALL ioconf_setatt('UNITS', '-')
308	CALL ioconf_setatt('LONG_NAME','Dry bare soil albedo')
309	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_dry, "gather", nbp_glo, index_g)
310
311	var_name= 'soilalbedo_wet'
312	CALL ioconf_setatt('UNITS', '-')
313	CALL ioconf_setatt('LONG_NAME','Wet bare soil albedo')
314	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_wet, "gather", nbp_glo, index_g)
315
316	var_name= 'soilalbedo_moy'
317	CALL ioconf_setatt('UNITS', '-')
318	CALL ioconf_setatt('LONG_NAME','Mean bare soil albedo')
319	CALL restget_p (rest_id, var_name, nbp_glo, 2, 1, kjit, .TRUE., soilalb_moy, "gather", nbp_glo, index_g)
320	!
321	IF ( MINVAL(soilalb_wet) .EQ. MAXVAL(soilalb_wet) .AND. MAXVAL(soilalb_wet) .EQ. val_exp .OR.&
322	& MINVAL(soilalb_dry) .EQ. MAXVAL(soilalb_dry) .AND. MAXVAL(soilalb_dry) .EQ. val_exp .OR.&
323	& MINVAL(soilalb_moy) .EQ. MAXVAL(soilalb_moy) .AND. MAXVAL(soilalb_moy) .EQ. val_exp) THEN
324	CALL condveg_soilalb(kjpindex, lalo, neighbours, resolution, contfrac, soilalb_dry,soilalb_wet)
325	WRITE(numout,*) '---> val_exp ', val_exp
326	WRITE(numout,*) '---> ALBEDO_wet VIS:', MINVAL(soilalb_wet(:,ivis)), MAXVAL(soilalb_wet(:,ivis))
327	WRITE(numout,*) '---> ALBEDO_wet NIR:', MINVAL(soilalb_wet(:,inir)), MAXVAL(soilalb_wet(:,inir))
328	WRITE(numout,*) '---> ALBEDO_dry VIS:', MINVAL(soilalb_dry(:,ivis)), MAXVAL(soilalb_dry(:,ivis))
329	WRITE(numout,*) '---> ALBEDO_dry NIR:', MINVAL(soilalb_dry(:,inir)), MAXVAL(soilalb_dry(:,inir))
330	WRITE(numout,*) '---> ALBEDO_moy VIS:', MINVAL(soilalb_moy(:,ivis)), MAXVAL(soilalb_moy(:,ivis))
331	WRITE(numout,*) '---> ALBEDO_moy NIR:', MINVAL(soilalb_moy(:,inir)), MAXVAL(soilalb_moy(:,inir))
332	ENDIF
333	!
334	ELSE
335	WRITE (numout,*) ' l_first_condveg false . we stop '
336	STOP 'condveg_init'
337	ENDIF
338	!! test de commentaires
339
340	IF (long_print) WRITE (numout,*) ' condveg_init done '
341
342	END SUBROUTINE condveg_init
343	!!
344	!!
345	SUBROUTINE condveg_clear ()
346
347	l_first_condveg=.TRUE.
348
349	IF (ALLOCATED (soilalb_dry)) DEALLOCATE (soilalb_dry)
350	IF (ALLOCATED(soilalb_wet)) DEALLOCATE (soilalb_wet)
351	! Ajout Nathalie
352	IF (ALLOCATED(soilalb_moy)) DEALLOCATE (soilalb_moy)
353	IF (ALLOCATED(albedo_snow)) DEALLOCATE (albedo_snow)
354	IF (ALLOCATED(albedo_glob)) DEALLOCATE (albedo_glob)
355	IF (ALLOCATED(alb_bare)) DEALLOCATE (alb_bare)
356	IF (ALLOCATED(alb_veget)) DEALLOCATE (alb_veget)
357	!
358	END SUBROUTINE condveg_clear
359
360	!! Algorithm:
361	!! - initialisation of local arry
362	!! - reads map for emissivity, albedo or roughness
363	!!
364	SUBROUTINE condveg_var_init (ldrestart_read, kjpindex, veget, veget_max, frac_nobio, totfrac_nobio,&
365	& drysoil_frac, zlev, height, deadleaf_cover, emis, albedo, z0, roughheight)
366
367	! interface description
368	! input scalar
369	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
370	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
371	! input fields
372	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
373	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
374	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
375	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+...
376	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
377	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
378	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
379	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: deadleaf_cover !! Fraction of soil covered by dead leaves
380	! output scalar
381	! output fields
382	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
383	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
384	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Roughness
385	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
386	!
387	! local declaration
388	INTEGER(i_std) :: ier !! Error code
389	INTEGER(i_std) :: jv
390	! initialisation of variables
391	!
392
393	!
394	!Config Key = IMPOSE_AZE
395	!Config Desc = Should the surface parameters be prescribed
396	!Config Def = n
397	!Config Help = This flag allows the user to impose the surface parameters
398	!Config (Albedo Roughness and Emissivity). It is espacially interesting for 0D
399	!Config simulations. On the globe it does not make too much sense as
400	!Config it imposes the same vegetation everywhere
401	!
402	impaze = .FALSE.
403	CALL getin_p('IMPOSE_AZE', impaze)
404
405	!!
406	!! calculs de emis
407	!!
408
409	IF ( impaze ) THEN
410	!
411	!Config Key = CONDVEG_EMIS
412	!Config Desc = Emissivity of the surface for LW radiation
413	!Config Def = 1.0
414	!Config If = IMPOSE_AZE
415	!Config Help = The surface emissivity used for compution the LE emission
416	!Config of the surface in a 0-dim version. Values range between
417	!Config 0.97 and 1.. The GCM uses 0.98.
418	!
419	emis_scal = un
420	CALL getin_p('CONDVEG_EMIS', emis_scal)
421	emis(:) = emis_scal
422
423	ELSE
424	! Some day it will be moisture dependent
425	emis_scal = un
426
427	emis(:) = emis_scal
428	ENDIF
429
430	!!
431	!! calculs de albedo
432	!!
433	!
434	IF ( impaze ) THEN
435	!
436	!Config Key = CONDVEG_ALBVIS
437	!Config Desc = SW visible albedo for the surface
438	!Config Def = 0.25
439	!Config If = IMPOSE_AZE
440	!Config Help = Surface albedo in visible wavelengths to be used
441	!Config on the point if a 0-dim version of SECHIBA is used.
442	!Config Look at the description of the forcing data for
443	!Config the correct value.
444	!
445	albedo_scal(ivis) = 0.25_r_std
446	CALL getin_p('CONDVEG_ALBVIS', albedo_scal(ivis))
447	albedo(:,ivis) = albedo_scal(ivis)
448	!
449	!Config Key = CONDVEG_ALBNIR
450	!Config Desc = SW near infrared albedo for the surface
451	!Config Def = 0.25
452	!Config If = IMPOSE_AZE
453	!Config Help = Surface albedo in near infrared wavelengths to be used
454	!Config on the point if a 0-dim version of SECHIBA is used.
455	!Config Look at the description of the forcing data for
456	!Config the correct value.
457	!
458	albedo_scal(inir) = 0.25_r_std
459	CALL getin_p('CONDVEG_ALBNIR', albedo_scal(inir))
460	albedo(:,inir) = albedo_scal(inir)
461	ELSE
462	!
463	CALL condveg_albcalc (kjpindex,deadleaf_cover,veget,drysoil_frac,albedo)
464	!
465	ENDIF
466	!!
467	!! calculs de z0
468	!!
469	!
470	!Config Key = Z0CDRAG_AVE
471	!Config Desc = Average method for z0
472	!Config Def = y
473	!Config Help = If this flag is set to true (y) then the neutral Cdrag
474	!Config is averaged instead of the log(z0). This should be
475	!Config the prefered option. We still wish to keep the other
476	!Config option so we can come back if needed. If this is
477	!Config desired then one should set Z0CDRAG_AVE=n
478	z0cdrag_ave = .TRUE.
479	CALL getin_p('Z0CDRAG_AVE', z0cdrag_ave)
480	!!
481	IF ( impaze ) THEN
482	!
483	!Config Key = CONDVEG_Z0
484	!Config Desc = Surface roughness (m)
485	!Config Def = 0.15
486	!Config If = IMPOSE_AZE
487	!Config Help = Surface rougness to be used on the point if a 0-dim version
488	!Config of SECHIBA is used. Look at the description of the forcing
489	!Config data for the correct value.
490	!
491	z0_scal = 0.15_r_std
492	CALL getin_p('CONDVEG_Z0', z0_scal)
493	z0(:) = z0_scal
494	!
495	!Config Key = ROUGHHEIGHT
496	!Config Desc = Height to be added to the height of the first level (m)
497	!Config Def = 0.0
498	!Config If = IMPOSE_AZE
499	!Config Help = ORCHIDEE assumes that the atmospheric level height is counted
500	!Config from the zero wind level. Thus to take into account the roughness
501	!Config of tall vegetation we need to correct this by a certain fraction
502	!Config of the vegetation height. This is called the roughness height in
503	!Config ORCHIDEE talk.
504	!
505	roughheight_scal = zero
506	CALL getin_p('ROUGHHEIGHT', roughheight_scal)
507	roughheight(:) = roughheight_scal
508	!
509	ELSE
510	!
511	IF ( z0cdrag_ave ) THEN
512	CALL condveg_z0cdrag(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, &
513	& height, z0, roughheight)
514	ELSE
515	CALL condveg_z0logz(kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
516	& z0, roughheight)
517	ENDIF
518	!
519	ENDIF
520	!!
521	!!
522	IF (long_print) WRITE (numout,*) ' condveg_var_init done '
523
524	END SUBROUTINE condveg_var_init
525	!!
526	!! Algorithm:
527	!! - Simply update the emissivity, albedo and roughness fields
528	!!
529	SUBROUTINE condveg_var_update (ldrestart_read, kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, &
530	& drysoil_frac, zlev, height, deadleaf_cover, emis, albedo, z0, roughheight)
531
532	! interface description
533	! input scalar
534	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
535	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
536	! input fields
537	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
538	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: veget_max !! Fraction of vegetation type
539	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
540	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+...
541	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
542	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
543	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: height !! Vegetation Height (m)
544	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: deadleaf_cover !! Fraction of soil covered by dead leaves
545	! output scalar
546	! output fields
547	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis !! Emissivity
548	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
549	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0 !! Roughness
550	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: roughheight !! Effective height for roughness
551	!
552	! local
553	INTEGER(i_std) :: ji, jv
554	!!
555	!! calculs de emis
556	!!
557
558	emis(:) = emis_scal
559
560	!!
561	!! calculs de albedo
562	!!
563	!
564	IF ( impaze ) THEN
565	!
566	albedo(:,ivis) = albedo_scal(ivis)
567	albedo(:,inir) = albedo_scal(inir)
568	!
569	ELSE
570	!
571	CALL condveg_albcalc (kjpindex,deadleaf_cover,veget,drysoil_frac,albedo)
572	!
573	ENDIF
574	!
575	!!
576	!! Calculs de la rugosite
577	!!
578
579	IF ( impaze ) THEN
580
581	DO ji = 1, kjpindex
582	z0(ji) = z0_scal
583	roughheight(ji) = roughheight_scal
584	ENDDO
585
586	ELSE
587	!
588	IF ( z0cdrag_ave ) THEN
589	CALL condveg_z0cdrag (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, zlev, height, &
590	& z0, roughheight)
591	ELSE
592	CALL condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
593	& z0, roughheight)
594	ENDIF
595	!
596	ENDIF
597
598	IF (long_print) WRITE (numout,*) ' condveg_var_update done '
599
600	END SUBROUTINE condveg_var_update
601
602	!! Algorithm:
603	!! - The snow albedo is updated by the snow within the mesh
604	!! This is done as a function of snow mass, snow age and vegetation type
605	!! The model is described in Chalita 1992
606	!!
607	SUBROUTINE condveg_snow (ldrestart_read, kjpindex, veget, frac_nobio, totfrac_nobio, &
608	snow, snow_age, snow_nobio, snow_nobio_age, albedo, albedo_snow, albedo_glob)
609
610	! interface description
611	! input scalar
612	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
613	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
614	! input fields
615	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation types
616	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
617	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
618	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow !! Snow mass [Kg/m^2] in vegetation
619	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow mass [Kg/m^2]on ice, lakes, ...
620	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: snow_age !! Snow age
621	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio_age !! Snow age on ice, lakes, ...
622	! output scalar
623	! output fields
624	REAL(r_std),DIMENSION (kjpindex,2), INTENT (inout) :: albedo !! Albedo
625	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_snow !! Albedo de la neige
626	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: albedo_glob !! Mean albedo
627	!
628	! local declaration
629	INTEGER(i_std) :: ji, jv, jb !! Loop index
630	REAL(r_std), DIMENSION(kjpindex) :: frac_snow_veg !! Fraction of snow on vegetation
631	REAL(r_std), DIMENSION(kjpindex,nnobio) :: frac_snow_nobio !! Fraction of snow on ice, lakes, ...
632	REAL(r_std), DIMENSION(kjpindex) :: snowa_veg !! Total albedo of snow covered area on vegetation
633	REAL(r_std), DIMENSION(kjpindex,nnobio) :: snowa_nobio !! albedo of snow covered area on ice, lakes, ...
634	REAL(r_std), DIMENSION(kjpindex) :: fraction_veg !! total vegetation fraction
635	REAL(r_std), DIMENSION(kjpindex) :: agefunc_veg !! age dependency of snow albedo on vegetation
636	REAL(r_std), DIMENSION(kjpindex,nnobio) :: agefunc_nobio !! age dependency of snow albedo on ice, lakes, ..;
637	REAL(r_std) :: alb_nobio
638	!
639
640	DO ji = 1, kjpindex
641	albedo_snow(ji) = zero
642	albedo_glob(ji) = zero
643	ENDDO
644
645	IF (ABS(fixed_snow_albedo - undef_sechiba) .GT. EPSILON(undef_sechiba)) THEN
646	snowa_veg(:) = fixed_snow_albedo
647	snowa_nobio(:,:) = fixed_snow_albedo
648	ELSE
649	!
650	! calculate first age dependence
651	!
652	DO ji = 1, kjpindex
653	agefunc_veg(ji) = EXP(-snow_age(ji)/tcst_snowa)
654	ENDDO
655	!
656	!
657	DO jv = 1, nnobio
658	DO ji = 1, kjpindex
659	agefunc_nobio(ji,jv) = EXP(-snow_nobio_age(ji,jv)/tcst_snowa)
660	ENDDO
661	ENDDO
662	!
663	! snow albedo on vegetated surfaces
664	!
665	fraction_veg(:) = 1. - totfrac_nobio(:)
666	snowa_veg(:) = 0.
667	DO jv = 1, nvm
668	DO ji = 1, kjpindex
669	IF ( fraction_veg(ji) .GT. min_sechiba ) THEN
670	snowa_veg(ji) = snowa_veg(ji) + &
671	veget(ji,jv)/fraction_veg(ji) * ( snowa_ini(jv)+snowa_dec(jv)*agefunc_veg(ji) )
672	ENDIF
673	ENDDO
674	ENDDO
675	!
676	! snow albedo on other surfaces
677	!
678	DO jv = 1, nnobio
679	DO ji = 1, kjpindex
680	snowa_nobio(ji,jv) = ( snowa_ini(1) + snowa_dec(1) * agefunc_nobio(ji,jv) )
681	ENDDO
682	ENDDO
683	ENDIF
684
685	frac_snow_veg(:) = MIN(MAX(snow(:),zero)/(MAX(snow(:),zero)+snowcri_alb),un)
686	DO jv = 1, nnobio
687	frac_snow_nobio(:,jv) = MIN(MAX(snow_nobio(:,jv),zero)/(MAX(snow_nobio(:,jv),zero)+snowcri_alb),un)
688	ENDDO
689
690
691	DO jb = 1, 2
692	!
693	albedo(:,jb) = ( fraction_veg(:) ) * &
694	( (un-frac_snow_veg(:)) * albedo(:,jb) + &
695	( frac_snow_veg(:) ) * snowa_veg(:) )
696	albedo_snow(:) = albedo_snow(:) + (fraction_veg(:)) * (frac_snow_veg(:)) * snowa_veg(:)
697	!
698	DO jv = 1, nnobio
699	!
700	IF ( jv .EQ. iice ) THEN
701	alb_nobio = alb_ice(jb)
702	ELSE
703	WRITE(numout,*) 'jv=',jv
704	STOP 'DO NOT KNOW ALBEDO OF THIS SURFACE TYPE'
705	ENDIF
706	!
707	albedo(:,jb) = albedo(:,jb) + &
708	( frac_nobio(:,jv) ) * &
709	( (un-frac_snow_nobio(:,jv)) * alb_nobio + &
710	( frac_snow_nobio(:,jv) ) * snowa_nobio(:,jv) )
711	albedo_snow(:) = albedo_snow(:) + &
712	( frac_nobio(:,jv) ) * ( frac_snow_nobio(:,jv) ) * &
713	snowa_nobio(:,jv)
714	albedo_glob(:) = albedo_glob(:) + albedo(:,jb)
715	!
716	ENDDO
717	!
718	END DO
719	!
720	DO ji = 1, kjpindex
721	albedo_snow(ji) = (albedo_snow(ji))/2.
722	albedo_glob(ji) = (albedo_glob(ji))/2.
723	ENDDO
724
725	IF (long_print) WRITE (numout,*) ' condveg_snow done '
726
727	END SUBROUTINE condveg_snow
728	!
729	!
730	!
731	SUBROUTINE condveg_soilalb(nbpt, lalo, neighbours, resolution, contfrac, soilalb_dry,soilalb_wet)
732	!
733	!
734	! This subroutine should read the soil color maps from the Henderson-Sellers & Wilson database. This data
735	! is then interpolated to the models resolution and transformed into dry and wet albedos. We have chosen to
736	! do both these operations at the same time as one can average the albedo but not the color types.
737	!
738	! We make the assumption in this code that the grid of the data is regular and that it covers the globe.
739	! For the model grid we base the calculation of the borders of the grid on the resolution.
740	!
741	!
742	! 0.1 INPUT
743	!
744	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
745	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
746	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
747	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
748	REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
749	REAL(r_std), INTENT(inout) :: soilalb_dry(nbpt,2) ! albedo for the dry bare soil
750	REAL(r_std), INTENT(inout) :: soilalb_wet(nbpt,2) ! albedo for the wet bare soil
751
752	!
753	!
754	! 0.3 LOCAL
755	!
756	INTEGER(i_std) :: nbvmax
757	!
758	CHARACTER(LEN=80) :: filename
759	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, ilf, lastjp, nbexp
760	REAL(r_std) :: lev(1), date, dt, coslat, sgn
761	INTEGER(i_std) :: itau(1)
762	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, soilcol
763	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
764	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
765	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
766	CHARACTER(LEN=30) :: callsign
767	!
768	LOGICAL :: ok_interpol = .FALSE. ! optionnal return of aggregate_2d
769	!
770	INTEGER :: ALLOC_ERR
771	!
772	!
773	! Needs to be a configurable variable
774	!
775	!
776	!Config Key = SOILALB_FILE
777	!Config Desc = Name of file from which the bare soil albedo
778	!Config Def = ../surfmap/soils_param.nc
779	!Config If = !IMPOSE_AZE
780	!Config Help = The name of the file to be opened to read the soil types from
781	!Config which we derive then the bare soil albedos. This file is 1x1
782	!Config deg and based on the soil colors defined by Wilson and Henderson-Seller.
783	!
784	filename = '../surfmap/soils_param.nc'
785	CALL getin_p('SOILALB_FILE',filename)
786	!
787	IF (is_root_prc) CALL flininfo(filename,iml, jml, lml, tml, fid)
788	CALL bcast(iml)
789	CALL bcast(jml)
790	CALL bcast(lml)
791	CALL bcast(tml)
792	!
793	! soils_param.nc file is 1° soit texture file.
794	!
795	ALLOC_ERR=-1
796	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
797	IF (ALLOC_ERR/=0) THEN
798	PRINT *,"ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
799	STOP
800	ENDIF
801	ALLOC_ERR=-1
802	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
803	IF (ALLOC_ERR/=0) THEN
804	PRINT *,"ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
805	STOP
806	ENDIF
807	ALLOC_ERR=-1
808	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
809	IF (ALLOC_ERR/=0) THEN
810	PRINT *,"ERROR IN ALLOCATION of mask : ",ALLOC_ERR
811	STOP
812	ENDIF
813	ALLOC_ERR=-1
814	ALLOCATE(soilcol(iml,jml), STAT=ALLOC_ERR)
815	IF (ALLOC_ERR/=0) THEN
816	PRINT *,"ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
817	STOP
818	ENDIF
819	!
820	IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
821	CALL bcast(lon_rel)
822	CALL bcast(lat_rel)
823	CALL bcast(lev)
824	CALL bcast(itau)
825	CALL bcast(date)
826	CALL bcast(dt)
827	!
828	IF (is_root_prc) CALL flinget(fid, 'soilcolor', iml, jml, lml, tml, 1, 1, soilcol)
829	CALL bcast(soilcol)
830	!
831	IF (is_root_prc) CALL flinclo(fid)
832	!
833	! Mask of permitted variables.
834	!
835	mask(:,:) = zero
836	DO ip=1,iml
837	DO jp=1,jml
838	IF (soilcol(ip,jp) > min_sechiba) THEN
839	mask(ip,jp) = un
840	ENDIF
841	ENDDO
842	ENDDO
843	!
844	nbvmax = 200
845	!
846	callsign = 'Soil color map'
847	!
848	DO WHILE ( .NOT. ok_interpol )
849	WRITE(numout,*) "Projection arrays for ",callsign," : "
850	WRITE(numout,*) "nbvmax = ",nbvmax
851	!
852	ALLOC_ERR=-1
853	ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
854	IF (ALLOC_ERR/=0) THEN
855	PRINT *,"ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
856	STOP
857	ENDIF
858	sub_area(:,:)=zero
859	ALLOC_ERR=-1
860	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
861	IF (ALLOC_ERR/=0) THEN
862	PRINT *,"ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
863	STOP
864	ENDIF
865	sub_index(:,:,:)=0
866	!
867	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
868	& iml, jml, lon_rel, lat_rel, mask, callsign, &
869	& nbvmax, sub_index, sub_area, ok_interpol)
870	!
871	IF ( .NOT. ok_interpol ) THEN
872	DEALLOCATE(sub_area)
873	DEALLOCATE(sub_index)
874	ENDIF
875	!
876	nbvmax = nbvmax * 2
877	ENDDO
878	!
879	nbexp = 0
880	!
881	! Check that we found some points
882	!
883	soilalb_dry(:,:) = zero
884	soilalb_wet(:,:) = zero
885	! Ajout Nathalie
886	soilalb_moy(:,:) = zero
887	!
888	DO ib=1,nbpt
889	!
890	! GO through the point we have found
891	!
892	!
893	fopt = COUNT(sub_area(ib,:) > zero)
894	!
895	IF ( fopt .EQ. 0) THEN
896	nbexp = nbexp + 1
897	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
898	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
899	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
900	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
901	! Ajout Nathalie
902	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
903	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
904	ELSE
905	sgn = zero
906	!
907	! Compute the average bare soil albedo parameters
908	!
909	DO ilf = 1,fopt
910	!
911	ip = sub_index(ib,ilf,1)
912	jp = sub_index(ib,ilf,2)
913	!
914	IF ( NINT(soilcol(ip,jp)) .LE. classnb) THEN
915	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
916	soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
917	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
918	soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
919	! Ajout Nathalie
920	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
921	soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
922	sgn = sgn + sub_area(ib,ilf)
923	ELSE
924	WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
925	STOP
926	ENDIF
927	!
928	ENDDO
929	!
930	! Normalize the surface
931	!
932	IF ( sgn .LT. min_sechiba) THEN
933	nbexp = nbexp + 1
934	soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
935	soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
936	soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
937	soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
938	! Ajout Nathalie
939	soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
940	soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
941	ELSE
942	soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis)/sgn
943	soilalb_dry(ib,inir) = soilalb_dry(ib,inir)/sgn
944	soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis)/sgn
945	soilalb_wet(ib,inir) = soilalb_wet(ib,inir)/sgn
946	! Ajout Nathalie
947	soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis)/sgn
948	soilalb_moy(ib,inir) = soilalb_moy(ib,inir)/sgn
949	ENDIF
950	!
951	ENDIF
952	!
953	ENDDO
954	!
955	IF ( nbexp .GT. 0 ) THEN
956	WRITE(numout,*) 'CONDVEG_soilalb : The interpolation of the bare soil albedo had ', nbexp
957	WRITE(numout,*) 'CONDVEG_soilalb : points without data. This are either coastal points or'
958	WRITE(numout,*) 'CONDVEG_soilalb : ice covered land.'
959	WRITE(numout,*) 'CONDVEG_soilalb : The problem was solved by using the average of all soils'
960	WRITE(numout,*) 'CONDVEG_soilalb : in dry and wet conditions'
961	ENDIF
962	!
963	DEALLOCATE (lat_rel)
964	DEALLOCATE (lon_rel)
965	DEALLOCATE (mask)
966	DEALLOCATE (sub_index)
967	DEALLOCATE (sub_area)
968	DEALLOCATE (soilcol)
969	!
970	!
971	RETURN
972	!
973	END SUBROUTINE condveg_soilalb
974	!
975	!
976	!
977	SUBROUTINE condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
978	& z0, roughheight)
979
980	!
981	! 0. Declarations
982	!
983
984	! 0.1 Input
985	INTEGER(i_std), INTENT(in) :: kjpindex
986	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget
987	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max
988	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio
989	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio
990	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height
991
992	! 0.2 Output
993	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0
994	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight
995
996	! 0.3 Local
997	INTEGER(i_std) :: jv
998	REAL(r_std), DIMENSION(kjpindex) :: sumveg, ave_height
999	REAL(r_std), DIMENSION(kjpindex) :: d_veg, zhdispl
1000	REAL(r_std) :: z0_nobio
1001
1002	z0(:) = veget(:,1) * LOG(z0_bare)
1003	sumveg(:) = veget(:,1)
1004	ave_height(:) = zero
1005	!
1006	DO jv = 2, nvm
1007	!
1008	IF ( is_tree(jv) ) THEN
1009	! tree trunks influence the atmosphere even when there are no leaves
1010	d_veg(:) = veget_max(:,jv)
1011	ELSE
1012	! grasses only have an influence if they are really there!
1013	d_veg(:) = veget(:,jv)
1014	ENDIF
1015	!
1016	z0(:) = z0(:) + d_veg(:) * &
1017	LOG( MAX(height(:,jv)*z0_over_height,z0_bare) )
1018	sumveg(:) = sumveg(:) + d_veg(:)
1019	!
1020	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1021	!
1022	ENDDO
1023	!
1024	WHERE ( sumveg(:) > zero ) z0(:) = z0(:) / sumveg(:)
1025	!
1026	z0(:) = (un - totfrac_nobio(:)) * z0(:)
1027	!
1028	DO jv = 1, nnobio
1029	!
1030	IF ( jv .EQ. iice ) THEN
1031	z0_nobio = z0_ice
1032	ELSE
1033	WRITE(numout,*) 'jv=',jv
1034	STOP 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1035	ENDIF
1036	!
1037	z0(:) = z0(:) + frac_nobio(:,jv) * LOG(z0_nobio)
1038	!
1039	ENDDO
1040	!
1041	z0(:) = EXP( z0(:) )
1042	!
1043	! Temporarily we compute the zero plane displacement height
1044	!
1045	zhdispl(:) = ave_height(:) * height_displacement
1046	!
1047	! In order to get a variable independent of the height of the
1048	! vegetation we compute what we call the effective roughness height.
1049	! This is the height over which the roughness acts. It combines the
1050	! zero plane displacement height and the vegetation height.
1051	!
1052	roughheight(:) = ave_height(:) - zhdispl(:)
1053	!
1054	END SUBROUTINE condveg_z0logz
1055	!
1056	!
1057	!
1058	SUBROUTINE condveg_z0cdrag (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, &
1059	& z0, roughheight)
1060
1061	!
1062	! 0. Declarations
1063	!
1064
1065	! 0.1 Input
1066	INTEGER(i_std), INTENT(in) :: kjpindex
1067	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget
1068	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max
1069	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio
1070	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: totfrac_nobio
1071	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
1072	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height
1073
1074	! 0.2 Output
1075	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: z0
1076	REAL(r_std), DIMENSION(kjpindex), INTENT(out) :: roughheight
1077
1078	! 0.3 Local
1079	INTEGER(i_std) :: jv
1080	REAL(r_std), DIMENSION(kjpindex) :: sumveg, ztmp, ave_height
1081	REAL(r_std), DIMENSION(kjpindex) :: d_veg, zhdispl
1082	REAL(r_std) :: z0_nobio
1083	!
1084	! The grid average z0 is computed by averaging the neutral drag coefficients.
1085	! This is pretty straight forward except for the reference level which needs
1086	! to be chosen.
1087	!
1088	! We need a reference lever high enough above the canopy else we get into
1089	! singularities of the LOG.
1090	!
1091	ztmp(:) = MAX(10., zlev(:))
1092	!
1093	z0(:) = veget(:,1) * (ct_karman/LOG(ztmp(:)/z0_bare))**2
1094	sumveg(:) = veget(:,1)
1095	ave_height(:) = zero
1096	!
1097	DO jv = 2, nvm
1098	!
1099	IF ( is_tree(jv) ) THEN
1100	! tree trunks influence the atmosphere even when there are no leaves
1101	d_veg(:) = veget_max(:,jv)
1102	ELSE
1103	! grasses only have an influence if they are really there!
1104	d_veg(:) = veget(:,jv)
1105	ENDIF
1106	!
1107	z0(:) = z0(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)z0_over_height,z0_bare)))*2
1108	sumveg(:) = sumveg(:) + d_veg(:)
1109	!
1110	ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
1111	!
1112	ENDDO
1113	!
1114	WHERE ( sumveg(:) .GT. 0.0 ) z0(:) = z0(:) / sumveg(:)
1115	!
1116	z0(:) = (un - totfrac_nobio(:)) * z0(:)
1117	!
1118	DO jv = 1, nnobio
1119	!
1120	IF ( jv .EQ. iice ) THEN
1121	z0_nobio = z0_ice
1122	ELSE
1123	WRITE(numout,*) 'jv=',jv
1124	STOP 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1125	ENDIF
1126	!
1127	z0(:) = z0(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1128	!
1129	ENDDO
1130	!
1131	z0(:) = ztmp(:) / EXP(ct_karman/SQRT(z0(:)))
1132	!
1133	! Temporarily we compute the zero plane displacement height
1134	!
1135	zhdispl(:) = ave_height(:) * height_displacement
1136	!
1137	! In order to get a variable independent of the height of the
1138	! vegetation we compute what we call the effective roughness height.
1139	! This is the height over which the roughness acts. It combines the
1140	! zero plane displacement height and the vegetation height.
1141	!
1142	roughheight(:) = ave_height(:) - zhdispl(:)
1143	!
1144	END SUBROUTINE condveg_z0cdrag
1145	!
1146	!
1147	!
1148	SUBROUTINE condveg_albcalc (kjpindex,deadleaf_cover,veget,drysoil_frac,albedo)
1149
1150	!
1151	! 0. Declarations
1152	!
1153
1154	! 0.1 Input
1155	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1156	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: deadleaf_cover !! Fraction of soil covered by dead leaves
1157	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation types
1158	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: drysoil_frac !! Fraction of visibly Dry soil(between 0 and 1)
1159
1160	! 0.2 Output
1161	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo !! Albedo, vis(1) and nir(2)
1162	! 0.3 Local
1163	! REAL(r_std),DIMENSION (kjpindex) :: alb_bare
1164	REAL(r_std),DIMENSION (nvm,2) :: alb_leaf_tmp
1165	INTEGER(i_std) :: ks, jv
1166	!
1167	!!$ DS :Correction 11/02/2011 : update 2D parameters
1168	!!$ before the components were updated but not the parameter itself!
1169	alb_leaf(1:nvm) = alb_leaf_vis(:)
1170	alb_leaf(nvm+1:2*nvm) = alb_leaf_nir(:)
1171	!!$ maybe we could use directly alb_leaf_vis and alb_leaf_nir in alb_leaf_temp
1172	!
1173	!!$ alb_leaf_tmp(:,1) = alb_leaf_vis(:)
1174	!!$ alb_leaf_tmp(:,2) = alb_leaf_nir(:)
1175	!
1176	alb_leaf_tmp(:,1) = alb_leaf(1:nvm)
1177	alb_leaf_tmp(:,2) = alb_leaf(nvm+1:2*nvm)
1178	!
1179	!
1180	DO ks = 1, 2
1181	!
1182	IF ( alb_bare_model ) THEN
1183	alb_bare(:,ks) = soilalb_wet(:,ks) + drysoil_frac(:) * (soilalb_dry(:,ks) - soilalb_wet(:,ks))
1184	ELSE
1185	! Nouvelle formulation Nathalie, sans dependance en drysoil_frac.
1186	alb_bare(:,ks) = soilalb_moy(:,ks)
1187	ENDIF
1188	!
1189	! Correction Nathalie le 12 Avril 2006 - suppression de la dependance en deadleaf_cover
1190	!albedo(:,ks) = veget(:,1) * ( (1.-deadleaf_cover(:))*alb_bare(:) + &
1191	! deadleaf_cover(:)*alb_deadleaf(ks) )
1192	albedo(:,ks) = veget(:,1) * alb_bare(:,ks)
1193
1194	! vegetation
1195	alb_veget(:,ks) = zero
1196	DO jv = 2, nvm
1197	albedo(:,ks) = albedo(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
1198	alb_veget(:,ks) = alb_veget(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
1199	ENDDO
1200	!
1201	ENDDO
1202
1203	END SUBROUTINE condveg_albcalc
1204
1205	END MODULE condveg

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

Download in other formats: