source: tags/ORCHIDEE/src_sechiba/condveg.f90 @ 6

Last change on this file since 6 was 6, checked in by orchidee, 14 years ago

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