source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/condveg.f90 @ 6268

Last change on this file since 6268 was 843, checked in by didier.solyga, 12 years ago

Improved parameters documentation and presentation. The presentation is the same for all paramters so a script can automatically build a orchidee.default file.
  • Property svn:keywords set to Revision Date HeadURL Date Author Revision
File size: 43.0 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$, $Date$
9!!
10!< $HeadURL$
11!< $Date$
12!< $Author$
13!< $Revision$
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   = soils_param.nc
665  !Config If    = NOT(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  !Config Units = [FILE]
670  !
671  filename = 'soils_param.nc'
672  CALL getin_p('SOILALB_FILE',filename)
673  !
674  IF (is_root_prc) CALL flininfo(filename,iml, jml, lml, tml, fid)
675  CALL bcast(iml)
676  CALL bcast(jml)
677  CALL bcast(lml)
678  CALL bcast(tml)
679  !
680  ! soils_param.nc file is 1° soit texture file.
681  !
682  ALLOC_ERR=-1
683  ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
684  IF (ALLOC_ERR/=0) THEN
685     WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
686     STOP
687  ENDIF
688  ALLOC_ERR=-1
689  ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
690  IF (ALLOC_ERR/=0) THEN
691     WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
692     STOP
693  ENDIF
694  ALLOC_ERR=-1
695  ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
696  IF (ALLOC_ERR/=0) THEN
697     WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
698     STOP
699  ENDIF
700  ALLOC_ERR=-1
701  ALLOCATE(soilcol(iml,jml), STAT=ALLOC_ERR)
702  IF (ALLOC_ERR/=0) THEN
703     WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
704     STOP
705  ENDIF
706  !
707  IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
708  CALL bcast(lon_rel)
709  CALL bcast(lat_rel)
710  CALL bcast(lev)
711  CALL bcast(itau)
712  CALL bcast(date)
713  CALL bcast(dt)
714  !
715  IF (is_root_prc) CALL flinget(fid, 'soilcolor', iml, jml, lml, tml, 1, 1, soilcol)
716  CALL bcast(soilcol)
717  !
718  IF (is_root_prc) CALL flinclo(fid)
719  !
720  ! Mask of permitted variables.
721  !
722  mask(:,:) = zero
723  DO ip=1,iml
724     DO jp=1,jml
725        IF (soilcol(ip,jp) > min_sechiba) THEN
726           mask(ip,jp) = un
727        ENDIF
728     ENDDO
729  ENDDO
730  !
731  nbvmax = 200
732  !
733  callsign = 'Soil color map'
734  !
735  DO WHILE ( .NOT. ok_interpol )
736     WRITE(numout,*) "Projection arrays for ",callsign," : "
737     WRITE(numout,*) "nbvmax = ",nbvmax
738     !
739     ALLOC_ERR=-1
740     ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
741     IF (ALLOC_ERR/=0) THEN
742         WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
743        STOP
744     ENDIF
745     sub_area(:,:)=zero
746     ALLOC_ERR=-1
747     ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
748     IF (ALLOC_ERR/=0) THEN
749        WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
750        STOP
751     ENDIF
752     sub_index(:,:,:)=0
753     !
754     CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
755          &                iml, jml, lon_rel, lat_rel, mask, callsign, &
756          &                nbvmax, sub_index, sub_area, ok_interpol)
757     !
758     IF ( .NOT. ok_interpol ) THEN
759        DEALLOCATE(sub_area)
760        DEALLOCATE(sub_index)
761     ENDIF
762     !
763     nbvmax = nbvmax * 2
764  ENDDO
765  !
766  nbexp = 0
767  !
768  !    Check that we found some points
769  !
770  soilalb_dry(:,:) = zero
771  soilalb_wet(:,:) = zero
772  ! Ajout Nathalie
773  soilalb_moy(:,:) = zero
774  !
775  DO ib=1,nbpt
776     !
777     ! GO through the point we have found
778     !
779     !
780     fopt =  COUNT(sub_area(ib,:) > zero)
781     !
782     IF ( fopt .EQ. 0) THEN
783        nbexp = nbexp + 1
784        soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
785        soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
786        soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
787        soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
788        ! Ajout Nathalie
789        soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
790        soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
791     ELSE
792        sgn = zero
793        !
794        !   Compute the average bare soil albedo parameters
795        !
796        DO ilf = 1,fopt
797           !
798           ip = sub_index(ib,ilf,1)
799           jp = sub_index(ib,ilf,2)
800           !
801           IF ( NINT(soilcol(ip,jp)) .LE. classnb) THEN
802              soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis) + vis_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
803              soilalb_dry(ib,inir) = soilalb_dry(ib,inir) + nir_dry(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
804              soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis) + vis_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
805              soilalb_wet(ib,inir) = soilalb_wet(ib,inir) + nir_wet(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
806              ! Ajout Nathalie
807              soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis) + albsoil_vis(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
808              soilalb_moy(ib,inir) = soilalb_moy(ib,inir) + albsoil_nir(NINT(soilcol(ip,jp))) * sub_area(ib,ilf)
809              sgn = sgn + sub_area(ib,ilf)
810           ELSE
811              WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
812              STOP
813           ENDIF
814           !
815        ENDDO
816        !
817        ! Normalize the surface
818        !
819        IF ( sgn .LT. min_sechiba) THEN
820           nbexp = nbexp + 1
821           soilalb_dry(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
822           soilalb_dry(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
823           soilalb_wet(ib,ivis) = (SUM(vis_dry)/classnb + SUM(vis_wet)/classnb)/deux
824           soilalb_wet(ib,inir) = (SUM(nir_dry)/classnb + SUM(nir_wet)/classnb)/deux
825           ! Ajout Nathalie
826           soilalb_moy(ib,ivis) = SUM(albsoil_vis)/classnb
827           soilalb_moy(ib,inir) = SUM(albsoil_nir)/classnb
828        ELSE
829           soilalb_dry(ib,ivis) = soilalb_dry(ib,ivis)/sgn
830           soilalb_dry(ib,inir) = soilalb_dry(ib,inir)/sgn
831           soilalb_wet(ib,ivis) = soilalb_wet(ib,ivis)/sgn
832           soilalb_wet(ib,inir) = soilalb_wet(ib,inir)/sgn
833           ! Ajout Nathalie
834           soilalb_moy(ib,ivis) = soilalb_moy(ib,ivis)/sgn
835           soilalb_moy(ib,inir) = soilalb_moy(ib,inir)/sgn           
836        ENDIF
837        !
838     ENDIF
839     !
840  ENDDO
841  !
842  IF ( nbexp .GT. 0 ) THEN
843     WRITE(numout,*) 'CONDVEG_soilalb : The interpolation of the bare soil albedo had ', nbexp
844     WRITE(numout,*) 'CONDVEG_soilalb : points without data. This are either coastal points or'
845     WRITE(numout,*) 'CONDVEG_soilalb : ice covered land.'
846     WRITE(numout,*) 'CONDVEG_soilalb : The problem was solved by using the average of all soils'
847     WRITE(numout,*) 'CONDVEG_soilalb : in dry and wet conditions'
848  ENDIF
849  !
850  DEALLOCATE (lat_rel)
851  DEALLOCATE (lon_rel)
852  DEALLOCATE (mask)
853  DEALLOCATE (sub_index)
854  DEALLOCATE (sub_area)
855  DEALLOCATE (soilcol)
856  !
857  !
858  RETURN
859  !
860  END SUBROUTINE condveg_soilalb
861  !
862  !
863  !
864  SUBROUTINE condveg_z0logz (kjpindex, veget, veget_max, frac_nobio, totfrac_nobio, height, &
865       &                     z0, roughheight)
866
867    !
868    ! 0. Declarations
869    !
870
871    ! 0.1 Input
872    INTEGER(i_std), INTENT(in)                       :: kjpindex
873    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget
874    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max
875    REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio
876    REAL(r_std), DIMENSION(kjpindex), INTENT(in)     :: totfrac_nobio
877    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height
878   
879    ! 0.2 Output
880    REAL(r_std), DIMENSION(kjpindex), INTENT(out)    :: z0
881    REAL(r_std), DIMENSION(kjpindex), INTENT(out)    :: roughheight
882   
883    ! 0.3 Local
884    INTEGER(i_std)                                 :: jv
885    REAL(r_std), DIMENSION(kjpindex)                 :: sumveg, ave_height
886    REAL(r_std), DIMENSION(kjpindex)                 :: d_veg, zhdispl
887    REAL(r_std)                                      :: z0_nobio
888   
889    z0(:) = veget(:,1) * LOG(z0_bare)
890    sumveg(:) = veget(:,1)
891    ave_height(:) = zero
892    !
893    DO jv = 2, nvm
894       !
895       IF ( is_tree(jv) ) THEN
896          ! tree trunks influence the atmosphere even when there are no leaves
897          d_veg(:) = veget_max(:,jv)
898       ELSE
899          ! grasses only have an influence if they are really there!
900          d_veg(:) = veget(:,jv)
901       ENDIF
902       !
903       z0(:) = z0(:) + d_veg(:) * &
904            LOG( MAX(height(:,jv)*z0_over_height,z0_bare) )
905       sumveg(:) = sumveg(:) + d_veg(:)
906       !
907       ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
908       !
909    ENDDO
910    !
911    WHERE ( sumveg(:) > zero ) z0(:) = z0(:) / sumveg(:)
912    !
913    z0(:) = (un - totfrac_nobio(:)) * z0(:)
914    !
915    DO jv = 1, nnobio
916       !
917       IF ( jv .EQ. iice ) THEN
918          z0_nobio = z0_ice
919       ELSE
920          WRITE(numout,*) 'jv=',jv
921          STOP 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
922       ENDIF
923       !
924       z0(:) = z0(:) + frac_nobio(:,jv) * LOG(z0_nobio)
925       !
926    ENDDO
927    !
928    z0(:) = EXP( z0(:) )
929    !
930    ! Temporarily we compute the zero plane displacement height
931    !
932    zhdispl(:) = ave_height(:) * height_displacement
933    !
934    ! In order to get a variable independent of the height of the
935    ! vegetation we compute what we call the effective roughness height.
936    ! This is the height over which the roughness acts. It combines the
937    ! zero plane displacement height and the vegetation height.
938    !
939    roughheight(:) = ave_height(:) - zhdispl(:)
940    !
941  END SUBROUTINE condveg_z0logz
942  !
943  !
944  !
945  SUBROUTINE condveg_z0cdrag (kjpindex,veget,veget_max,frac_nobio,totfrac_nobio,zlev, height, &
946       &                      z0, roughheight)
947   
948    !
949    ! 0. Declarations
950    !
951   
952    ! 0.1 Input
953    INTEGER(i_std), INTENT(in)                       :: kjpindex
954    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget
955    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget_max
956    REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(in) :: frac_nobio
957    REAL(r_std), DIMENSION(kjpindex), INTENT(in)     :: totfrac_nobio
958    REAL(r_std),DIMENSION (kjpindex), INTENT (in)    :: zlev             !! Height of first layer
959    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: height   
960   
961    ! 0.2 Output
962    REAL(r_std), DIMENSION(kjpindex), INTENT(out)    :: z0
963    REAL(r_std), DIMENSION(kjpindex), INTENT(out)    :: roughheight
964   
965    ! 0.3 Local
966    INTEGER(i_std)                                 :: jv
967    REAL(r_std), DIMENSION(kjpindex)                 :: sumveg, ztmp, ave_height
968    REAL(r_std), DIMENSION(kjpindex)                 :: d_veg, zhdispl
969    REAL(r_std)                                      :: z0_nobio
970    !
971    ! The grid average z0 is computed by averaging the neutral drag coefficients.
972    ! This is pretty straight forward except for the reference level which needs
973    ! to be chosen.
974    !
975    ! We need a reference lever high enough above the canopy else we get into
976    ! singularities of the LOG.
977    !
978    ztmp(:) = MAX(10., zlev(:))
979    !
980    z0(:) = veget(:,1) * (ct_karman/LOG(ztmp(:)/z0_bare))**2
981    sumveg(:) = veget(:,1)
982    ave_height(:) = zero
983    !
984    DO jv = 2, nvm
985       !
986       IF ( is_tree(jv) ) THEN
987          ! tree trunks influence the atmosphere even when there are no leaves
988          d_veg(:) = veget_max(:,jv)
989       ELSE
990          ! grasses only have an influence if they are really there!
991          d_veg(:) = veget(:,jv)
992       ENDIF
993       !
994       z0(:) = z0(:) + d_veg(:) * (ct_karman/LOG(ztmp(:)/MAX(height(:,jv)*z0_over_height,z0_bare)))**2
995       sumveg(:) = sumveg(:) + d_veg(:)
996       !
997       ave_height(:) = ave_height(:) + veget_max(:,jv)*height(:,jv)
998       !
999    ENDDO
1000    !
1001    WHERE ( sumveg(:) .GT. zero ) z0(:) = z0(:) / sumveg(:)
1002    !
1003    z0(:) = (un - totfrac_nobio(:)) * z0(:)
1004    !
1005    DO jv = 1, nnobio
1006       !
1007       IF ( jv .EQ. iice ) THEN
1008          z0_nobio = z0_ice
1009       ELSE
1010          WRITE(numout,*) 'jv=',jv
1011          STOP 'DO NOT KNOW ROUGHNESS OF THIS SURFACE TYPE'
1012       ENDIF
1013       !
1014       z0(:) = z0(:) + frac_nobio(:,jv) * (ct_karman/LOG(ztmp(:)/z0_nobio))**2
1015       !
1016    ENDDO
1017    !
1018    z0(:) = ztmp(:) / EXP(ct_karman/SQRT(z0(:)))
1019    !
1020    ! Temporarily we compute the zero plane displacement height
1021    !
1022    zhdispl(:) = ave_height(:) * height_displacement
1023    !
1024    ! In order to get a variable independent of the height of the
1025    ! vegetation we compute what we call the effective roughness height.
1026    ! This is the height over which the roughness acts. It combines the
1027    ! zero plane displacement height and the vegetation height.
1028    !
1029    roughheight(:) = ave_height(:) - zhdispl(:)
1030    !
1031  END SUBROUTINE condveg_z0cdrag
1032  !
1033  !
1034  !
1035  SUBROUTINE condveg_albcalc (kjpindex,deadleaf_cover,veget,drysoil_frac,albedo)
1036
1037    !
1038    ! 0. Declarations
1039    !
1040
1041    ! 0.1 Input
1042    INTEGER(i_std), INTENT(in)                       :: kjpindex        !! Domain size
1043    REAL(r_std), DIMENSION(kjpindex), INTENT(in)     :: deadleaf_cover  !! Fraction of soil covered by dead leaves
1044    REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: veget           !! Fraction of vegetation types
1045    REAL(r_std), DIMENSION(kjpindex), INTENT(in)     :: drysoil_frac    !! Fraction of visibly Dry soil(between 0 and 1)
1046   
1047    ! 0.2 Output
1048    REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo  !! Albedo, vis(1) and nir(2)
1049    ! 0.3 Local
1050!    REAL(r_std),DIMENSION (kjpindex)                 :: alb_bare
1051    REAL(r_std),DIMENSION (nvm,2)                    :: alb_leaf_tmp
1052    INTEGER(i_std)                                   :: ks, jv
1053    !
1054    alb_leaf_tmp(:,ivis) = alb_leaf_vis(:)
1055    alb_leaf_tmp(:,inir) = alb_leaf_nir(:)
1056    !
1057    !
1058    DO ks = 1, 2
1059       !
1060       IF ( alb_bare_model ) THEN
1061          alb_bare(:,ks) = soilalb_wet(:,ks) + drysoil_frac(:) * (soilalb_dry(:,ks) -  soilalb_wet(:,ks))
1062       ELSE
1063          ! Nouvelle formulation Nathalie, sans dependance en drysoil_frac.
1064          alb_bare(:,ks) = soilalb_moy(:,ks)
1065       ENDIF
1066       !
1067       ! Correction Nathalie le 12 Avril 2006 - suppression de la dependance en deadleaf_cover
1068       !albedo(:,ks) = veget(:,1) * ( (un-deadleaf_cover(:))*alb_bare(:) + &
1069       !                              deadleaf_cover(:)*alb_deadleaf(ks)    )
1070       albedo(:,ks) = veget(:,1) * alb_bare(:,ks)
1071
1072       ! vegetation
1073       alb_veget(:,ks) = zero
1074       DO jv = 2, nvm
1075          albedo(:,ks) = albedo(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
1076          alb_veget(:,ks) = alb_veget(:,ks) + veget(:,jv)*alb_leaf_tmp(jv,ks)
1077       ENDDO
1078       !
1079    ENDDO
1080   
1081  END SUBROUTINE condveg_albcalc
1082
1083END MODULE condveg
Note: See TracBrowser for help on using the repository browser.