source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE/src_sechiba/condveg.f90 @ 55

Last change on this file since 55 was 45, checked in by mmaipsl, 14 years ago

MM: Tests with lf95 compiler : correct f95 strict norm problems.

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