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

Last change on this file since 426 was 386, checked in by didier.solyga, 13 years ago

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